#  Axion mass in antiferromagnetic insulators

Koji Ishiwata  Institute for Theoretical Physics, Kanazawa University,
Kanazawa 920-1192, Japan

###### Abstract

We calculate axion potential in antiferromagnetic insulators from path
integral. It is derived from the effective potential for the order parameter
of the antiferromagnetic phase in insulators. Static and dynamical axions are
defined consistently from the potential. Consequently
antiferromagnetic/paramagnetic and topological/normal orders are classified.
The dynamical axion is predicted in all phases and its mass turns out to be
various values up to eV scale.

††preprint: KANAZAWA-21-02

##  I introduction

Axion is a hypothetical elementary particle that is considered to solve the
strong CP problem in QCD [1]. It is originally massless due to shift symmetry
of the Lagrangian, and it becomes massive after the QCD transition via a
nonperturbative effect. Other than the QCD axion, axion particles are known to
emerge from superstring theory. Motivated partly by that, axion-like particles
(ALPs) are studied in phenomenological, astrophysical, and cosmological
contexts. The mass of ALPs is usually considered as a free parameter, leading
to more rich phenomenology compared to the QCD axion (see, e.g., Ref. [2] for
review).

Besides in the field of particle physics, axion is discussed in condensed
matter physics. For example, the magnetoelectric response due to ${\bf\it
E}\cdot{\bf\it B}$ coupling (${\bf\it E}$ and ${\bf\it B}$ are electric and
magnetic fields, respectively) is already known, and the ${\bf\it
E}\cdot{\bf\it B}$ coupling is the interaction term of axion with the gauge
field. Therefore, such magnetoelectric response may be able to understood in
terms of “axion.” There have been several studies in that context so far.
Different from axion in elementary particle physics, there are two types of
“axions” in condensed matter physics: static axion and dynamical axion.
(Hereafter we use “axion” as an axion in condensed matter physics unless
otherwise mentioned.) In terms of Lagrangian, they are written as

$\displaystyle{\cal L}_{\theta}=\frac{\alpha}{\pi}(\theta_{0}+a){\bf\it
E}\dotproduct{\bf\it B}\,,$ (1)

where $\alpha=e^{2}/(4\pi)$ is the fine-structure constant and we have denoted
$\theta_{0}$ as the static axion (constant) and $a$ as the dynamical axion.
The former is pointed out by Refs. [3, 4, 5, 6, 7, 8] as the ${\bf\it
E}\cdot{\bf\it B}$ coupling. Nonzero $\theta_{0}$ breaks both time-reversal
${\cal T}$ and parity ${\cal P}$. However, when $\theta_{0}=\pi$ those
symmetries are intact. This nontrivial state is known as an topological state.
Ref. [9], on the other hand, studied the latter to predict that the dynamical
axion causes the total reflection of light irradiated to the topological
insulators under an antiferromagnetic (AFM) order. They solved equations of
motion of photon and the dynamical axion that interact with each other via the
${\bf\it E}\dotproduct{\bf\it B}$ interaction term and found an energy gap
emerged in the topological insulators. In the similar context, Ref. [10] shows
that an electric field is expected to be converted to a magnetic field due to
the dynamical axion when the electric field exceeds certain critical value.
(See also Refs. [11, 12, 13, 14, 15].) Such recent developments in both
particle physics and condensed matter physics have motivated particle axion
and ALPs search by using axions in insulators [16, 17, 18, 19].

In this paper we calculate the effective potential for the order parameter of
the AFM phase from the partition function of the system defined by the path
integral. Using the effective potential, axion potential is derived. It is
found that the effective potential and axion potential are intuitive and
useful in order to distinguish the AFM/paramagnetic (PM) order and the
topological/normal phase depending on the model parameters. We clarify the
static and dynamical axions in the AFM insulators and how they are related to
the AFM order and the topological state. Consequently, the mass of dynamical
axion is determined. While it is estimated to be $\order{\rm meV}$ in Ref.
[9], the axion mass turns out to be larger up to $\order{{\rm eV}}$ or a much
suppressed value depending on the property of insulators.

This paper is organized as follows. In Sec. II, the effective potential for
order parameter of the AFM phase is derived. Using the result, the AFM/PM
phases are determined. The effective potential is interpreted in terms of
axion field and axion potential is derived in Sec. III. Sec. IV contains
conclusions and discussion.

##  II potential of antiferromagnetic order parameter

The topological state in insulators is consider to be realized by a (quasi)
gapless state via the spin-orbit coupling. For modeling such a state, a
minimum setup is to consider 4-by-4 matrix, i.e., two energy states (or orbits
or sublattices) with up and down spins. Such a model is often described by the
gamma matrices (even in lower or higher dimension space [3]). In the momentum
space, it is given by

$\displaystyle H({\bf\it k})=\epsilon_{0}({\bf\it k}){\bf 1}_{4\times
4}+\sum_{a=1}^{5}d^{a}({\bf\it k})\Gamma^{a}\,,$ (2)

where ${\bf\it k}=(k_{x},\,k_{y},\,k_{z})$ is the wavenumber and $\Gamma^{a}$
are the gamma matrices satisfying

$\displaystyle\\{\Gamma^{a},\Gamma^{b}\\}=$ $\displaystyle\,2\delta^{ab}{\bf
1}_{4\times 4}\,,$ (3) $\displaystyle{\rm
tr}(\Gamma^{a}\Gamma^{b})=4\delta^{ab}\,,$ $\displaystyle~{}~{}~{}{\rm
tr}(\Gamma^{a})=0\,.$ (4)

The representation of the gamma matrices depends on the basis. In Refs. [9,
20], for example, the basis of the antibonding and bonding are chosen, i.e.,
($\ket{P1_{z}^{+},\sigma}$, $\ket{P2_{z}^{-},\sigma}$)
($\sigma=\uparrow,\,\downarrow$):

$\displaystyle\Gamma^{a}=(\sigma^{x}\otimes s^{x},\sigma^{x}\otimes
s^{y},\sigma^{y}\otimes{\bf 1},\sigma^{z}\otimes{\bf 1},\sigma^{x}\otimes
s^{z})\,,$ (5)

where $\sigma^{j}$ and $s^{j}$ ($j=x,y,z$) are the Pauli matrices. The
coefficients $d^{a}$, on the other hand, are given by theoretical model and/or
the first-principles computation. Among models, the so-called Dirac model is a
low energy effective model to describe the topological insulators expanding a
Hamiltonian around a Dirac point. In the Dirac model, $d^{a}$ are given as

$\displaystyle(d^{1},\,d^{2},\,d^{3},\,d^{4},\,d^{5})=(Ak_{x},\,Ak_{y},\,Ak_{z},\,m+Bk^{2},m_{5})\,,$
(6)

where $k^{2}=|{\bf\it k}|^{2}$ and the gamma matrices are given in Eq. (5).
$\epsilon_{0}$, $A$, $B$, $m$, and $m_{5}$ have a mass dimension of one.
($\epsilon_{0}$ is not important in our later discussion, which is explained
in Appendix C.) The Dirac model can be derived from a certain model, such as
an extended Fu-Kane-Mele-Hubbard model on a diamond lattice at half-filling,
discussed in Refs. [21, 19], or the first-principles calculation of the
layered, stoichiometric crystals Sb2Te3, Sb2Se3, Bi2Te3 and Bi2Se3 [20]. In
general, the coefficient $A$ corresponds to the spin-orbit coupling.111To be
precise, the spin-orbit interaction itself is local and it is constant in
momentum space. The momentum dependence comes from nonlocal interaction
coupled to spin-orbit coupling.  If it is zero, then $\theta$ is exactly zero
and no topological phase appears, which will be easily understood from Eq.
(28). The coefficient $m$ and $B$ describe the energy gap. They originate in,
for instance, the spin-orbit coupling [20], or hopping strength anisotropy due
to the lattice distortion in the diamond lattice [21]. $m$ (and $B$)
effectively determines whether the insulator is topological or normal. Namely,
$m/B<0$ ($m/B>0$) corresponds to the topological (normal) state for $B\neq 0$.
When $B=0$, negative (positive) $m$ gives the topological (normal) state.
Finally $m_{5}$ is an order parameter of the AFM phase and plays an important
role for describing axion in insulators, which will be discussed in detail
soon.

While in the numerical study we adopt the Dirac model, the following
analytical calculation is applied for generic model described by Eq. (2).
Besides, in Appendix A we give another numerical results in the effective
model for 3D topological insulators given in Eq. (37) [9].

The $m_{5}\Gamma^{5}$ term appears due to the Hubbard-Stratonovich
transformation or the mean field approximation of the Hubbard interaction
term, which causes the AFM order written in a sublattice basis ($\ket*{{\rm
A},\sigma}$, $\ket*{{\rm B},\sigma}$),

$\displaystyle{\cal H}_{{\rm int}}=\frac{UV}{N}\int d^{3}x~{}(n_{{\rm
A}\uparrow}n_{{\rm A}\downarrow}+n_{{\rm B}\uparrow}n_{{\rm B}\downarrow})\,,$
(7)

where $n_{{\rm A}\sigma}=\psi^{\dagger}_{{\rm A}\sigma}\psi_{{\rm A}\sigma}$
and $n_{{\rm B}\sigma}=\psi^{\dagger}_{{\rm B}\sigma}\psi_{{\rm B}\sigma}$ are
number density of electron with spin $\sigma=\uparrow,\,\downarrow$ at
sublattice A and B, respectively. $U$ is a parameter with mass dimension one,
$V$ is the volume of the material and $N$ is the number of sublattice $A$
($B$) in the material. This interaction term involves four electrons and it is
difficult to analyze. However, it can be rewritten by the Hubbard-Stratonovich
transformation with introducing a scalar field $\phi$ in the path integral.
(See Appendix B for details.) As a result, the interaction term becomes

$\displaystyle\exp\left[-i\int dt{\cal H}_{{\rm int}}\right]$
$\displaystyle=\int{\cal D}\phi\exp\left\\{-i\int
d^{4}x\left[M^{2}\phi^{2}+\phi(n_{{\rm A}-}-n_{{\rm
B}-})\right]+\cdots\right\\}\,,$ (8)

where $n_{{\rm I}-}=n_{{\rm I}\uparrow}-n_{{\rm I}\downarrow}$ (${\rm I}=$A,
B), and

$\displaystyle M^{2}=\int\frac{d^{3}k}{(2\pi)^{3}}\frac{2}{U}\,.$ (9)

Here the momentum integral is executed in the first Brillouin zone unless
otherwise mentioned, which is applied to other momentum integral in the later
discussion. The term proportional to $n_{{\rm A}-}-n_{{\rm B}-}$ gives
$m_{5}\Gamma^{5\prime}$ term in the sublattice basis. (See Appendix D for the
explicit expression of $\Gamma^{a\prime}$.) Here we have renamed $m_{5}$ as

$\displaystyle\phi\equiv m_{5}\,.$ (10)

When $m_{5}=0$ the Hamiltonian is ${\cal T}$ and ${\cal P}$ conserving while a
finite value of $m_{5}$ breaks both. In addition, it should be noticed that
the mass term for $\phi$ emerges, which must be included in the effective
Lagrangian for $\phi$.

For later convenience we separate the Hamiltonian as

$\displaystyle H({\bf\it k})=H_{0}({\bf\it k})+\delta H({\bf\it k})\,,$ (11)

where

$\displaystyle H_{0}({\bf\it k})$ $\displaystyle=\epsilon_{0}({\bf\it k}){\bf
1}_{4\times 4}+\sum_{a=1}^{4}d^{a}_{0}({\bf\it k})\Gamma^{a}\,,$ (12)
$\displaystyle\delta H({\bf\it k})$ $\displaystyle=m_{5}\Gamma^{5}\,.$ (13)

Here $d_{0}^{a}$ is defined by

$\displaystyle d_{0}^{a}=(d^{1},\,d^{2},\,d^{3},\,d^{4},\,0)\,.$ (14)

From the Hamiltonian, the partition function is defined by functional integral
of four-component electron field $\psi$ and $\phi$,

$\displaystyle Z=\int{\cal D}\psi{\cal D}\psi^{\dagger}{\cal
D}\phi~{}e^{iS+iS_{\phi}^{\rm mass}}\,,$ (15)

where

$\displaystyle
S=\int\differential[4]{x}\psi^{\dagger}(x)\left[i\partial_{t}-H\right]\psi(x)\,,$
(16) $\displaystyle S_{\phi}^{\rm
mass}=-\int\differential[4]{x}M^{2}\phi^{2}\,.$ (17)

Here $H$ is the Hamiltonian in the coordinate space. Now we integrate out
$\psi$ in Eq. (15) to derive the effective Lagrangian for $\phi$. It is
straightforward to get

$\displaystyle Z$ $\displaystyle=\int{\cal D}\phi~{}e^{iS_{\phi}^{\rm
mass}}\,{\rm det}\left[-i\partial_{t}+H\right]$ $\displaystyle=\int{\cal
D}\phi~{}{\rm exp}\left\\{iS_{\phi}^{\rm mass}+{\rm
Tr}\log\left[-i\partial_{t}+H\right]\right\\}\,,$ (18)

Here “Tr” includes that trace over matrices and spacetime integral. To further
rewrite the above expression, let us define Green’s function as

$\displaystyle G_{0}^{-1}\equiv i\partial_{t}-H_{0}\,,$ (19)

where $H_{0}$ is $H_{0}({\bf\it k})$ written in the coordinate space. Then,
separating $H$ as $H=H_{0}+\delta H$ and using

$\displaystyle\log\left[-i\partial_{t}+H\right]$
$\displaystyle=\log(-G_{0}^{-1})-\sum_{n=1}^{\infty}\frac{1}{n}(G_{0}\delta
H)^{n}\,,$ (20)

we get

$\displaystyle Z$ $\displaystyle=\int{\cal D}\phi\,{\rm
exp}\Bigl{\\{}iS_{\phi}^{\rm mass}+{\rm Tr}\left[\log(-G^{-1}_{0})\right]$
$\displaystyle~{}~{}~{}~{}~{}~{}-{\rm
Tr}\Bigl{[}\sum_{n=1}^{\infty}\frac{1}{n}(G_{0}\delta
H)^{n}\Bigr{]}\Bigr{\\}}\,.$ (21)

Figure 1: $V_{\phi}$ as function of $\phi$ for various values of $m$ [eV] in
Dirac model. $A/a=1$ eV, $B/a^{2}=0.5$ eV is taken. At each panel, $U$ [eV] is
taken to 1.5 to 2.5 from top line to bottom line. We put a cutoff to the
normalized momentum integral as $-1<\tilde{{\bf\it k}}<1$.

Figure 2: $V_{\theta}$ as function of $\theta/(2\pi)$ for various values of
$m$ [eV] in Dirac model. The parameters and line contents are the same as in
Fig. 1.

The effective potential for $\phi$ is derived as polynomial expansion from Eq.
(21). Using the Fourier transformation of the propagator,

$\displaystyle G_{0}(x)$
$\displaystyle=\int\frac{d^{4}k}{(2\pi)^{4}}\tilde{G}_{0}(k)e^{-ik\dotproduct
x}\,,$ (22) $\displaystyle\tilde{G}_{0}(k$
$\displaystyle)=\frac{k^{0}-\epsilon_{0}({\bf\it
k})+\sum_{a=1}^{4}d_{0}^{a}({\bf\it k})\Gamma^{a}}{(k^{0}-\epsilon_{0}({\bf\it
k}))^{2}-|d_{0}({\bf\it k})|^{2}}\,,$ (23)

where $k\dotproduct x=k^{0}x^{0}-{\bf\it k}\dotproduct{\bf\it x}$ and
$|d_{0}|^{2}\equiv\sum_{a=1}^{4}d^{a}_{0}d^{a}_{0}$, it is obtained as (see
Appendix C for detail)

$\displaystyle V_{\phi}$
$\displaystyle=-2\int\frac{d^{3}k}{(2\pi)^{3}}(\sqrt{|d_{0}|^{2}+\phi^{2}}-|d_{0}|)+M^{2}\phi^{2}\,.$
(24)

A similar calculation is done in Refs. [15, 18, 22]. However, the expansion is
truncated at $n=2$. On top of that, the mass term, the second term in Eq. (24)
is missing.222There are other typos in Ref. [22]. We thank A. Sekine for
confirming this.  In fact, both the negative term derived from the last term
in Eq. (21) and the mass term are crucial to understand the AFM/FM phase and
the topological/normal state, which will be shown below.

Fig. 1 shows $V_{\phi}(\phi)$ as function of $\phi$. In the calculation we
introduce dimensionless wavenumber $\tilde{k}\equiv ka$ where $a$ is the
lattice size (typically Å) and take $A/a=1$ eV and $B/a^{2}=0.5$ eV333 A
larger value of $B$ is obtained in first-principles calculation [20]. This may
indicate it is difficult to have the AFM order in the materials such as Fe-
doped Bi2Te3 (see also later discussion). In our analysis, we choose
$\order{1}$ eV for $A$ and $B$.  and the integral is executed in $[-1,\,1]$
region. This is because the Dirac model is low energy effective model and the
(normalized) momentum cannot take as large as $\pi$. (Instead if the effective
model for 3D topological insulators given in Ref. [9] is used, then the
momentum integral is done in $[-\pi,\,\pi]$. See Appendix A for the results.)
It is seen that the curvature at $\phi=0$ becomes negative when $U$ gets
larger. This situation can be understood by expanding $V_{\phi}$ around
$\phi=0$,

$\displaystyle
V_{\phi}=\int\frac{d^{3}k}{(2\pi)^{3}}\Bigl{[}\frac{2}{U}-\frac{1}{|d_{0}|}\Bigr{]}\
\phi^{2}+{\cal O}(\phi^{4})\,,$ (25)

Here we have used Eq. (9). It is seen that the coefficient of $\phi^{2}$ term
becomes negative when $U>U_{\rm crit}$ where $U_{\rm crit}$ satisfies

$\displaystyle\int\frac{d^{3}k}{(2\pi)^{3}}\left[\frac{2}{U_{\rm
crit}}-\frac{1}{|d_{0}|}\right]=0\,.$ (26)

Namely, for $U<U_{\rm crit}$ the minimum of the potential located at $\phi=0$,
corresponding to the PM phase. When $U>U_{\rm crit}$, on the other hand, the
phase transition occurs to get nonzero $\phi_{0}$, which satisfies

$\displaystyle\int\frac{d^{3}k}{(2\pi)^{3}}\Bigl{(}\frac{2}{U}-\frac{1}{\sqrt{|d_{0}|^{2}+\phi_{0}^{2}}}\Bigr{)}=0\,,$
(27)

and the AFM order appears. In the next section, we show how the situation is
described in terms of the axion field.

##  III Axion potential

So far we have derived the effective potential for $\phi$. It can be described
by axion field because $\phi$ is related to the axion field $\theta$. To be
concrete, $\theta$ is given by $d^{a}$ [9],

$\displaystyle\theta=\frac{1}{4\pi}\int
d^{3}k\frac{2|d|+d^{4}}{(|d|+d^{4})^{2}|d|^{3}}\epsilon^{ijkl}d^{i}\partial_{k_{x}}d^{j}\partial_{k_{y}}d^{k}\partial_{k_{z}}d^{l}\,,$
(28)

where $|d|^{2}=\sum_{a=1}^{5}d^{a}d^{a}$ and $\epsilon^{ijkl}$ is Levi-Civita
symbol with $i,j,k,l$ being 1, 2, 3, 5. It is checked that $\theta=0$ or
$\pm\pi$ when $d^{5}=\phi=0$ and the value depends on $d^{4}$. (We note that
$\theta=+\pi$ and $-\pi$ are identical and hereafter we focus on $\theta=\pi$
unless otherwise stated. See also Appendix D for another derivation of
$\theta$ using chiral anomaly. However, it is not equivalent to the above
expression, except for a specific case.) For instance, in the Dirac model
given in Eq. (6), $m/B>0$ ($<0$) gives $\theta=0$ ($\pi$) when $d^{5}=0$ [23].
We have confirmed that the results are consistent with Ref. [23] (see Figs. 4
and 5 in Appendix A). Using Eq. (28), we can derive $\phi$ as function of
$\theta$. Let us write

$\displaystyle\phi=\Phi(\theta)\,.$ (29)

Then, the potential as function of $\theta$ can be obtained by

$\displaystyle V_{\phi}(\Phi(\theta))\equiv V_{\theta}(\theta)\,.$ (30)

Figure 3: Axion mass on ($U$, $m$) plane in the Dirac model. $A/a=1$ eV,
$B/a^{2}=0.5$ eV is taken. Solid (yellow) line shows the critical points.
Upper and lower regions separated by the solid line correspond to
antiferromagnetic order (“AFM”) and paramagnetic order (“PM”), respectively.
Dotted (green) line is boundary of topological insulator phase (“TI”) and
normal insulator phase (“NI”).

Fig. 2 shows $V_{\theta}$ as function of $\theta/(2\pi)$ in the Dirac model.
It is seen the axion potential nontrivially depends on the sign of $m$. When
$U$ is smaller than $U_{\rm crit}$, the potential minimum is located at
$\theta=0$ $(\pi)$ for $m>0$ $(m<0)$. $V_{\theta}$ is multivalued function of
$\theta$ since $\phi$ is multivalued function of $\theta$ for $m>0$. For
$m<0$, on the other hand, $\theta$ is monovalent function, i.e.,
$\theta\to\pm\pi$ for $\phi\to 0\pm$. (See Fig. 4 in Appendix A.) In terms of
potential minimum of $V_{\phi}$, both $\theta=0$ and $\pi$ correspond to
$\phi$ (or $m_{5}$) $=0$, i.e., the PM phase. In terms of $\theta$, on the
other hand, they are different topologically. Namely, $\theta=0$ ($\pi$) is
classified as the normal insulator (topological insulator). When $U$ gets
larger than $U_{\rm crit}$, nonzero $\phi_{0}$ becomes a minimum of the
potential and insulator is in the AFM order. In terms of axion, both
$\theta=0$ and $\pi$ becomes unstable and a new stable point appears. It
should be noticed that the effective potential given in Eq. (24) is crucial to
derive this result, which consists of the negative potential term derived from
integrating out the electron field and the positive mass term that comes
consequently from the Hubbard-Stratonovich transformation of the Hubbard
interaction term (7).

Now the definition of static and dynamical axions are clear. The static axion
$\theta_{0}$ is defined as a stationary point of $V_{\theta}(\theta)$, and the
dynamical axion $a$ as a fluctuation around $\theta_{0}$. Denoting $\phi_{0}$
as corresponding value of $\theta_{0}$, we expand $\phi$ and the potential
$V_{\phi}$,

$\displaystyle\phi$
$\displaystyle=\phi_{0}+\derivative{\Phi(\theta)}{\theta}\evaluated{}_{\theta=\theta_{0}}a+\cdots\,,$
$\displaystyle\equiv\phi_{0}+ga+\cdots\,,$ (31) $\displaystyle V_{\phi}$
$\displaystyle=V_{\phi}(\phi_{0})+\frac{1}{2}g^{2}\left[\derivative[2]{V_{\phi}}{\phi}\right]_{\phi=\phi_{0}}a^{2}+\cdots$
$\displaystyle\equiv V_{\phi}(\phi_{0})+Jg^{2}m_{a}^{2}a^{2}+\cdots\,,$ (32)

where $J$ is called stiffness and obtained as (see Appendix C for derivation)

$\displaystyle
J=\int\frac{d^{3}k}{(2\pi)^{3}}\frac{|d_{0}|^{2}}{4(|d_{0}|^{2}+\phi_{0}^{2})^{5/2}}\,.$
(33)

Using Eqs. (III) and (32), the action for the dynamical axion $a$ is obtained,

$\displaystyle
S_{a}=Jg^{2}\int\differential[4]{x}a\left[-\partial_{t}^{2}+({\bf\it
v}\cdot\nabla)^{2}-m_{a}^{2}\right]a+{\cal O}(a^{4})\,.$ (34)

Therefore, the axion mass is obtained by

$\displaystyle
Jm_{a}^{2}=\int\frac{d^{3}k}{(2\pi)^{3}}\left[\frac{2}{U}-\frac{1}{|d_{0}|}\right]\,,$
(35)

for $U<U_{\rm crit}$ and

$\displaystyle
Jm_{a}^{2}=\int\frac{d^{3}k}{(2\pi)^{3}}\frac{\phi_{0}^{2}}{(|d_{0}|^{2}+\phi_{0}^{2})^{3/2}}\,,$
(36)

for $U>U_{\rm crit}$. Eqs. (33) and (36) are consistent with the results in
Ref. [9], except for a factor of 4 larger compared to ones in the literature.
In the literature, however, how the value of $\phi_{0}$ is determined is not
clearly stated. Rather it is treated as another parameter that is independent
of the model parameter in $d^{a}$. By using the effective potential, on the
contrary, $\phi_{0}$ is determined uniquely as a stationary point of the
potential when the model parameters are fixed.

To demonstrate the result, we plot the axion mass on phase plane in Fig. 3.
There are phases classified by the AFM/PM and the topological/normal orders,
and they are determined by $U$ and $m$ for the given parameters $A$ and $B$ of
the Dirac model. The AFM (PM) phase corresponds to $\phi_{0}=$nonzero ($0$),
and the topological (normal) phase corresponds to $m<0$ ($m>0$). In the PM
phase, $\theta_{0}=0$ ($\pi$) in the normal (topological) phase meanwhile
$\theta_{0}$ takes a value of $0<\theta_{0}<\pi$ in the AFM phase depending on
$\phi_{0}$. It is found that the axion mass is typically $\order{1}$ eV in the
four phases. (The same result is obtained in the effective model for 3D
topological insulators. See Fig. 8 in Appendix A.) This result seems to be
different from the value argued in the literature [9] where the $\order{\rm
meV}$ axion in the topological AFM insulators is discussed. Let us elaborate
on this issue closely.

First of all, it is clear that the existence of dynamical axion does not
depend on whether insulator is topological or not. This point is mentioned in
Ref. [10]. Fig. 3 shows that the axion mass can be much suppressed near the
critical point of the AFM phase in both the topological and normal insulators.
In Ref. [9], on the other hand, $m_{5}\,(=\phi_{0})\,=1$ meV is taken, which
leads to meV scale axion mass. At a glance this looks possible; a small value
of $U$ can give a suppressed $m_{5}$ since $m_{5}\propto UM_{z}^{-}$ in mean
field approximation ($M_{z}^{-}$ in Ref. [9] corresponds to $m^{z}$ in
Appendix B.3). However, it is no longer the AFM order if $U<U_{\rm crit}$.
Therefore, suppressed $U$ cannot be naively considered in this simple setting.
This relates to an intrinsic issue that ${\rm Bi}_{2}{\rm Se}_{3}$ does not
have an AFM phase. A possible way out suggested in Ref. [9] is doping ${\rm
Fe}$ to create the AFM order. In this way, the AFM order could be realized.
However, doping Fe may not be successful due to the first-principles
calculation [24]. A current candidate is Mn-doped topological insulator,
Mn2Bi2Te5 suggested in Ref. [23]. The first-principles calculation done by
Ref. [25] shows that Mn2Bi2Te5 has rich magnetic topological quantum states,
including the AFM, FM and other magnetic states. The energy gap of the AFM
phase indicates $m_{5}\sim{\cal O}(10\,\mathchar 45\relax\,100~{}{\rm meV})$.
The total energy of each magnetic states are also found to be the same order,
which might indicate that controlling the realization of the magnetic states
would be challenging. Therefore, the next step is to confirm these magnetic
states experimentally.444Another interesting candidate is MnBi2Te4-family
materials. Although they are not supposed to have dynamical axion, they are
found to have controllable AFM and ferromagnetic orders. Therefore it is worth
studying their properties for the future application. See, for example, Refs.
[26, 27, 28, 29, 30, 31].

In a nutshell, the axion mass is roughly $\order{\rm eV}$ but it can be
suppressed near the critical boundary between the AFM and PM orders. Such
dynamical axion with various mass range might be utilized in the future
particle axion or ALPs search. On the other hand, in order to precisely
predict the axion mass in such specific circumstances, the first-principles
calculation to determine the parameters is required. Besides, we need
experimental evidence of having the AFM order in Fe- or Mn-doped Bismuth
Selenide or Bismuth Telluride.

It is worth noticing that the dynamical axion should exist in the PM
phase.555Ref. [11] pointed out the existence of dynamical axion in topological
PM insulators. On the other hand, there would be other degree of freedom in
the PM phase. For example, $\phi$ is introduced in the circumstance where the
AFM order is assumed (see, e.g., Eq. (18)). If this assumption is relaxed,
another scalar field should appear. Therefore, it would be nontrivial to find
out dynamical axion mode under such contaminated circumstances. If the
possible obstacles are overcome, however, the dynamical axion in the PM phase
might be another tool for the future particle axion search.

##  IV conclusions and discussion

We have investigated axion in AFM insulators. In the analysis, the effective
potential for the order parameter of the AFM phase is derived from the path
integral. From the effective potential, the order parameter of the AFM is
determined as the stationary point of the potential, which distinguishes the
AFM or PM order. Axion potential is obtained from the effective potential. It
is shown that the axion field value $\theta=0$ or $\pi$ is the minimum of the
axion potential in the PM order, meanwhile in the AFM order different values
of $\theta$ are obtained as the minimum of the potential. The field value of
the stationary point of axion is the static axion, and the dynamical axion is
defined as a fluctuation around stationary value. The dynamical axion is
predicted in four types of phases distinguished by the AFM/PM and the
topological/normal states. The mass of axion turns out to be typically
$\order{\rm eV}$. On the other hand, it can be much suppressed near the phase
boundary between the AFM and PM orders. In either case, it is crucial to
determine the value of the order parameter in a possible candidate, e.g.,
Mn2Bi2Te5, to predict the axion mass.

Another possible direction to investigate the static/dynamical axion is to
utilize the magnetoelectric effect. The magnetoelectric couplings in generic
insulators are obtained by the first-principles calculation by Ref. [32]. For
example, $\theta=1.3\times 10^{-3}$ and $1.07\pi$ are obtained for Cr2O3 and
Bi2Se3. In addition, if Bi is magnetized, the value of $\theta$ is shifted to
$\theta=\pi\pm 0.55$. On the experimental side, quantized Faraday and Kerr
rotation due to Bi2Se3 film is observed precisely [33]. Recently, anisotropic
topological magnetoelectric effect is observed in axion insulators [34]. Those
are examples of “axion electrodynamics” as discussed in Refs. [9, 10]. Such
methods using optical responses would make it possible to directly probe the
order parameter of the magnetic states and the properties of the
static/dynamical axion. The formalism given here would be a help for finding
and modeling realistic materials for the future particle axion or ALPs search.
Searching a wide range of materials including ones that have no topological
states or have ferromagnetic/paramagnetic phases might open a new possibility
for particle axion search.

###### Acknowledgements.

We are grateful to Makoto Naka for valuable discussions in the early stage of
this project and careful reading of the manuscript. We also thank Akihiko
Sekine for useful discussions. This work was supported by JSPS KAKENHI Grant
Numbers JP17K14278, JP17H02875, JP18H05542, and JP20H01894.

##  Appendix A Effective model for 3D topological insulators

We give numerical results in the effective model for 3D topological insulators
studied in Ref. [9]. The model is parametrized as,

$\displaystyle(d^{1},\,d^{2},\,d^{3},\,d^{4},\,d^{5})$
$\displaystyle=(A_{2}\sin k_{x},\,A_{2}\sin k_{y},\,A_{1}\sin k_{z},\,{\cal
M},m_{5})\,,$ $\displaystyle{\cal M}$
$\displaystyle=M-2B_{1}-4B_{2}+2B_{1}\cos k_{z}+2B_{2}(\cos k_{x}+\cos
k_{y})\,,$ (37)

where the gamma matrices are the same as Eq. (5). If we take $A_{1}=A_{2}=A$,
$B_{1}=B_{2}=-B$, and $M=m$, and expand around ${\bf\it k}=0$, this model
reduces to the Dirac model given in Eq. (6). Fig. 4 (5) shows $\theta$ as
function of $\phi$ ($M$ or $m$). As comparison, the results in the Dirac model
are also plotted. It is seen that the results are qualitatively the same as
the Dirac model.

Similar to Figs. 1, 2, and 3, $V_{\phi}$, $V_{\theta}$, and the axion mass
computed in the effective model for 3D topological insulators are plotted in
Figs. 6, 7, and 8, respectively.

Figure 4: $\theta$ as function of $\phi$ for $M\,{[\rm eV]}=\pm 0.01,$ $\pm
0.1,$ and $\pm 0.5$ in the effective model for 3D topological insulators with
$A_{1}/a=A_{2}/a=1$ eV and $B_{1}/a^{2}=B_{2}/a^{2}=-0.5$ eV (left) and in the
Dirac model given with $m=M$, $A=1$ eV, and $B=0.5$ eV (right).

Figure 5: $\theta$ as function of $M$ (or $m$) for $m_{5}\,[{\rm eV}]=\pm
0.01,$ $\pm 0.1,$ and $\pm 0.5$ in the effective model for 3D topological
insulators (left) and in the Dirac model (right). The other parameters are the
same as Fig. 4.

Figure 6: $V_{\phi}$ as function of $\phi$ for various values of $M$ in the
effective model for 3D topological insulators. $A_{1}$, $A_{2}$, $B_{1}$,
$B_{2}$, and $M$ are taken the same as Fig. 4.

Figure 7: $V_{\theta}$ as function of $\theta/(2\pi)$ in the effective model
for 3D topological insulators. Parameters are taken the same as Fig. 6.
Figure 8: Axion mass on ($U$, $M$) plane calculated in the effective model for
3D topological insulators. $A_{1}/a^{2}=A_{2}/a^{2}=1$ eV and
$B_{1}/a^{2}=B_{2}/a^{2}=-0.5$ eV are taken. Solid (yellow) line shows the
critical points; upper and lower regions separated by the solid line
correspond to antiferromagnetic order (“AFM”) and paramagnetic order (“PM”),
respectively. Dotted (green) line is boundary of topological insulator phase
(“TI”) and normal insulator phase (“NI”).

##  Appendix B Origin of $d^{5}\Gamma^{5}$ term

In this section, we discuss a possible origin of $\delta H$ term given in Eq.
(13). For pedagogical purpose, we provide three types of derivations;
derivations in the continuum space (Appendix B.1), in a discrete space
(Appendix B.2) and by the mean field approximation (Appendix B.3).

###  B.1 In continuum space

To begin with, we consider one of sublattices since it is enough to get the
basic idea of emerging a scalar degree of freedom that corresponds to
$m_{5}\,(=\phi)$. The interaction term to cause the AFM order is given by
following Hamiltonian:666This term corresponds to the continuum limit of Eq.
(23) by making a replacement, $\frac{1}{N}\sum_{i}\to\frac{1}{V}\int d^{3}x$
and $n_{i\sigma}\to\frac{V}{N}n_{\sigma}({\bf\it x})$.

$\displaystyle{\cal H}_{\rm
int}=\frac{UV}{N}\int\differential[3]{x}n_{\uparrow}({\bf\it
x})n_{\downarrow}({\bf\it x})\,,$ (1)

where arrow stands for the spin, $V$ is the volume of the material, $N$ is the
number of site, and

$\displaystyle n_{\sigma}({\bf\it
x})=\psi^{\dagger}_{\sigma}\psi_{\sigma}~{}~{}~{}~{}~{}~{}(\sigma=\,\uparrow\,,\downarrow)~{}\,.$
(2)

It is convenient to define the following quantities:

$\displaystyle n_{\pm}({\bf\it x})\equiv n_{\uparrow}({\bf\it x})\pm
n_{\downarrow}({\bf\it x})\,,$ (3)

and take

$\displaystyle\expectationvalue{n_{+}({\bf\it x})}={\rm const.}\,,$ (4)

since the total number of electrons are unchanged. (This corresponds to Eq.
(46) in the mean field approximation.) Therefore, only $n_{-}$ is a dynamical
degree of freedom, and ${\cal H}_{\rm int}$ becomes777Strictly speaking, there
is always a fluctuation around $\expectationvalue{n_{+}}$. However, it is not
one we are interested in. See also the discussion when two sublattices are
introduced.

$\displaystyle{\cal H}_{\rm
int}=-\frac{UV}{4N}\int\differential[3]{x}n_{-}^{2}({\bf\it x})+{\rm const.}$
(5)

For later discussion we define the Fourier expansion of $\psi_{\sigma}$ and
$n_{\sigma}$ as

$\displaystyle\psi_{\sigma}({\bf\it x})$
$\displaystyle=\frac{1}{\sqrt{V}}\sum_{{\bf\it k}}c_{{\bf\it
k}\sigma}e^{i{\bf\it k}\dotproduct{\bf\it x}}\,,$ (6) $\displaystyle
c_{{\bf\it k}\sigma}$
$\displaystyle=\frac{1}{\sqrt{V}}\int\differential[3]{x}\psi_{\sigma}({\bf\it
x})e^{-i{\bf\it k}\dotproduct{\bf\it x}}\,,$ (7) $\displaystyle
n_{\sigma}({\bf\it x})$ $\displaystyle=\frac{1}{V}\sum_{{\bf\it
k}}\tilde{n}_{\sigma}({\bf\it k})e^{i{\bf\it k}\dotproduct{\bf\it x}}\,,$ (8)
$\displaystyle\tilde{n}_{\sigma}({\bf\it k})$
$\displaystyle=\int\differential[3]{x}n_{\sigma}({\bf\it x})e^{-i{\bf\it
k}\cdot{\bf\it x}}=\sum_{{\bf\it q}}c^{\dagger}_{{\bf\it q}\sigma}c_{{\bf\it
q}+{\bf\it k}\sigma}\,.$ (9)

Consequently, ${\cal H}_{\rm int}$ in the momentum space is given by

$\displaystyle{\cal H}_{\rm int}=\frac{U}{N}\sum_{{\bf\it
k}}\tilde{n}_{\uparrow}({\bf\it k})\tilde{n}_{\downarrow}(-{\bf\it
k})=-\frac{U}{4N}\sum_{{\bf\it k}}\tilde{n}_{-}({\bf\it
k})\tilde{n}_{-}(-{\bf\it k})+{\rm const.}$ (10)

In terms of Euclidean time $\tau=it$, action of the interaction term is
written by

$\displaystyle-S_{{\rm int}\,E}\equiv iS_{\rm
int}\equiv-\int\differential{\tau}{\cal H}_{\rm int}\,.$ (11)

Now we execute the Hubbard-Stratonovich transformation. It is also used in
Ref. [35], which studies topological superconductors and superfluids. (See
Ref. [36] for the renormalization group approach.) It is effectively the same
as the inverse procedure of integrating out the heavy scalar particles to give
the effective Lagrangian. (See Appendix E for details.) Introducing
$\tilde{\phi}({\bf\it k})$ and using the following identity,

$\displaystyle-\Bigl{\\{}\frac{N}{U}\tilde{\phi}({\bf\it
k})\tilde{\phi}(-{\bf\it k})+\frac{1}{2}\left[\tilde{\phi}({\bf\it
k})\tilde{n}_{-}(-{\bf\it k})+\tilde{\phi}_{I}(-{\bf\it
k})\tilde{n}_{-}({\bf\it k})\right]\Bigr{\\}}$
$\displaystyle=-\frac{N}{U}\left[\tilde{\phi}({\bf\it
k})+\frac{U}{2N}\tilde{n}_{-}({\bf\it k})\right]\left[\tilde{\phi}(-{\bf\it
k})+\frac{U}{2N}\tilde{n}_{-}(-{\bf\it
k})\right]+\frac{U}{4N}\tilde{n}_{-}({\bf\it k})\tilde{n}_{-}(-{\bf\it k})\,,$
(12)

the Euclidean action gives

$\displaystyle\exp\left[-\int\differential{\tau}{\cal H}_{\rm int}\right]$
$\displaystyle=\int{\cal
D}\tilde{\phi}\exp\left\\{-\int\differential{\tau}\left\\{\sum_{{\bf\it
k}}\frac{N}{U}\tilde{\phi}({\bf\it k})\tilde{\phi}(-{\bf\it
k})+\frac{1}{2}\sum_{{\bf\it k}}\left[\tilde{\phi}({\bf\it
k})\tilde{n}_{-}(-{\bf\it k})+\tilde{\phi}(-{\bf\it k})\tilde{n}_{-}({\bf\it
k})\right]\right\\}\right\\}$ $\displaystyle=\int{\cal
D}\phi\exp\left\\{-\int\differential{\tau}\differential[3]{x}\left[\frac{N}{VU}\phi^{2}({\bf\it
x})+\phi({\bf\it x})n_{-}({\bf\it x})\right]\right\\}\,.$ (13)

Here we have omitted irrelevant constant term in the last step and we have
used

$\displaystyle\tilde{\phi}({\bf\it
k})=\frac{1}{V}\int\differential[3]{x}\phi(x)e^{-i{\bf\it k}\cdot{\bf\it
x}}\,.$ (14)

The second term on the rhs of Eq. (13) is given in the momentum space as

$\displaystyle\int\differential[3]{x}\phi({\bf\it x})n_{-}({\bf\it x})$
$\displaystyle=\sum_{{\bf\it q},{\bf\it k}}(c^{\dagger}_{{\bf\it
k}\uparrow}c_{{\bf\it k}+{\bf\it q}\uparrow}-c^{\dagger}_{{\bf\it
k}\downarrow}c_{{\bf\it k}+{\bf\it q}\downarrow})\tilde{\phi}(-{\bf\it q})$
$\displaystyle=\sum_{{\bf\it k}}(c^{\dagger}_{{\bf\it k}\uparrow}c_{{\bf\it
k}\uparrow}-c^{\dagger}_{{\bf\it k}\downarrow}c_{{\bf\it
k}\downarrow})\tilde{\phi}(0)+{\cal O}({\bf\it q})$
$\displaystyle=\sum_{{\bf\it k},\alpha,\beta}c^{\dagger}_{{\bf\it
k}\alpha}(\tilde{\phi}(0)\sigma^{z})_{\alpha\beta}c_{{\bf\it k}\beta}+{\cal
O}({\bf\it q})\,.$ (15)

Now we apply the above result to two sublattices, A and B. The interaction
term to start with is

$\displaystyle{\cal H}_{\rm int}^{{\rm A}+{\rm
B}}=\frac{UV}{N}\int\differential[3]{x}(n_{{\rm A}\uparrow}({\bf\it x})n_{{\rm
A}\downarrow}({\bf\it x})+n_{{\rm B}\uparrow}({\bf\it x})n_{{\rm
B}\downarrow}({\bf\it x}))\,,$ (16)

where $N$ is the number of the sublattice A (B). As in the previous
discussion, we define $n_{{\rm I}\pm}({\bf\it x})$ (I=A, B),

$\displaystyle n_{{\rm I}\pm}({\bf\it x})\equiv n_{{\rm I}\uparrow}({\bf\it
x})\pm n_{{\rm I}\downarrow}({\bf\it x})\,.$ (17)

Since we are interested in the AFM order, $n_{{\rm A}-}$ and $n_{{\rm B}-}$
are not independent and related by888As discussed below Eq. (4), there are
fluctuations around $\expectationvalue{n_{{\rm A}+}}$ and
$\expectationvalue{n_{{\rm B}+}}$. In addition, if there is no the AFM order,
i.e., the PM order, then $\expectationvalue{n_{{\rm A}-}}$ and
$\expectationvalue{n_{{\rm B}-}}$ should be treated as independent degrees of
freedom.

$\displaystyle n_{{\rm A}-}+n_{{\rm B}-}\approx\expectationvalue{n_{{\rm
A}-}+n_{{\rm B}-}}=0\,.$ (18)

Therefore, we obtain

$\displaystyle\exp\left[-\int\differential{\tau}{\cal H}_{\rm int}^{{\rm
A}+{\rm B}}\right]$ $\displaystyle=\int{\cal
D}\tilde{\phi}\exp\left\\{-\int\differential{\tau}\left\\{\sum_{{\bf\it
k}}\frac{2N}{U}\tilde{\phi}({\bf\it k})\tilde{\phi}(-{\bf\it
k})+\frac{1}{2}\sum_{{\bf\it k}}\left[\tilde{\phi}({\bf\it k})(\tilde{n}_{{\rm
A}-}(-{\bf\it k})-\tilde{n}_{{\rm B}-}(-{\bf\it k}))+\tilde{\phi}(-{\bf\it
k})(\tilde{n}_{{\rm A}-}({\bf\it k})-\tilde{n}_{{\rm B}-}({\bf\it
k}))\right]\right\\}\right\\}$ $\displaystyle=\int{\cal
D}\phi\exp\left\\{-\int\differential{\tau}\differential[3]{x}\left[\frac{2N}{VU}\phi^{2}({\bf\it
x})+\phi({\bf\it x})(n_{{\rm A}-}({\bf\it x})-n_{{\rm B}-}({\bf\it
x}))\right]\right\\}\,.$ (19)

and

$\displaystyle\int\differential[3]{x}\phi({\bf\it x})(n_{{\rm A}-}({\bf\it
x})-n_{{\rm B}-}({\bf\it x}))=\sum_{{\bf\it
k}}\tilde{\phi}(0)\left\\{c^{\dagger}_{{\rm A}{\bf\it k}}\sigma^{z}c_{{\rm
A}{\bf\it k}}-c^{\dagger}_{{\rm B}{\bf\it k}}\sigma^{z}c_{{\rm B}{\bf\it
k}}\right\\}+{\cal O}({\bf\it q})\,.$ (20)

Here the summation over spin indices are implicit. This term corresponds to
$d^{5}\Gamma^{5\prime}\,(=\phi\Gamma^{5\prime})$ term in Eq. (2) after
changing the basis from ($\ket{P1_{z}^{+},\sigma}$, $\ket{P2_{z}^{-},\sigma}$)
to $(\ket{{\rm A},\sigma},\ket{{\rm B},\sigma})$ where
$\sigma=\,\uparrow,\downarrow$. In this basis, ${\Gamma}^{5\prime}={\rm
diag}\\{1,-1,-1,1\\}$ and Eq. (19) becomes

$\displaystyle\int{\cal
D}\phi\exp\left\\{-\int\differential{\tau}\differential[3]{x}\left[M^{2}\phi^{2}({\bf\it
x})-\phi({\bf\it x})(n_{{\rm A}-}({\bf\it x})-n_{{\rm B}-}({\bf\it
x}))\right]+\cdots\right\\}$ $\displaystyle=\int{\cal
D}\tilde{\phi}\exp\left\\{-\int\differential{\tau}\sum_{{\bf\it
k}}\left[VM^{2}\tilde{\phi}({\bf\it k})\tilde{\phi}(-{\bf\it k})-{\bf\it
c}^{\dagger}_{{\bf\it k}}\Gamma^{5\prime}\tilde{\phi}(0){\bf\it c}_{{\bf\it
k}}+\cdots\right]\right\\}\,,$ (21)

where ${\bf\it c}_{{\bf\it k}}=(c_{{\rm A}{\bf\it k}\uparrow},\,c_{{\rm
A}{\bf\it k}\downarrow},\,c_{{\rm B}{\bf\it k}\uparrow},\,c_{{\rm B}{\bf\it
k}\downarrow})$ and999We use $(V/N)^{-1}=\int\frac{d^{3}q}{(2\pi)^{3}}$ where
the integral is defined in the first Brillouin zone.

$\displaystyle
M^{2}\equiv\frac{2}{(V/N)U}=\left(\int\frac{d^{3}q}{(2\pi)^{3}}\right)\frac{2}{U}\,.$
(22)

###  B.2 In discrete space

As in the previous subsection, we first consider one of the sublattices and
extend the result to the other later. In the discrete space, we use
annihilation operator $c_{i\sigma}$ where $i$ is a label of site instead of
the wavefunction $\psi_{\sigma}({\bf\it x})$. The interaction Hamiltonian is

$\displaystyle{\cal H}_{\rm int}=U\sum_{i}n_{i\uparrow}n_{i\downarrow}\,,$
(23)

where

$\displaystyle n_{i\sigma}=c_{i\sigma}^{\dagger}c_{i\sigma}\,.$ (24)

While defining $n_{i\pm}$ as

$\displaystyle n_{i\pm}\equiv n_{i\uparrow}\pm n_{i\downarrow}\,,$ (25)

only $n_{i-}$ is dynamical since $n_{i+}=$const. as in Eq. (4). Consequently,
the interaction Hamiltonian is given by

$\displaystyle{\cal H}_{\rm int}=\frac{U}{4}\sum_{i}n_{i-}^{2}+{\rm
const.}\,.$ (26)

It is useful to give Fourier transformation of $c_{i\sigma}$ and $n_{i\sigma}$
,

$\displaystyle c_{i\sigma}$ $\displaystyle=\frac{1}{\sqrt{N}}\sum_{{\bf\it
k}}c_{{\bf\it k}\sigma}e^{i{\bf\it k}\dotproduct{\bf\it x}}\,,$ (27)
$\displaystyle c_{{\bf\it k}\sigma}$
$\displaystyle=\frac{1}{\sqrt{N}}\sum_{i}c_{i\sigma}e^{-i{\bf\it
k}\dotproduct{\bf\it x}}\,,$ (28) $\displaystyle n_{i\sigma}$
$\displaystyle=\frac{1}{N}\sum_{{\bf\it k}}\tilde{n}_{{\bf\it k}}e^{i{\bf\it
k}\dotproduct{\bf\it x}}\,,$ (29) $\displaystyle\tilde{n}_{{\bf\it k}\sigma}$
$\displaystyle=\sum_{i}n_{i\sigma}e^{-i{\bf\it k}\dotproduct{\bf\it
x}}=\sum_{{\bf\it q}}c^{\dagger}_{{\bf\it q}\sigma}c_{{\bf\it q}+{\bf\it
k}\sigma}\,,$ (30)

which leads to the Hamiltonian in the momentum space,

$\displaystyle{\cal H}_{\rm int}=\frac{U}{N}\sum_{{\bf\it
k}}\tilde{n}_{{\bf\it k}\uparrow}\tilde{n}_{-{\bf\it
k}\downarrow}=-\frac{U}{4N}\sum_{{\bf\it k}}\tilde{n}_{{\bf\it
k}-}\tilde{n}_{-{\bf\it k}-}+{\rm const.}\,$ (31)

As in the continuum case, the Hubbard-Stratonovich transformation is done by
using Eq.(12). Then the Euclidean action of the interaction term is given by

$\displaystyle\exp\left[-\int d\tau{\cal H}_{\rm int}\right]$
$\displaystyle=\int{\cal
D}\tilde{\phi}\exp\left\\{-\int\differential{\tau}\left\\{\sum_{{\bf\it
k}}\frac{N}{U}\tilde{\phi}_{{\bf\it k}}\tilde{\phi}_{-{\bf\it
k}}+\sum_{{\bf\it k}}\frac{1}{2}\left[\tilde{\phi}_{{\bf\it
k}}\tilde{n}_{-{\bf\it k}-}+\tilde{\phi}_{-{\bf\it k}}\tilde{n}_{{\bf\it
k}-}\right]\right\\}\right\\}$ $\displaystyle=\int{\cal
D}\tilde{\phi}\exp\left\\{-\int\differential{\tau}\sum_{i}\left[\frac{1}{U}\phi_{i}^{2}+\phi_{i}n_{i-}\right]\right\\}\,.$
(32)

Here we have defined the Fourier expansion of $\phi_{i}$ as

$\displaystyle\phi_{i}$ $\displaystyle=\sum_{{\bf\it k}}\tilde{\phi}_{{\bf\it
k}}e^{i{\bf\it k}\dotproduct{\bf\it x}}\,,$ (33)
$\displaystyle\tilde{\phi}_{{\bf\it k}}$
$\displaystyle=\frac{1}{N}\sum_{i}\phi_{i}e^{-i{\bf\it k}\dotproduct{\bf\it
x}}\,,$ (34)

and the second term on the rhs of Eq. (32) is written in the momentum space as

$\displaystyle\sum_{i}\phi_{i}n_{i-}$ $\displaystyle=\sum_{{\bf\it q},{\bf\it
k}}(c^{\dagger}_{{\bf\it k}\uparrow}c_{{\bf\it k}+{\bf\it
q}\uparrow}-c^{\dagger}_{{\bf\it k}\downarrow}c_{{\bf\it k}+{\bf\it
q}\downarrow})\tilde{\phi}_{-{\bf\it q}}$ $\displaystyle=\sum_{{\bf\it
k}}(c^{\dagger}_{{\bf\it k}\uparrow}c_{{\bf\it
k}\uparrow}-c^{\dagger}_{{\bf\it k}\downarrow}c_{{\bf\it
k}\downarrow})\tilde{\phi}_{\mathbf{0}}+{\cal O}({\bf\it q})$
$\displaystyle=\sum_{{\bf\it k},\alpha,\beta}c^{\dagger}_{{\bf\it
k}\alpha}(\tilde{\phi}_{\mathbf{0}}\sigma^{z})_{\alpha\beta}c_{{\bf\it
k}\beta}+{\cal O}({\bf\it q})\,.$ (35)

Extension to the other sublattice is trivial. As in the previous subsection,
let us call the two sublattice as A and B. Then the interaction Hamiltonian is

$\displaystyle{\cal H}_{\rm int}^{\rm A+B}=U\sum_{i}(n_{{\rm
A}i\uparrow}n_{{\rm A}i\downarrow}+n_{{\rm B}i\uparrow}n_{{\rm
B}i\downarrow})\,.$ (36)

Repeating the discussion so far, the Euclidean action is obtained as

$\displaystyle\exp\left[-\int\differential{\tau}{\cal H}_{\rm int}^{{\rm
A}+{\rm B}}\right]$ $\displaystyle=\int{\cal
D}\tilde{\phi}\exp\left\\{-\int\differential{\tau}\left\\{\sum_{{\bf\it
k}}\frac{2N}{U}\tilde{\phi}_{{\bf\it k}}\tilde{\phi}_{-{\bf\it
k}}+\frac{1}{2}\sum_{{\bf\it k}}\left[\tilde{\phi}_{{\bf\it
k}}(\tilde{n}_{{\rm A}-{\bf\it k}-}-\tilde{n}_{{\rm B}-{\bf\it
k}-})+\tilde{\phi}_{-{\bf\it k}}(\tilde{n}_{{\rm A}{\bf\it
k}-}-\tilde{n}_{{\rm B}{\bf\it k}-})\right]\right\\}\right\\}$
$\displaystyle=\int{\cal
D}\phi\exp\left\\{-\int\sum_{i}\left[\frac{2}{U}\phi^{2}_{i}+\phi_{i}(n_{i{\rm
A}-}-n_{i{\rm B}-})\right]\right\\}\,.$ (37)

The interaction part in momentum space is

$\displaystyle\sum_{i}\phi_{i}(n_{i\rm{A}-}-n_{i\rm{B}-})=\sum_{{\bf\it
k}}\tilde{\phi}_{\mathbf{0}}\left\\{c^{\dagger}_{{\rm A}{\bf\it
k}}\sigma^{z}c_{{\rm A}{\bf\it k}}-c^{\dagger}_{{\rm B}{\bf\it
k}}\sigma^{z}c_{{\rm B}{\bf\it k}}\right\\}+{\cal O}({\bf\it q})\,.$ (38)

Therefore we have obtained the same result as Eq.(20) but
$\tilde{\phi}_{2}(0)$ is replaced by $\tilde{\phi}_{\mathbf{0}}$. Finally Eq.
(37) becomes

$\displaystyle\int{\cal
D}\phi\exp\left\\{-\int\differential{\tau}\sum_{i}\left[\frac{2}{U}\phi_{i}^{2}+{\bf\it
c}_{i}^{\dagger}\Gamma^{5\prime}\phi_{i}{\bf\it
c}\right]+\cdots\right\\}=\int{\cal
D}\tilde{\phi}\exp\left\\{-\int\differential{\tau}\sum_{{\bf\it
k}}\left[\frac{2N}{U}\tilde{\phi}_{{\bf\it k}}\tilde{\phi}_{-{\bf\it
k}}+{\bf\it c}^{\dagger}_{{\bf\it k}}\Gamma^{5\prime}\tilde{\phi}_{{\bf\it
0}}{\bf\it c}_{{\bf\it k}}+\cdots\right]\right\\}\,,$ (39)

where ${\bf\it c}_{i}=(c_{{\rm A}i\uparrow},\,c_{{\rm A}i\downarrow},\,c_{{\rm
B}i\uparrow},\,c_{{\rm B}i\downarrow})$. When we make a continuum limit for
space, i.e.,

$\displaystyle\frac{1}{N}\sum_{i}$
$\displaystyle\to\frac{1}{V}\int\differential[3]{x}\,,$ (40)
$\displaystyle\phi_{i}$ $\displaystyle\to\phi({\bf\it x})\,,$ (41)

the mass term becomes

$\displaystyle\sum_{i}\frac{2}{U}\phi^{2}_{i}\to\int\differential[3]{x}\frac{2}{(V/N)U}\phi^{2}({\bf\it
x})\,,$ (42)

which agrees with Eq. (21).

###  B.3 Mean field approximation

Another way to see the appearance of $d^{5}\Gamma^{5}$ term is use the mean
field approximation (MFA). We refer to Refs. [21, 19]. In the MFA,

$\displaystyle n_{i\uparrow}n_{i\downarrow}\approx
n_{i\uparrow}\expectationvalue{n_{i\downarrow}}+\expectationvalue{n_{i\uparrow}}n_{i\downarrow}-\expectationvalue{n_{i\uparrow}}\expectationvalue{n_{i\downarrow}}-\expectationvalue*{c^{\dagger}_{i\uparrow}c_{i\downarrow}}c^{\dagger}_{i\downarrow}c_{i\uparrow}-\expectationvalue*{c^{\dagger}_{i\downarrow}c_{i\uparrow}}c^{\dagger}_{i\uparrow}c_{i\downarrow}+\expectationvalue*{c^{\dagger}_{i\uparrow}c_{i\downarrow}}\expectationvalue*{c^{\dagger}_{i\downarrow}c_{i\uparrow}}\,.$
(43)

The expectation values are written by introducing ${\bf\it m}$

$\displaystyle m^{j}$
$\displaystyle=\frac{1}{2}\sum_{\alpha,\beta}\expectationvalue*{c^{\dagger}_{i\alpha}\sigma^{j}_{\alpha\beta}c_{i\beta}}~{}~{}~{}~{}~{}~{}~{}~{}(j=x,y)\,,$
(44) $\displaystyle m^{z}$
$\displaystyle=\frac{1}{2}\sum_{\alpha,\beta}\expectationvalue*{c^{\dagger}_{i\alpha}\sigma^{z}_{\alpha\beta}c_{i\beta}}+\frac{1}{2}\expectationvalue*{c^{\dagger}_{i\uparrow}c_{i\uparrow}+c^{\dagger}_{i\downarrow}c_{i\downarrow}}\,,$
(45)

and we take

$\displaystyle\frac{1}{2}\expectationvalue*{c^{\dagger}_{i\uparrow}c_{i\uparrow}+c^{\dagger}_{i\downarrow}c_{i\downarrow}}=1\,,$
(46)

since the half-filling model is considered. Then, it is straightforward to get

$\displaystyle n_{i\uparrow}n_{i\downarrow}\approx{\bf\it
m}^{2}-\sum_{\alpha,\beta}c^{\dagger}_{i\alpha}{\bf\it
m}\dotproduct{\bf\it\sigma}_{\alpha\beta}c_{i\beta}\,.$ (47)

Extending the argument to two sublattices is trivial. Now taking

$\displaystyle\expectationvalue*{S_{{\rm A}}^{j}}$
$\displaystyle=m^{j}=\frac{1}{2}\sum_{\alpha,\beta}\expectationvalue*{c^{\dagger}_{{\rm
A}i\alpha}\sigma^{j}_{\alpha\beta}c_{{\rm
A}i\beta}}~{}~{}~{}~{}~{}~{}~{}~{}(j=x,y)\,,$ (48)
$\displaystyle\expectationvalue*{S_{{\rm A}}^{z}}$
$\displaystyle=m^{z}=\frac{1}{2}\sum_{\alpha,\beta}\expectationvalue*{c^{\dagger}_{{\rm
A}i\alpha}\sigma^{z}_{\alpha\beta}c_{{\rm
A}i\beta}}+\frac{1}{2}\expectationvalue*{c^{\dagger}_{{\rm A}i\uparrow}c_{{\rm
A}i\uparrow}+c^{\dagger}_{{\rm A}i\downarrow}c_{{\rm A}i\downarrow}}\,,$ (49)

and

$\displaystyle\expectationvalue*{{\bf\it S}_{\rm
A}}+\expectationvalue*{{\bf\it S}_{\rm B}}=0\,,$ (50)

since we are interested in the AFM order, we get

$\displaystyle U\sum_{i}(n_{{\rm A}i\uparrow}n_{{\rm A}i\downarrow}+n_{{\rm
B}i\uparrow}n_{{\rm B}i\downarrow})\approx 2U\sum_{i}{\bf\it
m}^{2}-U\sum_{i}{c}^{\dagger}_{{\rm A}i}{\bf\it
m}\dotproduct{\bf\it\sigma}c_{{\rm A}i}+U\sum_{i}{c}^{\dagger}_{{\rm
B}i}{\bf\it m}\dotproduct{\bf\it\sigma}c_{{\rm B}i}\,.$ (51)

In the current model, we consider that only $m^{z}$ is nonzero. Changing the
variable to $\phi$ by $\phi=-Um^{z}$, we finally obtain

$\displaystyle U\sum_{i}(n_{{\rm A}i\uparrow}n_{{\rm A}i\downarrow}+n_{{\rm
B}i\uparrow}n_{{\rm
B}i\downarrow})\approx\sum_{i}\frac{2}{U}\phi^{2}+\sum_{i}{\bf\it
c}_{i}^{\dagger}\phi\Gamma^{5\prime}{\bf\it c}_{i}\,.$ (52)

It is seen that both mass term for $\phi$ and $m_{5}\Gamma^{5\prime}$ term
appear.

##  Appendix C Derivation of effective Lagrangian for $\phi$

We show how to derive the effective potential (24) and stiffness $J$ (33) from
${\rm Tr}(G_{0}\delta H)^{n}$ term in Eq. (21). It corresponds to calculation
of one-loop diagram with external scalar fields those three-momenta are zero.
It is similar to the Coleman-Weinberg potential [37], but the the result turns
out to be non-logarithmic function. In the following calculation, the
continuum limit in the coordinate and momentum space is taken. We define the
Fourier expansion of $G_{0}$ and $\phi$ as

$\displaystyle G_{0}(x-y)$
$\displaystyle=\int\frac{d^{4}k}{(2\pi)^{4}}\tilde{G}_{0}(k)e^{-ik\cdot(x-y)}\,,$
(53) $\displaystyle\phi(x)$
$\displaystyle=\int\frac{d^{4}k}{(2\pi)^{4}}\tilde{\phi}(k)e^{-ik\cdot x}\,,$
(54)

where $k\cdot x=k^{0}x^{0}-{\bf\it k}\cdot{\bf\it x}$.

First we show how to derive the effective potential. Note that the trace
vanishes when $n$ is odd. Thus, let us begin with $n=2$:

$\displaystyle{\rm Tr}(G_{0}\delta H)^{2}$
$\displaystyle=\int\prod_{i=1}^{3}\differential[4]{x_{i}}~{}{\rm
tr}\left[G_{0}(x_{1}-x_{2})\phi(x_{2})G_{0}(x_{2}-x_{3})\phi(x_{3})\right]\delta^{(4)}(x_{3}-x_{1})$
$\displaystyle=\int\frac{d^{4}k}{(2\pi)^{4}}\frac{d^{4}q}{(2\pi)^{4}}~{}{\rm
tr}\left[\tilde{G}_{0}(q)\Gamma^{5}\tilde{G}_{0}(q+k)\Gamma^{5}\right]\tilde{\phi}(-k)\tilde{\phi}(k)\,.$
(55)

Here trace part gives

$\displaystyle{\rm
tr}\left[\tilde{G}_{0}(q)\Gamma^{5}\tilde{G}_{0}(q+k)\Gamma^{5}\right]=4\frac{\left\\{q^{0}-\epsilon_{0}({\bf\it
q})\right\\}\left\\{q^{0}+k^{0}-\epsilon_{0}({\bf\it q}+{\bf\it
k})\right\\}-d_{0}^{a}({\bf\it q})d_{0}^{a}({\bf\it q}+{\bf\it
k})}{((q^{0}-\epsilon_{0}({\bf\it q}))^{2}-|d_{0}({\bf\it
q})|^{2})((q^{0}+k^{0}-\epsilon_{0}({\bf\it q}+{\bf\it k}))^{2}-|d_{0}({\bf\it
q}+{\bf\it k})|^{2})}\,,$ (56)

where index $a$ is contracted. To give the effective potential we can take
$k=0$. As a result, it is given in a simple form,

$\displaystyle\left.{\rm
tr}\left[\tilde{G}_{0}(q)\Gamma^{5}\tilde{G}_{0}(q+k)\Gamma^{5}\right]\right|_{k=0}=4\frac{1}{(q^{0}-\epsilon_{0}({\bf\it
q}))^{2}-|d_{0}({\bf\it q})|^{2}}\,.$ (57)

In the similar manner, higher polynomials are calculated as

$\displaystyle{\rm Tr}(G_{0}\delta H)^{n}$
$\displaystyle=\int\prod_{i=1}^{n+1}\differential[4]{x_{i}}{\rm
tr}[G_{0}(x_{1}-x_{2})\phi(x_{2})G_{0}(x_{2}-x_{3})\phi(x_{3})\cdots
G_{0}(x_{n}-x_{n+1})\phi(x_{n+1})]\delta^{(4)}(x_{n+1}-x_{1})$
$\displaystyle=\int\prod_{i=1}^{n-1}\frac{d^{4}k_{i}}{(2\pi)^{4}}\frac{d^{4}q}{(2\pi)^{4}}~{}{\rm
tr}[\tilde{G}_{0}(q)\Gamma^{5}\tilde{G}_{0}(q+k_{1})\Gamma^{5}\tilde{G}_{0}(q+k_{1}+k_{2})\Gamma^{5}\cdots\tilde{G}_{0}(q+\sum_{i}^{n-1}k_{i})\Gamma^{5}]$
$\displaystyle~{}~{}~{}~{}~{}~{}\times\tilde{\phi}(-k_{1})\tilde{\phi}(-k_{2})\tilde{\phi}(-k_{3})\cdots\tilde{\phi}(-k_{n-1})\tilde{\phi}(\sum_{i}^{n-1}k_{i})\,,$
(58)

where $n\geq 2$. To get the $\phi^{n}$ terms, we can take $k_{i}=0$. Then the
trace part is

$\displaystyle\left.{\rm tr}[\cdots]\right|_{k_{i}=0}$ $\displaystyle={\rm
tr}[(\tilde{G}_{0}(q)\Gamma^{5})^{n}]$
$\displaystyle=4\frac{1}{\left[(q^{0}-\epsilon_{0}({\bf\it
q}))^{2}-|d_{0}({\bf\it q})|^{2}\right]^{m}}\,.$ (59)

where $n=2m$ $(m=1,2,3,\cdots)$. Using this formula, the last term in Eq. (21)
gives

$\displaystyle-{\rm Tr}\left[\sum_{n=1}^{\infty}\frac{1}{n}(G_{0}\delta
H)^{n}\right]\supset$
$\displaystyle-4i\int\differential[4]{x}\int\frac{d^{3}q}{(2\pi)^{3}}\frac{dq^{0}_{E}}{2\pi}\left[\sum_{m=1}^{\infty}\frac{1}{2m}\frac{(-1)^{m}\phi^{2m}}{\left\\{(q^{0}_{E})^{2}+|d_{0}|^{2}\right\\}^{m}}\right]$
$\displaystyle=2i\int\differential[4]{x}\int\frac{d^{3}q}{(2\pi)^{3}}\frac{dq^{0}_{E}}{2\pi}\log\left[1+\frac{\phi^{2}}{(q^{0}_{E})^{2}+|d_{0}|^{2}}\right]\,.$
(60)

Here we assume that $q^{0}$ integral is from $-\infty$ to $+\infty$, and we
use the Wick rotation of the the integral path as

$\displaystyle\int_{-\infty}^{\infty}\differential{q^{0}}=i\int_{-\infty}^{\infty}\differential{q_{E}^{0}}\,,$
(61)

Finally the effective potential for $\phi\,(=m_{5})$ is obtained as

$\displaystyle V_{\phi}$
$\displaystyle=-2\int\frac{d^{3}q}{(2\pi)^{3}}\frac{dq^{0}_{E}}{2\pi}\log\left[1+\frac{\phi^{2}}{(q^{0}_{E})^{2}+|d_{0}|^{2}}\right]+M^{2}\phi^{2}$
$\displaystyle=-2\int\frac{d^{3}q}{(2\pi)^{3}}(\sqrt{|d_{0}|^{2}+\phi^{2}}-|d_{0}|)+M^{2}\phi^{2}\,.$
(62)

The stiffness $J$ is obtained as an coefficient of the time-derivative term of
the dynamical field. It is determined after the potential minimum is found.
Therefore, instead of using Eqs .(12) and (13), we take

$\displaystyle H_{0}^{\prime}({\bf\it k})$ $\displaystyle=\epsilon_{0}({\bf\it
k}){\bf 1}_{4\times 4}+\sum_{a=1}^{5}d^{a}({\bf\it k})\Gamma^{a}\,,$ (63)
$\displaystyle\delta H^{\prime}({\bf\it k})$
$\displaystyle=\delta\phi\Gamma^{5}\,,$ (64)

where now $d^{5}=\phi_{0}$ (and the others are the same) and
$\delta\phi=\phi-\phi_{0}$. The propagator is given similarly to Eq. (23) as

$\displaystyle\tilde{G}(k)=\frac{k^{0}-\epsilon_{0}({\bf\it k})+d^{a}({\bf\it
k})\Gamma^{a}}{(k^{0}-\epsilon_{0}({\bf\it k}))^{2}-|d({\bf\it k})|^{2}}\,,$
(65)

where $|d|^{2}\equiv\sum_{a=1}^{4}d^{a}d^{a}+\phi_{0}^{2}$. To derive the
time-derivative kinetic term, the quadratic term of $\delta H^{\prime}$ should
be expanded in terms of $k^{0}$ (but ${\bf\it k}=0$ can be taken). The trace
part is given by

$\displaystyle{\rm
tr}\Bigl{[}\tilde{G}(q)\Gamma^{5}\tilde{G}(q+k)\Gamma^{5}\Bigr{]}\evaluated{}_{{\bf\it
k}=0}$ $\displaystyle=4\frac{(q^{0}-\epsilon_{0}({\bf\it
q}))^{2}+(q^{0}-\epsilon_{0}({\bf\it q}))k^{0}-|d_{0}({\bf\it
q})|^{2}+\phi_{0}^{2}}{((q^{0}-\epsilon_{0}({\bf\it q}))^{2}-|d({\bf\it
q})|^{2})((q^{0}+k^{0}-\epsilon_{0}({\bf\it q}))^{2}-|d({\bf\it q})|^{2})}\,.$
(66)

After expanding with respect to $k^{0}$, we get the kinetic term

$\displaystyle-{\rm Tr}\left[\frac{1}{n}(G\delta
H^{\prime})^{n}\right]\evaluated{}_{n=2}\supset~{}~{}$ $\displaystyle
i\int\frac{d^{4}k}{(2\pi)^{4}}\delta\tilde{\phi}(-k)\delta\tilde{\phi}(k)\int\frac{d^{3}q}{(2\pi)^{3}}\frac{|d_{0}|^{2}}{4(|d_{0}|^{2}+\phi_{0}^{2})^{5/2}}(k^{0})^{2}$
$\displaystyle=$
$\displaystyle\,J\int\differential[4]{x}\delta\phi(-\partial_{t}^{2})\delta\phi\,.$
(67)

where $\delta\tilde{\phi}$ is Fourier transformation of $\delta\phi$ and

$\displaystyle
J=\int\frac{d^{3}q}{(2\pi)^{3}}\frac{|d_{0}|^{2}}{4(|d_{0}|^{2}+\phi_{0}^{2})^{5/2}}\,.$
(68)

##  Appendix D $\theta$ term as chiral anomaly

In a special case of the Dirac model, $\theta$ can be derived as the path
integral measure due to the chiral rotation of the electron field $\psi$ [21,
19]. In this section we apply the technique to the present model and clarify
an issue claimed in Ref. [21] that $\theta$ derived in this technique can not
be applied in some circumstances.

When $B=0$ and $\epsilon_{0}=0$, the Hamiltonian is given by

$\displaystyle H({\bf\it
k})=Ak_{x}\Gamma^{1}+Ak_{y}\Gamma^{2}+Ak_{z}\Gamma^{3}+m\Gamma^{4}+\phi\Gamma^{5}\,.$
(69)

In the ($\ket*{{\rm A},\sigma}$, $\ket*{{\rm B},\sigma}$)
($\sigma=\uparrow,\downarrow$) basis, the Gamma matrix is given as

(70)

To be concrete,

$\Gamma^{a\prime}$ corresponds to $\tilde{\alpha}_{\mu}$ rotated by $U_{3}$ in
Ref. [19], except for an extra sign in $\tilde{\alpha}_{5}$. Or if one makes
another transformation,

$\displaystyle\tilde{\Gamma}^{a}$
$\displaystyle=\tilde{U}\Gamma^{a}\tilde{U}^{\dagger}\,,$ (71) $\displaystyle
U$ $\displaystyle=\frac{1}{\sqrt{2}}{\rm diag}(1+i,\,1-i,\,1-i,\,1+i)\,,$ (72)

the Dirac $\alpha$ matrices in the Dirac representation are obtained,

Here the Dirac $\alpha$ matrices corresponds to
$({\bf\it\alpha},\,\beta)=(\tilde{\Gamma}^{1},\,\tilde{\Gamma}^{2}\,,\tilde{\Gamma}^{3}\,,\tilde{\Gamma}^{4})$.
Therefore,

$\displaystyle\tilde{U}H({\bf\it k})\tilde{U}^{\dagger}$
$\displaystyle={\bf\it\alpha}\dotproduct{\bf\it k}^{\prime}+\beta
m+\phi\tilde{\Gamma}^{5}$
$\displaystyle=\beta({\bf\it\gamma}\dotproduct{\bf\it
k}^{\prime}+m+\phi\gamma_{5})\,,$ (73)

where we have changed variable ${\bf\it k}$ as $A{\bf\it k}={\bf\it
k}^{\prime}$, and introduced the Dirac gamma matrix
$\gamma^{\mu}=(\beta,\,\beta{\bf\it\alpha})$ and
$\gamma_{5}=i\gamma^{0}\gamma^{1}\gamma^{2}\gamma^{3}$. Consequently, this
Hamiltonian leads to the following action,

$\displaystyle
S=\int\differential[4]{x}\bar{\psi}[i\gamma^{\mu}(\partial_{\mu}-ieA_{\mu})-m-i\phi\gamma_{5}]\psi\,,$
(74)

where $\bar{\psi}=\psi^{\dagger}\gamma^{0}$ and we have added vector potential
$A_{\mu}$ of $U(1)_{\rm em}$. Then chiral rotation of $\psi$ gives the chiral
anomaly term in the QED [38, 39],

$\displaystyle\Delta S=-\frac{\alpha}{4\pi}\int\differential[4]{x}\Theta
F_{\mu\nu}\tilde{F}^{\mu\nu}=\frac{\alpha}{\pi}\int\differential[4]{x}\Theta{\bf\it
E}\dotproduct{\bf\it B}\,,$ (75)

where
$\tilde{F}^{\mu\nu}=\frac{1}{2}\epsilon^{\mu\nu\rho\sigma}F_{\rho\sigma}$
($\epsilon^{0123}=+1$) and

$\displaystyle\Theta=\frac{\pi}{2}[1-{\rm sgn}(m)]{\rm
sgn}(\phi)+\tan^{-1}\frac{\phi}{m}\,.$ (76)

Here we have set the domain of $\tan^{-1}$ as $[-\pi,\pi]$. This is consistent
with Ref. [21], except that the sign of $m$ and we have used single degree of
freedom for electron field $\psi$. In the model considered in Refs. [21, 19],
there are three Dirac points and they expand the Hamiltonian at each Dirac
points. That is why they have three degree of freedom for electron field. (The
sign of $M_{f}e^{i\kappa_{f}\gamma_{5}}$ in Eq. (30) of Ref. [21] is supposed
to be a plus. Comparing to the result in Ref. [19] (v1 on arXiv), their result
is a factor of two smaller. In addition, it is claimed that first term is 1/2,
instead of $\pi/2$, in the topological phase.101010Those are just typos. We
thank S. Chigusa for confirming it.) The authors of Ref. [21] claim that the
expression of $\theta$ cannot be applied when $Um_{f}/\delta t^{\prime}\gg 1$.
($Um_{f}/\delta t^{\prime}$ corresponds to $\phi/m$ in our notation.) This is
misleading because the derivation by chiral rotation in the path integral is
exact and there is no approximation. However, it should be noted that in this
calculation the cutoff of the momentum integral is taken to be infinity. To
put it more correctly, therefore, this result is reliable when $Um_{f}/\delta
t^{\prime}$ is smaller than the cutoff momentum. Although the Dirac model is
the low energy effective model, if the cutoff momentum is taken to be
infinity, then we get from Eq. (28)

$\displaystyle\theta={\rm
sgn}(\phi)\Bigl{[}\frac{\pi}{2}-\tan^{-1}\frac{m}{|\phi|}\Bigr{]}\,,$ (77)

which is the same result with Eq. (76).

##  Appendix E Another derivation of the Lagrangian for $\phi$

In the Dirac model with $B=0$ discussed in Appendix D, we can derive the
Lagrangian for $\phi$ in a way that is familiar to particle physicists. We
start with the following effective Lagrangian,

$\displaystyle{\cal L}_{\rm
eff}=\bar{\psi}[i\gamma^{\mu}(\partial_{\mu}-ieA_{\mu})-m]\psi-\frac{1}{2M_{\phi}^{2}}(\bar{\psi}\gamma_{5}\psi)^{2}\,.$
(78)

The Hubbard-Stratonovich transformation is the inverse procedure of
integrating out heavy scalar particles. In fact, above effective Lagrangian is
given by integrating out a scalar field $\phi$ that has a Yukawa interaction
with $\psi$,

$\displaystyle{\cal
L}=\bar{\psi}[i\gamma^{\mu}(\partial_{\mu}-ieA_{\mu})-m-i\phi\gamma_{5}]\psi-
M_{\phi}^{2}\phi^{2}\,.$ (79)

The last so-called four-Fermi interaction term in Eq. (78) corresponds to Eq.
(7). Taking into account the background of the AFM order $n_{{\rm
A}-}=-n_{{\rm B}-}$, $M^{2}_{\phi}$ is related to $U$ as

$\displaystyle-\frac{1}{2M_{\phi}^{2}}=-\frac{UV}{4N}\,,$ (80)

to obtain

$\displaystyle M_{\phi}^{2}=\frac{2}{(V/N)U}\,.$ (81)

This is the same result with Eq. (22). Therefore, we can use the Lagrangian
(79) as starting point, and compute the effective Lagrangian for $\phi$. As a
result, we get

$\displaystyle V_{\phi}^{(\rm
Dirac)}=-2\int\frac{d^{3}q}{(2\pi)^{3}}(\sqrt{E_{{\bf\it
q}}^{2}+\phi^{2}}-E_{{\bf\it q}})+M_{\phi}^{2}\phi^{2}\,,$ (82)

where $E_{{\bf\it q}}=\sqrt{{\bf\it q}^{2}+m^{2}}$. This agrees with Eq. (24)
but taking $A=1$ and $B=0$. As mentioned a couple of times, the momentum
integral should have a cutoff since the Dirac model is low energy effective
model.

## References

  * [1] R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440-1443 (1977) doi:10.1103/PhysRevLett.38.1440
  * [2] P.A. Zyla et al. [Particle Data Group], PTEP 2020, no.8, 083C01 (2020) doi:10.1093/ptep/ptaa104
  * [3] X. L. Qi, T. Hughes and S. C. Zhang, Phys. Rev. B 78, 195424 (2008) doi:10.1103/PhysRevB.78.195424 [arXiv:0802.3537 [cond-mat.mes-hall]].
  * [4] X. L. Qi, R. Li, J. Zang and S. C. Zhang, Science 323, 1184-1187 (2009) doi:10.1126/science.1167747 [arXiv:0811.1303 [cond-mat.mes-hall]].
  * [5] F. W. Hehl, Y. N. Obukhov, J. P. Rivera and H. Schmid, Phys. Rev. A 77, 022106 (2008) doi:10.1103/PhysRevA.77.022106 [arXiv:0707.4407 [cond-mat.other]].
  * [6] I. E. Dzyaloshinskii, Sov. Phys. JETP 10, 628 (1959); D. N. Astrov, Sov. Phys. JETP 11, 708 (1960); ibid. 13, 729 (1961).
  * [7] A. M. Essin, J. E. Moore, D. Vanderbilt, Phys. Rev. Lett. 102, 146805 (2009)
  * [8] L. Wu, M. Salehi, N. Koirala, J. Moon, S. Oh, N. P. Armitage, Science Vol. 354, Issue 6316, pp. 1124-1127 (2016)
  * [9] R. Li, J. Wang, X.-L. Qi, S.-C. Zhang, Nature Phys 6, 284 - 288 (2010) https://doi.org/10.1038/nphys1534 [arXiv:0908.1537 [cond-mat]].
  * [10] H. Ooguri and M. Oshikawa, Phys. Rev. Lett. 108, 161803 (2012) doi:10.1103/PhysRevLett.108.161803 [arXiv:1112.1414 [cond-mat.mes-hall]].
  * [11] J. Wang, B. Lian and S. C. Zhang, Phys. Rev. B 93, no.4, 045115 (2016) doi:10.1103/PhysRevB.93.045115 [arXiv:1512.00534 [cond-mat.mes-hall]].
  * [12] T .maeda, Y. Kawaguchi, Y. Tanaka, M. Sato, J. Phys. Soc. Jap. 88, no.2, 024402 (2019) doi:10.7566/JPSJ.88.024402 [arXiv:1808.10466 [cond-mat.mes-hall]].
  * [13] K. Taguchi, T. Imaeda, T. Hajiri, T. Shiraishi, Y. Tanaka, N. Kitajima and T. Naka, Phys. Rev. B 97, no.21, 214409 (2018) doi:10.1103/PhysRevB.97.214409 [arXiv:1806.01548 [cond-mat.mes-hall]].
  * [14] N. Varnava, D. Vanderbilt Phys. Rev. B 98, no.24, 245117 (2018) [ arXiv:1809.02853 [cond-mat.mtrl-sci]]
  * [15] A. Sekine and K. Nomura, [arXiv:2011.13601 [cond-mat.mes-hall]].
  * [16] D. J. E. Marsh, K. C. Fong, E. W. Lentz, L. Smejkal and M. N. Ali, Phys. Rev. Lett. 123, no.12, 121601 (2019) doi:10.1103/PhysRevLett.123.121601 [arXiv:1807.08810 [hep-ph]].
  * [17] S. Chigusa, T. Moroi and K. Nakayama, Phys. Rev. D 101, no.9, 096013 (2020) doi:10.1103/PhysRevD.101.096013 [arXiv:2001.10666 [hep-ph]].
  * [18] J. Schütte-Engel, D. J. E. Marsh, A. J. Millar, A. Sekine, F. Chadha-Day, S. Hoof, M. Ali, K. C. Fong, E. Hardy and L. Šmejkal, [arXiv:2102.05366 [hep-ph]].
  * [19] S. Chigusa, T. Moroi and K. Nakayama, [arXiv:2102.06179 [hep-ph]].
  * [20] H. Zhang, C. X. Liu, X. L. Qi, X. Dai, Z. Fang and S. C. Zhang, Nature Phys. 5, 438-442 (2009) doi:10.1038/nphys1270
  * [21] A. Sekine and K. Nomura, J. Phys. Soc. Jap. 83, no.10, 104709 (2014) doi:10.7566/JPSJ.83.104709 [arXiv:1401.4523 [cond-mat.str-el]].
  * [22] A. Sekine and K. Nomura, Phys. Rev. Lett. 116, no.9, 096401 (2016) doi:10.1103/PhysRevLett.116.096401 [arXiv:1508.04590 [cond-mat.str-el]].
  * [23] J. Zhang, D. Wang, M. Shi, T. Zhu, H. Zhang and J. Wang, Chin. Phys. Lett. 37, no.7, 077304 (2020) doi:10.1088/0256-307X/37/7/077304 [arXiv:1906.07891 [cond-mat.mes-hall]].
  * [24] J. M. Zhang, W. Ming, Z. Huang, G. B. Liu, X. Kou, Y. Fan, K. L. Wang, Y. Yao, Phys. Rev. B 88, 235131 (2013) doi:10.1103/PhysRevB.88.235131.
  * [25] Y. Li, Y. Jiang, J. Zhang, Z. Liu, Z. Yang and J. Wang, Phys. Rev. B 102, no.12, 121107 (2020) doi:10.1103/PhysRevB.102.121107 [arXiv:2001.06133 [cond-mat.mes-hall]].
  * [26] J. Li, C. Wang, Z. Zhang, B. L. Gu, W. Duan, and Y. Xu, Phys. Rev. B 100, 121103 (2019) doi:10.1103/PhysRevB.100.121103.
  * [27] H. Li et al., Phys. Rev. X 9, 041039 (2019) doi:10.1103/PhysRevX.9.041039.
  * [28] J. Li et al., Science Advances Vol. 5, no. 6, eaaw5685, doi:10.1126/sciadv.aaw5685.
  * [29] D. Zhang, M. Shi, T. Zhu, D. Xing, H. Zhang, J. Wang, Phys. Rev. Lett. 122, 206401 (2019), doi:10.1103/PhysRevLett.122.206401.
  * [30] C. Yue, Y. Xu, Z. Song, H. Weng, Y. M. Lu, C. Fang, X. Dai, Nature Physics Vol. 15, 577–581 (2019), doi:10.1038/s41567-019-0457-0.
  * [31] Y. J. Hao et al., Phys. Rev. X 9, 041038 (2019), doi:10.1103/PhysRevX.9.041038.
  * [32] S. Coh, D. Vanderbilt, A. Malashevich, I. Souza, Phys. Rev. B 83, 085108 (2011) doi:10.1103/PhysRevB.83.085108.
  * [33] L. Wu, M. Salehi, N. Koirala, J. Moon, S. Oh and N. P. Armitage, Science 354, 1124 (2016) doi:10.1126/science.aaf5541 [arXiv:1603.04317 [cond-mat.mes-hall]].
  * [34] Z. Liu and J. Wang, Phys. Rev. B 101, no.20, 205130 (2020) doi:10.1103/PhysRevB.101.205130 [arXiv:1912.07786 [cond-mat.mes-hall]].
  * [35] K. Shiozaki and S. Fujimoto, Phys. Rev. B 89, no.5, 054506 (2014) doi:10.1103/PhysRevB.89.054506 [arXiv:1310.4982 [cond-mat.supr-con]].
  * [36] B. Roy, P. Goswami and J. D. Sau, Phys. Rev. B 94, no.4, 041101 (2016) doi:10.1103/PhysRevB.94.041101 [arXiv:1507.00722 [cond-mat.mes-hall]].
  * [37] S. R. Coleman and E. J. Weinberg, Phys. Rev. D 7, 1888-1910 (1973) doi:10.1103/PhysRevD.7.1888
  * [38] K. Fujikawa, Phys. Rev. Lett. 42, 1195-1198 (1979) doi:10.1103/PhysRevLett.42.1195
  * [39] K. Fujikawa, Phys. Rev. D 21, 2848 (1980) [erratum: Phys. Rev. D 22, 1499 (1980)] doi:10.1103/PhysRevD.21.2848


# Higgs with Hadronic Signatures

W.J. Murray  on Behalf of the ATLAS and CMS collaborations  Department of
Physics, Warwick University, CV4 7AL, UK
& Rutherford Appleton Lab., Didcot, OX11 0QX, UK

###### Abstract

The decay of the Higgs boson to b quarks should be the dominant decay mode,
but it has not yet been conclusively established. The LHC run 1 results are
recalled and the current status of the LHC Run 2 studies is reviewed. The
analysis is approaching decisive sensitivity.

##  1 Introduction

The ATLAS and CMS combined Higgs coupling results[1] from LHC Run 1 included
$H\rightarrow\mbox{$b\overline{b}$}$ searches in association with vector
bosons and top quarks. When Standard Model production is assumed then
Br($H\rightarrow\mbox{$b\overline{b}$}$) is $\mu=0.70^{+0.29}_{-0.27}$, where
$\mu$ is measured over expected. If both production and decay rates are
fitted, then Br($H\rightarrow\mbox{$b\overline{b}$}$)/Br($H\rightarrow ZZ$)
has $\mu=0.20^{+0.20}_{-0.12}$.

This report discusses the results on $H\rightarrow\mbox{$b\overline{b}$}$ from
Run 2, which are currently a fragmentary picture but give indications of the
direction. Cross-sections rise as the beam energy goes from 8 TeV to 13 TeV,
especially in $t\overline{t}$$H$ where the (expected) number of Higgs events
analysed exceeds Run 1. The analyses have large backgrounds, so machine
learning techniques such as Boosted Decision Trees (BDTs)[2] are used to
identify events with enhanced signal to background.

VBF production is studied by CMS[4] and ATLAS[5], and is discussed in section
2. VH production from ATLAS is in section 3 and $t\overline{t}$$H$ from both
experiments[7, 8] is in section 4. CMS tH [3] is omitted for space reasons.
Finally section 5 contains some conclusions.

##  2 VBF Higgs production

The highest-rate Higgs production process used for
$H\rightarrow\mbox{$b\overline{b}$}$ is vector boson fusion, with a cross-
section of 3.8 pb. ATLAS and CMS have released results using different
approaches. CMS, following Run 1, and with 2.3 fb-1, used a dedicated trigger
selecting the topology of a central pair of b-jets and two forward jets, known
as tag jets. ATLAS instead, from 12.6 fb-1, use the much rarer topology where
an additional photon is produced. This simplifies the trigger and raises the
signal to background but with a large price in rate.

The CMS trigger requires 3 jets at L1, and four at the HLT. Two categories are
selected, either requiring a jet identified as a b-jet and a mass over 460 GeV
for the tag pair, or two b jets and only 200 GeV mass in the tag pair. The
signal efficiency is 6%. ATLAS used the photon as an L1 trigger, and then
require a 700 GeV tag jet pair mass at HLT with no b-tag requirement.

Both ATLAS and CMS rely on BDTs to enrich in signal purity. They use different
variables, but are sensitive to the same six areas: the dynamics of the VBF
system, the b-jet content of the event, the quark or gluon nature of the tag
jets, the suppression of jet activity between the tag jets, the H production
angles and the overall momentum of the system of four jets. The modelling in
the ATLAS case is corrected using side-bands. The resulting BDT output is used
to define several signal regions.

Figure 1: The most sensitive regions of the CMS (left) and ATLAS (right) VBF
$H\rightarrow\mbox{$b\overline{b}$}$ searches.

The BDT selections are designed not to bias the mass of the Higgs dijet,
allowing the signal rate to be extracted from a fit to a smooth background and
a signal peak, which reduces dependence upon the background simulation. In the
CMS case the background shape is constrained using the lower purity regions.
Figure 1 shows the most sensitive regions. The resulting estimates of the
signal strength are $\mu=-3.7^{+2.4}_{-2.5}$ for CMS and
$\mu=-3.9^{+2.8}_{-2.7}$ for ATLAS.

##  3 Vector-boson associated

The total cross-sections for ZH and WH are 0.88 pb and 1.37 pb respectively,
but only non-hadronic vector boson decays are used. Furthermore, the
$V$$b\overline{b}$ background in an inclusive analysis is very large, but a
boosted topology provides the sensitivity. ATLAS have released results[6]
using 13.2 fb-1 which are briefly reported here.

The analysis uses selections for 0, 1 or 2 leptons, (muons or electrons),
targeting $Z\rightarrow\nu\overline{\nu}$, $W\rightarrow l\nu$ and
$Z\rightarrow ll$ respectively. In each case two b-jets are required from the
Higgs, with possible extra jets. The backgrounds and triggers vary
considerably: two leptons provides the cleanest, but low rate, channel; one
lepton has the highest rate but large top backgrounds and the zero lepton
channel has the best overall sensitivity.

The trigger comes from either the charged leptons or
E${}_{\mathrm{T}}^{\mathrm{m}iss}$, and is most problematic in the zero lepton
channel, where an online threshold of 90 GeV was used in 2016. All channels
apply a pT threshold of 150 GeV to the Higgs candidate, with the exception of
the di-lepton, where an additional signal region of lower pT events is used. A
BDT is used to produce signal-enriched regions, but the mass is included,
implying that the background modelling is critical. Figure 2 shows the impact
on $\mu$ of various systematic errors; dominated by flavour tagging and the
normalisation of $Z$ plus heavy flavour background.

Figure 2: The VH analysis. Left is the leading constraints on nuisance
parameters of the fit, and right the combined BDTs of all the VH channels.

The resulting signal strength is $\mu=0.21\pm 0.36(stat)\pm 0.36(sys)$. A test
using the same methodology, but searching for $VZ$, yields $\mu=0.91\pm
0.17(stat)\pm 0.30(sys)$, giving confidence in the modelling.

##  4 Top-pair associated

Both CMS and ATLAS have released $t\overline{t}$$H$ analyses[7, 8] based on
about 13 fb-1 of data. The $t\overline{t}$$H$ channel cross-section rose by a
factor 3.8 from 8 TeV to 0.51 pb. The leptonic top decay provides an excellent
tag and trigger, but the complex final state means that $t\overline{t}$, and
especially $t\overline{t}$$b\overline{b}$, modelling issues dominate the
analysis.

Figure 3: The distribution of events in jet categories in ATLAS before (left)
and after (right) the fit.

Figure 3 shows the distribution of numbers of events in bins of the number of
jets and the number of b-tagged jets in ATLAS before and after the fit.
Significant mismodelling can be seen - much larger that the signal size. CMS
and ATLAS both used PowHeg v. 2 to simulate $t\overline{t}$$b\overline{b}$,
but CMS used a custom tune based on 8 TeV data and get a better pre-fit
agreement.

The fits are to BDT or Matrix Element Method (MEM) distributions.
Uncertainties on background shapes are taken from scale variations, PDFs etc.,
but both ATLAS and CMS add additional uncertainty to the rate of top pairs
with additional heavy-flavour quarks. The fit to the data constrains the
modelling so that finally there is good agreement across the distribution and
much reduced uncertainties. This is essential if the analysis is to be
sensitive.

, ,

Figure 4: The most sensitive regions of the CMS $t\overline{t}$$H$ search in
di-lepton and single-lepton (centre) channels. The plots are shown pre-fit,
and the single-lepton distribution is also shown post-fit (right).

Figure 4 shows the most signal-enriched regions of the CMS analysis, in one
and two lepton modes. The impact of the fit in reducing the uncertainties can
again be seen.

##  5 Summary and Conclusions

Table 1: The measurements of the $H\rightarrow\mbox{$b\overline{b}$}$ rate from LHC run 2 in units of the standard model expectation. Production mode | Expt. | Luminosity | Rate
---|---|---|---
$t\overline{t}$$H$ | ATLAS | 13.2 fb-1 | 2.1 ${}^{+1.0}_{-0.9}$
$t\overline{t}$$H$ | CMS | 12.9 fb-1 | $-0.19\pm 0.80$
$VH$ | ATLAS | 13.2 fb-1 | $0.21\pm 0.51$
VBF | CMS | 2.3 fb-1 | $-3.7^{+2.4}_{-2.5}$
VBF+$\gamma$ | ATLAS | 12.6 fb-1 | $-3.9^{+2.8}_{-2.7}$
Naive average | n.a. | 0.2$\pm 0.4$

The measured rates of $H\rightarrow\mbox{$b\overline{b}$}$ in the various
analyses discussed here are presented in table 1, together with an average
which symmetrises the error bars and combines all the estimates as
uncorrelated. The run 2 data do not yet show clear evidence for the
$H\rightarrow\mbox{$b\overline{b}$}$ process, but the precision is
interesting. Future measurements will have much larger samples, which will not
only improve the statistical precision but also allow better control of the
systematics and in many cases the identification of higher purity regions. The
current, murky, situation will improve dramatically.

## References

## References

  * [1] ATLAS and CMS collaborations, Journal of High Energy Physics 8, 45 (2016).
  * [2] L. Breiman, J.H. Friedman, R.A. Olshen and C.J. Stone, Classification and Regression Trees, Wadsworth, Stamford, 1984.
  * [3] CMS Collaboration, CMS-PAS-HIG-16-019, 2016.
  * [4] CMS Collaboration, CMS-PAS-HIG-16-03, 2016.
  * [5] ATLAS Collaboration, ATLAS-CONF-2016-063, 2016.
  * [6] ATLAS Collaboration, ATLAS-CONF-2016-091, 2016.
  * [7] CMS collaboration, CMS-PAS-HIG-16-038, 2016.
  * [8] ATLAS collaboration, ATLAS-CONF-2016-080, 2016.