# Collective excitations in magnetic topological insulators and axion dark
matter search

Koji Ishiwata  Institute for Theoretical Physics, Kanazawa University,
Kanazawa 920-1192, Japan  Kentaro Nomura  Department of Physics, Kyushu
University, Fukuoka 819-0395, Japan

###### Abstract

We investigate collective excitations in magnetic topological insulators (TIs)
and their impact on particle axion detection. In the three-dimensional TI
model with the Hubbard term, the effective action of magnons and amplitude
modes is formulated by dynamical susceptibility under the antiferromagnetic
and ferromagnetic states. One of the amplitude modes is identified as
“axionic” quasi-particle and its effective coupling to the electromagnetic
fields turns out to be enhanced by about four orders of magnitude larger than
the previous estimate, which may drastically change the sensitivity of the
particle axion search using “axion” in magnetic TIs.

††preprint: KANAZAWA-24-06

##  I Introduction

Particle  axion and axion-like particles are candidates for dark matter of the
universe. While many observations have been dedicated to the search, it is yet
to be found so far. Recently, the dynamical ‘axion’, as a quasi-particle in
condensed matter, has drawn attention in the context of the axion dark matter
search. The ‘axion’ couples to the electromagnetic field and Ref. [1] proposed
an axion search in the mass range of meV using the magnetic topological
insulators (TIs).111See, for examples, Refs. [2, 3, 4, 5], for recent reviews
of magnetic TIs.

The mass scale is based on an estimation of the mass of the dynamical ‘axion’
in the TIs [6]. Bi2Se3, Bi2Te3, and Sb2Se3 are the representative examples of
three-dimensional (3D) TIs. Since the TIs do not have magnetic order, it is
assumed in this estimation that magnetic impurities, such as Fe or Cr, are
doped to Bi2Se3 so that the materials have ferromagnetism (FM) or
antiferromagnetism (AFM). Then Ref. [6] proposes that the magnetism is
described by introducing the Hubbard interaction term. Namely, the total model
consists of the effective Hamiltonian for the 3D TIs and the Hubbard term. In
the literature the ‘axion’ mass is estimated to be about 1 meV. On the other
hand, it was shown in Ref. [7] that the typical mass of the ‘axion’ is of the
order of eV, independent of the topology of insulators, in the same model. The
study also indicates that it can be suppressed near the phase boundary of the
AFM order. Similar results are reported in a different model in which the spin
interaction term between the magnetic impurity and electron is introduced [8].
In the study the phase transitions of the AFM and FM are shown and the ‘axion’
mass is typically eV in both phases.

Determination of the ‘axion’ mass is crucial for the proposed axion search.
Furthermore, it is important to find a magnetic environment that is suitable
for the axion detection. In fact, the first-principles calculation shows that
Mn2Bi2Te5, which is recently synthesized experimentally [9], has a variety of
magnetic states [10]. In the literature the topological phases are tuned by
layer magnetization and the energy gap is computed in each magnetic state.
Tuning of magnetism by an electric field in magnetically doped TIs is
discussed in Ref. [11]. These studies inspire us to search for other possible
magnetic excitations that may be utilized for the axion search.

In this paper, we investigate possible collective excitations in the magnetic
TI model and study their impact on the particle axion detection. To this end,
we revisit the original model of the magnetic TIs, i.e., the effective model
for 3D TIs augmented by the Hubbard term. Starting from the partition
function, the effective action for the excitations is derived using the
dynamical susceptibility. The emergence of the AFM and FM orders is shown and
we analyze the gap and dispersion relation of eight magnetic excitations that
consist of the amplitude modes and magnons with various types.222The amplitude
mode is also called ‘Higgs amplitude mode’ and various experimental efforts
have been made to detect it, for instance, in low-dimensional magnetic
materials [12, 13, 14, 15, 16, 17] or frustrated magnetic materials [18, 19].
For the theory of the spin wave, see, for example, Refs. [20, 21, 22, 23, 24].
We identify one of the amplitude modes as ‘axion’ and one of the magnons also
can play the role of ‘axion’ in a specific case. The validity of the effective
description of the excitations is clarified and some excitations are found to
be unstable.

The paper is organized as follows. The effective action is derived from the
model in Sec. II. From the effective action the effective potential for the
order parameters of the AFM and FM is calculated in Sec. III, and the features
of the magnetic orders are discussed. In Sec. IV the effective actions of the
eight types of collective excitations are derived and the dispersion relation
and the stability of the excitations are described. Here we identify the
‘axion’ quasi-particle and investigate its properties. Based on the results,
we discuss a direction to realistic particle axion search in Sec. V. Sec. VI
is devoted to conclusion.

In our study, we adopt natural units where $\hbar=c=k_{B}=1$.

##  II The model

To discuss possible magnetic orders, we consider an effective Hamiltonian for
3D TIs [25, 6, 26] that is augmented by the Hubbard term. The Hamiltonian of
the model is given by

$\displaystyle H=H^{\rm TI}+H_{U}\,,$ (1)

where

$\displaystyle H^{\rm TI}$ $\displaystyle=\sum_{i}c^{\dagger}_{i}{\cal H}^{\rm
TI}c_{i}\,,$ (2) $\displaystyle H_{U}$
$\displaystyle=U\sum_{I=A,B}\sum_{i}n_{Ii\uparrow}n_{Ii\downarrow}\,.$ (3)

Here $U$ is a positive constant, $I=A,B$ and $i$ are the indices of the
sublattice and site, respectively.
$n_{Ii\sigma}=c^{\dagger}_{Ii\sigma}c_{Ii\sigma}$ is the number density of the
electrons at the site with spin $\sigma=\uparrow,\downarrow$ and
$c_{i}=(c_{Ai\uparrow},c_{Ai\downarrow},c_{Bi\uparrow},c_{Bi\downarrow})^{T}$.
${\cal H}^{\rm TI}$ will be given later. In our study, we consider $z$ axis as
the easy axis. By the Stratonovich-Hubbard transformation, the Hubbard term is
written by introducing two fields $\phi_{Ii}$:

$\displaystyle H_{U}\to$
$\displaystyle\frac{U}{4}\sum_{I,i}\left[\phi^{2}_{Ii}+2\phi_{Ii}(n_{Ii\uparrow}-n_{Ii\downarrow})\right]$
(4) $\displaystyle=$
$\displaystyle\frac{U}{2}\sum_{i}(\phi^{2}_{fi}+\phi^{2}_{ai})+\frac{U}{2}\sum_{i}c^{\dagger}_{i}\left[\phi_{ai}\Gamma^{5}+\phi_{fi}\Gamma^{12}\right]c_{i}\,,$
(5)

where $\Gamma^{12}$ and $\Gamma^{5}$ are 4 by 4 matrices defined as

Here $\vec{\sigma}=(\sigma^{1},\sigma^{2},\sigma^{3})$ are the Pauli matrices.
Appendix A.1 gives the list of the Gamma matrices that are used in this study.
$\phi_{Ii}$ is physically the spin moment in $z$ direction at site $i$ of
sublattice $I$ and we have introduced

$\displaystyle\phi_{ai}\equiv(\phi_{Ai}-\phi_{Bi})/2\,,$ (6)
$\displaystyle\phi_{fi}\equiv(\phi_{Ai}+\phi_{Bi})/2\,.$ (7)

When $\phi_{ai}$ and $\phi_{fi}$ are constants, then their constant values
correspond to the order parameters of the AFM and FM, respectively. The term
proportional to $\Gamma^{12}$ is considered in Refs. [27, 28, 29] in different
models. Then (Euclidean) action is given by

$\displaystyle
S=\int^{\beta}_{0}d\tau\sum_{i}\left[c^{\dagger}_{i}(\partial_{\tau}+{\cal
H})c_{i}+\frac{U}{2}(\phi^{2}_{ai}+\phi^{2}_{fi})\right]$ (8)

where $\tau=it$ is the imaginary time and

$\displaystyle{\cal H}={\cal H}^{\rm
TI}+(U/2)(\phi_{ai}\Gamma^{5}+\phi_{fi}\Gamma^{12})\,.$ (9)

$\beta=T^{-1}$ is the inverse temperature. The effective action $S_{\rm eff}$
is given by integrating out the electron field $c_{i}$ in the partition
function $Z$, i.e.,

$\displaystyle Z$ $\displaystyle=\int{\cal D}c{\cal D}c^{\dagger}{\cal
D}\phi_{a}{\cal D}\phi_{f}e^{-S}$ $\displaystyle=\int{\cal D}\phi_{a}{\cal
D}\phi_{f}e^{-S_{\rm eff}}\,,$ (10)

which leads to

$\displaystyle S_{\rm
eff}=\int^{\beta}_{0}d\tau\frac{U}{2}\sum_{i}(\phi^{2}_{ai}+\phi^{2}_{fi})-\ln\det(\partial_{\tau}+{\cal
H})\,.$ (11)

We will use $Z$ and $S_{\rm eff}$ to derive the effective potential for the
order parameters in the next section.

##  III Magnetic orders

Figure 1: Effective potential on $(\phi_{a},\phi_{f})$ plane at $T=0$. Each
panel corresponds to different values of $U$. The other parameters are
$A_{1}=A_{2}=1~{}{\rm eV}$, $B_{1}=B_{2}=-0.05~{}{\rm eV}$, and $M_{0}=0.01$
eV.

In this section we find possible global orders of the magnetism and identify
the ground state. For this purpose, we ignore the local quantum fluctuations
and compute the effective potential following the technique of Ref. [8].

###  III.1 The effective potential

The effective potential is derived from the grand potential. The grand
potential $\Omega$ is defined by the partition function $Z$:333Though we give
the formulae with finite temperature, we will focus on the zero temperature
limit in the later numerical study.

$\displaystyle\Omega=-\beta^{-1}\ln Z\,.$ (12)

The effective potential per site is then obtained by taking $\phi_{ai}$ and
$\phi_{fi}$ as constants, i.e., $\phi_{ai}=\phi_{a}$ and $\phi_{fi}=\phi_{f}$,

$\displaystyle V_{\rm eff}$ $\displaystyle\equiv-\beta^{-1}\ln e^{-S_{\rm
eff}}/N|_{(\phi_{ai},\phi_{fi})=(\phi_{a},\phi_{f})}$
$\displaystyle=\frac{U}{2}(\phi^{2}_{f}+\phi^{2}_{a})-\frac{1}{\beta
N}\sum_{h,{\bf\it k}}\ln(1+e^{-\beta(E_{h{\bf\it k}}-\mu)})\,,$ (13)

where $N$ is the number of the cite, $\mu$ is the chemical potential, $h$ is
the band index, and ${\bf\it k}$ is the wavenumber. $E_{h{\bf\it k}}$ will be
given soon. In the derivation, we have adopted the Matsubara formalism, which
gives

$\displaystyle\ln\det(\partial_{\tau}+{\cal H})=\sum_{n}\sum_{h,{\bf\it
k}}\ln(-i\omega_{n}+E_{h{\bf\it k}}-\mu)\,,$ (14)

where $\omega_{n}=(2n+1)\pi/\beta$ is the Matsubara frequency for fermion.
$E_{h{\bf\it k}}$ are the eigenvalues of ${\cal H}_{{\bf\it k}}$, which is
${\cal H}$ in the wavenumber space:

$\displaystyle{\cal H}_{\mathbf{k}}={\cal H}^{\rm
TI}_{\mathbf{k}}+d_{5}\Gamma^{5}+d_{12}\Gamma^{12}\,,$ (15)

where

$\displaystyle{\cal H}^{\rm TI}_{\mathbf{k}}$ $\displaystyle=\epsilon_{0}{\bf
1}+\sum_{a=1}^{4}d_{a}\Gamma^{a}\,,$ (16) $\displaystyle(d_{5},d_{12})$
$\displaystyle=(U/2)(\phi_{a},\phi_{f})\,.$ (17)

$\Gamma^{a}$ ($a=1$ \- 4) are the Gamma matrices given by

and we take $\phi_{a},\phi_{f}>0$ without loss of generality. $\epsilon_{0}$
is a constant and $d^{a}$ is parametrized as [25, 6, 26]

$\displaystyle(d_{1},d_{2},d_{3},d_{4})=$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}(A_{2}\sin k_{x}a_{x},\,A_{2}\sin
k_{y}a_{y},\,A_{1}\sin k_{z}a_{z},\,{\cal M})\,,$ (18)

where ${\cal M}=M_{0}-2B_{1}-4B_{2}+2B_{1}\cos k_{z}a_{z}+2B_{2}(\cos
k_{x}a_{x}+\cos k_{y}a_{y})$. In the following study, we consider a cubic
lattice and use dimensionless wavenumber for simplicity, i.e.,
$a_{x}=a_{y}=a_{z}=1$. The first and second terms of ${\cal H}_{{\bf\it k}}$
have the time-reversal invariance, which is one of the features of the TIs,
and it describes Bi2Se3 family of materials. The eigenvalues are obtained as
$(E_{1{\bf\it k}},E_{2{\bf\it k}},E_{3{\bf\it k}},E_{4{\bf\it
k}})=(\epsilon_{0}-\epsilon_{1{\bf\it k}},\epsilon_{0}-\epsilon_{2{\bf\it
k}},\epsilon_{0}+\epsilon_{1{\bf\it k}},\epsilon_{0}+\epsilon_{2{\bf\it k}})$,
where

$\displaystyle\epsilon_{1{\bf\it k}}$
$\displaystyle=\sqrt{d_{0}^{2}-d_{s}^{2}+[d_{12}+(d_{s}^{2}+d_{5}^{2})^{1/2}]^{2}}\,,$
(19) $\displaystyle\epsilon_{2{\bf\it k}}$
$\displaystyle=\sqrt{d_{0}^{2}-d_{s}^{2}+[d_{12}-(d_{s}^{2}+d_{5}^{2})^{1/2}]^{2}}\,,$
(20)

Here $d_{0}\equiv(\sum_{a=1}^{4}d_{a}^{2})^{1/2}$ and
$d_{s}\equiv(d_{3}^{2}+d_{4}^{2})^{1/2}$. Since we focus on the insulated
states, we consider that the electrons are half-filled. Namely, we take
$\epsilon_{0}-\mu\simeq 0$ and two lower bands $E_{1{\bf\it k}}$ and
$E_{2{\bf\it k}}$ are filled.

For later analysis, we derive the stationary conditions for the possible
magnetic orders. For the AFM, it is given by $\partialderivative*{V_{\rm
eff}}{\phi_{a}}|_{(\phi_{a},\phi_{f})=(\phi_{a0},0)}=0$, and the result is

$\displaystyle 1+\frac{U}{2N}\sum_{{\bf\it
k}}\\{-n_{F}(-d)+n_{F}(d)\\}\frac{1}{d}=0\,,$ (21)

where $\phi_{a0}$ is a nonzero stationary value,
$d\equiv(\sum_{a=1}^{5}d_{a}^{2})^{1/2}$ and $n_{F}$ is the Fermi-Dirac
distribution function. Similarly, $\partialderivative*{V_{\rm
eff}}{\phi_{f}}|_{(\phi_{a},\phi_{f})=(0,\phi_{f0})}=0$ gives a nonzero
$\phi_{f0}$ that satisfies

$\displaystyle 1+\frac{1}{2N\phi_{f0}}\sum_{{\bf\it k}}\Bigr{[}$
$\displaystyle\\{-n_{F}(-\epsilon_{1})+n_{F}(\epsilon_{1})\\}\frac{d_{12}+d_{s}}{\epsilon_{1}}$
$\displaystyle+$
$\displaystyle\\{-n_{F}(-\epsilon_{2})+n_{F}(\epsilon_{2})\\}\frac{d_{12}-d_{s}}{\epsilon_{2}}\Bigl{]}=0\,.$
(22)

For both the AFM and FM cases, it is easy to check that
$\phi_{a0},\phi_{f0}\to 1$ as $U\to\infty$, which corresponds to a case where
the magnetism is saturated.

Figure 2:  (a) Effective potential under $\phi_{f}=0$ as function of
$\phi_{a}$ with various values of $U$ indicated in the color bar. (b) The same
as (a) but under $\phi_{a}=0$ as function of $\phi_{f}$. For both panels, the
parameters $A_{1},A_{2},B_{1},B_{2}$, and $M_{0}$ are the same as Fig. 1.

Before computing the effective potential, we discuss conditions for the AFM
and FM orders. To this end we take $T\to 0$ limit, which leads to

$\displaystyle V_{\rm
eff}|_{T=0}=\frac{U}{2}(\phi^{2}_{f}+\phi^{2}_{a})-\frac{1}{N}\sum_{{\bf\it
k}}(\epsilon_{1{\bf\it k}}+\epsilon_{2{\bf\it k}})\,.$ (23)

Regarding the AFM, we expand the effective potential up to
$\order{\phi_{a}^{2}}$ while taking $\phi_{f}=0$ to obtain

$\displaystyle V_{\rm
eff}|_{T=0,\phi_{f}=0}=\frac{U}{2}\phi_{a}^{2}\left[1-\frac{U}{2N}\sum_{{\bf\it
k}}\frac{1}{d_{0}}\right]+\order{\phi_{a}^{4}}+{\rm const.}$ (24)

Therefore we expect to have the AFM order when

$\displaystyle 1-\frac{U}{2N}\sum_{{\bf\it k}}\frac{1}{d_{0}}<0\,.$ (25)

Similarly, one may get a condition to have the FM order as

$\displaystyle 1-\frac{U}{2N}\sum_{{\bf\it
k}}\frac{d_{1}^{2}+d_{2}^{2}}{d_{0}^{3}}<0\,.$ (26)

However, it is not an exact condition for the FM order; the FM order may arise
even when (26) is not satisfied. This can be understood by expanding the
effective potential at $\order{\phi_{f}^{4}}$:

$\displaystyle V_{\rm eff}|_{T=0,\phi_{a}=0}=$
$\displaystyle\frac{U}{2}\phi_{f}^{2}\left[1-\frac{U}{2N}\sum_{{\bf\it
k}}\frac{d_{1}^{2}+d_{2}^{2}}{d_{0}^{3}}\right]$ $\displaystyle+$
$\displaystyle\frac{U^{4}\phi_{f}^{4}}{64}\sum_{{\bf\it
k}}\frac{d_{0}^{4}-6d_{0}^{2}d_{s}^{2}+5d_{s}^{2}}{d_{0}^{7}}+\order{\phi_{f}^{6}}\,.$
(27)

The coefficient of a term proportional to $\phi_{f}^{4}$ can be negative,
which may lead to the FM order. We will confirm this numerically later. Even
though it is not an exact condition, (26) can be used to understand the
qualitative behavior of the emergence of the magnetic order. I.e., we expect
that the FM order will appear for relatively larger $d_{1}$ and $d_{2}$ than
$d_{3}$ and $d_{4}$, in other words, larger $A_{2}$ compared to $A_{1}$,
$B_{1}$, $B_{2}$, and $M_{0}$.

On the other hand, if $d_{1}=d_{2}=0$ is exactly satisfied, the effective
potential becomes

$\displaystyle V_{\rm
eff}|_{T=0,\phi_{a}=0}=\frac{U}{2}\phi_{f}(\phi_{f}-2)\,.$ (28)

Thus, the FM order appears without other conditions, and the stationary point
is given by $\phi_{f0}=1$. We will discuss this specific case later.

Both conditions (25) and (26) are satisfied in the limit $U\to\infty$. Thus,
the AFM and FM order should appear in the large $U$ limit.

###  III.2 Numerical results

Now we are ready to show the effective potential on $(\phi_{a},\phi_{f})$
plane. Fig. 1 shows the effective potential for a given set of parameters. We
take $T=0$ and the normal insulator (NI) phase is considered. The TI phase can
be calculated by taking $M_{0}<0$ with the others being unchanged and we found
the result merely changes. It is seen the potential minima appear as $U$
increases. As $U$ gets large, the AFM order forms first since the condition
for the AFM order, i.e., (25), is relatively easy to satisfy, and eventually
the FM order shows up.

Figure 3: (Left) Classification of the magnetic excitations. Excitation modes
are amplitude mode and magnon. The amplitude mode has ‘AFM-type’ and ‘FM-
type,’ defined in Eqs. (34) and (35), and the magnon has ‘$\alpha$-type’ and
‘$\beta$-type,’ defined in Eqs. (56) and (57). Images of the amplitude mode
and magnons are shown in Figs. 4 and 6, respectively. (Right) Diagram
corresponding to dynamical susceptibility $\tilde{\chi}^{\xi}_{M}(k)$ and the
(dimensionless) inverse propagator $\tilde{\Gamma}^{\xi}_{M}(k)$.
$k=(i\omega_{n},{\bf\it k})$ and $\ell=(i\omega_{\ell},{\bf\it\ell})$ are
external and loop momenta, respectively. Dashed line shows the magnetic
excitations, amplitude mode or magnon, and solid line shows electron. In the
definitions of $\tilde{\chi}_{M}^{\xi}$ and $\tilde{\Gamma}_{M}^{\xi}$, the
external lines of the excitations are excluded. Blobs are vertices and a set
of $\Gamma^{5}$, $\Gamma^{12}$, $\Gamma^{\pm}$, or $\Sigma^{\pm}$ enters
depending on the excitations. Crossed blob stands for the tree-level vertex or
tadpole.

To see the emergence of the magnetic orders, Fig. 2 shows the effective
potential as function $\phi_{a}$ or $\phi_{f}$. It should be noted that
$\phi_{a0}$ appears continuously from zero meanwhile the non-zero $\phi_{f0}$
is obtained discontinuously. For the FM order, it is seen that another vacuum
appears as the $U$ becomes larger while the coefficient of the quadratic term
is positive. This is due to the negative $\phi_{f}^{4}$ term in Eq. (27). Such
a situation never happens for the AFM since the quartic term is always
positive. In fact, the phase transition to the AFM order happens when the
value of $U$ exceeds a critical value $U_{c}^{\rm AFM}$ defined by

$\displaystyle U_{c}^{\rm AFM}=\left(\frac{1}{2N}\sum_{{\bf\it
k}}\frac{1}{d_{0}}\right)^{-1}\,.$ (29)

In the present parameter set, i.e., $A_{1}=A_{2}=1$ eV, $B_{1}=B_{2}=-0.05$
eV, and $M_{0}=0.01$ eV, we find $U_{c}^{\rm AFM}\sim 2~{}$eV.

Regarding to the ground state of the magnetic order, we have numerically
checked that the minimum of the effective potential is obtained when the
magnetic state is the AFM. This is true for other sets of parameters.
Therefore, the ground state of the system is the AFM order. On the other hand,
both minima for the AFM and FM orders get degenerate as $U\to\infty$. This can
be understood analytically. In the large $U$ limit, the stationary points
$(\phi_{a0},\phi_{f0})$ reduce to $(1,0)$ and $(0,1)$ and $\epsilon_{1{\bf\it
k}}$, $\epsilon_{2{\bf\it k}}\simeq U/2$. Thus both minima approach to

$\displaystyle V_{\rm eff}|_{T=0}\to\frac{U}{2}-\frac{1}{N}\sum_{{\bf\it
k}}2\times\frac{U}{2}=-\frac{U}{2}\,.$ (30)

Therefore, when the magnetization is saturated, i.e.,
$\phi_{a0},\phi_{f0}\simeq 1$, both AFM and FM orders are expected to be
realized as quasi-degenerate vacua.

To summarize, we have seen the emergence of the AFM and FM orders by computing
the effective potential from the grand potential and found that the AFM is the
ground state. On the other hand, the AFM and FM states are degenerate as $U$
gets larger.

##  IV The amplitude mode and magnon under magnetic order

In this section, we evaluate the gap and dispersion relation of the magnetic
excitations under the magnetic orders. As seen in the previous section, the
possible magnetic orders are the AFM and FM and the former is the ground
state. To consider the excitations under both magnetic environments, we take
$(\phi_{ai},\phi_{fi})=(\phi_{a0},0)$ or $(0,\phi_{f0})$, where $\phi_{a0}$
and $\phi_{f0}$ are constants that satisfy Eqs. (21) and (22), respectively.
We will see that there are two magnetic excitations, that are amplitude mode
and magnon, and each of them has two types. Fig. 3 (left panel) shows the
classification of the excitations. We will also discuss an excitation called
‘axion’ in the magnetic TIs.

###  IV.1 Dynamical susceptibility and inverse propagator

Figure 4: Image of the amplitude mode under the AFM (top) and FM (bottom)
orders. The AFM-type (left) and FM-type (right) is shown. The plot is an
example of the excitation with a wavenumber ${\bf\it k}=(\pi/4,0,0)$.

For the study we introduce a dynamical susceptibility, which is defined by

$\displaystyle\chi_{M}(k;O_{1},O_{2})$
$\displaystyle~{}~{}~{}~{}~{}=-\frac{1}{\beta
N}\sum_{i\omega_{\ell},{\bf\it\ell}}{\rm
Tr}[\tilde{G}_{M}(\ell)O_{1}\tilde{G}_{M}(\ell+k)O_{2}]\,,$ (31)

where $\tilde{G}_{M}$ is the Green’s function of the electron under the global
magnetic order, AFM or FM, in the wavenumber space, which will be defined
soon. The index $M=A$ and $F$ stands for the magnetic order, i.e., AFM and FM,
respectively. $O_{1}$ and $O_{2}$ are four by four matrices, such as Gamma
matrices. For the arguments of the susceptibility and the Green’s function, we
write $k$ and $\ell$ as $k=(i\omega_{n},{\bf\it k})$ and
$\ell=(i\omega_{\ell},{\bf\it\ell})$. Instead of $\chi_{M}$, we will use a
dynamical susceptibility that is symmetric with respect to $O_{1}$ and
$O_{2}$:

$\displaystyle\tilde{\chi}^{\xi}_{M}(k)\equiv\frac{1}{2}\left\\{\chi_{M}(k;O_{1},O_{2})+\chi_{M}(k;O_{2},O_{1})\right\\}\,,$
(32)

Here $\xi$ is a label for a specific set of operators $O_{1}$ and $O_{2}$,
which will be given explicitly in this section.

With the dynamical susceptibility, we will see that
$\tilde{\Gamma}^{\xi}_{M}(k)$ defined by

$\displaystyle\tilde{\Gamma}^{\xi}_{M}(k)\equiv
1-\frac{U}{4}\tilde{\chi}^{\xi}_{M}(k)\,,$ (33)

plays the role of the inverse propagator, which is normalized to be
dimensionless, of the excitations at the one-loop level. Fig. 3 (right panel)
shows diagrams that correspond to the dynamical susceptibility
$\tilde{\chi}_{M}^{\xi}$ and the inverse propagator
$\tilde{\Gamma}^{\xi}_{M}$.

###  IV.2 The amplitude mode

To discuss the amplitude mode, we take

$\displaystyle\phi_{ai}$ $\displaystyle=\phi_{a0}+\delta\phi_{ai}\,,$ (34)
$\displaystyle\phi_{fi}$ $\displaystyle=\phi_{f0}+\delta\phi_{fi}\,.$ (35)

$\delta\phi_{ai}$ and $\delta\phi_{fi}$ are two independent fluctuations of
the AFM and FM order parameters, respectively, and we call them ‘AFM-type’ and
‘FM-type.’ The example of the fluctuations are shown in Fig. 4. One can
imagine that there are four excitations by looking at the effective potential
in Fig. 1; the AFM- and FM-type amplitude modes under the AFM correspond to
the fluctuations around $(\phi_{a},\phi_{f})=(\phi_{a0},0)$ in $\phi_{a}$ and
$\phi_{f}$ directions, respectively, and a similar interpretation applies to
the fluctuations around $(\phi_{a},\phi_{f})=(0,\phi_{f0})$. From the figure
we naively expect all the amplitude modes to be stable. This intuition,
however, needs to be modified, which will be discussed below.

Figure 5: Dispersion relation for the AFM-type amplitude mode under the AFM.
(a) Contour of $\omega(k_{x},k_{y},0)$. (b) $\omega({\bf\it k})$ as function
of $k_{z}$ for various values of $k_{x},k_{y}$. We take $U=5$ eV and the other
parameters are the same as Fig. 1.

With Eqs. (34) and (35), we expand the second term of Eq. (11) as444Similar
expansion is done in Refs. [30, 5, 31] to derive ‘axion’ mass, meanwhile the
first term of Eq. (5), i.e., $U\phi^{2}_{ai}/2$ that will give rise to a term
$U\delta\phi^{2}_{ai}/2$ in Eq. (39), is missing at the beginning. When this
term is missing, then the resultant ‘mass term’ gets the opposite sign,
leading to the tachyonic instability. The authors, however, claim that they
obtain the same result as one in Ref. [6] though their calculation is
inconsistent. See also discussion below Eq. (45).

$\displaystyle-\ln\det(\partial_{\tau}+{\cal H})$ $\displaystyle=-{\rm
Tr}\ln(\partial_{\tau}+{\cal H})$ $\displaystyle=-{\rm
Tr}\ln(-G_{M}^{-1})+\sum_{n=1}^{\infty}{\rm Tr}(G_{M}\delta{\cal H})^{n}\,,$
(36)

by decomposing ${\cal H}={\cal H}^{\rm MTI}+\delta{\cal H}$ where

$\displaystyle{\cal H}^{\rm MTI}$ $\displaystyle={\cal H^{\rm
TI}}+(U/2)(\phi_{a0}\Gamma^{5}+\phi_{f0}\Gamma^{12})\,,$ (37)
$\displaystyle\delta{\cal H}$
$\displaystyle=(U/2)(\delta\phi_{ai}\Gamma^{5}+\delta\phi_{fi}\Gamma^{12})\,.$
(38)

Then the Green’s function is given by $G_{M}^{-1}=-\partial_{\tau}-{\cal
H^{\rm MTI}}$ where $\phi_{f0}=0$ and $\phi_{a0}=0$ for $M=A$ and $F$,
respectively. Using the definition of the Fourier expansion of the field and
the Green’s function given in Appendix A.2, the effective action at
$\order{\delta\phi_{ai}^{2},\delta\phi_{fi}^{2}}$ is

$\displaystyle S_{\rm eff}$
$\displaystyle\supset\int_{0}^{\beta}d\tau\sum_{i}\sum_{\xi=a,f}\frac{U}{2}\delta\phi_{\xi
i}^{2}+\frac{1}{2}{\rm Tr}(G\delta{\cal H})^{2}$ (39)
$\displaystyle=\frac{U}{2}\sum_{i\omega_{n},{\bf\it
k}}\sum_{\xi=a,f}\tilde{\delta\phi_{\xi}}(k)\tilde{\delta\phi_{\xi}}(-k)\tilde{\Gamma}^{\xi}_{M}(k)\,,$
(40)

where $\tilde{\Gamma}^{\xi}_{M}$ is given in Eq. (33) and

$\displaystyle\tilde{\chi}_{M}^{a}(k)$
$\displaystyle=\chi_{M}(k;\Gamma^{5},\Gamma^{5})\,,$ (41)
$\displaystyle\tilde{\chi}_{M}^{f}(k)$
$\displaystyle=\chi_{M}(k;\Gamma^{12},\Gamma^{12})\,.$ (42)

The first term of $\tilde{\Gamma}^{\xi}_{M}(k)$ comes from the first term in
Eq. (39), which is the tree-level contribution mass term and the second one
corresponds to the loop diagram. $\tilde{\Gamma}^{\xi}_{M}(k)$ plays a role of
an inverse propagator of the excitations at the one-loop level while it is
normalized to dimensionless in our notation. Thus
$1/\tilde{\Gamma}^{\xi}_{M}(k)$ plays a role of the connected two-point
Green’s function that is given by summing over the infinite series of the one-
loop diagram. $\tilde{\chi}^{\xi}_{M}(k)$ can be evaluated in the basis where
the Green’s functions are diagonalized. The calculation is shown in Appendix
B. From Eq. (40), the gap and dispersion relation $\omega=\omega({\bf\it k})$
of the amplitude mode is given by solving

$\displaystyle\tilde{\Gamma}^{\xi}_{M}(i\omega_{n}=\omega+i\delta,{\bf\it
k})=0\,,$ (43)

where $\delta\ll 1$.

Before the numerical study, we comment on the ‘axionic’ excitation in the
magnetic TIs. We claim that the AFM-type amplitude mode under the AFM order
corresponds to the quasi-particle ‘axion,’ regardless of the topological
nature of the insulators. Namely, the ‘axion’ corresponds to the AFM-type
amplitude mode in both TI and NI phases.555The existence of the dynamical
‘axion’ in NIs with the AFM state is already pointed out in Ref. [32]. The
dynamical ‘axion’ is also studied under the FM order in Ref. [29, 8]. First of
all, the gap of the AFM-type amplitude mode under the AFM order, obtained by
Eq. (43), corresponds to the mass of the ‘axion’. This fact is understood from
the fact that $\tilde{\Gamma}^{m}_{M}$ is the inverse propagator of the
excitation and that the pole of the propagator gives the physical mass. To
derive the effective action for the ‘axion,’ let us take the zero temperature
limit. First, the mass $m_{a}$ of ‘axion’ is given by

$\displaystyle\tilde{\Gamma}^{a}_{A}(i\omega_{n}=m_{a},{\bf\it 0})=0\,.$ (44)

It is easy to find a solution for the equation by using an analytic expression
for $\tilde{\chi}^{a}_{A}$ given Eq. (103) of Appendix C:

$\displaystyle m_{a}=U\phi_{a0}\,(=2d_{5})\,.$ (45)

Then, expanding $\tilde{\Gamma}^{a}_{A}$ at $i\omega_{n}=m_{a}$ and ${\bf\it
k}={\bf\it 0}$ and taking a continuum limit for the space, we get the action
of the ‘axion’ in the Minkowski spacetime as

$\displaystyle S_{a}=\left(\frac{U}{2}\right)^{2}J\int
d^{4}x~{}\delta\phi_{a}\left[-\partial_{t}^{2}+v_{i}^{2}\partial^{2}_{i}-m_{a}^{2}\right]\delta\phi_{a}\,,$
(46)

where an index $i$ is summed from 1 to 3 and

$\displaystyle
J=\int\frac{d^{3}\ell}{(2\pi)^{3}}\frac{d^{2}_{0}}{4d^{5}}\frac{1}{1-(d_{5}/d)^{2}}\,,$
(47)

Here $J$ and $v_{i}$ are the stiffness and velocity, respectively. See
Appendix D for the derivation of the action. The results of the mass and
stiffness given in Eqs. (45) and (47) are different from the previous results
[7]666$J$ and $Jm_{a}^{2}$ in Refs. [6, 31, 30, 5] are smaller than those in
Ref. [7] by a factor of 4, which was already pointed out in Ref. [7].

$\displaystyle J^{\rm old}$
$\displaystyle=\int\frac{d^{3}\ell}{(2\pi)^{3}}\frac{d^{2}_{0}}{4d^{5}}\,,$
(48) $\displaystyle J^{\rm old}(m_{a}^{\rm old})^{2}$
$\displaystyle=\int\frac{d^{3}\ell}{(2\pi)^{3}}\frac{d_{5}^{2}}{d^{3}}\,.$
(49)

The reason is that the previous results are given by expanding
$\tilde{\Gamma}^{a}_{A}(i\omega_{n},{\bf\it 0})$ at $i\omega_{n}=0$, not
$i\omega_{n}=m_{a}$. This is based on the assumption that the expansion around
$i\omega_{n}=0$ is a good approximation. However, the quantitative argument
regarding the validity of the expansion around $i\omega_{n}=0$ has not been
clarified explicitly in the previous works.

Compared to the previous calculation, we have two improvements. First, the
mass and stiffness given in the present work are more accurate. This is
because we do not use perturbative expansion to derive the mass and that the
stiffness is derived by the derivative of $\tilde{\Gamma}^{a}_{A}$ with
respect to $(i\omega_{n})^{2}$ at the mass. Second, it is possible to check
whether the gap (or mass) obtained at ${\bf\it k}=0$ is the minimum. Namely,
there is a possibility that a smaller gap is obtained for a nonzero ${\bf\it
k}$. In such a case, we must expand $\tilde{\Gamma}^{a}_{A}$ at the wavenumber
to derive the effective action for the lower-energy excitation. We will check
this possibility numerically below.

By the use of the effective action, we can evaluate the validity condition of
the perturbative expansion of the past works. As a result, we found that the
expansion is valid for $m_{a}\lesssim\order{{\rm eV}}$. The details are given
in Appendix D.

Figure 6: Image of the magnon under the AFM (top) and FM (bottom) orders. The
plot is an example of the excitation with a wavenumber ${\bf\it
k}=(\pi/4,0,0)$.

Fig. 5 shows the dispersion relation of the AFM-type amplitude mode under the
AFM order. We take $U=5$ and the other parameters are the same as Fig. 1. We
found that ${\bf\it k}={\bf\it 0}$ is the stable minimum. Here ‘stable’ means
that the curvature of the dispersion relation $\omega=\omega({\bf\it k})$ is
positive in all ${\bf\it k}$ directions around ${\bf\it k}={\bf\it 0}$. In
other word, $v_{i}^{2}$ are positive for all $i=1$ – 3. This can be seen
explicitly in the figure. We have checked that there are no smaller gaps, thus
the gap is surely given by $m_{a}=U\phi_{a0}\sim\order{\rm eV}$.

The situation is different for the other amplitude modes. As an example, let
us discuss the FM-type amplitude mode with the AFM order. The dynamical
susceptibility of the excitation is analytically obtained by (see Appendix C)

$\displaystyle\tilde{\chi}_{A}^{f}(\omega+i\delta,{\bf\it k})$
$\displaystyle=$
$\displaystyle\frac{2}{N}\sum_{{\bf\it\ell}}\frac{df-\sum_{a=1}^{5}d_{a}f_{a}+2(d_{1}f_{1}+d_{2}f_{2})}{df(d+f)}$
$\displaystyle\times\frac{1}{1-(\omega^{2}+i\delta)/(d+f)^{2}}\,,$ (50)

where $d_{a}$ and $f_{a}$ $(a=1$ – $4)$ depend on the wavenumber as
$d_{a}=d_{a}({\bf\it\ell})$, $f_{a}=d_{a}({\bf\it\ell}+{\bf\it k})$,
$f_{5}=d_{5}$, and $f=(\sum_{a=1}^{5}f_{a}f_{a})^{1/2}$. When ${\bf\it k}=0$,
for instance, we found that the gap is found to satisfy $\omega({\bf\it
0})>\omega_{c}$, where $\omega_{c}={\rm Min}(2d)$. In addition, the imaginary
part of $\tilde{\chi}^{f}_{A}$ is enhanced. This indicates that the amplitude
mode tends to dissipate to an electron-hole pair. One can understand the
interpretation of the imaginary part from a corresponding diagram to the
dynamical susceptibility, shown in Fig. 3. Similar behavior is found for
nonzero ${\bf\it k}$. Therefore, the excitation is not stable and dissipates
even if it is created. We have also found no stable configuration for both
AFM- and FM-type amplitude modes under the FM order. In addition, the
expansion for $\omega$ is inappropriate since $\omega/(d+f)$ can be larger
than unity. Therefore, the effective description of the form of Eq. (46) is
invalid for those excitations.

###  IV.3 The magnon

To describe the magnon, we need to rewrite the Hubbard term in an isotropic
form. Namely, the Stratonovich-Hubbard transformation gives

$\displaystyle H_{U}\to$
$\displaystyle\frac{U}{4}\sum_{I,i}\left[\phi^{2}_{I}+2\phi_{I}c^{\dagger}_{Ii}\vec{n}_{Ii}\dotproduct\vec{\sigma}c_{Ii}\right]\,,$
(51)

where $\vec{n}_{Ii}=(n_{Ii}^{x},n_{Ii}^{y},n_{Ii}^{z})$ that satisfies
$|\vec{n}_{Ii}|=1$. Under the AFM, $\phi_{Ai}=-\phi_{Bi}=\phi_{a0}/2$, while
$\phi_{Ai}=\phi_{Bi}=\phi_{f0}/2$ for the FM. Recalling that we are interested
in the case where $z$ direction is the easy axis, the fluctuation of
$n_{Ii}^{x},n_{Ii}^{y}$ corresponds to the magnon. For later analysis we
introduce $\vec{n}_{\alpha i}$ and $\vec{n}_{\beta i}$ as

$\displaystyle\vec{n}_{\alpha i}\equiv(\vec{n}_{Ai}+\vec{n}_{Bi})/2\,,$ (52)
$\displaystyle\vec{n}_{\beta i}\equiv(\vec{n}_{Ai}-\vec{n}_{Bi})/2\,.$ (53)

The labels $\alpha$ and $\beta$ stand for two different types of magnon
excitations, which we call ‘$\alpha$-type’ and ‘$\beta$-type,’ respectively.
We will formulate each excitation under both the AFM and FM orders below.

First, let us discuss the magnon under the AFM order. $\vec{n}_{A}$ and
$\vec{n}_{B}$ are parallel for the $\alpha$-type, meanwhile rotation phases of
$\vec{n}_{A}$ and $\vec{n}_{B}$ coincide for the $\beta$-type. See Fig. 6 (top
panel) for an example of the $\alpha$\- and $\beta$-types under the AFM. The
interaction Hamiltonian gives

$\displaystyle
H_{U}\supset\frac{U\phi_{a0}}{2}\sum_{i}c^{\dagger}_{i}\Bigl{[}$
$\displaystyle\Gamma^{5}n_{\alpha i}^{z}+\Gamma^{1}n_{\alpha
i}^{x}+\Gamma^{2}n_{\alpha i}^{y}$ $\displaystyle+$
$\displaystyle\Gamma^{12}n_{\beta i}^{z}+\Gamma^{25}n_{\beta
i}^{x}+\Gamma^{51}n_{\beta i}^{y}\Bigr{]}c_{i}\,,$ (54)

where

Assuming that $n^{x}_{\alpha i},n^{y}_{\alpha i},n^{x}_{\beta i},n^{y}_{\beta
i},\ll 1$ the Hamiltonian reduces to ${\cal H}={\cal H}^{\rm MTI}+\delta{\cal
H}$ where

$\displaystyle\delta{\cal H}=\frac{U\phi_{a0}}{2}\Bigl{[}$
$\displaystyle-\frac{1}{2}\Gamma^{5}\\{(n^{x}_{\alpha i})^{2}+(n^{y}_{\alpha
i})^{2}+(n^{x}_{\beta i})^{2}+(n^{y}_{\beta i})^{2}\\}$
$\displaystyle+\Gamma^{1}n^{x}_{\alpha i}+\Gamma^{2}n^{y}_{\alpha i}$
$\displaystyle+4\Gamma^{12}(n^{x}_{\alpha i}n^{x}_{\beta i}+n^{y}_{\alpha
i}n^{y}_{\beta i})$ $\displaystyle+\Gamma^{25}n^{x}_{\beta
i}+\Gamma^{51}n^{y}_{\beta i}\Bigr{]}\,.$ (55)

The term proportional to $\Gamma^{12}$ does not contribute to AFM order and we
will ignore it hereafter. Now it is convenient to introduce

$\displaystyle n^{\pm}_{\alpha i}\equiv n^{x}_{\alpha i}\pm in^{y}_{\alpha
i}\,,$ (56) $\displaystyle n^{\pm}_{\beta i}\equiv n^{x}_{\beta i}\pm
in^{y}_{\beta i}\,.$ (57)

Using them, $\delta{\cal H}$ is written as

$\displaystyle\delta{\cal H}=\frac{U\phi_{a0}}{2}\Bigl{[}$
$\displaystyle-\frac{1}{2}\Gamma^{5}(n^{+}_{\alpha i}n^{-}_{\alpha
i}+n^{+}_{\beta i}n^{-}_{\beta i})$ $\displaystyle+\Gamma^{+}n^{-}_{\alpha
i}+\Gamma^{-}n^{+}_{\alpha i}+\Sigma^{+}n^{-}_{\beta i}+\Sigma^{-}n^{+}_{\beta
i}\Bigl{]}\,,$ (58)

where

Now we derive the effective action in a similar way as in Sec. IV.2. Computing
the terms up to $\order{(n^{\pm}_{\alpha i})^{2},(n^{\pm}_{\beta i})^{2}}$,
the effective action for the excitations is

$\displaystyle S_{\rm
eff}\supset\frac{U\phi_{a0}^{2}}{2}\sum_{i\omega_{n},{\bf\it
k}}\sum_{\xi=\alpha,\beta}\tilde{n}_{\xi}^{+}(k)\tilde{n}_{\xi}^{-}(-k)\tilde{\Gamma}_{A}^{\xi}(k)\,,$
(59)

where $\tilde{\Gamma}^{\xi}_{M}$ is given in Eq. (33) and
$\tilde{\chi}_{A}^{\xi}(k)$ ($\xi=\alpha,\beta$) are given by

$\displaystyle\tilde{\chi}_{A}^{\alpha}(k)$
$\displaystyle=\frac{1}{2}\\{\chi_{A}(k;\Gamma^{+},\Gamma^{-})+\chi_{A}(k;\Gamma^{-},\Gamma^{+})\\}\,,$
(60) $\displaystyle\tilde{\chi}_{A}^{\beta}(k)$
$\displaystyle=\frac{1}{2}\\{\chi_{A}(k;\Sigma^{+},\Sigma^{-})+\chi_{A}(k;\Sigma^{-},\Sigma^{+})\\}\,.$
(61)

Here the first term of $\tilde{\Gamma}_{A}^{\xi}$ is given by the term
proportional to $\Gamma^{5}$ in Eq. (58). Diagrammatically it is a tadpole and
we have found that it can be rewritten by using Eq. (21). As a result, it
turns out to have the same structure as the amplitude modes.

Figure 7: Dispersion relation for the $\alpha$-type magnon under the AFM. (a)
Contour of $\omega(k_{x},k_{y},0)$. (b) $\omega({\bf\it k})$ as function of
$k_{z}$ for various values of $k_{x},k_{y}$. We take $U=5$ eV and the other
parameters are the same as Fig. 1.

We repeat the argument for the FM order. In this case, we take
$\phi_{Ai}=\phi_{Bi}=\phi_{f0}/2$ to give the Hamiltonian $H_{U}$

$\displaystyle
H_{U}\supset\frac{U\phi_{f0}}{2}\sum_{i}c^{\dagger}_{i}\Bigl{[}$
$\displaystyle\Gamma^{12}n_{\alpha i}^{z}+\Gamma^{25}n_{\alpha
i}^{x}+\Gamma^{51}n_{\alpha i}^{y}$ $\displaystyle+$
$\displaystyle\Gamma^{5}n_{\beta i}^{z}+\Gamma^{1}n_{\beta
i}^{x}+\Gamma^{2}n_{\beta i}^{y}\Bigr{]}c_{i}\,.$ (62)

As seen, the Hamiltonian is given by a replacement
$\Gamma^{5}\leftrightarrow\Gamma^{12}$,
$\Gamma^{1}\leftrightarrow\Gamma^{25}$,
$\Gamma^{2}\leftrightarrow\Gamma^{51}$, and $\phi_{a0}\to\phi_{f0}$ in Eq.
(54). Then the perturbative term $\delta{\cal H}$ in the Hamiltonian ${\cal
H}$ becomes,

$\displaystyle\delta{\cal H}=\frac{U\phi_{f0}}{2}\Bigl{[}$
$\displaystyle-\frac{1}{2}\Gamma^{12}(n^{+}_{\alpha i}n^{-}_{\alpha
i}+n^{+}_{\beta i}n^{-}_{\beta i})$ $\displaystyle+\Sigma^{+}n^{-}_{\alpha
i}+\Sigma^{-}n^{+}_{\alpha i}+\Gamma^{+}n^{-}_{\beta i}+\Gamma^{-}n^{+}_{\beta
i}\Bigl{]}\,.$ (63)

Here we have omitted a term proportional to $\Gamma^{5}$ that is ineffective
for the FM order. The effective action is then obtained by

$\displaystyle S_{\rm
eff}\supset\frac{U\phi_{f0}^{2}}{2}\sum_{i\omega_{n},{\bf\it
k}}\sum_{\xi=\alpha,\beta}\tilde{n}_{\xi}^{+}(k)\tilde{n}_{\xi}^{-}(-k)\tilde{\Gamma}^{\xi}_{F}(k)\,,$
(64)

and

$\displaystyle\tilde{\chi}_{F}^{\alpha}(k)$
$\displaystyle=\frac{1}{2}\\{\chi_{F}(k;\Sigma^{+},\Sigma^{-})+\chi_{F}(k;\Sigma^{-},\Sigma^{+})\\}\,,$
(65) $\displaystyle\tilde{\chi}_{F}^{\beta}(k)$
$\displaystyle=\frac{1}{2}\\{\chi_{F}(k;\Gamma^{+},\Gamma^{-})+\chi_{F}(k;\Gamma^{-},\Gamma^{+})\\}\,.$
(66)

The first term of $\tilde{\Gamma}^{\xi}_{F}$ comes from the term proportional
to $\Gamma^{12}$. It corresponds to the tadpole diagram and it can be
rewritten by using Eq. (22). Thus the structure is the same as
$\tilde{\Gamma}_{A}^{\xi}$. With the FM order, the $\alpha$-type and
$\beta$-type of magnon corresponds to $n_{\alpha i}^{\pm}$ and $n_{\beta
i}^{\pm}$, respectively, which are shown in Fig. 6 (bottom).

From the expressions of the effective action regarding magnons, i.e., Eqs.
(59) and (64), the gap and dispersion relation of the magnon is given by

$\displaystyle\tilde{\Gamma}^{\xi}_{M}(i\omega_{n}=\omega+i\delta,{\bf\it
k})=0~{}~{}~{}(\xi=\alpha,\beta,~{}~{}M=A,F)\,.$ (67)

Now we show the numerical results of the dispersion relation of the magnons.
Fig. 7 shows a plot of the dispersion relation regarding the $\alpha$-type
magnon under the AFM order. We have found the dispersion relation
$\omega({\bf\it k})$ is minimized at ${\bf\it k}={\bf\it 0}$. On the other
hand, there are additional quasi-degenerate values of $\omega({\bf\it k})$ at
${\bf\it k}=(\pi,0,0),(0,\pi,0),(\pi,\pi,0)$. The degeneracy is about 1%,
compared to $\omega({\bf\it 0})$ and the curvatures of the wavenumbers are
positive at the ${\bf\it k}$s. Therefore we expect four different stable
configurations of the magnon. The excitation energy is found to be around
eV.777When we derive the effective action for those configurations of magnon,
we need to expand $\tilde{\Gamma}^{\alpha}_{A}(\omega,{\bf\it k})$ around each
wavenumber, i.e., ${\bf\it k}=(\pi,0,0),(0,\pi,0),(0,\pi,\pi)$.

Figure 8: Same as Fig. 7 but for $\beta$-type magnon under the AFM. (a)
Contour of $\omega(k_{x},k_{y},\pi)$. (b) $\omega({\bf\it k})$ as function of
$k_{z}$ for various values of $k_{x},k_{y}$.

Fig. 8 shows the dispersion relation for the $\beta$-type magnon under the AFM
order. For the excitation similar result to the $\alpha$-type is obtained but
with a different configurations; we have found $\omega({\bf\it k})$ has the
minimum at ${\bf\it k}=(\pi,\pi,\pi)$, and cross values are found at ${\bf\it
k}=(0,0,\pi),(\pi,0,\pi),(0,\pi,\pi)$ with a degeneracy of around 1%. The
energy scale turns out to be the same as the $\alpha$-type one, i.e.,
$\order{\rm eV}$.

It is worth discussing the origin of the gap of the magnons. Let us consider
the $\alpha$-type magnon under the AFM at zero temperature. Assuming that the
expansion around $\omega=0$ is valid, we get the stiffness and mass of the
magnon as (see Appendix D for the derivation)

$\displaystyle J^{\prime}_{\alpha}$
$\displaystyle=\int\frac{d^{3}\ell}{(2\pi)^{3}}\frac{2d^{2}-d_{1}^{2}-d_{2}^{2}}{8d^{5}}\,,$
(68) $\displaystyle J^{\prime}_{\alpha}m_{\alpha}^{\prime\,2}$
$\displaystyle=\int\frac{d^{3}\ell}{(2\pi)^{3}}\frac{d_{1}^{2}+d_{2}^{2}}{2d^{3}}\,.$
(69)

It is clear that magnon is gapped at the zero temperature. From the
expression, we see that gap is typically $\order{\rm eV}$, which is consistent
with the numerical result shown in Fig. 7. The typical scale is the same as
the AFM-type amplitude mode under the AFM. However, the parameter dependence
is different; in the $U\to\infty$ limit, $J^{\prime}_{\alpha}\propto{\rm
eV}^{3}/U^{3}$ and $J^{\prime}_{\alpha}m_{\alpha}^{\prime\,2}\sim
A_{2}^{2}~{}{\rm eV}^{2}/U^{3}$ leads to $m^{\prime}_{\alpha}\sim A_{2}$.
Therefore, the gap remains to be $\order{A_{2}}$ even if $U$ is much larger.
We have also confirmed this numerically.

Based on the above argument, we expect that the gap of the magnon becomes zero
when $d_{1},d_{2}\to 0$. This can be confirmed without the expansion with
respect to $\omega$; considering the zero wavenumber and taking the limit at
the zero temperature, Eq. (67) with $M=A$ and $\xi=\alpha$ gives rise to

$\displaystyle
1-\frac{U}{2N}\sum_{{\bf\it\ell}}\frac{1}{d}\frac{1}{1-\omega^{2}/(4d^{2})}=0\,.$
(70)

Recalling Eq. (21), it is easy to find $\omega=0$ is the solution to satisfy
the gap equation. Therefore nonzero $d_{1}$ and $d_{2}$ (or $A_{2}$) are the
origin of the $\alpha$-type magnon gap under the AFM. This can be understood
by considering the symmetry breaking of the system. When $d_{1}=d_{2}=0$,
SO(3) symmetry in the $\vec{\phi}_{Ai}-\vec{\phi}_{Bi}\equiv 2\vec{\phi}_{ai}$
space is restored. Under the AFM order, the nonzero
$\vec{\phi}_{ai}=(0,0,\phi_{a0})$ breaks the SO(3) to SO(2). Consequently two
massless Nambu-Goldstone bosons appear, which correspond to $n^{\pm}_{\alpha
i}$. See Appendix E for more detail.

Since we already have identified Nambu-Goldstone bosons, another excitation,
$\beta$-type magnon under the AFM order, should be massive even if
$d_{1}=d_{2}=0$. In fact, by expanding
$\tilde{\Gamma}_{A}^{\beta}(\omega,{\bf\it 0})$ around $\omega=0$, the
stiffness and mass are given by

$\displaystyle J^{\prime}_{\beta}$
$\displaystyle=\int\frac{d^{3}\ell}{(2\pi)^{3}}\frac{2d_{5}^{2}+d_{1}^{2}+d_{2}^{2}}{8d^{5}}\,,$
(71) $\displaystyle J^{\prime}_{\beta}m_{\beta}^{\prime\,2}$
$\displaystyle=\int\frac{d^{3}\ell}{(2\pi)^{3}}\frac{2d_{0}^{2}-d_{1}^{2}-d_{2}^{2}}{2d^{3}}\,.$
(72)

Thus a finite mass (or gap) should appear even if $d_{1},d_{2}\to 0$.

The results for $\alpha$-type and $\beta$-type magnons under the FM state are
similarly computed. It is found that qualitative features of the $\alpha$-type
and $\beta$-type ones are the same as the $\beta$\- and $\alpha$-types under
the AFM order, respectively. Typical energy scale of the excitations are a bit
smaller but they are $\order{\rm eV}$. See their plots in Appendix F.

To summarize this section, we have formulated the effective actions of the the
amplitude modes and magnons under both the AFM and FM orders. Each mode has
two types; the amplitude mode has AFM- and FM-types and the magnon has the
$\alpha$\- and $\beta$-types. From the effective action, we have computed the
dispersion relations for all excitations. Among them, the AFM-type amplitude
mode is identified as the ‘axionic’ quasi-particle. We have found that the
mass of ‘axion’ coincides with the gap and it is derived more precisely
compared to the past works. The typical energy scale turns out to be eV and
the other excitations have the same scale.

##  V Implication to particle axion search

As seen in Sec. IV.2, the ‘axionic’ excitation corresponds to the AFM-type
amplitude mode (up to normalization of the field) under the AFM, and it is
shown its mass (or gap) $m_{a}$ is given by $U\phi_{a0}$. On the other hand,
the energy bands of the electrons of the magnetic TIs are computed by the
first-principles calculation. For ${\rm Mn}_{2}{\rm Bi}_{2}{\rm Te}_{5}$, as
an example, Ref. [33] shows that the gap of the electrons under the AFM order
is $\order{10\,\mathchar 45\relax\,10^{2}}$ meV. With this knowledge, we
discuss the possibility of the particle axion detection by using the
interaction of the ‘axion’ with the electromagnetic fields.

The ‘axion’ field $\theta$ has an ${\bf\it E}\dotproduct{\bf\it B}$ coupling
in the Lagrangian [6]:

$\displaystyle{\cal L}_{EB}=-\frac{\alpha}{4\pi}\theta
F_{\mu\nu}\tilde{F}^{\mu\nu}=\frac{\alpha}{\pi}\theta{\bf\it
E}\dotproduct{\bf\it B}\,,$ (73)

where $\alpha$ is the fine-structure constant and $F_{\mu\nu}$ is the field
strength of the electromagnetism.888$\mu,\nu,\rho,\sigma$ are the Lorentz
indices that are 0, 1, 2, 3.  $\tilde{F}^{\mu\nu}$ is the dual of
$F_{\mu\nu}$, defined as
$\tilde{F}^{\mu\nu}=\frac{1}{2}\epsilon^{\mu\nu\rho\sigma}F_{\rho\sigma}$
($\epsilon^{\mu\nu\rho\sigma}$ is the Levi-Civita symbol with
$\epsilon^{0123}=+1$) and $\theta$ is related to the model parameters of the
3D TIs as [6]

$\displaystyle\theta=\frac{1}{4\pi}\int
d^{3}k\frac{2d+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}\,,$
(74)

with $i,j,k,l$ being 1, 2, 3, and 5.

Figure 9: Effective decay constant $m_{a}/c^{\rm eff}_{a}$ as function of gap
$m_{a}$ of the amplitude mode. $M_{0}\,[{\rm meV}]=\pm 10$ (left) and $\pm
10^{2}$ (right) are taken, where positive (negative) $M_{0}$ corresponds to
the NI (TI) phase. The other parameters ($A_{1}$, $A_{2}$, $B_{1}$ and
$B_{2}$) are the same as Fig. 1.

Recalling that the effective action for the AFM-type amplitude mode under the
AFM is given in Eq. (46), we define a canonically normalized field
$\delta\hat{\phi}_{a}$:

$\displaystyle\frac{U}{2}\sqrt{J}\delta\phi_{a}\equiv\frac{1}{\sqrt{2}}\delta\hat{\phi}_{a}\,.$
(75)

Therefore, the interaction term can be written in a form,

$\displaystyle{\cal L}_{EB}\supset-\frac{\alpha}{4\pi}c_{a}^{\rm
eff}\frac{\delta\hat{\phi}_{a}}{m_{a}}F_{\mu\nu}\tilde{F}^{\mu\nu}\,,$ (76)

where

$\displaystyle c^{\rm
eff}_{a}=\partialderivative{\theta}{\ln\phi_{a}}\sqrt{\frac{2}{J}}~{}\biggr{|}_{\phi_{a}=\phi_{a0}}\,.$
(77)

Therefore, $m_{a}/c^{\rm eff}_{a}$ works as an effective“decay constant” for
axion.999Be aware that there is an additional factor of $\alpha/(4\pi)$ in the
$\delta\hat{\phi}_{a}F\tilde{F}$ coupling. In terms of a variable
$g^{-1}=\partialderivative*{\theta}{m_{5}}$ (or
$\partialderivative*{\theta}{d_{5}}$ in our notation) in the literature [6,
34], it is given by

$\displaystyle g^{2}J=\frac{1}{2}\left(\frac{m_{a}}{c^{\rm
eff}_{a}}\right)^{2}\,.$ (78)

Or it is equivalent to $f_{Q}$ used in Ref. [1], i.e.,

$\displaystyle f_{Q}=\frac{m_{a}}{c^{\rm eff}_{a}}\,.$ (79)

We plot the effective decay constant of the amplitude mode as function of
$m_{a}$ in Fig. 9. Here we have used the fact that the AFM order parameter
$\phi_{a0}$ becomes nonzero continuously as $U$ increases shown in Sec. III.
Therefore $m_{a}$ can take any value from zero to $U$ if the value of $U$ is
properly chosen. In the plot the results for the NI and TI phases are given by
taking $M_{0}=\pm 10,\,\pm 10^{2}$ meV for comparison.101010The enhancement
point corresponds to a root of $c^{\rm eff}_{a}$. While the results of the NI
and TI are similar for $M_{0}=\pm 10$ meV, those with $M_{0}=\pm 10^{2}$ meV
are different. This quantitative difference comes from the 3D TI model, which
will be discussed below. To put it simply, we found that the effective
coupling is $\order{10\,\mathchar 45\relax\,10^{2}}$ meV for $m_{a}\lesssim$
100 meV.

Ref. [1], which proposes the particle axion detection using ‘axion’ in the
magnetic TIs, claims that the effective decay constant $f_{Q}$ is about $190$
eV, based on Ref. [6]. However, this is larger than our result by
$\order{10^{3}\,\mathchar 45\relax\,10^{4}}$. Let us see this more closely. In
Ref. [6], they use a variable $b^{2}\equiv\alpha^{2}B_{0}^{2}/2\pi^{2}\epsilon
g^{2}J$,111111We have changed the expression in the literature to one in the
natural unit. where $B_{0}=2$ T and $\epsilon=100$, and estimate $b$ as
$b=0.5$ meV. Using this relation, $f_{Q}=\sqrt{2g^{2}J}$ can be calculated as
$180$ eV, which is close to the value given in Ref. [1]. Namely, the estimate
of $b=0.5$ meV in Ref. [6] indicates $\sqrt{g^{2}J}\simeq 130$ eV. One may
think that this discrepancy comes from another choice of the input parameters.
However, such a possibility is unlikely from an analytic estimation below.

Let us focus on the case where $d_{5}\,(=m_{a}/2)=\order{1~{}{\rm meV}}$ that
is considered in Ref. [6]. To take an analytic approach, we consider the Dirac
model instead of the 3D TI model. Since this scale is much smaller than the
energy scale of the electrons, $\theta$ in the Dirac model is approximately
given by [7]

$\displaystyle\theta\sim\frac{\pi}{2}[1-{\rm sgn}(M_{0})]{\rm
sgn}(d_{5})+\frac{d_{5}}{M_{0}}\,,$ (80)

which leads to

$\displaystyle g^{-1}\sim 1/M_{0}\,.$ (81)

In addition, $J$ is estimated as

$\displaystyle
J\simeq\int\frac{d^{3}\ell}{(2\pi)^{3}}\frac{1}{4d_{0}^{3}}\sim\order{0.1\,\mathchar
45\relax\,1}\,.$ (82)

The estimation above is also checked numerically. Therefore, we get

$\displaystyle\sqrt{g^{2}J}=\frac{m_{a}}{\sqrt{2}c^{\rm
eff}_{a}}\sim|M_{0}|\,,$ (83)

for $d_{5}\sim{\rm meV}$. This estimation is roughly consistent with the
results shown in Fig. 9.121212The result of the TI case with $M_{0}=-100$ meV
is much smaller than the estimation. This is due to the difference in the
model used in the calculation. We have computed the effective decay constant
with the Dirac model and checked this point. See Appendix F for the plots of
the results. We note that $|M_{0}|$ does not exceed the typical energy scale
of the electron, i.e., $\order{\rm eV}$. This is the reason that the value
$\sqrt{g^{2}J}\sim 10^{2}\,{\rm eV}$ is unlikely.131313Ref. [31] indicates
$f_{Q}\sim$ 30 eV, which is still much larger than the present estimation.

With the estimate (83) and $M_{0}=10$ meV, the signal strength of the axion
detection proposal given in Ref. [1] is enhanced by about four orders of
magnitude, which drastically improves the sensitivity of the detection if the
technique claimed in the literature is feasible.

To close this section, we point out a possible interaction of the magnon with
the electromagnetic fields. We define a canonically normalized field for the
$\alpha$-type magnon under the AFM order as

$\displaystyle\frac{U\phi_{a0}}{2}\sqrt{J_{\alpha}}n_{\alpha}^{\pm}\equiv\hat{n}^{\pm}_{\alpha}\,.$
(84)

Then the interaction with the electromagnetic field is

$\displaystyle{\cal L}_{EB}\supset-\frac{\alpha}{4\pi}c_{\alpha}^{\rm
eff}\frac{\hat{n}^{+}_{\alpha}\hat{n}^{-}_{\alpha}}{m^{2}_{a}}F_{\mu\nu}\tilde{F}^{\mu\nu}\,,$
(85)

where

$\displaystyle c^{\rm
eff}_{\alpha}=\partialderivative{\theta}{\ln\phi_{a}}\frac{2}{J_{\alpha}}~{}\biggr{|}_{\phi_{a}=\phi_{a0}}\,.$
(86)

The interaction indicates that two magnons are excited under the
electromagnetic fields, which can be another interesting signal for the
particle axion detection. However, we have found that the magnon tends to
dissipate when its mass is below $\order{{\rm eV}}$. Therefore, it would be
challenging to detect such excitations. Instead, there may be a possibility of
finding a relic of the dissipation of the two magnons, which is induced by the
axion. We leave it for future investigation.

##  VI Conclusion

We have investigated possible excitations in the effective model of 3D
magnetic topological insulators and discussed their impact on the particle
axion detection. The model consists of the 3D effective model of TIs and the
Hubbard term. In the current study we focus on the zero temperature case.
Computing the effective potential from the grand potential, the AFM and FM
orders are found. Both orders are triggered by a large value of $U$, which is
a dimensionful coefficient of the Hubbard term. It turns out that the AFM
state is lower energy than the FM, i.e., the AFM is the global minimum of the
system. The order parameter $\phi_{a0}$ of the AFM appears continuously from
zero, meanwhile that of the FM state turns out to emerge discontinuously.

Under the magnetic orders, we have derived the effective action for possible
magnetic excitations. The magnetic excitations are classified by amplitude
mode and magnon. Furthermore, the amplitude mode and magnon have two types:
AFM/FM-types and $\alpha$/$\beta$-types. We consider those four excitations
under both the AFM and FM states, that are in total eight magnetic
excitations. We found that the effective actions for all excitations are given
by the inverse propagator that is composed of the dynamical susceptibility.
The formulation provides not only the dispersion relation and the
correspondence between the gap and mass of the excitations but also a
criterion for the stability and validity of the effective description of the
excitations.

With the formalism, we have found the AFM-type amplitude mode under the AFM is
stable and the other amplitude modes tend to dissipate. The gap (or mass) of
the AFM-type amplitude mode under the AFM is given by $m_{a}=U\phi_{a0}$.
Namely, the gap is typically $U\,(\sim$eV) for a saturated magnetization
$\phi_{a0}\sim 1$. On the other hand, it can be suppressed when $\phi_{a0}\ll
1$.

The four magnons turn out to be stable and their typical mass scale is
$\order{\rm eV}$. For stable magnons, we have found there are several quasi-
degenerate configurations. Taking the $\alpha$-type magnon under the AFM as an
example, the states with wavenumber $(\pi,0,0)$, $(0,\pi,0)$, and
$(\pi,\pi,0)$ are found to be stable in addition to one with $(0,0,0)$. On the
other hand, we discovered the magnon tends to dissipate for $\phi_{a0}\lesssim
1$.

We are especially interested in the AFM-type amplitude modes because they
relate to ‘axionic’ quasi-particle and particle axion detection. First of all,
we have identified ‘axion’ as the AFM-type amplitude mode under the AFM.
Besides, we have determined its mass and the effective coupling to the
electromagnetic fields more accurately than the past works. Since the mass
ranges from zero to $U$, it is possible to consider a situation where the mass
is less than $\order{10^{2}\,{\rm meV}}$. In the scale, the effective decay
constant of ‘axion,’ determined by the effective coupling and the mass, turns
out to be $\order{10\,\mathchar 45\relax\,10^{2}}$ meV. This value is about
three to four orders of magnitude smaller than the previous estimate, which
has a significant impact on the proposal of the particle axion detection using
‘axion’ in the magnetic TIs. We also point out that the nature of ‘axion’ is
insensitive to the topology of the magnetic insulators, i.e., topological
phase or normal phase.

The suppression of the effective decay constant of ‘axion’ is good news for
the partcile axion detection. In addition, the topological nature of the
materials is unnecessary to utilize ‘axion’ for the search. The $\alpha$-type
magnon under the AFM state is another possibility to access to particle axion.
Therefore, it is worth pursuing materials from broader candidates by studying
their magnetic states and collective excitations that are suitable for the
axion search. On the theoretical side, extensions of the model and simulations
by the first-principles’ calculation would be the next step to further
investigation of the magnetic excitations. We will leave it for future study.

###### Acknowledgements.

We are grateful to Makoto Naka for valuable discussions in the early stage of
this project. This work was supported by JSPS KAKENHI Grant Numbers
JP18H05542, JP20H01894, and JSPS Core-to-Core Program Grant No.
JPJSCCA20200002 (KI), and JST CREST, Grant No. JPMJCR18T2 (KN).

##  Appendix A Notation

We summarize the Gamma matrices and the notation for the wavefunction and
Green’s function we use in this paper.

###  A.1 Gamma matrices

In this study we mainly use the Gamma matrices in the so-called sublattice
basis. In the basis the Gamma matrices $\Gamma^{a}$ ($a=1,\cdots,5$) are
defined as

$\Gamma^{5}$ can also be written by
$\Gamma^{5}=-\Gamma^{1}\Gamma^{2}\Gamma^{3}\Gamma^{4}$. In addition, we define

$\displaystyle\Gamma^{ab}=\frac{1}{2i}\,[\Gamma^{a},\Gamma^{b}]\,.$ (87)

The explicit form of $\Gamma^{25}$, $\Gamma^{51}$, and $\Gamma^{12}$ are The
operators $\Gamma^{a}$ ($a=1,2,5$) and $\Gamma^{ab}$ ($ab=25,51,12$) represent
the AFM and FM order parameters of $(x,y,z)$ directions, respectively.
Especially $\Gamma^{ab}$ ($ab=25,51,12$) correspond to the spin operator in
the Dirac Gamma matrices in particle physics.

Recalling the Dirac Gamma matrices, one may wonder that $\Gamma^{ab}$
($ab=23,31,12$) should be used instead of Eq. (LABEL:eq:Sigma_a). This
confusion is due to different representations of the Gamma matrices. Namely,
$\Gamma^{a}$ ($a=1,2,3,5$) can be replaced each other by a unitary
transformation. For instance, a unitary transformation
$\Gamma^{a}\to\Gamma^{a\,\prime}=U\Gamma^{a}U^{\dagger}$ by a unitary matrix

which is $U_{3}$ in Ref. [35], gives $\Gamma^{3}\to\Gamma^{5}$,
$\Gamma^{5}\to-\Gamma^{3}$, and the others remain. This indicates that we can
construct a representation of the Lorentz group from three matrices out of
four ones. In the sublattice basis, we can use $\Gamma^{1}$, $\Gamma^{2}$, and
$\Gamma^{5}$. Let us see this explicitly below.

To begin with, we define $J_{i}$ ($i=1,2,3$) and $K_{i}$ ($i=1,2,3$) as

$\displaystyle
J_{i}=\frac{1}{4}\epsilon_{iab}\Gamma^{ab}\,,~{}K_{i}=\frac{i}{4}\Gamma^{4i}\,,$
(88)

where $\epsilon_{iab}$ is the Levi-Civita symbol and indices $a$ and $b$ take
1, 2 or 5, such as $J_{1}=\Gamma^{25}/2,J_{2}=\Gamma^{51}/2$ etc. Then they
obey the following commutation relation:

$\displaystyle[J_{i},J_{j}]=i\epsilon_{ijk}J_{k}\,,~{}[K_{i},K_{j}]=-i\epsilon_{ijk}J_{k}\,,~{}[J_{i},K_{j}]=i\epsilon_{ijk}K_{k}\,,$
(89)

which shows that $J_{i}$ and $K_{i}$ are the generators of the Lorentz
transformation. Especially $J_{i}$ are the operators for the rotational
transformations. Let us focus on $J_{1}$ as an example. An operator for a
finite $\theta$ rotation along $i=1$, i.e., $x$ axis, is given by

$\displaystyle
S=e^{i\frac{\theta}{2}\Gamma^{25}}=\cos\frac{\theta}{2}+i\Gamma^{25}\sin\frac{\theta}{2}\,,$
(90)

Then we obtain

$\displaystyle S^{-1}\Gamma^{1}S$ $\displaystyle=\Gamma^{1}\,,$ (91)
$\displaystyle S^{-1}\Gamma^{2}S$
$\displaystyle=\Gamma^{2}\cos\theta+\Gamma^{5}\sin\theta\,,$ (92)
$\displaystyle S^{-1}\Gamma^{5}S$
$\displaystyle=-\Gamma^{2}\sin\theta+\Gamma^{5}\cos\theta\,.$ (93)

Therefore, it is confirmed that $J_{i}$ are the generators for rotation in
three dimensional space with axes $i=1,2,5$. Since the generators of the
rotational transformation are the spin operators up to a factor,
$\Gamma^{25}$, $\Gamma^{51}$, and $\Gamma^{12}$ are introduced for the FM
order. Accordingly, $\Gamma^{1}$, $\Gamma^{2}$, and $\Gamma^{5}$ are used for
the AFM order.

###  A.2 Wavefunction and Green’s function

In a finite temperature the wavefunction $\Phi_{i}$ and the Green’s function
$G$ is Fourier expanded as

$\displaystyle\Phi_{i}(\tau)$ $\displaystyle=\frac{1}{\sqrt{\beta
N}}\sum_{i\omega_{n},{\bf\it k}}\tilde{\Phi}(i\omega_{n},{\bf\it
k})e^{-i\omega_{n}\tau+i{\bf\it k}\dotproduct{\bf\it x}_{i}}\,,$ (94)
$\displaystyle G(x_{i}-x_{j})$ $\displaystyle=\frac{1}{\beta
N}\sum_{i\omega_{n},{\bf\it k}}\tilde{G}(i\omega_{n},{\bf\it
k})e^{-i\omega_{n}(\tau_{i}-\tau_{j})+i{\bf\it k}\dotproduct({\bf\it
x}_{i}-{\bf\it x}_{j})}\,,$ (95)

where ${\bf\it x}_{i}$ shows the location of the cite $i$ and
$x_{i}=(\tau_{i},{\bf\it x}_{i})$. $\Phi_{i}$ includes $c_{i}$,
$\delta\phi_{ai}$, $\delta\phi_{fi}$, and $n^{\pm}_{\alpha i}$,
$n^{\pm}_{\beta i}$. Regarding the Green’s function, we change the definition
of $G$ according to the magnetic order of the system. For instance, if there
is no magnetic order, then we choose $G^{-1}=-\partial_{\tau}-{\cal H}^{\rm
TI}$, leading to $\tilde{G}=(-i\omega_{n}+d_{0})^{-1}$.

##  Appendix B Dynamical susceptibility

In this section, we give the expression for the dynamical susceptibility used
in the numerical calculation. We have introduced eight kinds of
susceptibilities and they can be written in a generic form as

$\displaystyle\chi_{M}(k;O_{1},O_{2})$ $\displaystyle=-\frac{1}{\beta
N}\sum_{i\omega_{\ell},{\bf\it\ell}}{\rm
Tr}[\tilde{G}_{M}(\ell)O_{1}\tilde{G}_{M}(\ell+k)O_{2}]\,,$ (96)

where $O_{1}$ and $O_{2}$ stand for the Gamma matrices, such as $\Gamma^{a}$,
$\Gamma^{\pm}$, or $\Sigma^{\pm}$. Recalling that
$\tilde{G}_{M}(\ell)=(-i\omega_{\ell}+{\cal H}_{{\bf\it\ell}})^{-1}$,
$\tilde{G}_{M}(\ell)$ can be diagonalized by a unitary matrix as

$\displaystyle
U_{{\bf\it\ell}}^{\dagger}\tilde{G}_{M}(\ell)U_{{\bf\it\ell}}=\tilde{G}_{M}(\ell)_{\rm
diag}\,,$ (97)

where $\tilde{G}_{M}(\ell)_{\rm diag}$ is a diagonal matrix and
$\tilde{G}^{-1}_{M}(\ell)={\rm
diag}(-i\omega_{\ell}+E_{1{\bf\it\ell}},-i\omega_{\ell}+E_{2{\bf\it\ell}},-i\omega_{\ell}+E_{3{\bf\it\ell}},-i\omega_{\ell}+E_{4{\bf\it\ell}})$.
Similarly, $\tilde{G}_{M}(\ell+k)$ is diagonalized as

$\displaystyle U_{{\bf\it\ell}+{\bf\it
k}}^{\dagger}\tilde{G}_{M}(\ell+k)U_{{\bf\it\ell}+{\bf\it
k}}=\tilde{G}_{M}(\ell+k)_{\rm diag}\,.$ (98)

Then the trace part is

$\displaystyle{\rm Tr}[\tilde{G}_{M}(\ell)O_{1}\tilde{G}_{M}(\ell+k)O_{2}]$
$\displaystyle={\rm Tr}[\tilde{G}_{M}(\ell)_{\rm
diag}O^{\prime}_{1}\tilde{G}_{M}(\ell+k)_{\rm diag}O^{\prime}_{2}]$
$\displaystyle=\sum_{i,j}\frac{1}{i\omega_{n}-E_{i{\bf\it
k}}}\frac{1}{i\omega_{n}+i\omega_{\ell}-E_{j{\bf\it\ell}+{\bf\it
k}}}(O^{\prime}_{1})_{ij}(O^{\prime}_{2})_{ji}$ (99)

where

$\displaystyle O^{\prime}_{1}$
$\displaystyle=U_{{\bf\it\ell}}^{\dagger}O_{1}U_{{\bf\it\ell}+{\bf\it k}}\,,$
(100) $\displaystyle O^{\prime}_{2}$ $\displaystyle=U_{{\bf\it\ell}+{\bf\it
k}}^{\dagger}O_{2}U_{{\bf\it\ell}}\,.$ (101)

Summing over $i\omega_{\ell}$, we finally get

$\displaystyle\chi_{M}(i\omega_{n},{\bf\it k};O_{1},O_{2})$
$\displaystyle=-\sum_{{\bf\it\ell}}\sum_{i,j}\frac{(O^{\prime}_{1})_{ij}(O^{\prime}_{2})_{ji}}{i\omega_{n}+E_{j{\bf\it\ell}+{\bf\it
k}}-E_{i{\bf\it k}}}(n_{F}(E_{i{\bf\it\ell}})-n_{F}(E_{j{\bf\it\ell}+{\bf\it
k}}))\,.$ (102)

##  Appendix C Dynamical susceptibility under AFM at zero temperature

With the AFM order at zero temperature, the susceptibility can be derived
analytically. In the limit $\beta^{-1}\sum_{i\omega_{\ell}}$ is replaced by
$\int d\ell_{E}^{0}/(2\pi)$ where $\omega_{\ell}=\ell^{0}_{E}$. We list the
results below:

$\displaystyle\tilde{\chi}_{A}^{a}(i\omega_{n}=\omega,{\bf\it k})$
$\displaystyle=\frac{2}{N}\sum_{{\bf\it\ell}}\frac{df+\sum_{a=1}^{5}d_{a}f_{a}-2d_{5}^{2}}{P}\,,$
(103) $\displaystyle\tilde{\chi}_{A}^{f}(i\omega_{n}=\omega,{\bf\it k})$
$\displaystyle=\frac{2}{N}\sum_{{\bf\it\ell}}\frac{df-\sum_{a=1}^{5}d_{a}f_{a}+2(d_{1}f_{1}+d_{2}f_{2})}{P}\,,$
(104) $\displaystyle\tilde{\chi}_{A}^{\alpha}(i\omega_{n}=\omega,{\bf\it k})$
$\displaystyle=\frac{2}{N}\sum_{{\bf\it\ell}}\frac{df+\sum_{a=1}^{5}d_{a}f_{a}-(d_{1}f_{1}+d_{2}f_{2})}{P}\,,$
(105) $\displaystyle\tilde{\chi}_{A}^{\beta}(i\omega_{n}=\omega,{\bf\it k})$
$\displaystyle=\frac{2}{N}\sum_{{\bf\it\ell}}\frac{df-\sum_{a=1}^{5}d_{a}f_{a}+(2d_{5}^{2}+d_{1}f_{1}+d_{2}f_{2})}{P}\,,$
(106)

where

$\displaystyle P=df(d+f)\\{1-\omega^{2}/(d+f)^{2}\\}\,.$ (107)

$d_{a}$ and $f_{a}$ $(a=1$ – $4)$ depend on the wavenumber as
$d_{a}=d_{a}({\bf\it\ell})$ and $f_{a}=d_{a}({\bf\it\ell}+{\bf\it k})$ and
$f_{5}=d_{5}$. $f$ is defined by $f=(\sum_{a=1}^{5}f_{a}f_{a})^{1/2}$.

##  Appendix D How to derive the effective action of the excitation at zero
temperature

In this section we derive the effective action for the excitation in the
continuum and zero temperature limit. Using the dynamical susceptibility, the
effective (Euclidean) action for the excitation has a form (see Eqs. (40),
(59) and (64)),

$\displaystyle S_{E}=\frac{U\phi_{0}}{2}\sum_{i\omega,{\bf\it
k}}\tilde{\eta}^{*}(k)\tilde{\eta}(-k)\tilde{\Gamma}(k)\,,$ (108)

where $\eta$ is the excitation ($\tilde{\eta}$ is the Fourier coefficient) and

$\displaystyle\tilde{\Gamma}(k)\equiv
1-\frac{U}{4}\tilde{\chi}_{M}(i\omega_{n},{\bf\it k})\,.$ (109)

$\phi_{0}$ is unity, $\phi^{2}_{a0}$, or $\phi^{2}_{f0}$, depending on the
excitation. Since $\tilde{\Gamma}$ is the inverse of the two-point Green’s
function, the effective action can be derived by expanding it at the mass $m$,

$\displaystyle\tilde{\Gamma}(i\omega_{n},{\bf\it k})$
$\displaystyle=\partialderivative{\tilde{\Gamma}}{(i\omega_{n})^{2}}\Bigl{|}_{i\omega_{n}=m,{\bf\it
k}={\bf\it 0}}((i\omega_{n})^{2}-m^{2})$
$\displaystyle+\sum_{i}\partialderivative{\tilde{\Gamma}}{k_{i}^{2}}\Bigl{|}_{i\omega_{n}=m,{\bf\it
k}={\bf\it 0}}k_{i}^{2}+\cdots\,.$ (110)

Here we have used the definition of the mass:

$\displaystyle\tilde{\Gamma}(i\omega_{n}=m,{\bf\it 0})=0\,.$ (111)

In the continuum coordinate space it can be written as

$\displaystyle S_{E}$
$\displaystyle=\frac{U\phi_{0}}{2}\int^{\beta}d\tau\sum_{i}\eta^{*}\frac{U}{2}a^{3}J(-\partial_{\tau}^{2}+v_{i}^{2}k_{i}^{2}+m^{2})\eta$
$\displaystyle\to\frac{U\phi_{0}}{2}\int^{\beta}d\tau\int
d^{3}x\eta^{*}\frac{U}{2}J(-\partial_{\tau}^{2}+v_{i}^{2}k_{i}^{2}+m^{2})\eta\,,$
(112)

where we have retrieved the lattice size $a$ to make the dimension of the
action intact. $J$ and $v_{i}$ are the stiffness and velocity defined as

$\displaystyle\frac{U}{2}a^{3}J=-\partialderivative{\tilde{\Gamma}}{(i\omega_{n})^{2}}\Bigl{|}_{i\omega_{n}=m,{\bf\it
k}={\bf\it 0}}\,,$ (113)
$\displaystyle\frac{U}{2}a^{3}Jv_{i}^{2}=\partialderivative{\tilde{\Gamma}}{k_{i}^{2}}\Bigl{|}_{i\omega_{n}=m,{\bf\it
k}={\bf\it 0}}\,.$ (114)

Using the real time ($t=-i\tau$), we get the (Minkowski) action as

$\displaystyle iS_{E}\to S=\left(\frac{U}{2}\right)^{2}\phi_{0}J\int
d^{4}x\eta^{*}(-\partial_{t}^{2}+v_{i}^{2}\partial_{i}^{2}-m^{2})\eta\,.$
(115)

Let us calculate the stiffness for the AFM-type amplitude mode in the AFM
order. In this case, the mass is given by $U\phi_{0}$. Then

$\displaystyle\frac{U}{2}a^{3}J$
$\displaystyle=\frac{U}{2N}\sum_{{\bf\it\ell}}\frac{d_{0}^{2}}{4d^{5}}\frac{1}{(1-d_{5}^{2}/d^{2})^{2}}$
$\displaystyle\to~{}~{}$ $\displaystyle
J=\int\frac{d^{3}\ell}{(2\pi)^{3}}\frac{1}{4dd_{0}^{2}}\,.$ (116)

In the previous works [6, 7], $\tilde{\Gamma}$ is expanded around
$i\omega_{n}=0$. Namely,

$\displaystyle\tilde{\Gamma}^{a}_{A}(i\omega_{n}=\omega,{\bf\it 0})$
$\displaystyle=\frac{U}{2N}\sum_{{\bf\it\ell}}\left[\frac{d^{2}_{5}}{d^{3}}-\frac{d_{0}^{2}}{4d^{5}}\omega^{2}+\order{\omega^{4}}\right]\,.$
(117)

As a result, we get

$\displaystyle J^{\rm old}$
$\displaystyle=\int\frac{d^{3}\ell}{(2\pi)^{3}}\frac{d^{2}_{0}}{4d^{5}}\,,$
(118) $\displaystyle J^{\rm old}(m_{a}^{\rm old})^{2}$
$\displaystyle=\int\frac{d^{3}\ell}{(2\pi)^{3}}\frac{d_{5}^{2}}{d^{3}}\,.$
(119)

It is clear that the expansion around $\omega=0$ in Eq. (117) is guaranteed if
$m_{a}^{\rm old}$ is smaller than $2d$. If this condition is satisfied, then
the expansion is valid and we expect $m_{a}^{\rm old}\simeq m_{a}$. The
discrepancy between the old results and the new ones is shown in Fig. 10. We
found that they agree at around $(10^{-3}\,\mathchar 45\relax\,10)$% for
$m_{a}\sim$ (meV - eV). Therefore the old results are a good approximation for
$m_{a}\lesssim\order{\rm eV}$.

Figure 10: Comparison of the mass and stiffness of axionic excitation.
$r_{m}=-(m_{a}-m_{a}^{\rm old})/m_{a}$ and $r_{J}=(J-J^{\rm old})/J$ are
plotted as function of $m_{a}$. The other parameters are the same as Fig. 1.

Similarly, the stiffness for the $\alpha$-type magnon with the AFM is

$\displaystyle
J_{\alpha}=\int\frac{d^{3}\ell}{(2\pi)^{3}}\frac{2d^{2}-d_{1}^{2}-d_{2}^{2}}{8d^{3}}\frac{1}{(1-m_{\alpha}^{2}/(4d^{2}))^{2}}\,,$
(120)

where the mass $m_{\alpha}$ is given by
$\tilde{\Gamma}_{A}^{\alpha}(i\omega_{n}=m_{\alpha},{\bf\it 0})=0$. On the
other hand, if we expand $\tilde{\Gamma}^{\alpha}_{A}$ around $i\omega_{n}=0$,
i.e.,

$\displaystyle\tilde{\Gamma}^{\alpha}_{A}(i\omega_{n}=\omega,{\bf\it 0})$
$\displaystyle=\frac{U}{2N}\sum_{{\bf\it\ell}}\left[\frac{d^{2}_{1}+d^{2}_{2}}{2d^{3}}-\frac{2d^{2}-d_{1}^{2}-d_{2}^{2}}{8d^{5}}\omega^{2}+\order{\omega^{4}}\right]\,,$
(121)

we obtain

$\displaystyle J^{\prime}_{\alpha}$
$\displaystyle=\int\frac{d^{3}\ell}{(2\pi)^{3}}\frac{2d^{2}-d_{1}^{2}-d_{2}^{2}}{8d^{5}}\,,$
(122) $\displaystyle J^{\prime}_{\alpha}m_{\alpha}^{\prime 2}$
$\displaystyle=\int\frac{d^{3}\ell}{(2\pi)^{3}}\frac{d_{1}^{2}+d_{2}^{2}}{2d^{3}}\,.$
(123)

##  Appendix E Symmetry of the Hamiltonian

To understand the appearance of the gapless modes, we use

$\displaystyle\vec{\phi}_{ai}$
$\displaystyle=(\vec{\phi}_{Ai}-\vec{\phi}_{Bi})/2\,,$ (124)
$\displaystyle\vec{\phi}_{fi}$
$\displaystyle=(\vec{\phi}_{Ai}+\vec{\phi}_{Bi})/2\,,$ (125)

in this section. The Hubbard term after the Stratonovich-Hubbard
transformation gives,

$\displaystyle
H_{U}\to\frac{U}{2}\sum_{i}(\vec{\phi}_{ai}^{\,2}+\vec{\phi}_{fi}^{\,2})+\frac{U}{2}\sum_{i}c^{\dagger}_{i}\Bigl{[}\vec{\phi}_{ai}\dotproduct\vec{\Gamma}+\vec{\phi}_{fi}\dotproduct\vec{\Sigma}\Bigr{]}c_{i}\,,$
(126)

where we have introduced $\vec{\Gamma}=(\Gamma^{1},\Gamma^{2},\Gamma^{5})$ and
$\vec{\Sigma}=(\Gamma^{25},\Gamma^{51},\Gamma^{12})$. It can be seen that the
rotational invariance of the Hubbard term is intact in the parameter space of
magnetic orders, $\vec{\phi}_{ai}$ and $\vec{\phi}_{fi}$, which corresponds to
$(\Gamma^{1},\Gamma^{2},\Gamma^{5})$ space. See also Appendix A.1. This
symmetry is broken by ${\cal H}^{\rm TI}$, i.e., $\Gamma^{1}$ and $\Gamma^{2}$
terms that couple to the wavenumbers. To see the structure of the symmetry, we
take $d_{1}=d_{2}=0$ hereafter. In addition we consider a uniform magnetic
configuration to study the magnons of the zero mode.

The Hamiltonian ${\cal H}$ is given by

$\displaystyle{\cal H}$ $\displaystyle={\cal H}^{\rm MTI}+\delta{\cal H}$
$\displaystyle=\sum_{a=3}^{4}d_{a}\Gamma^{a}+\frac{U}{2}(\vec{\phi}_{a}\dotproduct\vec{\Gamma}+\vec{\phi}_{f}\dotproduct\vec{\Sigma})\,.$
(127)

and the energy eigenvalues are

$\displaystyle\pm\sqrt{d_{s}^{2}+\vec{\phi}_{a}^{2}+\vec{\phi}_{f}^{2}+2\sqrt{\vec{\phi}_{a}\dotproduct\vec{\phi}_{f}+d_{s}^{2}(\vec{\phi}_{a}^{2}+\vec{\phi}_{f}^{2})}}\,,$
(128)
$\displaystyle\pm\sqrt{d_{s}^{2}+\vec{\phi}_{a}^{2}+\vec{\phi}_{f}^{2}-2\sqrt{\vec{\phi}_{a}\dotproduct\vec{\phi}_{f}+d_{s}^{2}(\vec{\phi}_{a}^{2}+\vec{\phi}_{f}^{2})}}\,.$
(129)

The expression clearly shows that the system has SO(3) symmetry. Under the AFM
state, with $(\phi_{az},\phi_{fz})=(\phi_{a0},0)$, you can see that SO(3) is
broken to SO(2). Therefore, massless Nambu-Goldstone bosons are expected,
which corresponds to the $\alpha$-type magnon. The symmetry breaking by the FM
order is the same. In this case, the $\alpha$-type magnon is the Nambu-
Goldstone bosons.

Figure 11: Dispersion relation for $\alpha$-type magnon under the FM. (a)
Contour of $\omega(k_{x},k_{y},\pi)$. (b) $\omega({\bf\it k})$ as function of
$k_{z}$ for various values of $k_{x},k_{y}$.

Figure 12: Same as Fig. 11 but for $\beta$-type magnon under the FM. (a)
Contour of $\omega(k_{x},k_{y},0)$. (b) $\omega({\bf\it k})$ as function of
$k_{z}$ for various values of $k_{x},k_{y}$. The other parameters are the same
as Fig. 1.

Figure 13: Same as Fig. 9 computed in the Dirac model. $M_{0}\,[{\rm meV}]=\pm
10$ (left) and $\pm 10^{2}$ (right) are taken, where positive (negative)
$M_{0}$ corresponds to the NI (TI) phases. The other parameters are taken as
$A_{1}=A_{2}=1~{}{\rm eV}$.

##  Appendix F Additional figures

The dispersion relations for $\alpha$\- and $\beta$-type magnons under the FM
order are shown in Figs. 11 and 12, respectively.

As the AFM order, we expect a gapless mode under the FM in the case of
$d_{1}=d_{2}=0$. In fact, we found the $\alpha$-type magnon becomes gapless at
${\bf\it k}={\bf\it 0}$. While it is gapless, the mode turns out to be
unstable. Here ‘unstable’ means

$\displaystyle\tilde{\Gamma}^{\alpha}_{F}(0,{\bf\it k})<0\,,$ (130)

around ${\bf\it k}={\bf\it 0}$. Namely, the spatial kinetic term has the
opposite sign. Thus, such a mode is regarded as unphysical.

Fig. 13 shows the effective decay constant computed in the Dirac model. We
parametrize the Dirac model as

$\displaystyle(d_{1},d_{2},d_{3},d_{4})=(A_{2}k_{x},\,A_{2}k_{y},\,A_{1}k_{z},\,M_{0})\,.$
(131)

## References

  * [1] D. J. E. Marsh, K.-C. Fong, E. W. Lentz, L. Smejkal, and M. N. Ali,  Phys. Rev. Lett. 123, 121601 (2019), arXiv:1807.08810.
  * [2] Y. Tokura, K. Yasuda, and A. Tsukazaki,  Nature Reviews Physics 1, 126 (2019).
  * [3] B. A. Bernevig, C. Felser, and H. Beidenkopf,  Nature 603, 41 (2022).
  * [4] J. Liu and T. Hesjedal,  Advanced Materials 35, 2102427 (2023).
  * [5] A. Sekine and K. Nomura,  Journal of Applied Physics 129 (2021).
  * [6] R. Li, J. Wang, X. Qi, and S.-C. Zhang,  Nature Phys. 6, 284 (2010), arXiv:0908.1537.
  * [7] K. Ishiwata,  Phys. Rev. D 104, 016004 (2021), arXiv:2103.02848.
  * [8] K. Ishiwata,  Phys. Rev. B 106, 195157 (2022), arXiv:2206.00841.
  * [9] L. Cao et al.,  Phys. Rev. B 104, 054421 (2021).
  * [10] Y. Li et al.,  Phys. Rev. B 102, 121107 (2020), arXiv:2001.06133.
  * [11] J. Wang, B. Lian, and S.-C. Zhang,  Physical Review Letters 115 (2015).
  * [12] B. Lake, D. A. Tennant, and S. E. Nagler,  Phys. Rev. Lett. 85, 832 (2000).
  * [13] A. Zheludev, K. Kakurai, T. Masuda, K. Uchinokura, and K. Nakajima,  Phys. Rev. Lett. 89, 197205 (2002).
  * [14] C. Rüegg et al.,  Phys. Rev. Lett. 93, 257201 (2004).
  * [15] C. Rüegg et al.,  Phys. Rev. Lett. 100, 205701 (2008).
  * [16] A. Jain et al.,  Nature Physics 13, 633–637 (2017).
  * [17] T. Hong et al.,  Nature Physics 13, 638–642 (2017).
  * [18] S. Hayashida et al.,  Phys. Rev. B 97, 140405 (2018).
  * [19] S. Hayashida et al.,  Science Advances 5 (2019).
  * [20] C. Herring and C. Kittel,  Phys. Rev. 81, 869 (1951).
  * [21] C. Herring,  Phys. Rev. 85, 1003 (1952).
  * [22] C. Herring,  Phys. Rev. 87, 60 (1952).
  * [23] J. König, T. Jungwirth, and A. H. MacDonald,  Phys. Rev. B 64, 184423 (2001).
  * [24] Y. Araki and K. Nomura,  Phys. Rev. B 93, 094438 (2016).
  * [25] H. Zhang et al.,  Nature Phys. 5, 438 (2009).
  * [26] C.-X. Liu et al.,  Physical Review B 82 (2010).
  * [27] G. Rosenberg and M. Franz,  Physical Review B 85 (2012).
  * [28] D. Kurebayashi and K. Nomura,  Journal of the Physical Society of Japan 83, 063709 (2014).
  * [29] J. Wang, B. Lian, and S.-C. Zhang,  Phys. Rev. B 93, 045115 (2016), arXiv:1512.00534.
  * [30] A. Sekine and K. Nomura,  Phys. Rev. Lett. 116, 096401 (2016), arXiv:1508.04590.
  * [31] J. Schütte-Engel et al.,  JCAP 08, 066 (2021), arXiv:2102.05366.
  * [32] H. Ooguri and M. Oshikawa,  Phys. Rev. Lett. 108, 161803 (2012), arXiv:1112.1414.
  * [33] Y. Li et al.,  Physical Review B 102 (2020).
  * [34] J. Zhang et al.,  Chinese Physics Letters 37, 077304 (2020).
  * [35] S. Chigusa, T. Moroi, and K. Nakayama,  JHEP 08, 074 (2021), arXiv:2102.06179.