\section{Introduction}

The Helmholtz boundary value problems with high wavenumber arise from many practical fields, such as, seismology, electromagnetics, underwater acoustics, medical imaging, and so on.  It is well known that traditional numerical methods do not work well for this problem and exhibit the so-called pollution effect \cite{BI1995}.
To reduce or avoid the pollution effect, many researchers have done some various attempts, for example, \cite{ALCR08,WL2011,GM2011,FX2013,LS2017}.

In this paper, we propose the HRTMDG method \cite{ES2010,ZZ2017,ZV2019} to solve Helmholtz equation, in which, the subproblems for the flux and unknown function are solved at element level and these variables are eliminated in favor of the Lagrange multiplier, identified as  the unknown function trace at the element interfaces, and the global system involves only the degrees of freedom associated with the multiplier, significantly reducing the computational cost. By the similar technique as in \cite{ES2010,ZZ2017},  with a new energy norm, we give the convergence analysis. The corresponding error estimate shows that the HRTMDG method has an optimal $L^2$-norm convergence accuracy which is independent of  wavenumber.

\section{Formulation of HRTMDG method}
\setcounter{equation}{0}
Here we consider the following Helmholtz equation
\begin{equation}\label{eq1}
\begin{split}
\triangle  u+\kappa^2u=\tilde{f}\quad &\textrm{in}\ \Omega,\\
u=g\quad &\textrm{on}\ \partial\Omega,
\end{split}
\end{equation}
where $\Omega\subset\mathbb{R}^d$ $(d=2,3)$, a Lipschitz polyhedral domain with boundary $\partial\Omega$; $\kappa>0$ is the wave number, $\tilde{f}\in L^2(\Omega)$ is the source term and boundary value $g\in H^{1/2}(\Omega)$.

By introducing an unknown variable $\bm\sigma=i \nabla u/\kappa$ and $f=i\tilde{f}/\kappa$,
 we can rewrite our problem into  an equivalent first order formulation:
\begin{equation}\label{eq2}
\begin{split}
i\kappa\bm\sigma+\nabla u=0\quad &\textrm{in}\ \Omega,\\
i\kappa u+\nabla\cdot\bm\sigma =f\quad &\textrm{in}\ \Omega,\\
u=g\quad &\textrm{on}\ \partial\Omega.
\end{split}
\end{equation}
If $\kappa^2$ is not an eigenvalue for the above problem, there exists a unique solution $(\bm\sigma,u)\in H(\textup{div};\Omega)\times H^1(\Omega)$ of \eqref{eq2} and this solution satisfies the standard elliptic regularity (see \cite{GM2011}):
\begin{equation}\label{eq3}
\|\bm\sigma\|_{H(\textup{div};\Omega)}+\|u\|_{H^1(\Omega)}\leq C\{\|f\|_{L^2(\Omega)}+\|g\|_{H^{1/2}(\partial(\Omega)}\}.
\end{equation}

In order to construct our procedure, we first give a partitioning of the domain $\Omega$. Let  $\mathcal{T}_{h}$  be a  regular division into simplices $K$. Let $\mathcal{E}_{h}$ $
=\{e:e$ is an edge of $K$ for all $K\in \mathcal{T}_{h}\}$, $\mathcal{E}_{h}^{i}$ $=\{e:e$ is an interior edge of $K\}$, $\mathcal{E}_{h}^{o }=$ $\mathcal{E }_{h}\cap \partial
\Omega $, and denote $\mathbf{n}$ the unit outward normal vector to the element boundary $e$.

Define the inner products
\[
\begin{split}
(u,v)_K:=\int_K u\bar{v}dx,\quad (u,v) _{\mathcal{T}_h}:=\sum_{K\in\mathcal{T}_h}( u, v )_K,\\
\left\langle \lambda,\mu\right\rangle_e:=\int_e  \lambda\bar{\mu} ds,\quad \left\langle  \lambda,\mu\right\rangle_{\mathcal{E}_h}:=\sum_{e\in\mathcal{E}_h}\left\langle  \lambda,\mu\right\rangle_e,
\end{split}
\]
where the bar denotes complex conjugation. For vector-valued functions, some modification is required obviously. The corresponding norms are denoted by
\[
\|u\|^2_{\mathcal{T}_h}:=( u,u)_{\mathcal{T}_h},\quad
|\mu|^2_{\mathcal{E}_h}:=\left\langle  \mu,\mu\right\rangle_{\mathcal{E}_h}.
\]

Define the finite element spaces:
\[
\mathcal{V}_{h}=\{v_h\in L^2(\Omega): v_h|_{K}\in \mathbb{P}_k(K), \forall K\in \mathcal{T}%
_{h}\},\quad
\mathcal{W}_{h}=\{{\bm\tau}_h\in H(\textup{div};\Omega): {\bm\tau}_h|_{K}\in RT_k(K), \forall
K\in \mathcal{T}_{h}\}
\]
where $\mathbb{P}_k(K)$ is the space of complex polynomials of degree $\leq k$ on element $K$, $RT_k(K)=  [\mathbb{P}_k(K)]^d\oplus  x\mathbb{P}_k(K)$ denotes the Raviart-Thomas mixed finite element space.
 So the traditional mixed variational formulation can be read as
 \begin{equation}\label{eq4}
\begin{array}{rcll}
(i\kappa \bm\sigma_h, {\bm\tau}_h)_{\mathcal{T}_h}-(u_h,\nabla\cdot{\bm\tau}_h)_{\mathcal{T}_h}&=&-(g,{\bm\tau}_h\cdot\mathbf{n})_{\mathcal{E}^o_h} \quad&{\bm\tau}_h\in \mathcal{W}_h,\\[.1in]
(i\kappa u_h,v_h)_{\mathcal{T}_h}+(\nabla\cdot\bm\sigma_h,v_h)_{\mathcal{T}_h}&=&(f,v_h)_{\mathcal{T}_h}\quad &v_h\in \mathcal{V}_h.
\end{array}
\end{equation}

 As we know,  the classical mixed finite element formulation \eqref{eq4} can keep the mass conservation on the discrete level, but it leads to a saddle-point problem and involves considerably more degree of freedoms than a standard $H^1$-conforming method. The HRTMDG method \cite{ES2010,ZZ2017,ZV2019} can overcome these problems.  Its main ideas are: by adding appropriate constraints, one can use completely discontinuous piecewise polynomial functions and ensure the continuity of the normal fluxes over element interfaces.

To give the HRTMDG procedure of the system \eqref{eq1}, we first introduce the piecewise Sobolev spaces given by
\[
\mathcal{H}^s(\mathcal{T}_{h})=\left\{z|_K\in {H}^s(K),\ \forall\ K\in \mathcal{T}_{h}\right\},\quad s\geq 0,\quad L^2(\mathcal{E}_h)=\{\mu\in L^2(e), \forall\ e\in \mathcal{E}_h\}.
\]
Redefine the finite element approximate spaces $\mathcal{W}_h$, $\mathcal{V}_h:$
\[
\mathcal{V}_h=\{v_h\in L^2(\mathcal{T}_h):\ v_h|_K\in  \mathbb{P}_k(K),\quad\forall\ K\in\mathcal{T}_{h}\},\quad \mathcal{W}_h=\{{\bm\tau}_h\in [\mathcal{H}^k(\mathcal{T}_h)]^d:\ {\bm\tau}_h|_K\in  RT_k(K),\quad \forall\ K\in\mathcal{T}_{h}\}
\]
and define the space for Lagrange multiplier
\[
\mathcal{M}_{h}=\{ \mu_h\in L^2(\mathcal{E }_{h}): \mu_h|_{e}\in \mathbb{P}_k(e), \forall e\in \mathcal{E}^i_{h},\ \mu_h|_{e}=0,\forall e\in \mathcal{E}^0_{h}\}
\]
where $\mathbb{P}_k(e)$ is the space of complex polynomials of degree $\leq k$ on edge $e$.

The corresponding HRTMDG finite element problem can be written as follows:
 \begin{equation}\label{eq5}
\begin{array}{rcll}
(i\kappa \bm\sigma_h, {\bm\tau}_h)_{\mathcal{T}_h}-(u_h,\nabla\cdot{\bm\tau}_h)_{\mathcal{T}_h}+\langle \lambda_h,{\bm\tau}_h\cdot\mathbf{n}\rangle_{\mathcal{E}_h}&=& -(g,{\bm\tau}_h\cdot\mathbf{n})_{\mathcal{E}^o_h} \quad&{\bm\tau}_h\in \mathcal{W}_h,\\[.1in]
-(i\kappa u_h,v_h)_{\mathcal{T}_h}-(\nabla\cdot\bm\sigma_h,v_h)_{\mathcal{T}_h}&=&-(f,v_h)_{\mathcal{T}_h}\quad &v_h\in \mathcal{V}_h,\\[.1in]
\langle  \bm\sigma_h\cdot\mathbf{n}, \mu_h\rangle_{\mathcal{E}_h}&=&0\quad&\mu_h\ \in \mathcal{M}_h.
\end{array}
\end{equation}

By integrating by parts, we have

\noindent
{\bf HRTMDG method.} Find $(\bm\sigma_h,u_h,\lambda_h)\in\mathcal{W}_h\times\mathcal{V}_h\times \mathcal{M}_h$ such that
  \begin{equation}\label{eq6}
\mathcal{A}(\bm\sigma_h,u_h,\lambda_h;{\bm\tau}_h,v_h,\mu_h)=\mathcal{F}({\bm\tau}_h,v_h,\mu_h),\quad ({\bm\tau}_h,v_h,\mu_h)\in \mathcal{W}_h\times\mathcal{V}_h\times \mathcal{M}_h,
\end{equation}
where $\mathcal{A}$ and $\mathcal{F}$ are defined by
\[
\begin{split}
\mathcal{A}(\bm\sigma_h,u_h,\lambda_h;{\bm\tau}_h,v_h,\mu_h)
:=&(i\kappa \bm\sigma_h, {\bm\tau}_h)_{\mathcal{T}_h}
-(i\kappa u_h,v_h)_{\mathcal{T}_h}
+(\bm\sigma_h,\nabla v_h)_{\mathcal{T}_h}
+(\nabla u_h,{\bm\tau}_h)_{\mathcal{T}_h}
\\&+\langle \lambda_h-u_h,{\bm\tau}_h\cdot\mathbf{n}\rangle_{\mathcal{E}_h}
+\langle  \bm\sigma_h\cdot\mathbf{n}, \mu_h-v_h\rangle_{\mathcal{E}_h}
\end{split}
\]
and
\[
\mathcal{F}({\bm\tau}_h,v_h,\mu_h)= -(f,v_h)_{\mathcal{T}_h}-(g,{\bm\tau}_h\cdot\mathbf{n})_{\mathcal{E}^o_h}.
\]

On every element, we can write (\ref{eq6}) into a matrix
equation of the form
\[
\left(
\begin{array}{ccc}
\mathbf{A} & \mathbf{B} & \mathbf{D} \\
\mathbf{B}^t & \mathbf{E} & 0 \\
\mathbf{D}^t & 0 & 0 \\
\end{array}
\right) \left(
\begin{array}{c}
\bm{\sigma}_h \\
u_h \\
{\lambda}_h \\
\end{array}
\right) = \left(
\begin{array}{c}
\mathbf{F}_1 \\
\mathbf{F}_2 \\
0 \\
\end{array}
\right),
\]
where $\mathbf{F}_1$ incorporates the Dirichlet boundary data and $\mathbf{F}_2$ is the vector with respect to the right-hand side $f$. The local matrix $\mathbf{A,B,D,E}$ can be computed by their corresponding finite
element basis. $\mathbf{D}^t$ and $\mathbf{B}^t$  denote the complex conjugate transpose matrices of $\mathbf{D}$ and $\mathbf{B}$, respectively. Both the vectors of $\bm{\sigma}_h$ and $u_h$ can now be easily eliminated to obtain an equation for the multiplier only, namely,
\[
\mathbf{D}^t\mathbf{M}^{-1}\mathbf{D}\bm{\lambda_h}=\mathbf{G},
\]
where $\mathbf{M}$ and $\mathbf{G}$ are given by
\[
\mathbf{M}=\mathbf{A}-\mathbf{B}\mathbf{E}^{-1}\mathbf{B}^t,  \quad
\mathbf{G}=\mathbf{D}^t\mathbf{M}^{-1}(\mathbf{F}_1-\mathbf{B}\mathbf{E}^{-1}\textbf{F}_2).
\]
Assembling the above equation on every elements, we can get the global
systems for solving the multiplier ${\lambda}_h$.

From the above, we can see that the HRTMDG method has several
advantages: (I) compared with the discontinuous Galerkin finite element method, the number of degrees of freedom of multiplier is remarkably
small; (II) once the multiplier ${\lambda}_h$ has been obtained, $\bm{\sigma}_h$ and $u_h$ can be efficiently computed element by element; (III) the matrix $\mathbf{D}^t\mathbf{M}^{-1}\mathbf{D}$ is symmetric and positive definite,
so we can solve the systems by using the conjugate gradient method.

\begin{thm}\label{thm1}
(Consistency)
HRTMDG method (\ref{eq6})  is consistent. That is,  let $u$ be the solution of \eqref{eq1},  $\bm\sigma=i\nabla u/\kappa$, and $\lambda=u$. Then the variational equation \eqref{eq6} holds if $\bm\sigma_h$, $u_h$ and $\lambda_h$ are replaced by  $\bm\sigma$, $u$ and $\lambda$.
\end{thm}
 \begin{proof}
 Let $u$ denote the solution of \eqref{eq1}, and make substitutions as mentioned in Theorem \ref{thm1}. Taking $({\bm\tau}_h,v_h,\mu_h)=({\bm\tau}_h,0,0)$ in \eqref{eq6}, we can get
\[
\mathcal{A}(\frac{i}{\kappa}\nabla u,u,u;{\bm\tau}_h,0,0)
=-(\nabla u, {\bm\tau}_h)_{\mathcal{T}_h}+(\nabla u,{\bm\tau}_h)_{\mathcal{T}_h}
-\langle u,{\bm\tau}_h\cdot\mathbf{n}\rangle_{\mathcal{E}^o_h}
=
-\langle u,{\bm\tau}_h\cdot\mathbf{n}\rangle_{\mathcal{E}^o_h}=-(g,{\bm\tau}_h\cdot\mathbf{n})_{\mathcal{E}^o_h}.
\]

Next, choosing $({\bm\tau}_h,v_h,\mu_h)=(0,v_h,0)$ in \eqref{eq6} and using Green's formula, we have
\[
\mathcal{A}(\frac{i}{\kappa}\nabla u,u,u;0,v_h,0)
=-(i\kappa u,v_h)_{\mathcal{T}_h}
+(\frac{i}{\kappa}\nabla u,\nabla v_h)_{\mathcal{T}_h}
-\langle \frac{i}{\kappa}\nabla u\cdot\mathbf{n}, v_h\rangle_{\mathcal{E}_h}
=-\frac{i}{\kappa}(\triangle  u+\kappa^2 u,v_h)_{\mathcal{T}_h}
=-\frac{i}{\kappa}(\tilde{f},v_h)_{\mathcal{T}_h},
\]
where we have used the fact that $u$ is the solution of \eqref{eq1} in the last equation.

Finally, testing with $ ({\bm\tau}_h,v_h,\mu_h)=(0,0,\mu_h)$, we can obtain
\[
\mathcal{A}(\frac{i}{\kappa}\nabla u,u,u;0,0,\mu_h)
=\langle  \frac{i}{\kappa}\nabla u\cdot\mathbf{n}, \mu_h\rangle_{\mathcal{E}_h}=0,
\]
which implies that the normal flux $\nabla u\cdot\mathbf{n}$ is continuous across element interfaces.
 \end{proof}

 \begin{thm}\label{thm2}
(Conservation)
HRTMDG method (\ref{eq6}) is locally and globally conservative.
\end{thm}
\begin{proof}
Let $\chi_K$ denote the characteristic function of a set $K\subset{\Omega}$. Taking $ (\omega_h,v_h,\mu_h)=(0,\chi_K,0)$ in \eqref{eq6}, we can get
\[
\begin{split}
-(i\kappa u_h,1)_{K}
-\sum_{e\in\partial K}\langle  \bm\sigma_h\cdot\mathbf{n}, 1\rangle_{e}=-(f,1)_K,
\end{split}
\]
the above equation implies that HRTMDG method (\ref{eq6}) keeps local mass balance, and hence it also keeps globally mass balance.
\end{proof}

\section{Existence and uniqueness}

In order to give the existence and uniqueness of HRTMDG method (\ref{eq6}), we first define the following mesh-dependent energy norm
\[
|\|({\bm\tau}, v,\mu)\||_{A}:=\left(\kappa\|{\bm\tau}\|^2_{\mathcal{T}_h}+\kappa \|u\|^2_{\mathcal{T}_h} +\frac{1}{\kappa}\|\nabla u\|^2_{\mathcal{T}_h}+\frac{1}{\kappa h}|\mu-v|^2_{\mathcal{E}_h}\right)^{\frac{1}{2}}.
\]

Next, we show the stability and boundedness of the bilinear form $\mathcal{A}(\cdot;\cdot)$ in the sense of energy norms $|\|\cdot\||_{A}$. The following result ( see Lemma 3.1 in \cite{ES2010} ) will be used.
\begin{lem}\label{lem1}
For given $v_h\in \mathcal{V}_h$ and $\mu_h\in\mathcal{M}_h$, there exists a unique solution $\tilde{{\bm\tau}}\in \mathcal{W}_h$ such that
\begin{equation}\label{eq7}
\begin{split}
&(\tilde{{\bm\tau}},p)_K=(\nabla v_h, p)_K \quad\forall p\in [\mathbb{P}_{k-1}(K)]^d,\\
&\langle \tilde{{\bm\tau}}\cdot \mathbf{n}, q\rangle_{\partial K}= \langle \mu_h, q\rangle_{\partial K} \quad\forall q\in\mathbb{P}_{k}(\partial K),
\end{split}
\end{equation}
and the following estimate holds
\begin{equation}\label{eq8}
\|\tilde{{\bm\tau}}\|_{\mathcal{T}_h}\leq c_I\left(\|\nabla v_h\|^2_{\mathcal{T}_h}+h|\mu_h|^2_{\mathcal{E}_h}\right)^{\frac{1}{2}},
\end{equation}
where  $c_I$ is  a constant independent of mesh size $h$ and $\kappa$.
\end{lem}

\begin{lem}\label{thm3}
(Stability) There exists a positive constant $c_A$  that is independent of the mesh size $h$  and $\kappa$ such that
\begin{equation}\label{eq9}
\sup_{({\bm\tau}_h,v_h,\mu_h)}
\frac{\left|\mathcal{A}(\bm\sigma_h,u_h,\lambda_h;{\bm\tau}_h,v_h,\mu_h)\right|}
{|\|({\bm\tau}_h,v_h,\mu_h)\||_{A}}
\geq c_A|\|(\bm\sigma_h,u_h,\lambda_h)\||_{A}
\end{equation}
holds for all $(\bm\sigma_h,u_h,\lambda_h)\in\mathcal{W}_h\times\mathcal{V}_h\times\mathcal{M}_h$.
\end{lem}
\begin{proof}
Taking $v_h=\frac{1}{\kappa}u_h$ and $\mu_h=\frac{1}{\kappa h}(\lambda_h-u_h)$ in \eqref{eq7}, we can obtain

\begin{equation}\label{eq10}
\begin{split}
(\tilde{{\bm\tau}},p)_K&=( \frac{1}{\kappa}\nabla u_h, p)_K,\\
\langle \tilde{{\bm\tau}}\cdot \mathbf{n}, q\rangle_{\partial K}&= \langle \frac{1}{\kappa h}(\lambda_h-u_h), q\rangle_{\partial K},
\end{split}
\end{equation}
and
\begin{equation}\label{eq11}
\|\tilde{{\bm\tau}}\|_{\mathcal{T}_h}\leq c_I\left(\frac{1}{\kappa^2}\|\nabla u_h\|^2_{\mathcal{T}_h}+\frac{1}{\kappa^2 h}|\lambda_h-u_h|^2_{\mathcal{E}_h}\right)^{\frac{1}{2}}.
\end{equation}
For any $\epsilon>0$, we have
\[
\begin{split}
\mathcal{A}(\bm\sigma_h,u_h,\lambda_h;\epsilon\tilde{{\bm\tau}},0,0)
=&\epsilon(i\kappa \bm\sigma_h, \tilde{{\bm\tau}})_{\mathcal{T}_h}
+\epsilon(\nabla u_h,\tilde{{\bm\tau}})_{\mathcal{T}_h}
+\epsilon\langle \lambda_h-u_h,\tilde{{\bm\tau}}\cdot\mathbf{n}\rangle_{\mathcal{E}_h}\\
=&\epsilon(i\kappa \bm\sigma_h, \tilde{{\bm\tau}})_{\mathcal{T}_h}
+\epsilon\left(\frac{1}{\kappa}\|\nabla u_h\|^2_{\mathcal{T}_h}
+\frac{1}{\kappa h}| \lambda_h-u_h|^2_{\mathcal{E}_h})\right)\\
\geq &-\frac{\kappa}{2}\|\bm\sigma_h\|^2_{\mathcal{T}_h}-\frac{\kappa\epsilon^2}{2}\|\tilde{{\bm\tau}}\|^2_{\mathcal{T}_h}+\epsilon\left(\frac{1}{\kappa}\|\nabla u_h\|^2_{\mathcal{T}_h}
+\frac{1}{\kappa h}| \lambda_h-u_h|^2_{\mathcal{E}_h})\right)
\\\geq&-\frac{\kappa}{2}\|\bm\sigma_h\|^2_{\mathcal{T}_h}+\left(\epsilon-\frac{c_I^2\epsilon^2}{2}\right)\left(\frac{1}{\kappa}\|\nabla u_h\|^2_{\mathcal{T}_h}
+\frac{1}{\kappa h}| \lambda_h-u_h|^2_{\mathcal{E}_h})\right)
\end{split}
\]
where we have used \eqref{eq10} in the second equation and \eqref{eq11} in the second inequality. Set $\epsilon=1/c_I$. From \eqref{eq11} we have
\begin{equation}\label{eq12}
\mathcal{A}(\bm\sigma_h,u_h,\lambda_h;\epsilon\tilde{{\bm\tau}},0,0)
\geq-\frac{\kappa}{2}\|\bm\sigma_h\|^2_{\mathcal{T}_h}+\frac{c_I}{2}\left(\frac{1}{\kappa}\|\nabla u_h\|^2_{\mathcal{T}_h}
+\frac{1}{\kappa h}| \lambda_h-u_h|^2_{\mathcal{E}_h})\right).
\end{equation}
Note that
\begin{equation}\label{eq13}
\mathcal{A}(\bm\sigma_h,u_h,\lambda_h;i\bm\sigma_h,-iu_h,-i\lambda_h)
=\kappa \|\bm\sigma_h\|^2_{\mathcal{T}_h}
+\kappa\| u_h\|^2_{\mathcal{T}_h}.
\end{equation}
Combining the above two inequalities for the two choices of test functions, we complete our proof of Lemma \ref{thm3}.
\end{proof}

\begin{lem}\label{thm4}
(Boundedness) There exists a constant $C_A$ independent of $h$  and $\kappa$, such that,
\begin{equation}\label{eq14}
|\mathcal{A}(\bm\sigma_h,u_h,\lambda_h;{\bm\tau}_h,v_h,\mu_h)|\leq C_A|\|(\bm\sigma_h,u_h,\lambda_h)\||_{A}|\|({\bm\tau}_h,v_h,\mu_h)\||_{A}
\end{equation}
holds for all $(\bm\sigma_h,u_h,\lambda_h), ({\bm\tau}_h,v_h,\mu_h)\in\mathcal{W}_h\times\mathcal{V}_h\times\mathcal{M}_h$.
\end{lem}
\begin{proof}
Using the standard arguments, we can easily get the above estiamte.
\end{proof}

By the stability and boundedness of the bilinear form $\mathcal{A}(\cdot,\cdot)$, and Lax-Milgram theorem, we obtain the following result.
\begin{thm}
HRTMDG method \eqref{eq6} has a unique solution.
\end{thm}

\section{A uniform error estimate}
For $K\in\mathcal{T}_h$, $e\in\mathcal{E}_h$ and functions $u\in L^2(K)$ and $\lambda\in L^2(e)$, we define the local $L^2$-projection operaters $\Pi^K$ and $\Pi_e$ by
\[
\left\langle u-\Pi^Ku,v_h\right\rangle_{K}=0,\quad\forall v_h\in \mathbb{P}_{k}(K)
\]
and
\[
\left\langle \lambda-\Pi^e\lambda,\mu_h\right\rangle_{e}=0,\quad\forall \mu_h\in \mathbb{P}_{k}(e).
\]
The following error estimates hold:
\begin{equation}\label{eq15}
\begin{array}{rcll}
  \|u-\Pi^Ku\|_{K}& \leq  &  C^{*}h^s|u|_{s,K}, & 0\leq s\leq k+1,\\
  \|\nabla(u-\Pi^Ku)\|_{K}& \leq  & C^{*}h^s|u|_{s+1,K}, & 0\leq s\leq k,  \\
   \|u-\Pi^Ku\|_{\partial K}+ \|u-\Pi^e u\|_{\partial K}& \leq &   C^{*}h^{s+1/2}|u|_{s+1,K}, & 0\leq s\leq k
\end{array}
\end{equation}
where $C^{*}$ is a constant independent of $h$ and $\kappa$.

Similarly, the interpolation operators for functions on $\mathcal{T}_h$ and $\mathcal{E}_h$ are defined element-wise and are denoted by the same symbols.

We also introduce the Raviart-Thomas projection operator (see \cite{BF1991}) such that
\begin{equation}\label{eq16}
\begin{array}{rcll}
  (\bm\sigma-\Pi^{RT}\bm\sigma,p_h)_K& =  &0, &\forall p_h\in[\mathbb{P}_{k-1}(K)]^d ,  \\
 ((\bm\sigma-\Pi^{RT}\bm\sigma)\cdot\mathbf{n},\mu_h)_e&=   & 0, &\forall\mu_h\in \mathbb{P}_{k}(e),\quad e\in\partial K.
\end{array}
\end{equation}
It is well known that, there exists a constant $C^{*}$ independent of $h$ and $\kappa$ such that the following approximate properties:
\begin{equation}\label{eq17}
\begin{array}{rcll}
  \|\bm\sigma-\Pi^{RT}\bm\sigma\|_{K}+ h^{1/2}\|\bm\sigma-\Pi^{RT}\bm\sigma\|_{\partial K}& \leq  &  C^{*}h^s|\bm\sigma|_{s,K}, & 1/2\leq s\leq k+1,\\
  \|\nabla\cdot(\bm\sigma-\Pi^{RT}\bm\sigma)\|_{K}& \leq  & C^{*}h^s|\nabla\cdot\bm\sigma|_{s,K}, & 1\leq s\leq k+1.
\end{array}
\end{equation}

As we know, the error of HRTMDG method (\ref{eq6}) can be divided into two parts:  an approximate error and a discrete error. We need to estimate the discrete error.

\begin{thm}\label{thm5}
Let $(u,\bm\sigma)$  and $(\bm\sigma_h,u_h,\lambda_h)$ be the solutions of \eqref{eq2} and \eqref{eq6}, respectively. Then there exists a constant $C$ independent of $h$ and $\kappa$  such that
\begin{equation}\label{eq18}
|\|(\Pi^{RT}\bm\sigma-\bm\sigma_h,\Pi^Ku-u_h,\Pi^e u-\lambda_h)\||_{A}
\leq C\sqrt{\kappa} \|\Pi^{RT}\bm\sigma-\bm\sigma\|_{\mathcal{T}_h}.
\end{equation}
\end{thm}
\begin{proof}
Using the stability of the bilinear form and Galerkin orthogonality, we have \begin{equation}\label{eq19}
\begin{split}
c_A|\|(\Pi^{RT}\bm\sigma-\bm\sigma_h,\Pi^Ku-u_h,\Pi^e u-\lambda_h)\||_{A}
\leq &\sup_{({\bm\tau}_h,v_h,\mu_h)}
\frac{\left|\mathcal{A}(\Pi^{RT}\bm\sigma-\bm\sigma_h,\Pi^Ku-u_h,\Pi^e u-\lambda_h;{\bm\tau}_h,v_h,\mu_h)\right|}
{|\|({\bm\tau}_h,v_h,\mu_h)\||_{A}}
\\= &\sup_{({\bm\tau}_h,v_h,\mu_h)}
\frac{\left|\mathcal{A}(\Pi^{RT}\bm\sigma-\bm\sigma,\Pi^Ku-u,\Pi^eu- u;{\bm\tau}_h,v_h,\mu_h)\right|}{|\|({\bm\tau}_h,v_h,\mu_h)\||_{A}}.
\end{split}
\end{equation}

In fact, utilizing the definitions of the bilinear form $\mathcal{A}$ and projection operators, we can get
\begin{equation}\label{eq21}
\begin{split}
& \mathcal{A}(\Pi^{RT}\bm\sigma-\bm\sigma,\Pi^Ku-u,\Pi^eu- u;{\bm\tau}_h,v_h,\mu_h)
\\&\quad=(i\kappa (\Pi^{RT}\bm\sigma-\bm\sigma), {\bm\tau}_h)_{\mathcal{T}_h}
-(i\kappa (\Pi^Ku-u),v_h)_{\mathcal{T}_h}
+(\Pi^{RT}\bm\sigma-\bm\sigma,\nabla v_h)_{\mathcal{T}_h}
+(\nabla (\Pi^Ku-u),{\bm\tau}_h)_{\mathcal{T}_h}
\\&\quad\quad+\langle (\Pi^eu- u)- (\Pi^Ku-u),{\bm\tau}_h\cdot\mathbf{n}\rangle_{\mathcal{E}_h}
+\langle(\Pi^{RT}\bm\sigma-\bm\sigma)\cdot\mathbf{n}, \mu_h-v_h\rangle_{\mathcal{E}_h}
\\&\quad=(i\kappa (\Pi^{RT}\bm\sigma-\bm\sigma), {\bm\tau}_h)_{\mathcal{T}_h}
-(i\kappa (\Pi^Ku-u),v_h)_{\mathcal{T}_h}
-(\nabla\cdot(\Pi^{RT}\bm\sigma-\bm\sigma), v_h)_{\mathcal{T}_h}-(\Pi^Ku-u,\nabla\cdot{\bm\tau}_h)_{\mathcal{T}_h}
\\&\quad\quad+\langle \Pi^eu- u,{\bm\tau}_h\cdot\mathbf{n}\rangle_{\mathcal{E}_h}
+\langle(\Pi^{RT}\bm\sigma-\bm\sigma)\cdot\mathbf{n}, \mu_h\rangle_{\mathcal{E}_h}
=(i\kappa (\Pi^{RT}\bm\sigma-\bm\sigma), {\bm\tau}_h)_{\mathcal{T}_h}.
\end{split}
\end{equation}
Substituting \eqref{eq21} into \eqref{eq19}, we  get the estimate \eqref{eq18}.
\end{proof}

So, we obtain the following optimal a priori error estimate.
\begin{thm}\label{thm7}
Let $(u,\bm\sigma)$  and $(\bm\sigma_h,u_h,\lambda_h)$ be the solutions of \eqref{eq2} and \eqref{eq6}, respectively. Then there exists a constant $C$ independent of the mesh size $h$ and $\kappa$ such that
\begin{equation}\label{eq23}
\|\bm\sigma-\bm\sigma_h\|_{\mathcal{T}_h}+\|u-u_h\|_{\mathcal{T}_h}
\leq Ch^s\{|\bm\sigma|_{s,\mathcal{T}_h}+ |u|_{s,\mathcal{T}_h}\},\quad 1/2\leq s\leq k+1.
\end{equation}
\end{thm}
\begin{proof}
Using the definition of the norm $|\|\cdot\||_{A}$ and the estimate \eqref{eq18}, we get
\[
\|\Pi^{RT}\bm\sigma-\bm\sigma_h\|_{\mathcal{T}_h}+\|\Pi^Ku-u_h\|_{\mathcal{T}_h}
\leq C\|\Pi^{RT}\bm\sigma-\bm\sigma\|_{\mathcal{T}_h}.
\]
Combining the above estimate with the approximate properties \eqref{eq15} and \eqref{eq17}, we obtain the estimate \eqref{eq23}.
\end{proof}

\section*{Acknowledgments}
J. Zhang's work was supported partially by the Major Scientific and Technological Projects of CNPC under Grant (ZD2019-183-008), the Natural Science Foundation of Shandong Province (ZR2019MA015) and the Fundamental Research Funds for the Central Universities  (20CX05011A). J. Zhu's work was supported partially by the National Council for Scientific and Technological Development (CNPq).

\section*{References}


\section{Introduction}

	In the 1970s, 't Hooft proposed that QCD has a nontrivial vacuum structure and solved the $U(1)_A$ problem~\cite{Hooft1976,Hooft1976a}. The nontrivial vacuum structure suggests that there is an additional topological term which violates $CP$ symmetry in the QCD Lagrangian. The $CP$ violating term brings a free parameter denoted as $\theta$. Due to the axial anomaly, a chiral rotation of quarks changes $\theta$ and the argument of Yukawa couplings of quarks, while the sum of them is invariant. Thus, the observable strong $CP$ phase in the Standard Model (SM) is
	\begin{equation} \label{eq: SM_theta}
	\overline{\theta}=\theta+\arg \det (Y_{u} Y_{d}) ,
	\end{equation}
	where $u,d$ represent up and down type quarks respectively. The measurement of neutron electric dipole moment (EDM) suggests that $|\bar\theta| < 10^{-10}$~\cite{Baker2006}. How to understand the extreme smallness of $\overline{\theta}$ is the well-known strong $CP$ problem.

	There are several ways to solve the strong $CP$ problem. One solution is the Nelson-Barr mechanism~\cite{Nelson1984NWCV,Barr1984StSwtPQS}, which assumes $CP$ to be conserved at high energy scale, i.e. the $CP$-violating phases of the CKM matrix and $\overline{\theta}$ are both zero at the high energy scale. While at the low energy scale the $CP$-violating phase $\delta$ of the CKM matrix is reproduced and $\overline{\theta}$ is still fixed to be zero.
	An alternative way is utilizing the chiral rotation of fermions to absorb $\overline{\theta}$, when a massless quark or an additional global chiral $U(1)$ symmetry exists. However, the solution with massless quark has been disfavored by the results of Lattice QCD with $m_u= 2.2 ^{+0.5}_{-0.4}~MeV$~\cite{Gasser1982}. The additional global chiral $U(1)$ solution was first proposed by Peccei and Quinn~\cite{Peccei1977,Peccei1977a}. Later, Weinberg~\cite{Weinberg1978ANLB} and Wilczek~\cite{Wilczek1978PoSitPoI} predicted a pseudoscalar as the Goldstone boson called axion from the spontaneous breaking of this additional $U(1)_{PQ}$ symmetry. Although the original Weinberg-Wilczek axion model has been ruled out by experiments~\cite{Bardeen1987}, derived models still survive. The KSVZ~\cite{Kim1979a,Shifman1980a} and the DFSZ~\cite{Dine1981a,zhitnitskii1980possible} are the most typical models among them.

	Although the KSVZ and the DFSZ survive from experimental constrains, there exists an additional theoretical issue, the quality problem. Constraints from Astrophysics give the lower bound of the decay constant of axion to be $f_{a}\gtrsim  10^{8} \mathrm{GeV}$~\cite{RaffeltAAB}. There is a general consensus that gravitational effects generate operators suppressed by the Planck scale $M_{\mathrm{Pl}}$, which  explicitly break global symmetries. As for the axion model, the explicit breaking effect is estimated by the operators~\cite{Barr1992,Kamionkowski1992}
	\begin{equation}\label{eq: gravitational violate}
	V(\phi)=g\frac{|\phi|^{2 m} \phi^{n}}{M_{\mathrm{Pl}}^{2 m+n-4}}+h . c . .
	\end{equation}

	Although these operators are suppressed by the Planck scale, the smallness of strong $CP$ phase and high scale of $\langle\phi\rangle\sim f_{a}$ cause non-negligible effects from the lower dimension operators. Consequentially, the minimum of the scalar potential is shifted and the $CP$ phase $\overline{\theta}$ reappears.

	In order to solve this problem, one can impose a discrete symmetry $Z_N$ with large $N$ on $\phi$ so that operators with dimension less than $N$ are forbidden. For $f_a=10^{12}$ $\mathrm{GeV}$, a scale that the axion can serve as dark matter, operators with dimension $2m+n < 14$ are forbidden~\cite{Hook2018TLotSCPaA}. Another idea is to reduce the VEV $\langle\phi\rangle$ while preserving a large $f_a$. The implementation is accomplished in the multiple axion models using the alignment mechanism ~\cite{Kim2005,Choi2014a,Higaki2016} or the Clockwork mechanism~\cite{Bonnefoy2018}. On the other hand, very heavy aixon with $m_{a} \gtrsim O(100) \mathrm{MeV}$ evades the astrophysical constraints~\cite{Fukuda2015AMoVQA}, thus a small $f_a$ is allowed to relax the quality problem. Furthermore, in some composite axion models, gauge invariance forbids operators of high dimension by arranging fermions suitably ~\cite{Randall1992CAMaPSP}.

	Additional non-Abelian gauge group is often introduced in various new physics, such as Hidden Valley~\cite{Strassler2006EoaHVaHC}, Vectorlike  Confinement~\cite{Kilic2009VCatL} and Twin Higgs~\cite{Chacko2006NEBfaMS}. These extensions of the SM are usually motivated by the Hierarchy problem or dark matter. Furthermore, the $U(1)$ global symmetry to solve the strong $CP$ problem could be the remnant of multi-$U(1)$ symmetry from the hidden strong dynamics. We call this kind of non-Abelian gauge groups as ``hidden QCD" in this paper.

	Similar to the SM QCD, the hidden QCD may also contain the $CP$ violation source from the $\theta'G'\tilde{G'}$ term and the relevant fermion Yukawa couplings, though it may or may not influcence the observed neutron EDM. Besides the QCD theta term, there are three types of effective operators directly contributing to the neutron EDM, including the quark EDM, the quark chromo EDM and the three-gluon Weinberg operator~\cite{Weinberg1989LHBETitNEDM,Pospelov2005EDMAPoNP}:
	\begin{equation}\begin{aligned}
	O_d&=-\frac{i}{2}  d  \bar{q} \sigma^{\mu \nu} \gamma_{5} q F_{\mu \nu}, \\
	O_{\tilde{d}}&=-\frac{i}{2}  \tilde{d}  \bar{q} \sigma^{\mu \nu} t^{a} \gamma_{5} q G_{\mu \nu}^{a}, \\
	O_w&=\frac{1}{3} w f^{a b c} G_{\mu \nu}^{a} \widetilde{G}_{\nu \beta}^{b} G_{\beta \mu}^{c},
	\end{aligned}\end{equation}
	where $d$ denotes EDM, $\tilde{d}$ denotes chromo EDM, $t^a$ is the generator of the QCD group and $f^{abc}$ denotes the QCD structure constant. The $\theta'G'\tilde{G'}$ term in hidden QCD is not directly related to these operators. However, two typical scenarios will result in the sensitivity of the neutorn EDM to the CP violation in the hidden sector:
	\begin{itemize}
		\item  Two strong $CP$ angels can be associated with each other by introducing a pseudo-scalar. The interaction between the pseudo-scalar and gauged fermions and the chiral rotation of the pseudo-scalar are
		\begin{equation}\label{first way}
		\begin{aligned}
		\mathcal{L}&\sim y_q e^{ina/f_a} q\bar{q} + y_Q e^{ima/f_a} Q\bar{Q}+h.c. , \\
		a/f_a&\rightarrow a/f_a+\alpha,
		\qquad \theta \rightarrow \theta-n\alpha,
		\qquad \theta' \rightarrow \theta'-m\alpha,
		\end{aligned}
		\end{equation}
		where $q$ and $Q$ notate fermions charged under QCD and hidden QCD respectively. $a$ is the pseudo-scalar. $n$ and $m$ are constants.
		\item Fermions charged under both QCD and hidden QCD can link these two $CP$ phases. The chiral rotation of these fermions can change both two phases as
		\begin{equation}\label{second way}
		\psi_{L}(N_c,N_h)\rightarrow \psi_{L}(N_c,N_h) e^{i\beta} , \qquad \theta \rightarrow \theta-N_h\beta ,
		\qquad \theta' \rightarrow \theta'-N_c\beta,
		\end{equation}
		where $N_c$ and $N_h$ are representation of fermions in QCD and hidden QCD respectively.
	\end{itemize}
	Under such circumstance, these two angles are needed to be small simultaneously to explain the strong $CP$ problem. An additional global chiral $U(1)$ symmetries is not enough.
    Therefore, two global chiral $U(1)$ symmetries are needed, or there is one chiral $U(1)$ symmetry with a $Z_2$ among the two sectors. In this work, we will focus on the first scenario, while the later one is discussed in Ref.~\cite{Rubakov1997,Berezhiani2001a,Albaid2015SCaS,Fukuda2015AMoVQA,Hook2014ASttSCP,Chiang2016}.

	Chiral $U(1)$ symmetries are classified into two kinds, according to whether the corresponding axion is elementary or composite.
	We refer to the one generating the elementary axion as the $U(1)_{PQ}$ arinsing from the phase of a complex scalar and the one generating composite axion as (also called dynamical axion ~\cite{Choi1985DA}) the $U(1)_A$ arising when massless fermions of the hidden QCD condensate at the high energy scale.

	In both scenarios, axions are the Goldstones of these chiral $U(1)$ symmetries. Instanton effects explicitly break these symmetries and determine the properties of axion such as mass and axion-photon coupling. Therefore, new instanton effects of the hidden QCD could enlarge axion mass in some models with $Z_2$ symmetry~\cite{Rubakov1997,Fukuda2015AMoVQA,Gherghetta2016a}
	and  it could also enhance axion-photon coupling in~\cite{Agrawal2017ETfPCotQA}. Moreover, an additional strong dynamics also provides solutions to quality problem and domain wall problem~\cite{Sikivie1982ADWatEU}.

	With hidden QCD introduced, we find new mixed $U(1)$ solutions which contain two sorts of chiral $U(1)$ symmetries. In these cases, spontaneous symmetry breaking of $U(1)_{PQ}$ could be triggered by the dynamical symmetry breaking in the hidden sector, besides the conventional Ginzburg-Landau potential method. The cancellation of $CP$ phases is always viable, except alignment situation. The cancellation of $CP$ phases in the QCD and the hidden QCD sector is shown in section~\ref{section3}. The lightest axion in these models are similar to the QCD axion. Some phenomena of these models are discussed in section~\ref{section4}.

	In addition to the new mixed $U(1)$ solution mentinoed above, we also propose a ``moose-like'' diagram method to visualize the cumbersome relations among gauge groups, new fields and CP phases in multi-axion models, from which the $U(1)$ charges for the fields and the potential of the relevant axions and CP phases can be easily read. Moreover, the diagram method helps to construct models containing more $U(1)$ and gauge groups.

	This paper is organized as follows. In section~\ref{section2}, we discuss two typical patterns to realize Chiral $U(1)$. Next, in section~\ref{section3}, we give solutions of the strong $CP$ problem, and propose a diagram method to present these solutions.
	In section~\ref{section4}, we study the axion mass, axion-photon coupling and axion deacy constant of one-axion solutions and two typical mixed two-axion solutions.
	Finally, we summarize the result in section~\ref{section5}.

	\section{Pecci-Quinn or Dynamical Solution}
	\label{section2}
	A large class of models solve the strong $CP$ problem by adding additional global $U(1)$ symmetries. There are mainly two kinds of $U(1)$s: the one in the Pecci-Quinn mechanism~\cite{Peccei1977,Peccei1977a}, denoted as $U(1)_{PQ}$, is associated with a elementary scalar $\phi$; the other is the axial $U(1)_A$ of some massless fermions charged under gauge groups, which induces the dynamical solution~\cite{Choi1985DA}.

	Although $U(1)_{PQ}$ and $U(1)_A$ symmetries are different, they solve the strong $CP$ problem with the same philosophy: to make the $CP$ phase dynamically cancelled by introducing a $U(1)$ pseudo-Goldstone boson (PGB) with an anomaly-induced potential. For both scenarios, there are fermions axially charged under the global $U(1)$,
	\begin{equation}\label{anormalous transformation}
	 f_L \rightarrow e^{i\alpha}f_L, \qquad f_R \rightarrow e^{-i\alpha}f_R.
\end{equation}
	However, the $U(1)$ is broken in different ways: for Pecci-Quinn mechanism, it is broken by the scalar potential $V(\phi)$, and the PGB turns out to be $a \sim\arg\phi$, the axion; in the dynamical solution, it is broken by condensation of the fermions, and the PGB is a composite of the fermion, also known as the dynamical axion. The common infrared behavior is that both PGB correspond to some axial currents, whose conservation is broken only by anomaly. It inspires a general description of these models, as presented in the following.

	In this section, we propose a ``moose-like'' diagram method to uniformly illustrate the structures of the models solving the strong $CP$ problem. The notation of all possible additional fermions introduced are listed in TABLE~\ref{tab: fermion classify}. As an example, we show the structure of the aligned axion model(left)~\cite{Kim2005} and the dynamical axion model(right)~\cite{Choi1985DA} in Fig.~\ref{fig:figure1}.  At the left end, we list all the axions, each representing a global $U(1)$, normalized as their proper contributions to the CP phase. Dashed lines link the axions to the fermions, and the numbers on  them indicate the corresponding $U(1)$ charges of the fermions. The solid lines with tags show the non-trivial representations of the fermions under the linked gauge groups at the right end, like the color $SU(3)_c$ and the hidden $SU(3)_h$. Circles ($\circ$) on the vertices represent massive fermions; Crosses ($\times$) on the vertexes represent massless fermions.

	\begin{table}[htbp]
		\caption{Fermion charge and representation}
		\centering
		\begin{tabular}{c|ccc}
			\hline
			& {$SU(3)_c$} & {$SU(3)_h$} & $U(1)_{PQ}$ \\
			\hline
			\hline
			$\psi$ & 3     & 3     & 0 \\
			$\chi$ & 1     & 3     & 0 \\
			\hline
			\hline
			$Q$   & 3     & 1     & m \\
			$Q'$  & 1     & 3     & n \\
			\hline
		\end{tabular}%

		\label{tab: fermion classify}%
	\end{table}%

	\begin{figure}[h]
		\centering
		\includegraphics[width=1\linewidth]{fig/figure1}
		\caption{ In diagram (a), the solid line between Q and $SU(3)_c$ shows that Q is under representation 3 of QCD and the solid line between $Q^\prime$ and $SU(3)_h$ shows that Q is under representation 3 of hidden QCD. The red dashed line between $a_1$ and Q($Q^\prime$) shows that the charge of Q($Q^\prime$) under $U(1)_{PQ}$ is $m_1$($n_1$). Similarly, blue dashed line shows that the charge of Q($Q^\prime$) under $U(1)_{PQ^\prime}$ is $m_2$($n_2$). $1/f_{a_1}$ and $1/f_{a_2}$ are scale factors for two $U(1)$s. In diagram (b), solid lines show that $\psi$ is under representation 3 for both QCD and hidden QCD, and $\chi$ is a singlet for QCD but under representation 3 for hidden QCD. Red dashed lines show that the ratio of $\psi$ and $\chi$ in $\eta_h^\prime$ current is 1:1, which is similar to $U(1)$ charge in the left one at infrared region. The only difference is one more normalization factor $\sqrt{2}/2$ in front of the scale factor $1 / f_{\pi_h}$. Similarly, blue lines show that the ratio of $\psi$ and $\chi$ in $\pi_h$ current is 1:-3. The normalization factor is $\sqrt{6}/6$ and the scale factor is $1 / f_{\pi_h}$.}
		\label{fig:figure1}
	\end{figure}

	The alignment axion model has been widely discussed in the literature~\cite{Kim2005,Choi2014a,Kappl2014a}. $SU(3)_h$ is introduced with a new free parameter of $CP$-violation $\theta'$. As shown in Fig.~\ref{fig:figure1}(a), massive fermions $Q(3,1)$ and $Q'(1,3)$ are charged under both $U(1)_{PQ}$ symmetry and $U(1)_{PQ'}$ symmetry. When $U(1)_{PQ}$ is broken, axion is left as a pseudo-Goldstone boson. The corresponding current is
	\begin{equation}\label{eq: Q current}
	\begin{aligned}
	J^\mu_{PQ}&=f_{a_1} \partial_{\mu}a_1+m_1\overline{Q}\gamma^\mu\gamma_5Q+n_1\overline{Q'}\gamma^\mu\gamma_5Q', \\
	J^\mu_{PQ'}&=f_{a_2} \partial_{\mu}a_2+m_2\overline{Q}\gamma^\mu\gamma_5Q+n_2\overline{Q'}\gamma^\mu\gamma_5Q', \\
	\partial_{\mu}J^\mu_{PQ}&=\frac{m_1g^{2}}{16 \pi^{2}} G_{a}^{\mu \nu} \tilde{G}_{a \mu \nu}+\frac{n_1{g'}^{2}}{16 \pi^{2}} {G'}_{A}^{\mu \nu} \tilde{G'}_{A \mu \nu}, \\
	\partial_{\mu}J^\mu_{PQ'}&=\frac{m_2g^{2}}{16 \pi^{2}} G_{a}^{\mu \nu} \tilde{G}_{a \mu \nu}+\frac{n_2{g'}^{2}}{16 \pi^{2}} {G'}_{A}^{\mu \nu} \tilde{G'}_{A \mu \nu}.
	\end{aligned}
	\end{equation}

	where $G_{a}^{\mu \nu}$ and ${G'}_{A}^{\mu \nu}$ are strength tensors of gauge fields for QCD and hidden QCD respectively. $g$ and $g'$ are couplings of gauge interactions. In some papers, the hidden sector is extended to a ``mirrored SM'' with massive fermions $u'$ and $d'$ that do not carry $U(1)_{PQ}$ charge like u and d in the SM~\cite{Gherghetta2016a,Hook2014ASttSCP}.

	The main structure of the dynamical solution~\cite{Choi1985DA} is shown in Fig~\ref{fig:figure1}(b).  Similar to the aligned axion model, $SU(3)_h$ is introduced with a new $CP$ parameter $\theta'$. Two massless fermions $\psi$ and $\chi$ are introduced to absorb $CP$-violating phases $\theta$ and $\theta'$ through a $U(1)_A$ transformation. It's common to assume that $SU(3)_h$ will confine just like QCD at a scale $f_{\pi_h}$ which is much higher than QCD confining scale $f_\pi$. The related currents are
	\begin{equation}\label{eq: anomalous current1}
	\begin{aligned}
	J^\mu_{A}(\psi)&=\overline{\psi}\gamma^\mu\gamma_5\psi, \quad J^\mu_{A}(\chi)=\overline{\chi}\gamma^\mu\gamma_5\chi  \\
	\partial_{\mu}J^\mu_{A}(\psi)&=\frac{3g^{2}}{16 \pi^{2}} F_{a}^{\mu \nu} \tilde{F}_{a \mu \nu}+\frac{3g'^{2}}{16 \pi^{2}} G_{A}^{\mu \nu} \tilde{G}_{A \mu \nu}, \\
	\partial_{\mu}J^\mu_{A}(\chi)&=\frac{g'^{2}}{16 \pi^{2}} G_{A}^{\mu \nu} \tilde{G}_{A \mu \nu},
	\end{aligned}
	\end{equation}

Below the scale $f_{\pi_h}$, we assume the hidden sector has a dynamical chiral symmetry breaking caused by the fermion condensate with
	\begin{equation}\label{condensate}
	\langle\bar{\psi} \psi\rangle \approx - c_\psi f_{\pi_h}^{3}, \qquad \langle\bar{\chi} \chi\rangle \approx - c_\chi f_{\pi_h}^{3},
	\end{equation}
	where $c_\psi$,$c_\chi$ are constants and $f_{\pi_h} \gg f_\pi$.  Considering QCD as an additional ``flavor symmetry'', there is a $SU(4)_L\times SU(4)_R$ symmetry between $\psi$ and $\chi$. After the condensate of hidden QCD, $\pi_h^a$ are Goldstones corresponding to   coset generators of $SU(4)_L\times SU(4)_R/SU(4)_V$. The decomposition of $\pi_h^a$ into $SU(3)_c$ is
	\begin{equation}
	15=8+3+\overline{3}+1.
	\end{equation}
	The QCD color-singlet scalar is denoted as $\pi_h\equiv \pi_h^{15}$ for short. And the rest 14 colored scalars are supposed to be heavy because of QCD condensate. However, there is one more color-singlet scalar $\eta_{h}^{\prime}$ related to the $U(1)_A$ symmetry. The corresponding currents related to these two fields are
	\begin{equation}\label{current 12}
	\begin{aligned}
	J_\mu(\pi_h)&=\dfrac{1}{\sqrt{6}}(\overline{\psi}^c \gamma_\mu\gamma_{5} \psi_c-3\overline{\chi} \gamma_\mu\gamma_{5}\chi), \\
	J_\mu(\eta_{h}^{\prime})&=\dfrac{1}{\sqrt{2}}(\overline{\psi}^c \gamma_\mu\gamma_{5} \psi_c+\overline{\chi} \gamma_\mu\gamma_{5}\chi),
	\end{aligned}
	\end{equation}
	where `c' is the color index for $\psi$. These currents have been normalized and we have chosen $tr(T^aT^b)=2\delta^{ab}$ as normalization of $SU(4)$ generators. The specific matrices of $SU(4)$ generators and details of derivation are shown in Appendix~\ref{A1}. As a result of condensate, we can use $\pi_h$ and $\eta_{h}^{\prime}$ standing for the $U(1)_A$ symmetry at low energy scale. $\theta$ and $\theta'$ can be offseted when $\pi_h$ and $\eta_{h}^{\prime}$ take VEVs.

	\section{Cancellation of the Strong $CP$ Phase}
	\label{section3}

	When considering the extension of the SM, new contributions to the strong $CP$ violation are introduced. Therefore, the strong $CP$ phase $\overline{\theta}$ needs to be modified. For distinction, the modified strong $CP$ phase is notated as $\theta_{phy}$. In this section, we discuss these new contributions and methods to solve the strong $CP$ problem. The difficulty is that $\overline{\theta}$ is at order $ 1$ but $\theta_{phy}$ is at order $  10^{-10}$. Axion models and dynamic solutions both solve this difficulty by canceling $\theta_{phy}$ at the minimal point of goldstone potential, which links $\theta_{phy}$ with $\overline{\theta}$.

	To calculate $\theta_{phy}$ which describes $CP$ violation effects at the scale lower than $\Lambda_{QCD}$, Chiral Perturbation Theory (ChPT) is used to match fields of quarks with fields of hadrons. Indeed, it is sufficient to consider ChPT with two light quarks $u$ and $d$. The Lagrangian for QCD with an axion at hadron level is~\cite{Kuster2008A}
	\begin{equation}
	\mathcal{L}_2=\dfrac{f_\pi^2}{4} \mathrm{Tr} \left( \partial_{\mu}U\partial^\mu U^\dagger\right) +A f_\pi^3 \mathrm{Tr} \left(MU^\dagger+UM^\dagger \right)+B f_\pi^4\left(\theta+ \frac{a}{f_a}+\frac{i}{2} \mathrm{Tr}\left(\log U-\log U^{\dagger}\right)\right)^2  ,
	\end{equation}
	with
	\begin{equation}
	U=\exp\left[  \dfrac{ i\left( \pi^a\tau^a+I \eta \right) }{f_\pi} \right], \qquad M=
	\begin{pmatrix}
	m_u & 0 \\
	0 & m_d
	\end{pmatrix}  ,
	\end{equation}
	where $\pi^a$, $\eta$ are mesons and $a$ is the aixon. $A$ and $B$ are dimensionless parameters matched by meson masses. Here $\tau^a$ are Pauli matrices and $I$ is the identity matrix. For convenience, we set the matrix of quark mass to be real in the rest of this paper, leading to $\theta=\overline{\theta}$. The potential part is shown as
	\begin{equation}\label{eq: QCD potential}
	V = -2Af_\pi^3 \left[ m_u\cos(\frac{\pi^0+\eta}{f_\pi})+m_d\cos(\frac{\pi^0-\eta}{f_\pi})\right] +Bf_\pi ^4\left( \frac{2\eta}{f_\pi}+\theta+ \frac{a}{f_a}\right) ^2.
	\end{equation}
	For convenience, $\left\langle\pi^0\right\rangle/f_\pi$, $\left\langle\eta\right\rangle/f_\pi$ and $\left\langle a\right\rangle/f_a$ are defined as phases $\phi_1$, $\phi_2$ and $\phi_a$ respectively.
	Considering $CP$ transformation, $\pi^0$, $\eta$ and $a$ are changed into $-\pi^0$, $-\eta$ and $-a$. Once $CP$ is conserved in eq.~\eqref{eq: QCD potential}, these phases need to meet
	\begin{equation}\label{$CP$ conserved1}
	\phi_1+\phi_2=0, \qquad \phi_1-\phi_2=0, \qquad 2\phi_2+\theta+\phi_a=0,
	\end{equation}
	which is equivalent to
	\begin{equation}\label{$CP$ conserved2}
	\phi_1=\phi_2=0, \qquad \theta+\phi_a=0.
	\end{equation}
	Considering the derivative of the effective potential in eq.~\eqref{eq: QCD potential} with respect to $\pi^0$, $\eta$ and $a$, eq.~\eqref{$CP$ conserved2} is exactly the solution at the minimal point of the potential. New contributions from $U(1)$ symmetries are described by phase $\phi_a$. Finally, the $CP$-violating observable $\theta_{phy}$ is gotten as
	\begin{equation}\label{eq:theta}
	\theta_{phy}=\theta+\phi_a=\theta+\frac{\left\langle a \right\rangle}{f_a}.
	\end{equation}
	The strong $CP$ problem is solved by the offset of axion VEV and the phase $\theta$ to ensure $\theta_{phy}=0$. Axion mass can also be derived from eq.~\eqref{eq: QCD potential} as
	\begin{equation}\label{axion mass}
	m_{a}^2 =m_{\pi}^2 \frac{f_\pi^2}{f_a^2}\frac{m_u m_d}{(m_u+m_d)^2}.
	\end{equation}
	From this QCD axion model, it can be found that $\theta_{phy}$ needs to be modified due to the new contribution from axion. Similar absorption happens in other axion models and dynamical solutions.

	For a general hidden QCD model with new fermions added, two UV free parameters, $ \theta$ and $ \theta'$, have to be introduced. In general cases, the potential term induced by the instanton effect is
	\begin{equation}\label{eq: theta_thetaprime potential}
	V_\theta = B_1f_\pi ^4\left( \sum_{i=1}^2 m_i \frac{a^i}{f_{a_i}}+\theta\right) ^2+B_2f_{\pi_h} ^4\left(  \sum_{i=1}^2 n_i \frac{a^i}{f_{a_i}}+\theta' \right) ^2,
	\end{equation}
	where $m$ and $n$ stand for different anomalous charges.
	$\footnote {It's more complex for dynamical solutions, because normalization factor also has to be included in this $m$ and $n$.}$

	Strong $CP$ problems for both QCD and hidden QCD can be solved only satisfying
	\begin{equation}\label{eq: theta_hid_QCD}
	\begin{aligned}
	\theta_{phy}&=\theta+  \frac{\left\langle \eta ' \right\rangle}{f_\pi}  + \theta_{U(1)} =0,\\
	\theta'_{phy}&=\theta'+\theta'_{U(1)} =0.
	\end{aligned}
	\end{equation}
	In the preceding equation, $ \left\langle \eta ' \right\rangle$ means the contribution from VEV of $\eta'$, which can be absorbed by redefinition of $\theta$. We are not interested in it, so we will absorb it with $\theta$ in the following discussion. $ \theta_{U(1)}$ and $ \theta'_{U(1)}$ come from axion VEVs.
	In following subsections, we discuss the $\theta$ cancellation of eq.~\eqref{eq: theta_hid_QCD} in some specific models with different $U(1)$ symmetries.

	\subsection{One U(1)}
	As shown in eq.~\eqref{eq: theta_hid_QCD}, there are two degrees of freedom needed to absorb both $\theta$ and $ \theta'$. When only one $U(1)$ is introduced, there is only one new degree of freedom. That means $\theta$ and $ \theta'$ cannot be absorbed simultaneously and the strong $CP$ problem cannot be solved.
	To solve this problem, a mirrored $Z_2$ symmetry has to be introduced between $SU(3)_{c}$ in the SM and $SU(3)_h$ in the hidden sector.
	In UV theories, some models can achieve this by embedding these two $SU(3)$ into a larger gauge symmetry group $SU(6)$~\cite{Gherghetta2016a}. Consequently, two $CP$ phases are forced to be identical, which means
	\begin{equation}\label{Z_2}
	\theta= \theta' .
	\end{equation}

	As shown in Fig.~\ref{fig:oneu1}, the additional $U(1)$ could be achieved by two methods. One method introduces an axion with $U(1)_{PQ}$(left) and the other uses massless $\psi(3,3)$ with $U(1)_{A}$(right). There are detailed discussions about both two methods in literature~\cite{Rubakov1997,Berezhiani2001a,Hook2014ASttSCP,Fukuda2015AMoVQA,Albaid2015SCaS,Chiang2016}. For the left one, cancellation eq.~\eqref{eq: theta_hid_QCD} becomes
	\begin{equation}
	\theta_{phy}^{\prime}= \theta_{phy}=\theta+\frac{\left\langle a\right\rangle}{f_a}=0 .
	\end{equation}
	For the right one, as shown in Appendix~\ref{A1}, the cancellation eq.~\eqref{eq: theta_hid_QCD} is deduced to
	\begin{equation}
	\theta_{phy}^{\prime}= \theta_{phy}= \theta+\sqrt{6}\dfrac{\left\langle \eta'_{h} \right\rangle }{f_{\pi_h}} =0 .
	\end{equation}
	In both cases, the strong $CP$ problem will be solved by only one $U(1)$ symmetry. And the $CP$-violating effects from $\theta$ and $\theta'$ will be counteracted by $ \left\langle a\right\rangle$ or $ \left\langle \eta'_{h} \right\rangle $.

	\begin{figure}[h]
		\centering
		\includegraphics[width=1\linewidth]{fig/figure2}
		\caption{(a): Structure of $U(1)_{PQ}$ models with mirror symmetry. (b): Structure of $U(1)_{A}$ dynamic solutions with mirror symmetry.}
		\label{fig:oneu1}
	\end{figure}

	\subsection{Two Similar U(1)s}

	As in the previous subsection, models with two $U(1)$ freedoms can also be sorted by the different mechanism of $U(1)$ symmetries. As shown in Fig.~\ref{fig:figure1}, models could contain two $U(1)_{PQ}$ or two $U(1)_{A}$ symmetries. For these kinds of models, axial currents are determined by representations of fermions in $SU(3)_c$ and $SU(3)_h$. The cancellation equation of $\theta$ and $\theta'$ will be determined by these currents and different $U(1)$ charge for fermions.

	Here shows another advantage of our diagram method that the cancellation equation can be read easily. Fig.~\ref{fig:figure1}(a) shows models with two $U(1)_{PQ}$ symmetries. Choosing $SU(3)_c$ as the starting point responsible for $\theta$, we can find all possible links end at left side. Notice these links can only go from right to left. The cancelling phase of each link is the product of the charge number on the dashed line and the scale factor at the end point on the left. Finally we can sum all offset terms from different links together. For $\theta'$, the starting point will be changed into $SU(3)_h$, then repeat these steps. The cancellation equations from our diagram read
	\begin{equation}\label{eq: theta_axion_}
	\begin{aligned}
	\theta_{phy}&=\theta+ m_1\dfrac{\left\langle a_1\right\rangle }{f_{a_1}} +m_2\dfrac{\left\langle a_2\right\rangle }{f_{a_2}}=0, \\
	\theta_{phy}^{\prime}&=\theta'+ n_1\dfrac{\left\langle a_1\right\rangle }{f_{a_1}}+n_2\dfrac{\left\langle a_2\right\rangle }{f_{a_2}}=0 .
	\end{aligned}
	\end{equation}
	This two-axion model is widely discussed. When the angle between vectors $(m_1,m_2)$ and $(n_1,n_2)$ is near zero (but not zero) on the contrast to the large m and n, it is exactly the `alignment axion' model~\cite{Kim2005,Choi2014a}.

	Models with two $U(1)_A$ symmetries provided by massless fermions have similar expressions for $CP$-violating angles. In Fig.~\ref{fig:figure1}(b), $M_\psi=0$ and $M\chi=0$ are supposed. The corresponding currents has been given in eq.~\eqref{current 12}. It's a little more complex to read cancellation equations for this model from the figure. There will be two more things to notice. First, the normalization factor has been already written together with scale factor. Second, the degeneracy should be multiplied for each link, which comes from the color index of the other $SU(3)$ outside the link. The degeneracy can be read from the number for representation on other solid lines attached to the fermion field in the link. For example, the degeneracy of the link, ``$SU(3)_c \rightarrow \psi \rightarrow \pi_h$", will be the number `3' on the solid line attached $\psi$ and $SU(3)_h$. When sum all links together, this link should be multiplied by 3. The final cancellation equations are

	\begin{equation}
	\begin{aligned}
	\theta_{phy}&=\theta+ \frac{\sqrt{6}}{2}\dfrac{\left\langle \pi_h\right\rangle }{f_{h}} +\frac{3\sqrt{2}}{2}\dfrac{\left\langle \eta'_{h}\right\rangle }{f_{h}}=0, \\
	\theta_{phy}^{\prime}&=\theta'+ 2\sqrt{2}\dfrac{\left\langle \eta'_{h}\right\rangle }{f_{h}}=0 .
	\end{aligned}
	\end{equation}
	The cancellation equations of other figures can also be read in this way.

	\subsection{Two Different U(1)s}

	In these two-$U(1)$ solutions, $U(1)$s need not to be the same, which means models of one $U(1)_{PQ}$ and one $U(1)_A$ is possible. Here we show two simple examples of that new kind of models in Fig.~\ref{fig: two U_1}.

	\begin{figure}[h]
		\centering
		\includegraphics[width=1\linewidth]{fig/figure3}
		\caption{Structure of models with mixed $U(1)$s. (a) Model A: the spontaneous breaking of $U(1)_{PQ}$ is independent with hidden QCD. (b) Model B: the dynamical symmetry breaking of hidden QCD induces $U(1)_{PQ}$ breaking.}
		\label{fig: two U_1}
	\end{figure}

	The mixed solution of $U(1)_{PQ}$ and $U(1)_A$ would be more complex, for the reason that two independent scales, $f_{h}$ and $f_a$, are simultaneously involved. $f_{h}$ is the condensate scale of hidden QCD fermions, while $f_a$ is the spontaneously breaking scale of the Peccei-Quinn symmetry. Fig.~\ref{fig: two U_1}(a) shows model A in which $U(1)_{PQ}$ is spontaneous broken independently with the hidden QCD. The result of cancellation equations in this model is

	\begin{equation}
	\begin{aligned}
	\theta_{phy}&=\theta+\sqrt{6}\dfrac{\left\langle \eta'_{h} \right\rangle }{f_{h}}+ m\dfrac{\left\langle a\right\rangle }{f_a}=0, \\
	\theta_{phy}^{\prime}&=\theta'+\sqrt{6}\dfrac{\left\langle \eta'_{h} \right\rangle }{f_{h}}+ n\dfrac{\left\langle a\right\rangle }{f_a}=0 .
	\end{aligned}
	\end{equation}
	Fig.~\ref{fig: two U_1}(b) shows model B in which the spontaneous breaking of $U(1)_{PQ}$ is induced by the hidden QCD . The cancellation equations are
	\begin{equation}
	\begin{aligned}
	\theta_{phy}&=\theta+\dfrac{3}{\sqrt{2}}\dfrac{\left\langle \eta'_{h} \right\rangle }{f_{h}}+\dfrac{3}{\sqrt{6}}\dfrac{\left\langle \pi_h \right\rangle }{f_{h}}=0, \\
	\theta_{phy}^{\prime}&=\theta'+2\sqrt{2}\dfrac{\left\langle \eta'_{h} \right\rangle }{f_{h}}+ n\dfrac{\left\langle a\right\rangle }{f_a}=0 .
	\end{aligned}
	\end{equation}
	It seems that redundancy exists in these equations, as there are three VEVs to absorb two $\theta$ parameters. However, there is an additional equation constraining these three VEVs from the Yukawa term of $Q'$, as shown in~\eqref{L_mass} in the next section.

	\section{Phenomenology}
	\label{section4}
	In this section, we study phenomena in the one-axion solution and the mixed two-axion solution, which can not directly be acquired from diagrams in the last section. We compare the axion masses in these two kinds of solutions. We also discuss the relation between the physical axion decay constant $F_a$ and  the axion-photon coupling $g_{a\gamma\gamma}$ in the mixed two-axion solution.

	\subsection{One-axion Solution}
	The one axion solution is exactly the implementation of  the visible axion model in which axion mass is enhanced by the condensate scale of the hidden QCD. We illustrate this with a specific model shown in  Fig.~\ref{fig:oneu1}(a). In this model, there is only one $U(1)_{PQ}$ chiral symmetry. The corresponding current and its divergence are
	\begin{equation}
	J^\mu_{PQ}=f_a \partial_{\mu}a+m\overline{Q}\gamma^\mu\gamma_5Q+m\overline{Q'}\gamma^\mu\gamma_5Q', \qquad
	\partial_{\mu}J^\mu_{PQ}=\frac{mg^{2}}{16 \pi^{2}} F_{a}^{\mu \nu} \tilde{F}_{a \mu \nu}+\frac{mg'^{2}}{16 \pi^{2}} G_{A}^{\mu \nu} \tilde{G}_{A \mu \nu}.
	\end{equation}
	At the low energy scale, the axion mass term stems from the instanton effects of QCD and the hidden QCD. Therefore, the mass terms are
	\begin{equation}
	-\mathcal{L}_{\text{mass}}=\cdots + \Lambda_c^4 \left(\frac{ma}{f_{a}}\right)^2+\Lambda_h^4 \left(\frac{ma}{f_{a}}\right)^{2},
	\end{equation}
	where $\Lambda_c$ and $\Lambda_h$ denote the confinement scale of QCD and the hidden QCD respectively. Due to the contribution from the high energy scale $\Lambda_h$, the axion mass could be much heavier than KSVZ or DFSZ axion. Different phenomena of the heavy axion may appear in collider experiments or rare decay of mesons~\cite{Gherghetta2016a}.
	\subsection{Two-axion Solution}
	In two-axion solutions with extra hidden QCD, the light mass eigenstate always exists and can be treated as an invisible axion, even if $\Lambda_h$ is much higher than QCD scale $\Lambda_c$. We illustrate this with the mass matrix of axions by a perturbative method. The Lagrangian of the mass part can be expressed as
	\begin{equation}\label{L_two_axions_mass}
	-\mathcal{L}_{\text{mass}}=\cdots + \Lambda_c^4 \left(\frac{m_1a_1}{f_{a_1}}+\frac{m_2a_2}{f_{a_2}}\right)^2+\Lambda_h^4 \left(\frac{n_1a_1}{f_{a_1}}+\frac{n_2a_2}{f_{a_2}}\right)^{2},
	\end{equation}
	where $f_{a_1}$ and $ f_{a_2}$ are scale factors. Anomaly  coefficients of axions are denoted as $m_1,\ m_2,\ n_1$ and $n_2$. The corresponding mass matrix can be expressed as
	\begin{equation}
	\Lambda_c^4\begin{pmatrix}
	m_1^2\frac{1}{f_{a_1}^2}	&m_1m_2\frac{1}{f_{a_1}f_{a_2}} \\
	m_1m_2\frac{1}{f_{a_1}f_{a_2}}	& m_1^2\frac{1}{f_{a_2}^2}
	\end{pmatrix}
	+
	\Lambda_h^4\begin{pmatrix}
	n_1^2\frac{1}{f_{a_1}^2}	&n_1n_2\frac{1}{f_{a_1}f_{a_2}} \\
	n_1n_2\frac{1}{f_{a_1}f_{a_2}}	& n_1^2\frac{1}{f_{a_2}^2}
	\end{pmatrix} .
	\end{equation}
	Since $\Lambda_h \gg \Lambda_c$, the first matrix with respect to $\Lambda_c $ can be treated as perturbation. Considering zero order approximation, the heavy degree of freedom is denoted as
	\begin{equation}\label{A}
	A \propto \frac{n_1a_1}{f_{a_1}}+\frac{n_2a_2}{f_{a_2}},
	\end{equation}
	while the massless one is in the orthogonal direction with $A$ in field space, which is
	\begin{equation}\label{a}
	a \propto \frac{n_2a_1}{f_{a_2}}-\frac{n_1a_2}{f_{a_1}}.
	\end{equation}
	Replacing $(a_1,\ a_2)$ with $(A,\ a)$, eq.~\eqref{L_two_axions_mass} becomes
	\begin{equation}\label{diagonalized}
	-\mathcal{L}_{\text{mass}}=\cdots + \Lambda_c^4 \left(\cos\alpha\frac{a}{F_a}+\sin\alpha\frac{A}{F_A}\right)^2+\Lambda_h^4 \left(\frac{A}{F_A}\right)^{2},
	\end{equation}
	where $F_A$ ($F_a$) is the decay constant of heavy axion ( light axion) after diagonalization and normalization. $\alpha$ is the mixing angle. In general, $F_a$ and $F_A$ are comparable. Therefore, the heavier one has a mass of $m_A^2 \sim \Lambda_h^4/F_A^2$ and the lighter one has a mass of $m_a \sim \Lambda_c^4/F_a^2$.
	The lighter axion is exactly the invisible QCD axion.

	In the following part, we discuss two specific models of different breaking mechanisms in mixed two-axion solutions. The model with two independently breaking $U(1)$ symmetries (model A) is shown in Fig.~\ref{fig: two U_1}(a)  and  the model with one breaking $U(1)$ symmetry induced by the other (model B) is shown in Fig.~\ref{fig: two U_1}(b). The main difference between these two models is whether the symmetry breaking of $U(1)_{PQ}$ is induced by the breaking of $U(1)_A$.
	Chiral symmetries are same in both cases. The corresponding currents and its divergence of these models are
	\begin{equation}\label{eq: anomalous current}
	\begin{aligned}
	J^\mu_{A}&=\dfrac{\sqrt{6}}{3}\overline{\psi}\gamma^\mu\gamma_5\psi, \qquad J^\mu_{PQ}=f_a \partial_{\mu}a+(m\overline{Q}\gamma^\mu\gamma_5Q)+n\overline{Q'}\gamma^\mu\gamma_5Q', \\
	\partial_{\mu}J^\mu_{A}&=\frac{\sqrt{6}g^{2}}{16 \pi^{2}} G_{a}^{\mu \nu} \tilde{G}_{a \mu \nu}+\frac{\sqrt{6}g'^{2}}{16 \pi^{2}} {G'}_{A}^{\mu \nu} \tilde{G'}_{A \mu \nu}, \\
	\partial_{\mu}J^\mu_{PQ}&=\frac{mg^{2}}{16 \pi^{2}} G_{a}^{\mu \nu} \tilde{G}_{a \mu \nu}+\frac{ng'^{2}}{16 \pi^{2}} {G'}_{A}^{\mu \nu} \tilde{G'}_{A \mu \nu},
	\end{aligned}
	\end{equation}
	where we set the anomaly coefficient $m=0$ in the second model for simplicity.
	\subsubsection{Model A}
	\label{indenpent}
	In model A, the spontaneous breaking of two chiral $U(1)$ symmetries comes from different mechanisms. The $U(1)_{PQ}$ symmetry is broken by the effective potential of the complex scalar $\phi$, while $U(1)_{A}$ is broken by the condensate of massless fermion $\psi$. Therefore, these two scales are independent. The only source of explicit breaking comes from the instanton effects contributing to masses of both axions. The mass-relevant Lagrangian is
	\begin{equation}\label{eq:L_mass}
	-\mathcal{L}_{\text{mass}}=\cdots + bf_\pi^4 \left(\frac{\sqrt{6}\eta^{\prime}_h}{f_{\pi_h}}+\frac{ma}{f_{a}}\right)^2+b'f_{\pi_h}^4 \left(\frac{\sqrt{6}\eta^{\prime}_h}{f_{\pi_h}}+\frac{na}{f_{a}}\right)^{2},
	\end{equation}
	where the first(second) term stems from QCD(hidden QCD) instanton effects with $f_{\pi_h} \gg f_\pi$.  After diagonalizing the mass matrix, we find the lighter mass eigenstate $a_{phy}$ to be
	\begin{equation}\label{axion mass 1}
	a_{\text{phy}}=\dfrac{nf_{\pi_h} \eta^{\prime}_h-\sqrt{6} f_{a} a}{\sqrt{{6f^2_a+n^2}f_{\pi_h}^2}},
	\end{equation}
    with its mass and decay constant are
	\begin{equation}\label{F_a}
	m_{a_\text{phy}}^2F_{a_\text{phy}}^2 =m_{\pi}^2 f_\pi^2\frac{m_u m_d}{(m_u+m_d)^2}, \qquad F_{a_\text{phy}}=\dfrac{\sqrt{{6f^2_a+n^2}f_{\pi_h}^2}}{\sqrt{6}|m-n|}.
	\end{equation}
	Furthermore, the axion-photon coupling $g_{{a_{phy}} \gamma \gamma}$ can be defined as
	\begin{equation}\label{eq:L_arr}
	\begin{aligned}
	\mathcal{L}_{{a_{phy}} \gamma \gamma}=\frac{1}{4}g_{{a_\text{phy}} \gamma \gamma}a F_{\mu \nu} \tilde{F}^{\mu \nu}=\frac{1}{4} (g_{{a_\text{phy}} \gamma \gamma}^\text{IR}+g_{{a_\text{phy}} \gamma \gamma}^\text{UV})a F_{\mu \nu} \tilde{F}^{\mu \nu} ,
	\end{aligned}
	\end{equation}
	where $g_{{a_\text{phy}}\gamma \gamma}^\text{IR}$ comes from the mixing between axion and $\pi^0$,  The latest result is accurate to NLO~\cite{Cortona2015}, which is
	\begin{equation}
 g_{{a_\text{phy}} \gamma \gamma}^\text{IR}=-1.92(4) \frac{\alpha_{\mathrm{em}}}{2 \pi F_{a_\text{phy}}}.
	\end{equation}
	The UV part is model-dependent. In this model, $g_{a \gamma \gamma}^\text{UV}$ has contribution from the mixing of $\eta'_h$ and $a$. the interactions by anomaly effect between  $\eta'_h$, $a$ and photon can be written as
	\begin{equation}\label{garr1}
	\mathcal{L}_{a\gamma \gamma}=\left(\sqrt 6q^2_\psi \frac{\eta^{\prime}_h}{f_{\pi_h}}+6(q_Q^2+q_{Q'}^2)\frac{a}{f_{a}} \right) \frac{\alpha_{\mathrm{em}}}{4\pi}F_{\mu \nu} \tilde{F}^{\mu \nu},
	\end{equation}
	where $q_\psi$, $q_Q$ and $q_{Q'}$ denote the electric charge of each particles. Therefore,
	\begin{equation}\label{coupling}
g_{{a_\text{phy}}\gamma\gamma}^\text{UV}=\dfrac{nq^2_\psi-6(q_Q^2+q_{Q'}^2)}{2|m-n|}\frac{\alpha_{\mathrm{em}}}{2 \pi F_{a_\text{phy}}}.
	\end{equation}
	When setting $n\gg|m-n|$, the UV contribution of the axion-photon coupling is enlarged. Meanwhile the physical axion decay constant meets the condition $F_{a_\text{phy}} \gg f_a $ as shown in eq.~\eqref{F_a}. This result is actually a special case of `alignment axion' models~\cite{Kim2005,Choi2014a}.

	\subsubsection{Model B}
	\label{induced}
	In model B, the spontaneous breaking of $U(1)_{PQ}$ is triggered by the condensate of the light fermion $Q'$ in the hidden sector. Thus two breaking scales are related. We assume the potential of the complex scalar $\phi$ containing only a quadratic term with $\mu^2 > 0$. The Yukawa term of $Q'$ will induce a tadpole term of $\phi$ after the chiral symmetry breaking in the hidden sector. The Lagrangian after the chiral symmery breaking is
	\begin{equation}
	\begin{aligned}\label{Lag}
	\mathcal{L}_\text{CSB}&=\partial_{\mu}\phi\partial^{\mu}\phi^*-\mu^2\phi^2
	+\dfrac{f_{\pi_h}^{2}}{4} \text{Tr} \partial_{\mu} \Sigma \partial^{\mu} \Sigma^{\dagger}+a' f_{\pi_h}^{3} \text{Tr}\left(H\Sigma^{\dagger}+H^{\dagger} \Sigma\right)+\cdots, \\
	\end{aligned}
	\end{equation}
	where
	\begin{equation}
	\Sigma=\exp\left[  \dfrac{ i\left( \pi_h^aT^a+\frac{1}{\sqrt{2}}I_{4\times 4} \eta'_h \right) }{f_\pi} \right] ,
	\qquad
	H=
	\begin{pmatrix}
	0 &  &  &  \\
	& 0 &  &  \\
	&  & 0 &  \\
	&  &  & y\phi
	\end{pmatrix} .
	\end{equation}
	Here $T^a$ are generators of the $SU(4)$ group meeting the trace condition $\text{Tr}(T^aT^b)=2\delta^{ab} $. The anomalous term is omitted in eq.~\eqref{Lag}. Matrix $H$, the Yukawa terms of $\psi$ and $Q'$, is treated as a spurion field transforming as an adjoint representation under the $SU(4)_V$. The third term in eq.~\eqref{Lag} is the linear term of $\phi$ which triggers $U(1)_{PQ}$ breaking. Specifically, the complex scalar $\phi$ is parameterized as $\phi=\rho e^{i\left( na/{f_a} \right)} $. The effective potential of $\rho$ and $a$ is
	\begin{equation}\label{potential phi}
	V(\rho,a)=\mu^{2} \rho^{2}-2a' y\rho f_{\pi_h}^{3} \cos \left(n\frac{a}{f_a}+\cdots \right),
	\end{equation}
	where $``\cdots"$ omits contributions from other composite scalars. Then we can estimate the VEV of $\rho$ to be $\left\langle \rho \right\rangle=f_a/\sqrt{2}  \sim 2y\dfrac{a'f_{\pi_h}^3}{\mu^2} $. Utilizing the  anomalous currents in eq.~\eqref{eq: anomalous current} and currents of $\pi_h$ and $\eta'_{h}$ in eq.~\eqref{current 12}, the mass terms of scalars in Lagrangian eq.~\eqref{Lag} becomes
	\begin{equation}\label{L_mass}
	-\mathcal{L}_{\text{mass}}=bf_\pi^4\left[ \sqrt{6}\left(\dfrac{1}{2}\dfrac{\pi_h}{f_{\pi_h}}
	+\dfrac{\sqrt{3}}{2}\dfrac{\eta^{\prime}_h}{f_{\pi_h}}\right)\right]^2+
	a'y\dfrac{f_af_{\pi_h}^{3}}{\sqrt{2}}\left(\dfrac{na}{f_a}-\dfrac{3}{\sqrt{6}}\frac{\pi_h}{f_{\pi_h}}+\dfrac{1}{\sqrt{2}} \frac{\eta_{h}^{\prime}}{f_{\pi_h}} \right)^2
	+
	b'f_{\pi_h}^4 \left(\frac{2\sqrt{2}\eta^{\prime}_h}{f_{\pi_h}}\right)^{2},
	\end{equation}
	where the second term stems from the Yukawa term of $Q'$.
	The $na/f_a$ in second term could move to the third term by a chiral rotation $\theta' \rightarrow \theta'+na/f_a $ and generates additional differential interaction of $a$.
	By means of the same steps in the former subsection, the lightest mass eigenstate is
	\begin{equation}\label{axion mass 2}
	a_{phy}=\dfrac{nf_{\pi_h} \pi_h+3/\sqrt{6} f_{a} a}{\sqrt{{3/2f^2_a+n^2}f_{\pi_h}^2}}, \qquad F_{a_\text{phy}}=\dfrac{2\sqrt{3/2f^2_a+n^2f_{\pi_h}^2}}{\sqrt{6}|n|}.
	\end{equation}
	The $g_{a\gamma\gamma}^{\text{UV}}$ coming from $\pi_h$ to $\gamma\gamma$ is
	\begin{equation}\label{garr2}
	\mathcal{L}_{\pi_h\gamma\gamma}=\dfrac{\sqrt{6}}{2}(q^2_\psi-q^2_{Q^{\prime}}) \frac{\pi_h}{f_{\pi_h}}  \frac{\alpha_{\mathrm{em}}}{4\pi}F^{\mu \nu} \tilde{F}_{\mu \nu}.
	\end{equation}
	Therefore, the UV part contribution of light axion is
	\begin{equation}\label{coupling2}
g_{{a_\text{phy}}\gamma\gamma}^\text{UV}=2(q^2_\psi-q^2_{Q^{\prime}})\frac{\alpha_{\mathrm{em}}}{2 \pi F_{a_\text{phy}}}.
	\end{equation}

	In this model, the axion-photon coupling $g_{{a_\text{phy}} \gamma \gamma}$ is independent of the anomaly coefficient $n$, while $F_{a_\text{phy}}$ is inversely proportional to $|n|$. Furthermore, there is another colorless scalar $\pi_h$ with $m^2_{\pi_h} \sim f_af_{\pi_h}$.

	\section{Conclusion and Discussion}
	\label{section5}
	We investigate all possible solutions to the strong $CP$ problem with anomalous $U(1)$ global symmetries when the hidden QCD is introduced, which include both the Peccei-Quinn mechanism and dynamical solutions. We find these two mechanisms can exist together, and a new class of solutions for the strong CP problem are proposed. Thus we classify two $U(1)$ models in an unified way, and discuss both one-$U(1)$ and two-$U(1)$ solutions in this work.

	A ``moose-like'' diagram method is used to illustrate how the strong CP problem is solved via cancellation between CP phases, which could be read  directly. When the gauge interaction and new particles are given, a diagram can be determined with a breaking mode of $U(1)$ symmetries. From these diagrams, cancellation between $CP$ phases can be easily read without detailed analysis of the model. The information is important for solving the strong $CP$ problem. This method could be used to find new types of solutions and easily extended to cases of multi-$U(1)$ and multi-gauge interaction.

	In one-axion solutions, additional $Z_2$ symmetry keeps $\theta=\theta'$ to solve the strong $CP$ problem.

	In two-axion solutions, two chiral $U(1)$s (either $U(1)_{PQ}$ or $U(1)_A$) have to be introduced. As a consequence, there always exists an invisible light axion and a heavy axion.

	Using the diagram method, we find two new mixed solutions. In Model A, the spontaneous breaking of $U(1)_{PQ}$ symmetry is triggered by the scalar potential, independent of the confinement of the hidden QCD, and the axion-photon couplings $g_{{a_\text{phy}}\gamma\gamma}^{UV}$ and the light axion decay constant $F_{a_\text{phy}}$ are both inversely proportional to the difference of anomalous charges $|m-n|$. In Model B, the  $U(1)_{PQ}$ is induced by the confinement of the hidden QCD, and  $g_{{a_\text{phy}}\gamma\gamma}^{UV}$ is independent of the anomaly coefficient $n$ of $Q'$ while $F_{a_\text{phy}}$ is inversely proportional to $|n|$.	At low energy scales, it is hard to distinguish this light axion from the QCD axion, even extending to multi-$U(1)$ symmetries and multi-gauge interactions. However, this kind of solutions could be probed by detecting particles at the hidden QCD scale. Those heavy particles may leave some hints in cosmology and at the colliders.

	\acknowledgments
We would like to thank Qing-Hong Cao and Shou-hua Zhu for valuable discussions and their support.
D.H. is supported in part by the National Science Foundation of China under Grants No. 11635001, 11875072.
H.R.J is supported in part by the National Science Foundation of China under Grant Nos. 11725520, 11675002, 11635001.
H.L.L and J.H.Y. are supported by the National Science Foundation of China (NSFC) under Grants No. 11875003.
M.L.X. is supported by the National Natural Science Foundation of China (NSFC) under grant No.2019M650856 and the 2019 International Postdoctoral Exchange Fellowship Program.
J.H.Y. is also supported by the National Science Foundation of China (NSFC) under Grants No. 11947302.


\subsection{Experimental system}

\begin{figure*}[!htb]
\centering
\includegraphics[width=140mm]{FigureS1.pdf}
\caption{Spectroscopy of spin exchange in double wells.
(a) Illustration of the spin-dependent potentials. For two spin configurations in a double well, $\ket{\downarrow,\uparrow}$ and $\ket{\uparrow,\downarrow}$, their energies can be well controlled by tuning the offset value $\delta_s$. For instance, at negative side, the state $\ket{\uparrow,\downarrow}$ has lower energy than $\ket{\downarrow,\uparrow}$.
(b) Spin superexchange under different staggered magnetic potentials.
The states are initialized to $\ket{\uparrow,\downarrow}$ in double wells, corresponding to the staggered magnetization of $M_z=-0.5$. Then, we quench the lattice depth to allow spin exchange dynamics at the coupling strength of $J=26.8(5)$ Hz, holding for 20 ms. The spectroscopy shows a symmetric spin-dependent effect. The solid curve is a guide for the eyes. Error bars are standard deviations throughout this supplementary material.}
\label{figS1}
\end{figure*}

Our experiment begins with a nearly pure Bose-Einstein condensate of $^{87}$Rb staying on the hyperfine state of  $5S_{1/2}\left|F=1,m_F=-1\right>$.
Then, we compress the atoms along the $z-$axis by increasing the confinement, followed by which atoms are loaded  into a single node of a 3-$\mu m$ spacing blue-detuned lattice.
Thus, we create a quasi-2D quantum gas.
The confinement in the $x-y$ plane is provided by a red-detuned optical dipole trap propagating along the $z$ direction.
To prepare a low-entropy Mott insulator, we implement a newly-developed cooling technique in optical lattices.
The cooling concept is that the entropy of these environment-isolated atoms is transferred into contacting superfluid reservoirs, which are removed from the system later.
In our region of interest (36 copies of 70-site chains), the probability of unity-filling for the cooled Mott insulator is 0.992(1).

To initialize the N\'eel order along the $x$-direction, we employ a site-selective spin addressing technique to overcome the optical diffraction limit.
The atoms in the Mott insulator state are transferred into an spin-dependent optical superlattice \cite{Yang:2017a}, which breaks the odd-even symmetry of the lattice sites and thereby separates the atoms into two subsystems.
For odd and even subsystems, the lattice potentials are controlled to split their transition frequencies between $\ket{\downarrow}$ and $\ket{\uparrow}$ hyperfine levels.
Applying an edited microwave pulse, we flip the atoms of odd sites into $\ket{\uparrow}$ with an efficiency of 0.995(3) through the rapid adiabatic passage.
For each site, the preparation fidelity of the N\'eel state can be estimated as 0.987(3).
Therefore, in the 70-site chain, a defect-free N\'eel order is achieved with a probability of $\sim$ 40\%.

\subsection{Staggered potentials}

The superlattice naturally forms staggered potentials for our ground-band ultracold atoms.
Here, we make use of two distinct effects of the superlattice, which are the spin-dependent and the  spin-independent potentials.
Such optical potential can be expressed as,

\begin{equation}
V(x) = V_{s,\sigma}\cos^2(k_s x+\theta_\sigma) + V^{\text{offset}}_{s,\sigma} - V_l\cos^2(k_s x/2 + \theta),
\end{equation}
where $k_s=2\pi/\lambda_s$ is the wave vector of the short-lattice laser.
Here, $\theta_\sigma$ is the spin-dependent phase shift of the short-lattice, where $\sigma$ denotes the two spin states.
The original trap depth of short-lattice $V_s$ is also altered to two spin-dependent terms $V_{s,\sigma}$, and $V^{\text{offset}}_{s,\sigma}$.
At $\theta_{\sigma} = 0$, the superlattice is simplified to the spin-independent form, $V^{\text{offset}}_{s,\sigma} =0$.
By tuning the applying voltage of a electro-optical modulator, we can control the spin-dependent energy with a bandwidth larger than 20 kHz \cite{Yang:2017a}.

The spin-dependent energy shift can be precisely controlled in our system.
To calibrate the zero-energy $\delta_s =0$, we perform a spectroscopy measurement of spin-exchange dynamics in isolated double-wells.
We prepare the atoms in $\ket{\uparrow,\downarrow}$ state and then quench the superlattice to $V_s=16.8(1) E_r$, $V_l=10.3(1)$ Er and $\theta=0$.
Here $E_r = h^2/(2 m \lambda_s^2)$ is the recoil energy with $h$ the Planck constant and $m$ the atomic mass.
In this balanced structure, the resonating superexchange frequency within double wells is $26.8(5)$ Hz.
The inter-well coupling is forbidden by the long-lattice barrier.
The degeneracy of the states $\ket{\uparrow,\downarrow}$ and $\ket{\downarrow,\uparrow}$ is lifted by a nonzero spin-dependent effect.
As we scan the $\delta_s$ and monitor the staggered magnetization $M_z$, a clear spectroscopy of the spin exchange dynamics is observed, as shown in Fig. \ref{figS1}.
The position of the central peak can be determined with a precision of $\sim$0.4 Hz.
In addition, the energy offsets of the side peaks agree well with our calculations based on the lattice structure.

\subsection{Establishing the isotropic Heisenberg Hamiltonian}

We design a special 1D superlattice along the $x$-direction for realizing the isotropic antiferromagnetic Heisenberg model.
When $V_l \ll V_s$, the inter-well coupling becomes comparably strong as the intra-well coupling.
In this 1D chain, there are two types of links connecting the neighboring sites, which can be denoted as odd-even (type I) and even-odd (type II) links.
If $\theta = \pi/4$, these two types of links are identical to each other.
While the spin-dependent effect will disappear in this configuration.
To enable the spin-dependent term, we set the superlattice phase to $\theta = \pi/3$ and keep the long-lattice at $V_l = 0.44(1) E_r$.
The tunneling strength on the type II links $t_2$ is roughly 10 percent larger than $t_1$ on the type I links \cite{Walters:2013s}.

In the second-order perturbation theory, the spin exchange coupling along these two types of links are,

\begin{equation}
J_{1}=-\frac{4t_{1}^2U}{U^2-(\delta_0-\Delta)^2}, \ \
J_{2}=-\frac{4t_{2}^2U}{U^2-(\delta_0+\Delta)^2}.
\end{equation}
Here, $\delta_0$ and $\Delta$ represent energy offsets between neighboring sites due to the spin-independent term and a linear potential, respectively.
$t_1$ and $t_2$ are tunneling strengths along the two types of nearest-neighbor links.
Since the lattice depths are almost the same for the two spin states, the tunneling strengths are considered as spin-independent.
$U$ corresponds to the onsite-interaction for all the spin state.
Here we neglect the slight spin-dependence of scattering lengths for $^{87}$Rb and adopt $U_{\uparrow\downarrow}=U_{\uparrow\uparrow}=U_{\downarrow\downarrow}=U$.
When $U < \delta_0 - \Delta \ \text{and} \ \delta_0 + \Delta$, we reverse the sign of $J_1$ and $J_2$ to positive, establishing an antiferromagnetic spin-exchange coupling \cite{Trotzky:2008}.

The difference between the spin exchange couplings $J_1$ and $J_2$ is compensated by the state-independent linear potential.
Since our $x$-axis lies at an angle of 4 degree related to the horizontal plane, the gravity leads to this constant 57-Hz/site gradient along the 1D chain.
Consequently, the discrepancy between these two types of exchange coupling becomes roughly 2 percent [see Figure \ref{figS2}(c)], one order of magnitude lower than the case without this compensation ($J_1$ and $J_2$ would have 20 percent difference in the absence of the linear gradient).

In our experimental conditions, the Hubbard parameters satisfy such relations $t_1, t_2 \ll \delta_0, U$, as shown in Fig. \ref{figS2}.
This Bose-Hubbard Hamiltonian can be mapped into a spin model as,

\begin{equation}
\hat{H}=\sum_{i}(J_{1}\hat{S}_{2i-1}\cdot\hat{S}_{2i} + J_{2}\hat{S}_{2i}\cdot\hat{S}_{2i+1}) + \delta_s\sum_{j}(S_{2j}^z-S_{2j+1}^z).
\end{equation}
Neglecting the slight difference between $J_1$ and $J_2$, we can adopt an isotropic Heisenberg model with the spin-exchange strength $J$, which is the geometric average of $J_1$ and $J_2$, as $J=\sqrt{J_{1}J_{2}}$.

We can understand the model in two limits.
In the limit of $\delta_s/J \rightarrow \pm \infty$, the two type N\'eel orders are the ground states, which are product states.
We should stress that even though atoms tend to occupy the low-energy site of the staggered potential in the ground state of the one-dimensional Hubbard model, the atoms cannot overcome the staggered potential via tunneling in our state-preparation time and the system can be faithfully described by the spin model.
In another limit, $\delta_s/J =0$, unlike the Ising-type spin model, this Heisenberg model owns a single energetic ground state $\ket{\psi_{\text{AF}}}$.
The SU(2) symmetry of this isotropic Hamiltonian leads to the same rotational symmetry of the ground state.

\begin{figure*}
\centering
\includegraphics[width=150mm]{FigureS2.pdf}
\caption{ From Bose-Hubbard model to Heisenberg spin model.
(a) The blue and black curves represent the tunneling strengths and spin-exchange coupling, respectively.
The lighter colors correspond to the type I links, representing $t_1$ and $J_1$.
While the darker colors correspond to the type II links , representing $t_2$ and $J_2$.
The discrepancy between the exchange coupling strengths has been mostly compensated by the gravity potential.
(b) The blue and black curves shows the on-site interaction $U$ and the spin-independent staggered potential $\delta_0$ under different short lattice depths.
(c) The ramping curves of the spin-exchange strength $J$ and the staggered magnetic potential $\delta_s$ in the adiabatic passage.}
\label{figS2}
\end{figure*}

\subsection{Adiabatic passage}

Starting from the initial N\'eel-ordered state $\left|\psi_0\right>=\left|\uparrow\downarrow\uparrow\downarrow\cdots\right>$, we experimentally optimize the ramping curves to minimize the non-adiabaticity to approach the Heisenberg antiferromagnet $\ket{\psi_{\text{AF}}}$.
The $y-$lattice with 42.0(3) $E_r$ isolates the system into 36 copies of 70-sites 1D chains.
The superlattice phases are switched to $\theta = \pi/3$ and $\theta_{\ket{\uparrow}} = -2 \theta_{\ket{\downarrow}} = 0.005 \pi$ before ramping down the short-lattice potential.

Then, we slowly ramp the lattice potentials in terms of the Hubbard parameters.
At 0 ms, the staggered magnetic gradient is $\delta_s/h=-17.6(3)$ Hz, and the coupling strength is $J/h= 0.24(1)$ Hz.
During the 60 ms ramping, we separate the curve into two stages with 30 ms interval, as shown in Fig. \ref{figS2}(c).
The exchange coupling is turned on gradually, and simultaneously the staggered gradient is also increased.
At 60 ms, $\delta_s$ is slightly larger than 0, but the zero value of the staggered magnetization $M_z = 0.00(2)$ indicates the achievement of a Heisenberg antiferromagnet.
In another 60 ms, we verify the coherence of the system by slowly ramping back the lattice potentials and increasing  $\delta_s$ oppositely.
To detect the property of the state at each time, we freeze the state by quenching the lattice barrier to disable the atom dynamics.

\subsection{Nearest-neighbor spin correlation}

The $in\ situ$ atomic densities are detected via an absorption imaging with the optical resolution of $\sim$1 $\mu$m.
Even though it is one kind of destructive imaging techniques, we can take two successive images for recording both the spin states, the first image for $\ket{\uparrow}$ and the second image $\ket{\downarrow}$.
We choose a region of interest (ROI) containing 36 copies of 70-site 1D chains in the center of the atom cloud for the statistical analysis.

The nearest-neighbor correlation function is defined as
\begin{equation}
C(d=1)= \frac{4}{N}\sum_{j}\left(\left<\hat{S}_j^z\hat{S}_{j+1}^z\right>-\left<\hat{S}_j^z\right>\left<\hat{S}_{j+1}^z\right>\right),
\end{equation}
where $N$ represents the number of sites.
In our spin-balanced system, the second term in the sum is equal to the square of the staggered magnetization, as $-1/N\sum_{j}\langle\hat{S}_j^z\rangle\langle\hat{S}_{j+1}^
z\rangle= M_z^2$.
For obtaining the first term in the sum, we perform correlated measurements on these two type of links respectively, as

\begin{equation}
\begin{split}
&\text{Type I :} \ \ \\
&-\frac{4}{N}\sum_{j}\langle\hat{S}_{2j-1}^z\hat{S}_{2j}^z\rangle =\frac{1}{2} \left(P_{\ket{\uparrow,\downarrow}}^{I} + P_{\ket{\downarrow,\uparrow}}^{I} -P_{\ket{\uparrow,\uparrow}}^{I}- P_{\ket{\downarrow,\downarrow}}^{I}\right),\\
\\
&\text{Type II :}  \ \ \\
&-\frac{4}{N}\sum_{j}\langle\hat{S}_{2j}^z\hat{S}_{2j+1}^z\rangle =\frac{1}{2} \left(P_{\ket{\uparrow,\downarrow}}^{II} + P_{\ket{\downarrow,\uparrow}}^{II} -P_{\ket{\uparrow,\uparrow}}^{II}- P_{\ket{\downarrow,\downarrow}}^{II}\right).
\end{split}
\end{equation}
Here, $P^{I(II)}$ represent the probabilities of spin states on the corresponding links.
Making use of the superlattice, we can combine the atoms on neighboring sites into a single lattice well.
Depending on the atom number on the combined site, we can derive the desired probabilities.
These two types of links are measured at two kinds of superlattice configurations, one has $\theta=0$ and the other has $\theta=\pi/2$.

\begin{figure*}
\centering
\includegraphics[width=140mm]{FigureS3.pdf}
\caption{Detection of the nearest-neighbor spin correlation.
(a) The detection scheme for the spin components on two types of nearest-neighbor links. In step (i), all $\ket{\uparrow}$ atoms are removed by a resonating imaging light pulse.
Then in (ii), the atoms on type I or type II links are respectively loaded into superlattice double-well units with $\theta=0$  and $\theta=\pi/2$.
In step (iii), the atoms within each double-well are combined and doublons in the long-lattice site get lost after the photoassociation collision.
Finally, in step (iv), we count the remaining atoms and deduce the probabilities of the state $\ket{\downarrow,\downarrow}$.
(b) The upper graph shows the state probabilities on type I links corresponding to the tunneling strength $t_1$, and for type II links the results are depicted in the lower graph.
In these two graphs, the red circles, the blue squares, the green triangles and the yellow inverted triangles represents the probabilities of $\ket{\uparrow,\downarrow},\ket{\downarrow,\uparrow},\ket{\uparrow,\uparrow}$ and $\ket{\downarrow,\downarrow}$, respectively.
The solid curves are numerical caculations of a 70-site spin chain.}
\label{figS3}
\end{figure*}

Not limited by our optical resolution, we probe the parity of atom number on lattice sites with a photoassociation (PA) laser.
The laser is red-detuned to the D2 line of $^{87}$Rb by 13.6 $\rm{cm}^{-1}$.
Assisted by this laser, atom pairs in hyperfine state $\ket{\downarrow}$ can be excited to the $\nu=17$ vibrational state of the $0_g^-$ long-range molecular channel \cite{Fioretti:2001s}.
This unstable molecule would decay back to free atoms and gains kinetic energy to escape from the trap.
Therefore, the ratio of atom loss is equivalent to the probability of the state $\ket{\downarrow,\downarrow}$ before the combining operation.

The detection scheme is illustrated in Fig. \ref{figS3}(a), where two sets of measurements are carried out for these two types of links.
There are several steps in the detection procedure.
The main concern is removing the atoms on $\ket{\uparrow}$ state before applying the atom combination.
The PA light has an intensity of 0.56 $\rm{W/cm}^2$ and lasts for 20 ms.
In this way, we obtain the probabilities of $P_{\ket{\downarrow,\downarrow}}^{I}$ and $P_{\ket{\downarrow,\downarrow}}^{II}$, as shown in Fig. \ref{figS3}(b).

For the other three spin-states $\ket{\uparrow,\uparrow},\ \ket{\uparrow,\downarrow},\ \ket{\downarrow,\uparrow}$, we employ the site-selective spin addressing technique to transfer them into the  $\ket{\downarrow,\downarrow}$ state before implementing the detection procedures \cite{Yang:2017a}.
For instance, we flip all of the atoms on odd sites for probing the probabilities $P_{\ket{\uparrow,\downarrow}}^{I}$ and $P_{\ket{\downarrow,\uparrow}}^{II}$.
The probabilities of these four states allow us to calculate the spin correlations $C(d=1)$.
However, all of the measurement errors would propagate into the final spin correlations.
We use the basic normalization condition $P_{\ket{\uparrow}}+P_{\ket{\downarrow}} =1$ to reduce the errors by only taking the residual atoms in the last step of Fig. \ref{figS3}(a) into calculations.

In Fig. \ref{figS3}(b), the probabilities of the four spin state for type I and type II nearest-neighbor links are depicted.
The deviation from the numerical calculations can be caused by the motional excitations in Hubbard model and the lattice heating during the adiabatic passage.
Besides, after steps of superlattice manipulation in the detection method, the accumulative inefficiency of the detection  also contributes to the reduction of the measured spin-state probabilities,

\subsection{Motional excitations}

\begin{figure*}[!htb]
\centering
\includegraphics[width=70mm]{FigureS4.pdf}
\caption{Motional excitations during the adiabatic passage. The blue (red) circles represent the ratios of the doublon excitations with identical (opposite) spin components. The solid curves are linear fits to the experimental data.}
\label{figS4}
\end{figure*}

Besides spin excitations in the state evolution, we also measure the ratios of motional excitations in these lattice sites, which are absent in the spin model.
In the Hubbard model, one particle-hole excitation emerges when an atom tunnels to its neighboring site.

The particle-hole excitations with doublons $\ket{\downarrow\downarrow}$ are detected by the same parity-projection technique making use of the PA laser.
After a global spin-flip, we can also acquire the doublons on the state $\ket{\uparrow\uparrow}$.
The particle excitations on the state $\ket{\uparrow\uparrow}$ and $\ket{\downarrow\downarrow}$ generated in the adiabatic passage are rare, since the atoms need to tunnel to its next-nearest neighboring sites.
Fig. \ref{figS4} shows that the total ratio of these excitations keeps constant at 0.01(1) in the state evolution.

However, the more probable particle-hole excitations result in the $\ket{\uparrow\downarrow}$ state, owing to the tunnelling of atoms to their neighboring sites.
From the normalization condition, the ratio of these excitations can be inferred as $1-\sum_{\mu,\nu}P_{\ket{\mu,\nu}}$, with $\mu,\nu=\ \uparrow,\downarrow$.
As depicted in Fig. \ref{figS4}, the ratio of these doublons excitations grows gradually during the 120 ms state evolution, which is 0.04(3) for our Heisenberg antiferromagnet.
The propagation of excitations is strongly suppressed by the spin-independent staggered offset $\delta_0$.
Besides, the linear potential also suppress the long-range motional transport \cite{Yang:2020s}.

\subsection{Next-nearest-neighboring spin correlation}

\begin{figure*}[!htb]
\centering
\includegraphics[width=150mm]{FigureS5.pdf}
\caption{Detection of next-nearest-neighbor spin correlation.
(a) The detection scheme for resolving the probabilities of state $\ket{\downarrow;\downarrow}$ on two types of links.
For probing the correlations on odd (even) sites, in step (i), we remove all the $\ket{\uparrow}$ atoms as well as the atoms residing on even (odd) sites.
In step (ii), the remaining atoms are transferred into the long-lattice potential.
Then in step (iii), the atoms are split into adjacent short-lattice sites with equal possibilities.
Finally, in step (iv), we combine the atoms and count the loss rate of doublons by removing them in the long-lattice potential.
(b) Probabilities of the states $\ket{\uparrow;\uparrow}$ and $\ket{\downarrow;\downarrow}$.
The upper and lower graph show the evolution of spin components for neighboring odd and even sites, respectively.
The blue circles represent the probabilities of state $\ket{\uparrow;\uparrow}$, and the red circles correspond to the state $\ket{\downarrow;\downarrow}$.
The solid curves are the numerical calculations of a 70-site spin chain.}
\label{figS5}
\end{figure*}

According to the Mermin-Wagner-Hohenberg theorem \cite{Mermin:1966s}, the one-dimensional isotropic Heisenberg Hamiltonian cannot support a long-range spin order \cite{Yang:2017s}.
As the distance grows, the magnetic correlation decays rapidly in the 1D Heisenberg antiferromagnet.
So the occurrence of nonzero next-nearest-neighbor spin correlation would be a distinctive signature of approaching a low-entropy antiferromagnet.

However, without the single-site imaging, the superlattice structure normally only provides the access to the nearest-neighbor spin correlation because of its $Z_2$ symmetry \cite{Trotzky:2010s,Dai:2016s}.
Here, owing to our precise control of the superlattice phase, we can extend the capability of the correlated spin measurement to next-nearest-neighbor links.

Figure \ref{figS5}(a) shows the procedure to detect the probability of spin state $\ket{\downarrow;\downarrow}$ along next-nearest-neighboring links, implemented on odd sites and even sites separately.
After a selective atom removal, the superlattice serves as an array of atomic splitters to redistribute the remaining atoms into their adjacent sites with equal probabilities.
Subsequently, changing the superlattice phase $\theta$ by $\pi/2$, we can transfer the $\ket{\downarrow}$ atoms originally lying two sites apart to the same long-lattice unit in a probabilistic way.
The atom loss in a PA collision can be used to deduce the probability of spin state $\ket{\downarrow;\downarrow}$.
Similarly, the $\ket{\uparrow; \uparrow}$ state on next-nearest-neighbor sites can be detected with an additional global spin-flip pulse before the step (i).

We define a next-nearest-neighbor correlator as,

\begin{equation}
C'(d=2)=P_{\ket{\uparrow;\uparrow}}+P_{\ket{\downarrow;\downarrow}},
\end{equation}
Here $P_{\ket{\uparrow;\uparrow}}$ and $P_{\ket{\downarrow;\downarrow}}$ represent the probabilities of spin state $\left|\uparrow;\uparrow\right>$ and $\left|\downarrow;\downarrow\right>$ along next-nearest-neighbor links, which equal to four times the corresponding atom density loss in the photoassociation collision.
For obtaining accurate probabilities, we subtract the contribution of the motional excitations in the PA collision.

Figure \ref{figS5}(b) shows the evolution of the probabilities $P_{\ket{\uparrow;\uparrow}}$ and $P_{\ket{\downarrow;\downarrow}}$ for odd sites and even sites separately.
Since only a quarter of atoms contribute to the signal, the probabilities measured here have significantly larger errors compared to their counterparts on nearest-neighbor links.

\subsection{Numerical simulation}

\begin{figure*}[!htb]
\centering
\includegraphics[width=140mm]{FigureS6.pdf}
\caption{Numerical calculations of the Heisenberg spin chains. (a) spin relaxation dynamics after quenching the N\'eel order to the Heisenberg spin model.
(b)(c) and (d) are respectively the numerical results of the evolution of $M_z$, $C(d=1)$ and $C'(d=2)$ in the adiabatic passage.
Different blue colors represent the chain length of 10, 12, 16, 20, 30, 40, 50 and 70, respectively.
The inset graphs show the discrepancies between numerics for different system lengths and the curve with the largest length (N = 70). The finite-size convergence is well below our detection errors.}
\label{figS6}
\end{figure*}

\begin{figure*}[!htb]
\centering
\includegraphics[width=154mm]{FigureS7.pdf}
\caption{A 70-site Heisenberg spin chain is simulated with the finite-temperature Monte Carlo method.
The nearest-neighbor spin correlation $C(d=1)$ and the entropy per particle versus temperatures are depicted in (a) and (b), respectively.
The red circles represent the case of our prepared antiferromagnet.
(c) Probabilities of the spin states in the nearest-neighbor and next-nearest-neighbor links.
The circles represent the measured probabilities of spin components, while the blue bars are the theoretical predictions of a thermalized Heisenberg antiferromagnet at $C(d=1)=-0.44$.}
\label{figS7}
\end{figure*}

The numerical data for the post-quench dynamics and the adiabatic passage of the Heisenberg models are obtained using the the $t$-DMRG algorithm in the framework of matrix product states.
In the numerical simulations, we neglect the slight difference between the two types of nearest-neighbor links and adopt the geometric average $J=\sqrt{J_{1}J_{2}}$ as the spin-exchange strength throughout the chain.
For the implementation, we use the OpenMPS library \cite{Jaschke:2018s} and fix the maximal bond dimension to $D=1500$.
Convergence are achieved at a time step of $3\times10^{-4}$ s and $10^{-2}\ h/J$ for the adiabatic passage and the quench dynamics, respectively, with truncation threshold of $10^{-6}$ per time step.
Finite-size effects are also investigated for several chain lengths.
As show in Fig. \ref{figS6}, the results for system size from $L=40$ to $70$ reveal similar post-quench dynamics, and converge with our experimental results in statistical errors.

For estimating the temperature and the entropy of our state, we numerical calculate the Heisenberg model using a finite-temperature quantum Monte Carlo method (QMC) \cite{Bauer:2011s}.
All the simulations are carried out in the grand-canonical ensemble on open chains of size with $N=70$, which is checked to be large enough to eliminate the finite-size effect.
After $10^5$ times of thermalization, we continue to perform the QMC sweeps for at least $5\times10^6$ times to get the results under each setting.
We choose the temperature $k_B T/J$ from 0.2 to 0.8 with a 0.03 interval.
From the nearest-neighbor correlation of interest [$C(d=1)=-0.44(3)$], we derive the temperature of our system to be 0.53(7)$J$, see Fig. \ref{figS7}(a).
The entropy per particle is obtained using another QMC method, which is called quantum Wang-Landau algorithm \cite{Wang:2001s}.
The cutoff of series expansion is set to be 1000.
As shown in Fig. \ref{figS7}(b), the entropy per particle is 0.42(5)$k_B$, which is deduced from the temperature of the spin chains.

Based on the comparison between the measured probabilities of the spin state and the prediction of the thermalized state, our system is inferred to be a thermalized state, see Fig. \ref{figS7}(c).