# Axion production in the $\eta\to\pi\pi a$ decay within $SU(3)$ chiral
perturbation theory

Jin-Bao Wanga,b, Zhi-Hui Guoa111<EMAIL_ADDRESS>Zhun Lub222
<EMAIL_ADDRESS>, Hai-Qing Zhoub333<EMAIL_ADDRESS>
a Department of Physics and Hebei Key Laboratory of Photophysics Research and
Application,
Hebei Normal University, Shijiazhuang 050024, China
b School of Physics, Southeast University, Nanjing 211189, China

###### Abstract

We study the axion and axion-like particle production from the $\eta\to\pi\pi
a$ decay within the $SU(3)$ chiral perturbation theory up to the one-loop
level. The conventional $SU(3)$ chiral low energy constants are found to be
able to reabsorb all the divergences from the chiral loops in the
$\eta\to\pi\pi a$ decay amplitude, and hence render the amplitude independent
of the renormalization scale. The unitarized $\eta\to\pi\pi a$ decay
amplitudes are constructed to take into account the $\pi\pi$ final-state
interactions and also properly reproduce the perturbative results from the
chiral perturbation theory. Detailed analyses between the perturbative
amplitudes and the unitarized ones are given in the phenomenological
discussions. By taking the values of the chiral low energy constants in
literature, we predict the Dalitz distributions, the spectra of the $\pi\pi$
and $a\pi$ systems, and also the branching ratios of the $\eta\to\pi\pi a$
process by varying $m_{a}$ from 0 to $m_{\eta}-2m_{\pi}$.

##  1 Introduction

The hunting of the hypothetical axion ($a$) and axion-like particles (ALPs)
constitutes one of the most active research subjects in particle physics
nowadays [1, 2, 3, 4, 5, 6]. As an elegant resolution of the strong CP problem
in QCD, the interaction of the axion and gluons, i.e.
$aG^{\mu\nu}\tilde{G}_{\mu\nu}\alpha_{s}/(8\pi f_{a})$, where $G^{\mu\nu}$ is
the gluon field strength tensor, $\tilde{G}_{\mu\nu}$ corresponds to its dual
and $f_{a}$ denotes the axion decay constant, is deemed to be the model
independent ingredient in the numerous constructions of axion models. In this
case, also dubbed as the QCD axion case, not only all the axion interactions
but also the axion mass are governed by the single parameter $f_{a}$ [7],
leading to a stringent constraint on the former parameter. The visible QCD
axion scenario, i.e., when $f_{a}$ lies at the electroweak symmetry breaking
scale $v_{EW}\simeq 256$ GeV, proposed in the early pioneer works of Refs. [8,
9, 7, 10], has been ruled out by various experiments [11]. As a result, the
invisible axion scenario with $f_{a}\gg v_{EW}$ becomes the focus in various
studies. Axion can only feebly interact with particles of the Standard Model
(SM), since the interaction strengths of the axion are generally suppressed by
$(1/f_{a})^{n}$, being $n$ the number of axions involved.

In this work, we stick to the model independent term
$aG^{\mu\nu}\tilde{G}_{\mu\nu}\alpha_{s}/(8\pi f_{a})$ to account for the
axion interaction, and this inevitably introduces the axion couplings with the
QCD hadrons in the low energy region. Regarding the axion mass, the
$aG\tilde{G}$ term gives a contribution behaving like $m_{a}^{2}\propto
m_{\pi}^{2}F_{\pi}^{2}/f_{a}^{2}$ [7], which becomes rather tiny in the
invisible axion scenario. In the ALP case, it is possible to introduce a bare
mass term $m_{a,0}$ to the axion from a UV theory, see the very recent
discussions in Refs. [4, 12, 13, 14, 15, 16], and the axion mass is then
disentangled from its decay constant $f_{a}$. As long as $m_{a,0}\ll f_{a}$,
the presence of the bare mass term $m_{a,0}$ can be considered as a soft
breaking of the Peccei-Quinn symmetry. In the phenomenological exploration,
$m_{a,0}$ can be taken as a free parameter in the region much below $f_{a}$.
We will concentrate on the axion and ALP production from the $\eta\to\pi\pi a$
process by varying the ALP mass from zero to $m_{\eta}-2m_{\pi}$.

The leading order (LO) amplitudes of $\eta^{(^{\prime})}\to\pi\pi a$ in a
$U(3)$ chiral framework have been calculated in Ref. [17], and it turns out
that apart from the tiny isospin breaking corrections there are only constant
terms, which are proportional to $m_{\pi}^{2}$, remaining in the amplitudes. A
variant chiral calculation in Ref. [18] gives comparable expressions as the LO
results in Ref. [17] that only contain the constant terms, and the energy
dependent parts are absent at this order. Therefore this implies that the
higher order effects could give important contributions, comparing with the LO
ones. The $\eta^{(^{\prime})}\to\pi\pi a$ decay widths are calculated in Ref.
[19] by including the bare resonance exchanges at tree level. Very recently,
the $\eta^{(^{\prime})}\to\pi\pi a$ amplitudes are calculated by incorporating
the $\pi\pi$ final state interaction into the LO $\eta^{(^{\prime})}\to\pi\pi
a$ amplitudes in Ref. [20]. In the current study, we pursue the $\eta\to\pi\pi
a$ calculation by working out the complete one-loop expressions in the $SU(3)$
chiral perturbation theory ($\chi$PT) and study the unitarization effects,
which are then employed to explore the phenomenological consequences, such as
the Dalitz decay distributions, the spectra of the $\pi\pi$ and $a\pi$
systems, and the $\eta\to\pi\pi a$ decay widths as a function of the ALP mass.

This paper is organized as follows. In Sec. 2, the axion $SU(3)$ chiral
Lagrangian up to $\mathcal{O}(p^{4})$ is elaborated. In Sec. 3, we discuss the
ALP-$\eta$ mixing and the ALP mass up to one-loop level. The decay amplitude
of the $\eta\to\pi\pi a$ process is calculated up to $\mathcal{O}(p^{4})$ in
Sec. 4, where the unitarization procedure is also discussed to include the
$\pi\pi$ final-state interaction. In Sec. 5, we study the phenomenologies in
the $\eta\to\pi\pi a$ decay by separately analyzing the perturbative
amplitudes and the unitarized ones. We give a short summary and conclusions in
Sec. 6.

##  2 Axion $SU(3)$ chiral Lagrangian

Although it is straightforward to include additional model-dependent unknown
axion couplings with quarks, we carry out the calculations by taking the
model-independent axion interaction, namely the couplings exclusively derived
from the $aG^{\mu\nu}\tilde{G}_{\mu\nu}\alpha_{s}/(8\pi f_{a})$ operator. In
this way it allows us to proceed the phenomenological discussions later with
more definite results, since all the unknown parameters except $f_{a}$ have
been determined in the studies with light-flavor mesons [21].

A practical way to handle the $aG^{\mu\nu}\tilde{G}_{\mu\nu}\alpha_{s}/(8\pi
f_{a})$ operator is to perform an axial transformation of the quark fields

$q\to\exp\left(i\frac{a}{2f_{a}}Q_{a}\gamma_{5}\right)q\,,$ (1)

where $Q_{a}$ is a $3\times 3$ hermitian matrix in light-quark flavor space
with its trace $\langle Q_{a}\rangle=1$. In such a manner one can eliminate
the $aG\widetilde{G}$ term and meanwhile introduces two new modifications,
namely the axion-quark derivative interaction and the modified quark mass
term. The resulting effective axion Lagrangian after the transformation of Eq.
(1) reads

$\mathcal{L}_{\text{ALP}}=\frac{1}{2}\partial_{\mu}a\partial^{\mu}a-\frac{1}{2}m_{a,0}^{2}a^{2}-\frac{\partial_{\mu}a}{2f_{a}}\bar{q}\gamma^{\mu}\gamma_{5}Q_{a}q-\bar{q}\exp\left(i\frac{a}{2f_{a}}Q_{a}\gamma_{5}\right)M\exp\left(i\frac{a}{2f_{a}}Q_{a}\gamma_{5}\right)q\,,$
(2)

where $m_{a,0}$ denotes the bare ALP mass and the light-flavor quark mass
matrix is given by $M={\rm diag}(m_{u},m_{d},m_{s})$. It is verified that the
strong isospin-breaking correction from the $\epsilon_{I}\equiv m_{d}-m_{u}$
factor in the $\eta\to\pi\pi a$ decay begins at
$\mathcal{O}(\epsilon_{I}^{2})$, which therefore can be safely neglected. As a
result we will work in isospin limit and take
$M=\text{diag}(\hat{m},\,\hat{m},\,m_{s})$ in later discussions.

Next we match the effective axion Lagrangian in Eq. (2) to the $SU(3)~{}\rm$.
The LO chiral Lagrangian with axion can be constructed as

$\mathcal{L}_{2}=\frac{F^{2}}{4}\langle\partial_{\mu}U\partial^{\mu}U^{\dagger}+\chi_{a}U^{\dagger}+U\chi_{a}^{\dagger}\rangle+\frac{\partial_{\mu}a}{2f_{a}}J_{A}^{\mu}\big{|}_{\mathrm{LO}}+\frac{1}{2}\partial_{\mu}a\partial^{\mu}a-\frac{1}{2}m_{a,0}^{2}a^{2}\,,$
(3)

where the octet of pseudo Nambu-Goldstone bosons (pNGBs) is collected in
$U={\exp}\big{(}i\Phi/F\big{)}$ with

$\Phi\,=\,\left(\begin{matrix}\pi^{0}+\frac{1}{\sqrt{3}}\eta_{8}&\sqrt{2}\pi^{+}&\sqrt{2}K^{+}\\\
\sqrt{2}\pi^{-}&-\pi^{0}+\frac{1}{\sqrt{3}}\eta_{8}&\sqrt{2}K^{0}\\\
\sqrt{2}K^{-}&\sqrt{2}\bar{K}^{0}&-\frac{2}{\sqrt{3}}\eta_{8}\end{matrix}\right)\,.$
(4)

It is noted that the physical $\eta$ meson is approximated by the $SU(3)$
octet $\eta_{8}$ field in the $SU(3)$ . The axion dressed quark mass matrix is
given by $\chi_{a}=2B_{0}M(a)$ with

$M(a)\equiv\exp\left(-i\frac{a}{2f_{a}}Q_{a}\right)M\exp\left(-i\frac{a}{2f_{a}}Q_{a}\right)\,.$
(5)

There is arbitrariness for the choice of $Q_{a}$ in Eq. (1) with the
constraint $\langle Q_{a}\rangle=1$, though different choices do not change
the physical results. It is customary to take
$Q_{a}=M_{q}^{-1}/\langle{M_{q}^{-1}}\rangle$ [22] which eliminates the mass
mixing between ALP and octet pNGBs at LO. The
$J_{A}^{\mu}\big{|}_{\mathrm{LO}}$ term in Eq. (3) corresponds to the hadronic
parameterization of the quark axial-vector current
$\bar{q}\gamma^{\mu}\gamma_{5}(-Q_{a})q$ in Eq. (2) at LO. The explicit
expression of $J_{A}^{\mu}\big{|}_{\mathrm{LO}}$ is given by

$J_{A}^{\mu}\big{|}_{\text{LO}}=-i\frac{F^{2}}{2}\langle
Q_{a}\left\\{\partial^{\mu}U,\,U^{\dagger}\right\\}\rangle\,.$ (6)

$Q_{a}$ contains both the flavor-octet and flavor-singlet components. It is
easy to verify that the the singlet component of the LO hadronic axial-vector
current in Eq. (6) vanishes and only the octet current contributes to the
derivative axion couplings at this order.

At next-to-leading order (NLO) in the chiral expansion, a physical amplitude
can receive contributions both from the one-loop Feynman diagrams and the
local operators of the $\mathcal{O}(p^{4})$ chiral Lagrangian [23]

$\displaystyle\mathcal{L}_{4}\supset$ $\displaystyle
L_{1}\langle\partial_{\mu}U\partial^{\mu}U^{\dagger}\rangle\langle\partial_{\nu}U\partial^{\nu}U^{\dagger}\rangle+L_{2}\langle\partial_{\mu}U\partial_{\nu}U^{\dagger}\rangle\langle\partial^{\mu}U\partial^{\nu}U^{\dagger}\rangle+L_{3}\langle\partial_{\mu}U\partial^{\mu}U^{\dagger}\partial_{\nu}U\partial^{\nu}U^{\dagger}\rangle$
(7)
$\displaystyle+L_{4}\langle\partial_{\mu}U\partial^{\mu}U^{\dagger}\rangle\langle\chi_{a}U^{\dagger}+U\chi_{a}^{\dagger}\rangle+L_{5}\langle\partial_{\mu}U\partial^{\mu}U^{\dagger}(\chi_{a}U^{\dagger}+U\chi_{a}^{\dagger})\rangle+L_{6}\left(\langle\chi_{a}U^{\dagger}+U\chi_{a}^{\dagger}\rangle\right)^{2}$
$\displaystyle+L_{7}\left(\langle\chi_{a}U^{\dagger}-U\chi_{a}^{\dagger}\rangle\right)^{2}+L_{8}\langle\chi_{a}U^{\dagger}\chi_{a}U^{\dagger}+U\chi_{a}^{\dagger}U\chi_{a}^{\dagger}\rangle+\frac{\partial_{\mu}a}{2f_{a}}J^{\mu}_{A}\big{|}_{\text{NLO}}\,,$

with the NLO hadronic axial-vector current

$\displaystyle J^{\mu}_{A}\big{|}_{\text{NLO}}=$
$\displaystyle-4iL_{1}\langle\bar{Q}_{a}\\{U^{\dagger},\,\partial^{\mu}U\\}\rangle\langle\partial_{\nu}U\partial^{\nu}U^{\dagger}\rangle$
(8)
$\displaystyle-2iL_{2}\langle\bar{Q}_{a}\\{U^{\dagger},\,\partial_{\nu}U\\}\rangle\langle\partial^{\mu}U\partial^{\nu}U^{\dagger}+\partial^{\nu}U\partial^{\mu}U^{\dagger}\rangle$
$\displaystyle-2iL_{3}\langle\partial^{\mu}U\\{\bar{Q}_{a},\,U^{\dagger}\\}\partial_{\nu}U\partial^{\nu}U^{\dagger}\rangle+2iL_{3}\langle\\{\bar{Q}_{a},\,U\\}\partial^{\mu}U^{\dagger}\partial_{\nu}U\partial^{\nu}U^{\dagger}\rangle$
$\displaystyle-2iL_{4}\langle\bar{Q}_{a}\\{U^{\dagger},\,\partial^{\mu}U\\}\rangle\langle\chi
U^{\dagger}+U\chi^{\dagger}\rangle$ $\displaystyle-
iL_{5}\langle\partial^{\mu}U\\{\bar{Q}_{a},\,U^{\dagger}\\}(\chi
U^{\dagger}+U\chi^{\dagger})\rangle+iL_{5}\langle\\{\bar{Q}_{a},\,U\\}\partial^{\mu}U^{\dagger}(\chi
U^{\dagger}+U\chi^{\dagger})\rangle\,,$

where $\bar{Q}_{a}$ denotes the octet component of $Q_{a}$. We have only shown
the relevant operators to this work in Eqs. (7) and (8). Although the singlet
component of the hadronic axial current at LO in Eq. (6) vanishes, it does
contribute at NLO in Eq. (8). Nevertheless, a more complete framework to
investigate the singlet axial-vector current would be the $U(3)$ that
explicitly includes the QCD $U_{A}(1)$ anomaly effect and the iso-singlet
$\eta_{0}$ state, and this is beyond the scope of the present study that
relies on the $SU(3)$ . In fact the inclusion of the singlet axial current in
the $U(3)$ with the explicit $\eta_{0}$ field will change the renormalization
of the $SU(3)$ LECs [24, 25]. Since we stick to the $SU(3)$ framework to
calculate the $\eta\to\pi\pi a$ process, the singlet component for the NLO
hadronic axial currents in Eq. (8) will be omitted as an approximation. In
this case, the $\eta\to\pi\pi a$ amplitude can be renormalized with the
conventional renormalization conditions for the standard $SU(3)$ LECs [23], as
will be demonstrated later.

##  3 ALP-$\eta$ mixing and ALP mass up to NLO

For the perturbative calculation of the amplitudes, we need to first solve the
mixing of the ALP and $\pi^{0}$, $\eta_{8}$ fields. The key difference between
the ALP-$\pi^{0}$ and ALP-$\eta_{8}$ mixing is that the former vanishes in the
isospin limit and the latter remains finite. It is pointed out that the
ALP-$\pi^{0}$ mixing will contribute to the $\eta\to\pi\pi a$ amplitude at the
level of $\mathcal{O}(\epsilon_{I}^{2})$, and we will neglect the strong
isospin breaking effects throughout. In this section, we deal with
ALP-$\eta_{8}$ mixing up to NLO in the $SU(3)$ and calculate the NLO chiral
correction to the ALP mass as well. The chiral correction to the ALP mass
comes from the couplings between ALP and pNGBs and thus is of
$\mathcal{O}(1/f_{a}^{2})$. Keeping the terms up to
$\mathcal{O}(1/f_{a}^{2})$, the LO quadratic terms of ALP and $\eta_{8}$
fields from LO Lagrangian Eq. (3) take the form

$\mathcal{L}_{\mathrm{mix}}^{\mathrm{LO}}=\frac{1}{2}\partial_{\mu}a\partial^{\mu}a+\frac{1}{2}\partial_{\mu}\eta_{8}\partial^{\mu}\eta_{8}+\frac{F}{f_{a}}C_{k}^{a\eta}\partial_{\mu}a\partial^{\mu}\eta_{8}-\frac{1}{2}\left(m_{a,0}^{2}+\frac{F^{2}}{f_{a}^{2}}C_{m}^{a}\right)a^{2}-\frac{1}{2}\bar{m}_{\eta_{8}}^{2}\eta_{8}^{2}\,,$
(9)

with

$C_{k}^{a\eta}=\frac{\bar{m}_{\eta_{8}}^{2}-\bar{m}_{\pi}^{2}}{2\sqrt{3}\bar{m}_{\eta_{8}}^{2}}\,,\quad
C_{m}^{a}=\frac{\bar{m}_{\pi}^{2}(3\bar{m}_{\eta_{8}}^{2}-\bar{m}_{\pi}^{2})}{12\bar{m}_{\eta_{8}}^{2}}\,.$
(10)

$\bar{m}_{\pi}$ and $\bar{m}_{\eta_{8}}$ are the LO masses of $\pi$ and
$\eta_{8}$ contributed by light-quark masses derived from Eq. (3)

$\bar{m}_{\pi}^{2}=2B_{0}\hat{m}\,,\quad\bar{m}_{\eta_{8}}^{2}=\frac{2B_{0}}{3}(\hat{m}+2m_{s})\,,$
(11)

and the LO kaon mass is $\bar{m}_{K}^{2}=B_{0}(\hat{m}+m_{s})$, which gives
the celebrated Gell-Mann-Okubo relation
$4\bar{m}_{K}^{2}=3\bar{m}_{\eta_{8}}^{2}+\bar{m}_{\pi}^{2}.$ Since we are
interested in $\eta\to\pi\pi a$ decay, for convenience, we will express the
light-quark masses $\hat{m}$ and $m_{s}$ in terms of $\bar{m}_{\pi}^{2}$ and
$\bar{m}_{\eta_{8}}^{2}$ in this work.

By performing the following field redefinitions for ALP and $\eta_{8}$ fields

$\displaystyle
a=\left[1+\frac{1}{2}{(C_{k}^{a\eta})}^{2}\frac{m_{a,0}^{2}(m_{a,0}^{2}-2\bar{m}_{\eta_{8}}^{2})}{(m_{a,0}^{2}-\bar{m}_{\eta_{8}}^{2})^{2}}\frac{F^{2}}{f_{a}^{2}}\right]\widetilde{a}-\frac{F}{f_{a}}C_{k}^{a\eta}\frac{\bar{m}_{\eta_{8}}^{2}}{\bar{m}_{\eta_{8}}^{2}-m_{a,0}^{2}}\widetilde{\eta}\,,$
(12)
$\displaystyle\eta_{8}=\left[1+\frac{1}{2}{(C_{k}^{a\eta})}^{2}\frac{\bar{m}_{\eta_{8}}^{2}(\bar{m}_{\eta_{8}}^{2}-2m_{a,0}^{2})}{(m_{a,0}^{2}-\bar{m}_{\eta_{8}}^{2})^{2}}\frac{F^{2}}{f_{a}^{2}}\right]\widetilde{\eta}-\frac{F}{f_{a}}C_{k}^{a\eta}\frac{m_{a,0}^{2}}{m_{a,0}^{2}-\bar{m}_{\eta_{8}}^{2}}\widetilde{a}\,,$
(13)

the LO quadratic terms Eq. (9) can be diagonalized to the canonical form up to
$\mathcal{O}(1/f_{a}^{2})$

$\mathcal{L}_{\mathrm{mix}}^{\mathrm{LO}}=\frac{1}{2}\partial_{\mu}\widetilde{a}\partial^{\mu}\widetilde{a}-\frac{1}{2}\bar{m}_{a}^{2}\widetilde{a}^{2}+\frac{1}{2}\partial_{\mu}\widetilde{\eta}\partial^{\mu}\widetilde{\eta}-\frac{1}{2}\bar{m}_{\eta}^{2}\widetilde{\eta}^{2}\,.$
(14)

The resulting LO masses of ALP and $\eta$ after the elimination of the mixing
read

$\displaystyle\bar{m}_{a}^{2}$ $\displaystyle=$ $\displaystyle
m_{a,0}^{2}+\frac{F^{2}}{f_{a}^{2}}\left(C_{m}^{a}+\frac{{(C_{k}^{a\eta})}^{2}m_{a,0}^{4}}{m_{a,0}^{2}-\bar{m}_{\eta_{8}}^{2}}\right)\,,$
(15) $\displaystyle\bar{m}_{\eta}^{2}$ $\displaystyle=$
$\displaystyle\bar{m}_{\eta_{8}}^{2}+\frac{F^{2}}{f_{a}^{2}}\frac{{(C_{k}^{a\eta})}^{2}\bar{m}_{\eta_{8}}^{4}}{\bar{m}_{\eta_{8}}^{2}-m_{a,0}^{2}}\,.$
(16)

By keeping the terms up to $\mathcal{O}(1/f_{a}^{2})$, we can replace
$m_{a,0}$ in Eqs. (12) and (13) by $\bar{m}_{a}$

$\displaystyle
a=\left[1+\frac{1}{2}{(C_{k}^{a\eta})}^{2}\frac{\bar{m}_{a}^{2}(\bar{m}_{a}^{2}-2\bar{m}_{\eta_{8}}^{2})}{(\bar{m}_{a}^{2}-\bar{m}_{\eta_{8}}^{2})^{2}}\frac{F^{2}}{f_{a}^{2}}\right]\widetilde{a}-\frac{F}{f_{a}}C_{k}^{a\eta}\frac{\bar{m}_{\eta_{8}}^{2}}{\bar{m}_{\eta_{8}}^{2}-\bar{m}_{a}^{2}}\widetilde{\eta}+\mathcal{O}(1/f_{a}^{3})\,,$
(17)
$\displaystyle\eta_{8}=\left[1+\frac{1}{2}{(C_{k}^{a\eta})}^{2}\frac{\bar{m}_{\eta_{8}}^{2}(\bar{m}_{\eta_{8}}^{2}-2\bar{m}_{a}^{2})}{(\bar{m}_{a}^{2}-\bar{m}_{\eta_{8}}^{2})^{2}}\frac{F^{2}}{f_{a}^{2}}\right]\widetilde{\eta}-\frac{F}{f_{a}}C_{k}^{a\eta}\frac{\bar{m}_{a}^{2}}{\bar{m}_{a}^{2}-\bar{m}_{\eta_{8}}^{2}}\widetilde{a}+\mathcal{O}(1/f_{a}^{3})\,,$
(18)

due to the fact $\bar{m}_{a}^{2}=m_{a,0}^{2}+\mathcal{O}(1/f_{a}^{2})$.

The $\widetilde{a}$ and $\widetilde{\eta}$ fields in Eq. (14) are now in the
canonical form at LO. However, the $\widetilde{a}$ and $\widetilde{\eta}$
fields will get mixed again at NLO. Define the two-point one-particle-
irreducible (1PI) amplitudes for ALP and $\eta$ fields as
$\Sigma_{aa}(p^{2})$, $\Sigma_{a\eta}(p^{2})$, $\Sigma_{\eta\eta}(p^{2})$,
whose explicit expressions at the one-loop level are provided in Appendix A.
Up to $\mathcal{O}(p^{4})$ they can be decomposed as
$\Sigma_{ij}^{(4)}(p^{2})=p^{2}\Sigma_{ij}^{(4)^{\prime}}(0)+\Sigma_{ij}^{(4)}(0)$
with
$\Sigma_{ij}^{(4)^{\prime}}(0)=\mathrm{d}\Sigma_{ij}^{(4)}(p^{2})/\mathrm{d}p^{2}|_{p^{2}=0}$,
with $i,j=a,\eta$. In fact, since
$\mathrm{d}\Sigma_{ij}^{(4)}(p^{2})/\mathrm{d}p^{2}$ does not depend on
$p^{2}$ at the $\mathcal{O}(p^{4})$ level, we will simply denote
$\Sigma_{ij}^{(4)^{\prime}}(0)$ as $\Sigma_{ij}^{(4)^{\prime}}$ in the
following discussions. The number $n$ inside the bracket in the superscript
denotes the chiral order of the amplitude. The quadric form of the ALP and
$\eta$ fields up to NLO can be written as

$\displaystyle\mathcal{L}_{\mathrm{mix}}^{\mathrm{NLO}}=$
$\displaystyle\frac{1}{2}(1-\Sigma_{aa}^{(4)^{\prime}})\partial_{\mu}\widetilde{a}\partial^{\mu}\widetilde{a}+\frac{1}{2}(1-\Sigma_{\eta\eta}^{(4)^{\prime}})\partial_{\mu}\widetilde{\eta}\partial^{\mu}\widetilde{\eta}-\Sigma_{a\eta}^{(4)^{\prime}}\partial_{\mu}\widetilde{a}\partial^{\mu}\widetilde{\eta}$
(19)
$\displaystyle-\frac{1}{2}\left[\bar{m}_{a}^{2}+\Sigma_{aa}^{(4)}(0)\right]\widetilde{a}^{2}-\frac{1}{2}\left[\bar{m}_{\eta}^{2}+\Sigma_{\eta\eta}^{(4)}(0)\right]\widetilde{\eta}^{2}-\Sigma_{a\eta}^{(4)}(0)\widetilde{a}\widetilde{\eta}\,.$

We first eliminate the kinetic mixing up to NLO with the field redefinitions

$\widetilde{a}=\left(1+\frac{1}{2}\Sigma_{aa}^{(4)^{\prime}}\right)\hat{a}+\frac{1}{2}\Sigma_{a\eta}^{(4)^{\prime}}\hat{\eta}\,,\quad\widetilde{\eta}=\left(1+\frac{1}{2}\Sigma_{\eta\eta}^{(4)^{\prime}}\right)\hat{\eta}+\frac{1}{2}\Sigma_{a\eta}^{(4)^{\prime}}\hat{a}\,.$
(20)

By using the field variables $\hat{a}$ and $\hat{\eta}$, the kinetic terms can
be diagonalized to the canonical form

$\mathcal{L}_{\mathrm{mix}}^{\mathrm{NLO}}=\frac{1}{2}\partial_{\mu}\hat{a}\partial^{\mu}\hat{a}+\frac{1}{2}\partial_{\mu}\hat{\eta}\partial^{\mu}\hat{\eta}-\frac{1}{2}\hat{m}_{a}^{2}\hat{a}^{2}-\frac{1}{2}\hat{m}_{\eta}^{2}\hat{\eta}^{2}-\left[\frac{1}{2}(\bar{m}_{a}^{2}+\bar{m}_{\eta}^{2})\Sigma_{a\eta}^{(4)^{\prime}}+\Sigma_{a\eta}^{(4)}(0)\right]\hat{a}\hat{\eta}\,,$
(21)

with

$\displaystyle\hat{m}_{a}^{2}=\bar{m}_{a}^{2}+\Sigma_{aa}^{(4)^{\prime}}\bar{m}_{a}^{2}+\Sigma_{aa}^{(4)}(0)=\bar{m}_{a}^{2}+\Sigma_{aa}^{(4)}(\bar{m}_{a}^{2})\,,$
(22)
$\displaystyle\hat{m}_{\eta}^{2}=\bar{m}_{\eta}^{2}+\Sigma_{\eta\eta}^{(4)^{\prime}}\bar{m}_{\eta}^{2}+\Sigma_{\eta\eta}^{(4)}(0)=\bar{m}_{\eta}^{2}+\Sigma_{\eta\eta}^{(4)}(\bar{m}_{\eta}^{2})\,.$
(23)

Next we deal with the $\hat{a}$-$\hat{\eta}$ mass mixing via a perturbative
real orthogonal transformation

$\hat{a}=a_{\mathrm{phy}}+\delta_{m}\,\eta_{\mathrm{phy}}\,,\quad\hat{\eta}=\eta_{\mathrm{phy}}-\delta_{m}\,a_{\mathrm{phy}}\,,$
(24)

with

$\delta_{m}=-\frac{\Sigma_{a\eta}^{(4)^{\prime}}\frac{1}{2}(\bar{m}_{a}^{2}+\bar{m}_{\eta}^{2})+\Sigma_{a\eta}^{(4)}(0)}{\hat{m}_{a}^{2}-\hat{m}_{\eta}^{2}}\,.$
(25)

Combining Eqs. (20) and (25) one can get the relations between field variables
$\widetilde{a}$, $\widetilde{\eta}$ and the NLO renormalized fields
$a_{\mathrm{phy}}$ and $\eta_{\mathrm{phy}}$,

$\widetilde{a}=\left(1+\frac{1}{2}\Sigma_{aa}^{(4)^{\prime}}\right)a_{\mathrm{phy}}+\frac{\Sigma_{a\eta}^{(4)}(m_{\eta}^{2})}{m_{\eta}^{2}-m_{a}^{2}}\eta_{\mathrm{phy}}\,,\quad\widetilde{\eta}=\left(1+\frac{1}{2}\Sigma_{\eta\eta}^{(4)^{\prime}}\right)\eta_{\mathrm{phy}}+\frac{\Sigma_{a\eta}^{(4)^{\prime}}(m_{a}^{2})}{m_{a}^{2}-m_{\eta}^{2}}a_{\mathrm{phy}}\,.$
(26)

Up to NLO $a_{\mathrm{phy}}$ and $\eta_{\mathrm{phy}}$ become the canonical
fields

$\mathcal{L}_{\mathrm{mix}}^{\mathrm{NLO}}=\frac{1}{2}\partial_{\mu}a_{\mathrm{phy}}\partial^{\mu}a_{\mathrm{phy}}-\frac{1}{2}m_{a}^{2}a_{\mathrm{phy}}^{2}+\frac{1}{2}\partial_{\mu}\eta_{\mathrm{phy}}\partial^{\mu}\eta_{\mathrm{phy}}-\frac{1}{2}m_{\eta}^{2}\eta_{\mathrm{phy}}^{2}\,,$
(27)

where $m_{a}$ and $m_{\eta}$ denote the physical masses of the ALP and $\eta$.
It is straightforward to verify that the mixing terms in the field
redefinitions (24) only affect the masses of ALP and $\eta$ at the level of
$\mathcal{O}(p^{6})$. As a result, the physical masses $m_{a}$ and $m_{\eta}$
in Eq. (27) simply equal to $\hat{m}_{a}$ and $\hat{m}_{\eta}$ in Eq. (21),
respectively. The explicit expression of the physical ALP mass up to NLO takes
the form

$\displaystyle
m_{a}^{2}=\bar{m}_{a}^{2}+\Sigma_{aa}^{(4)}(\bar{m}_{a}^{2})\,,$ (28)

where the LO ALP mass $\bar{m}_{a}^{2}$ is given in Eq. (15) and the NLO
chiral correction reads

$\displaystyle\Sigma_{aa}^{(4)}(\bar{m}_{a}^{2})=$
$\displaystyle\frac{1}{f_{a}^{2}}\Bigg{\\{}\frac{2L_{6}^{r}m_{\pi}^{2}(3m_{\eta}^{2}-m_{\pi}^{2})(m_{\eta}^{2}+m_{\pi}^{2})}{m_{\eta}^{2}}+\frac{2(3L_{7}^{r}+L_{8}^{r})m_{\pi}^{4}(3m_{\eta}^{2}-m_{\pi}^{2})^{2}}{3m_{\eta}^{4}}$
(29)
$\displaystyle-\frac{m_{\pi}^{2}(3m_{\eta}^{2}-m_{\pi}^{2})^{2}}{24m_{\eta}^{4}}\frac{m_{\pi}^{2}}{(4\pi)^{2}}\log\frac{m_{\pi}^{2}}{\mu^{2}}-\frac{m_{\pi}^{2}(9m_{\eta}^{4}-m_{\pi}^{4})}{36m_{\eta}^{4}}\frac{m_{K}^{2}}{(4\pi)^{2}}\log\frac{m_{K}^{2}}{\mu^{2}}$
$\displaystyle-\frac{m_{\pi}^{2}(3m_{\eta}^{2}-m_{\pi}^{2})(m_{\eta}^{2}+m_{\pi}^{2})}{72m_{\eta}^{4}}\frac{m_{\eta}^{2}}{(4\pi)^{2}}\log\frac{m_{\eta}^{2}}{\mu^{2}}\Bigg{\\}}$
$\displaystyle+\frac{X}{f_{a}^{2}}\Bigg{\\{}\frac{4\sqrt{3}L_{4}^{r}\bar{m}_{a}^{2}(m_{\eta}^{4}-m_{\pi}^{4})}{m_{\eta}^{2}}+\frac{8L_{5}^{r}\bar{m}_{a}^{2}(m_{\eta}^{2}-m_{\pi}^{2})}{\sqrt{3}}$
$\displaystyle+\frac{8\sqrt{3}(3L_{7}^{r}+L_{8}^{r})m_{\pi}^{2}(3m_{\eta}^{4}-4m_{\eta}^{2}m_{\pi}^{2}+m_{\pi}^{4})}{3m_{\eta}^{2}}$
$\displaystyle+\frac{m_{\pi}^{2}(3m_{\eta}^{2}-m_{\pi}^{2})}{2\sqrt{3}m_{\eta}^{2}}\frac{m_{\pi}^{2}}{(4\pi)^{2}}\log\frac{m_{\pi}^{2}}{\mu^{2}}-\frac{3m_{\eta}^{2}m_{\pi}^{2}-m_{\pi}^{4}}{6\sqrt{3}m_{\eta}^{2}}\frac{m_{\eta}^{2}}{(4\pi)^{2}}\log\frac{m_{\eta}^{2}}{\mu^{2}}$
$\displaystyle-\frac{\left[3m_{\eta}^{2}m_{\pi}^{2}-m_{\pi}^{4}+6\bar{m}_{a}^{2}(m_{\eta}^{2}-m_{\pi}^{2})\right]}{3\sqrt{3}m_{\eta}^{2}}\frac{m_{K}^{2}}{(4\pi)^{2}}\log\frac{m_{K}^{2}}{\mu^{2}}\Bigg{\\}}$
$\displaystyle+\frac{X^{2}}{f_{a}^{2}}\Bigg{\\{}-12L_{4}^{r}\bar{m}_{a}^{2}(m_{\eta}^{2}+m_{\pi}^{2})-8L_{5}^{r}\bar{m}_{a}^{2}m_{\eta}^{2}+24L_{6}^{r}m_{\eta}^{2}(m_{\eta}^{2}+m_{\pi}^{2})$
$\displaystyle+24L_{7}^{r}(m_{\eta}^{2}-m_{\pi}^{2})^{2}+8L_{8}^{r}(3m_{\eta}^{4}-2m_{\eta}^{2}m_{\pi}^{2}+m_{\pi}^{4})-\frac{m_{\pi}^{2}}{2}\frac{m_{\pi}^{2}}{(4\pi)^{2}}\log\frac{m_{\pi}^{2}}{\mu^{2}}$
$\displaystyle+\frac{3\bar{m}_{a}^{2}+m_{\pi}^{2}}{3}\frac{m_{K}^{2}}{(4\pi)^{2}}\log\frac{m_{K}^{2}}{\mu^{2}}-\frac{(4m_{\eta}^{2}-m_{\pi}^{2})}{6}\frac{m_{\eta}^{2}}{(4\pi)^{2}}\log\frac{m_{\eta}^{2}}{\mu^{2}}\Bigg{\\}}\,,$

with
$X=C_{k}^{a\eta}\frac{\bar{m}_{a}^{2}}{\bar{m}_{a}^{2}-\bar{m}_{\eta_{8}}^{2}}$.
In the QCD axion scenario with $m_{a,0}^{2}=0$, the LO mass of QCD axion is
given by $\bar{m}_{a}^{2}=(F^{2}/f_{a}^{2})C_{m}^{a}$. Then, the last two
terms proportional to $X$ and $X^{2}$ in the expression of
$\Sigma_{aa}^{(4)}(\bar{m}_{a}^{2})$ in Eq. (29) can be ignored up to
$\mathcal{O}(1/f_{a}^{2})$, and thus one can just keep the first term in Eq.
(29) when considering the QCD axion case. The NLO correction to the QCD axion
mass given here is in consistent with the result given in Ref. [26] in the
isospin limit. In the following, we will use the NLO renormalized canonical
fields $a_{\mathrm{phy}}$ and $\eta_{\mathrm{phy}}$ to perform the
perturbative calculation for the $\eta\to\pi\pi a$ decay.

##  4 $\eta\to\pi\pi a$ decay amplitude up to NLO and its unitarization

We are interested in the $\eta\to\pi\pi a$ decay with the ALP mass $m_{a}$
varying from $0$ to $m_{\eta}-2m_{\pi}$. The Mandelstam variables to describe
the kinematics of the $\eta(P_{\eta})\to\pi(p_{1})\pi(p_{2})a(p_{a})$ decay
are defined as

$\displaystyle s=(P_{\eta}-p_{a})^{2}=(p_{1}+p_{2})^{2}=m_{\pi\pi}^{2}\,,$
(30) $\displaystyle t=(P_{\eta}-p_{1})^{2}=(p_{a}+p_{2})^{2}=m_{a\pi}^{2}\,,$
$\displaystyle u=(P_{\eta}-p_{2})^{2}=(p_{a}+p_{1})^{2}\,,$

which satisfy the constraint $s+t+u=m_{\eta}^{2}+m_{a}^{2}+2m_{\pi}^{2}$. We
consider the $\eta\to\pi\pi a$ decay at the level of $\mathcal{O}(1/f_{a})$
and calculate the chiral corrections of the decay amplitudes up to
$\mathcal{O}(p^{4})$. According to Eq. (28), the difference between the
physical ALP mass $m_{a}$ and the bare one $m_{a,0}$ lies at the level of
$\mathcal{O}(1/f_{a}^{2})$, i.e.
$m_{a}^{2}=m_{a,0}^{2}+\mathcal{O}(1/f_{a}^{2})$, which indicates that the
replacement of $m_{a,0}$ with $m_{a}$ in the $\eta\to\pi\pi a$ amplitude will
only give effects at $\mathcal{O}(1/f_{a}^{3})$, beyond the
$\mathcal{O}(1/f_{a})$ precision we are considering in this work. Therefore,
in the following we do not need to distinguish the bare mass $m_{a,0}$ and the
physical mass $m_{a}$ for the ALP, both of which will be simply denoted as
$m_{a}$.

Up to $\mathcal{O}(p^{4})$ the decay amplitudes receive the contributions from
the Feynman diagrams shown in Fig. 1. The decay amplitudes for
$\eta\to\pi^{+}\pi^{-}a$ and $\eta\to\pi^{0}\pi^{0}a$ turn out to be the same
up to $\mathcal{O}(p^{4})$ in isospin limit. The LO part of the decay
amplitude has no energy dependence and its explicit expression is

$\mathcal{M}_{\eta;\pi\pi
a}^{(2)}=\left[\frac{2}{\sqrt{3}}+\left(1-\frac{m_{a}^{2}}{m_{a}^{2}-m_{\eta}^{2}}\right)\frac{m_{\eta}^{2}-m_{\pi}^{2}}{\sqrt{3}m_{\eta}^{2}}\right]\frac{m_{\pi}^{2}}{6f_{a}F_{\pi}}\,.$
(31)

We have replaced the LO quantities $F$, $\bar{m}_{\pi}^{2}$ and
$\bar{m}_{\eta_{8}}^{2}$ by the physical ones in the above equation, which is
legitimate at LO accuracy. While, up to NLO, this will bring an extra
$\mathcal{O}(p^{4})$ correction into the decay amplitude

$\displaystyle\mathcal{M}_{\eta;\pi\pi
a}^{(2)}\left(\delta_{F_{\pi}}^{(4)}-\frac{\delta_{m_{\pi}^{2}}^{(4)}}{m_{\pi}^{2}}\right)+\left(1-\frac{m_{a}^{2}}{m_{a}^{2}-m_{\eta}^{2}}\right)\left(\frac{m_{\eta}^{2}-m_{\pi}^{2}}{\sqrt{3}m_{\eta}^{2}}\frac{\delta_{m_{\eta}^{2}}^{(4)}}{m_{\eta}^{2}}-\frac{\delta_{m_{\eta}^{2}}^{(4)}-\delta_{m_{\pi}^{2}}^{(4)}}{\sqrt{3}m_{\eta}^{2}}\right)\frac{m_{\pi}^{2}}{6f_{a}F_{\pi}}$
(32)
$\displaystyle+\frac{m_{a}^{2}}{m_{a}^{2}-m_{\eta}^{2}}\frac{\delta_{m_{\eta}^{2}}^{(4)}}{m_{a}^{2}-m_{\eta}^{2}}\frac{m_{\eta}^{2}-m_{\pi}^{2}}{\sqrt{3}m_{\eta}^{2}}\frac{m_{\pi}^{2}}{6f_{a}F_{\pi}}\,,$

where $\delta_{F_{\pi}}^{(4)}$, $\delta_{m_{\eta}^{2}}^{(4)}$ and
$\delta_{m_{\pi}^{2}}^{(4)}$ denote the $\mathcal{O}(p^{4})$ corrections to
the pion mass, decay constant $F_{\pi}$ and the $\eta$ mass, respectively.
They have been worked out long time ago in the seminal paper of $SU(3)$ [23]
and their expressions are provided here for the sake of completeness

$\displaystyle F_{\pi}=F\left(1+\delta_{F_{\pi}}^{(4)}\right)\,,\quad
m_{\pi}^{2}=\bar{m}_{\pi}^{2}+\delta_{m_{\pi}^{2}}^{(4)}\,,\quad
m_{\eta}^{2}=\bar{m}_{\eta_{8}}^{2}+\delta_{m_{\eta}^{2}}^{(4)}\,,$ (33)
$\displaystyle\delta_{F_{\pi}}^{(4)}=\frac{2}{F_{\pi}^{2}}\left[3L_{4}^{r}(m_{\pi}^{2}+m_{\eta}^{2})+2L_{5}^{r}m_{\pi}^{2}\right]-\frac{m_{\pi}^{2}}{(4\pi
F_{\pi})^{2}}\log\frac{m_{\pi}^{2}}{\mu^{2}}-\frac{m_{K}^{2}}{2(4\pi
F_{\pi})^{2}}\log\frac{m_{K}^{2}}{\mu^{2}}\,,$ (34)
$\displaystyle\begin{aligned}
\delta_{m_{\pi}^{2}}^{(4)}=&m_{\pi}^{2}\Big{\\{}\frac{m_{\pi}^{2}}{2(4\pi
F_{\pi})^{2}}\log\frac{m_{\pi}^{2}}{\mu^{2}}-\frac{m_{\eta}^{2}}{6(4\pi
F_{\pi})^{2}}\log\frac{m_{\eta}^{2}}{\mu^{2}}-\frac{4}{F_{\pi}^{2}}\big{[}2(L_{5}^{r}-2L_{8}^{r})m_{\pi}^{2}\\\
&+3(L_{4}^{r}-2L_{6}^{r})(m_{\pi}^{2}+m_{\eta}^{2})\big{]}\Big{\\}}\,,\end{aligned}$
(35) $\displaystyle\begin{aligned}
\delta_{m_{\eta}^{2}}^{(4)}=&m_{\eta}^{2}\Big{\\{}\frac{m_{K}^{2}}{(4\pi
F_{\pi})^{2}}\log\frac{m_{K}^{2}}{\mu^{2}}-\frac{2}{3}\frac{m_{\eta}^{2}}{(4\pi
F_{\pi})^{2}}\log\frac{m_{\eta}^{2}}{\mu^{2}}-\frac{4}{F_{\pi}^{2}}\big{[}2(L_{5}^{r}-2L_{8}^{r})m_{\eta}^{2}\\\
&+3(L_{4}^{r}-2L_{6}^{r})(m_{\pi}^{2}+m_{\eta}^{2})\big{]}\Big{\\}}-m_{\pi}^{2}\Big{\\{}\frac{m_{\pi}^{2}}{2(4\pi
F_{\pi})^{2}}\log\frac{m_{\pi}^{2}}{\mu^{2}}-\frac{m_{K}^{2}}{3(4\pi
F_{\pi})^{2}}\log\frac{m_{K}^{2}}{\mu^{2}}\\\ &-\frac{m_{\eta}^{2}}{6(4\pi
F_{\pi})^{2}}\log\frac{m_{\eta}^{2}}{\mu^{2}}\Big{\\}}+\frac{8}{F_{\pi}^{2}}(3L_{7}^{r}+L_{8}^{r})\left(m_{\eta}^{2}-m_{\pi}^{2}\right)^{2}\,.\end{aligned}$
(36)

(a) (b) (c) (d) (e) (f)

Figure 1: Feynman diagrams for the $\eta\to\pi\pi a$ decay up to NLO. The
numbers in the vertices denote their chiral orders.

The decay amplitude acquires energy dependence at NLO. As shown in Fig. 1, the
NLO part of the decay amplitude receives the contributions from contact
diagram in Fig. 1(b) by taking the $\mathcal{O}(p^{4})$ vertices, and the loop
diagrams in Figs. 1(c)-1(f) by taking the $\mathcal{O}(p^{2})$ vertices, where
the intermediate states can be $\pi\pi,\,K\bar{K},\,\eta\eta$ for the
$s$-channel diagram in Fig. 1(d), and $\pi\eta,\,K\bar{K}$ for the $t$\- and
$u$-channel diagrams in Figs. 1(e) and 1(f). By plugging the field
redefinitions of Eq. (26) into the $\mathcal{O}(p^{2})$ Lagrangian, one gets
extra NLO contributions for the decay amplitude, which can be recognized as
the corrections from the wave function renormalization of $\eta$ and also the
NLO ALP-$\eta$ mixing. In addition, the wave function renormalization constant
for the pion is also needed in the calculation, while the ALP wave function
renormalization is neglected since it behaves as
$Z_{a}=1+\mathcal{O}(1/f_{a}^{2})$. In Appendix A, we give the full
$\mathcal{O}(p^{4})$ part of the $\eta\to\pi\pi a$ amplitude, which
incorporates the correction of Eq. (32). We have explicitly checked that the
divergences originated from the chiral loops in the $\eta\to\pi\pi a$
amplitude can be completely absorbed by the renormalization of the $SU(3)$
LECs determined in the seminal paper [23]. This also provides a way to cross
check about the correctness of the decay amplitudes.

Being the momentum expansion feature of , its description of the $\pi\pi$
system up to the $\eta$ mass in the $\eta\to\pi\pi a$ decay could be
insufficient, due to the strong $\pi\pi$ final-state interaction. The
unitarization procedure offers an efficient way to extend the applicable
energy range of the perturbative results. While for the $a\pi$ subsystem, the
result up $\mathcal{O}(p^{4})$, as discussed in detail later, indicates that
its final-state interaction could be perturbatively treated as a good
approximation, at least in the energy region up to $m_{\eta}-m_{\pi}$. As a
result, we will focus on the unitarization of the $\pi\pi$ system in the
$\eta\to\pi\pi a$ process. The unitarization can be conveniently implemented
through the partial-wave (PW) amplitude, which is obtained from the PW
projection

$\mathcal{M}_{\eta;\pi\pi
a}^{IJ}(s)=\frac{1}{2(\sqrt{2})^{N}}\int_{-1}^{+1}\mathrm{d}\cos\theta
P_{J}(\cos\theta)\mathcal{M}_{\eta;\pi\pi a}^{I}(s,t,u)\,,$ (37)

where $I$ and $J$ are the quantum numbers of the isospin and angular momentum
of the $\pi\pi$ system respectively, $\theta$ denotes the scattering angle in
the center of mass (CM) frame and $P_{J}(\cos\theta)$ corresponds to the
Legendre polynomial. The term $(\sqrt{2})^{N}$ stands for the symmetry factor,
with $N$ the number of the pairs of identical particles. $N=1$ in the
$\eta\to\pi\pi a$ case. The Mandelstam variables $t$ and $u$ in the CM frame
are given by

$\displaystyle
t=\frac{1}{2}(m_{\eta}^{2}+2m_{\pi}^{2}+m_{a}^{2}-s)-\frac{\sqrt{\lambda(m_{\eta}^{2},s,m_{a}^{2})\lambda(s,m_{\pi}^{2},m_{\pi}^{2})}}{2s}\cos\theta\,,$
(38) $\displaystyle
u=\frac{1}{2}(m_{\eta}^{2}+2m_{\pi}^{2}+m_{a}^{2}-s)+\frac{\sqrt{\lambda(m_{\eta}^{2},s,m_{a}^{2})\lambda(s,m_{\pi}^{2},m_{\pi}^{2})}}{2s}\cos\theta\,,$

with $\lambda(x,y,z)=x^{2}+y^{2}+z^{2}-2xy-2xz-2yz$. The $P$-wave $\pi\pi$
amplitude vanishes in the $\eta\to\pi\pi a$ process. The predominant
contribution comes from the $S$-wave amplitude, and the $D$-wave effect is
expected to be small in the kinematical region of the $\eta\to\pi\pi a$ decay.
Consequently, we will perform the unitarization of the $S$-wave $\pi\pi$
amplitude, and directly include the $D$-wave projection from the perturbative
results. Since we work in the isospin limit, the only possible isospin of the
two pions in $\eta\to\pi\pi a$ process is $I=0$.

Up to the $\eta$ mass, the final-state $\pi\pi$ interaction is well
approximated by including the single-channel unitarity effect and we propose
the following recipe [27] to perform the unitarization for the $S$-wave
$\pi\pi$ system

$\mathcal{M}_{\eta;\pi\pi a}^{00,{\rm Uni}}(s)=\frac{\mathcal{M}_{\eta;\pi\pi
a}^{00,{\rm L}}(s)}{1-G_{\pi\pi}(s)T_{\pi\pi\to\pi\pi}^{00,(2)}(s)}\,,$ (39)

where the quantity $\mathcal{M}_{\eta;\pi\pi a}^{00,{\rm L}}(s)$ excludes the
$s$-channel unitarity cut from the $\pi\pi$ case and contains the left hand
cuts and contact terms of the perturbative $\eta\to\pi\pi a$ amplitude,
$T_{\pi\pi\to\pi\pi}^{00,(2)}(s)$ stands for the LO $\pi\pi$ scattering
amplitude with $IJ=00$, and $G_{\pi\pi}(s)$ is the two-point one-loop function
evaluated in the dimensional regularization with $\overline{MS}-1$ subtraction
scheme. The explicit expression of the latter quantity is

$G_{\pi\pi}(s)=-\frac{1}{(4\pi)^{2}}\left(\log\frac{m_{\pi}^{2}}{\mu^{2}}-\sigma_{\pi}(s)\log\frac{\sigma_{\pi}(s)-1}{\sigma_{\pi}(s)+1}-1\right)\,,$
(40)

whose imaginary part is

$\mathrm{Im}G_{\pi\pi}(s)=\theta(s-4m_{\pi}^{2})\rho_{\pi\pi}(s)\,,$ (41)

with

$\rho_{\pi\pi}(s)=\frac{\sigma_{\pi}(s)}{16\pi}\,,\quad\sigma_{\pi}(s)=\sqrt{1-\frac{4m_{\pi}^{2}}{s}}\,.$
(42)

The LO PW $\pi\pi$ scattering amplitude in the convention of Eq. (37) with
$IJ=00$ reads

$T_{\pi\pi\to\pi\pi}^{00,(2)}(s)=\frac{2s-m_{\pi}^{2}}{2F_{\pi}^{2}}\,.$ (43)

The quantity $\mathcal{M}_{\eta;\pi\pi a}^{00,{\rm L}}(s)$ in Eq. (39) can be
determined by matching the right-hand side of this equation with the
perturbative result, and the explicit expression is given by

$\mathcal{M}_{\eta;\pi\pi a}^{00,{\rm L}}(s)=\mathcal{M}_{\eta;\pi\pi
a}^{00,(2)}(s)+\mathcal{M}_{\eta;\pi\pi
a}^{00,(4)}(s)-G_{\pi\pi}(s)\,\mathcal{M}_{\eta;\pi\pi
a}^{00,(2)}(s)\,T_{\pi\pi\to\pi\pi}^{00,(2)}(s)\,,$ (44)

which satisfies

${\rm Im}\mathcal{M}_{\eta;\pi\pi a}^{00,{\rm
L}}(s)=0\,,\qquad(2m_{\pi}<\sqrt{s}<2m_{K})\,,$ (45)

since the perturbative amplitude $\mathcal{M}_{\eta;\pi\pi a}^{00,(4)}(s)$
obeys

${\rm Im}\mathcal{M}_{\eta;\pi\pi
a}^{00,(4)}(s)=\rho_{\pi\pi}(s)\mathcal{M}_{\eta;\pi\pi
a}^{00,(2)}(s)\,T_{\pi\pi\to\pi\pi}^{00,(2)}(s)\,,\qquad(2m_{\pi}<\sqrt{s}<2m_{K})\,.$
(46)

It is easy to verify that the unitarized decay amplitude in Eq. (39) strictly
fulfills the unitarity relation

${\rm Im}\mathcal{M}_{\eta;\pi\pi a}^{00,{\rm
Uni}}(s)=\rho_{\pi\pi}\mathcal{M}_{\eta;\pi\pi a}^{00,{\rm
Uni}}(s)\left(T_{\pi\pi\to\pi\pi}^{00,{\rm
Uni}}(s)\right)^{*}\,,\qquad(2m_{\pi}<\sqrt{s}<2m_{K})\,,$ (47)

when taking the unitarized $\pi\pi$ amplitude as

$T_{\pi\pi\to\pi\pi}^{00,{\rm
Uni}}(s)=\frac{T_{\pi\pi\to\pi\pi}^{00,(2)}(s)}{1-G_{\pi\pi}(s)T_{\pi\pi\to\pi\pi}^{00,(2)}(s)}\,,$
(48)

which is a generalized $K$-matrix amplitude. The common quantity of
$1/[1-G_{\pi\pi}(s)T_{\pi\pi\to\pi\pi}^{00,(2)}(s)]$ in Eqs. (39) and (48)
practically plays the role of parameterizing the strong $\pi\pi$ final-state
interactions. It will be demonstrated later that the unitarized $\pi\pi$
amplitude can describe the $\pi\pi$ phase shifts with $IJ=00$ very well up to
around 0.7 GeV, implying that the neat forms in Eqs. (39) and (48) can capture
the key the unitarization effects in the $\pi\pi$ channel.

For the $D$-wave component, we will simply take the result from the PW
projection of the perturbative amplitude in Eq. (37). To combine the
unitarized $S$-wave amplitude and the perturbative $D$-wave amplitude, the
$\eta\to\pi\pi a$ decay amplitude can be reconstructed as

$\mathcal{M}_{\eta;\pi^{+}\pi^{-}a}^{I=0,{\rm
CPW}}(s,t,u)=\sqrt{2}\left[\mathcal{M}_{\eta;\pi\pi a}^{00,{\rm
Uni}}(s)+5P_{2}(\cos\theta)\mathcal{M}_{\eta;\pi\pi a}^{02}(s)\right]\,,$ (49)

where the factor $\sqrt{2}$ is introduced to recover the physical amplitude
from the PW convention in Eq. (37). From the difference of the perturbative
full amplitude $\mathcal{M}_{\eta;\pi\pi a}^{I=0}(s,t,u)$ given in the
Appendix and the combined partial-wave (CPW) decay amplitude
$\mathcal{M}_{\eta;\pi\pi a}^{I=0,{\rm CPW}}(s,t,u)$ in Eq. (49), one could
study how the strong final-state $\pi\pi$ interactions can affect the
$\eta\to\pi\pi a$ process.

##  5 Phenomenological discussions

As mentioned previously, the isospin breaking effects enter the $\eta\to\pi\pi
a$ process at the level of $\mathcal{O}(\epsilon_{I}^{2})$ and hence are
neglected in this work. In the isospin limit, the decay amplitudes for
$\eta\to\pi^{+}\pi^{-}a$ and $\eta\to\pi^{0}\pi^{0}a$ turn out to be the same
up to NLO. The decay rate for the $\pi^{0}\pi^{0}$ final state only differs
with that of $\pi^{+}\pi^{-}$ by an identical-particle factor $1/2$. So we
only show the results for the $\eta\to\pi^{+}\pi^{-}a$ decay in the following.

In the phenomenological discussions, we take $m_{\pi}=137~{}$MeV,
$m_{K}=496~{}$MeV, and $m_{\eta}=548~{}$MeV. The value of the pion decay
constant is taken as $F_{\pi}=92.1~{}$MeV [28]. Regarding the renormalized
LECs $L_{i}^{r}(\mu)$ at the renormalization scale $\mu=770~{}$MeV, we take
the values fitted at NLO given in Table 6 of Ref. [21], which are
$L_{1}^{r}=1.0(1)\times 10^{-3}$, $L_{2}^{r}=1.6(2)\times 10^{-3}$,
$L_{3}^{r}=-3.8(3)\times 10^{-3}$, $L_{4}^{r}=0.0(3)\times 10^{-3}$,
$L_{5}^{r}=1.2(1)\times 10^{-3}$, $L_{6}^{r}=0.0(4)\times 10^{-3}$,
$L_{7}^{r}=-0.3(2)\times 10^{-3}$ and $L_{8}^{r}=0.5(2)\times 10^{-3}$. To
estimate the uncertainty of the $\eta\to\pi\pi a$ amplitude, we will generate
a large amount of parameter configurations by randomly sampling all the
$L_{i=1,\cdots,8}^{r}$ with the assumptions that they are Gaussianly
distributed and not correlated with each other. Since the latter assumption
clearly introduces unrealistically large uncertainty, in another scenario we
will simply neglect the uncertainties of the LECs of $L_{4,6,7}^{r}$ and only
consider the error bars of $L_{1,2,3,5,8}^{r}$ in the random sampling
procedure. In later discussions, the uncertainty analysis by taking the error
bars of all the $L_{i=1,\cdots,8}^{r}$ in the random sampling will be named as
Set-I, while the uncertainty analysis by keeping the error bars of
$L_{1,2,3,5,8}^{r}$ in the random sampling procedure will be denoted as Set-
II.

The differential partial decay rate for the
$\eta(P_{\eta})\to\pi(p_{1})\pi(p_{2})a(p_{a})$ process reads

$\frac{\mathrm{d}^{2}\Gamma_{\eta;\pi\pi{a}}}{\mathrm{d}s\mathrm{d}t}=\frac{1}{32(2\pi)^{3}\mathcal{N}m_{\eta}^{3}}\left|\mathcal{M}_{\eta;\pi\pi{a}}\right|^{2}\,,$
(50)

with $\mathcal{N}=1$ for $\pi^{+}\pi^{-}$ channel and $\mathcal{N}=2$ for
$\pi^{0}\pi^{0}$ channel. The boundaries of the physical regions for $s$ and
$t$ are

$\displaystyle s_{\text{min}}=4m_{\pi}^{2}\,,\quad
s_{\text{max}}=(m_{\eta}-m_{a})^{2}\,,$ (51) $\displaystyle
t_{\text{max/min}}(s)=\frac{1}{2}\Bigg{[}m_{\eta}^{2}+m_{a}^{2}+2m_{\pi}^{2}-s\pm\frac{\sqrt{\lambda(s,m_{\eta}^{2},m_{a}^{2})\lambda(s,m_{\pi}^{2},m_{\pi}^{2})}}{s}\Bigg{]}\,,$
(52)

or

$\displaystyle t_{\mathrm{min}}=(m_{\pi}+m_{a})^{2}\,,\quad
t_{\mathrm{max}}=(m_{\eta}-m_{\pi})^{2}\,,$ (53) $\displaystyle\begin{aligned}
s_{\mathrm{max/min}}(t)=&\frac{m_{\eta}^{2}+m_{a}^{2}+2m_{\pi}^{2}-t}{2}+\frac{1}{2t}(m_{\eta}^{2}-m_{\pi}^{2})(m_{\pi}^{2}-m_{a}^{2})\pm\frac{\sqrt{\lambda(t,m_{\eta}^{2},m_{\pi}^{2})\lambda(t,m_{\pi}^{2},m_{a}^{2})}}{2t}\,.\end{aligned}$
(54)

Figure 2: $10^{6}\times f_{a}^{2}\mathrm{d}^{2}\Gamma/\mathrm{d}s\mathrm{d}t$
(see Eq. (50)) in unit of $\mathrm{GeV}^{-1}$. The results correspond to
taking all the three terms in the right hand side of Eq. (55), i.e. the
$\mathrm{NLO^{\prime}}$ case. See the text for details.

###  5.1 Analyses by taking the perturbative $\chi$PT amplitudes

In this part, we focus on the situation by taking the full perturbative
$\chi$PT amplitudes $\mathcal{M}_{\eta;\pi\pi a}^{(2+4)}(s,t,u)$, as given in
Eq. (65). When perturbatively including the NLO corrections in the decay
amplitude, we can write the amplitude squared as

$|\mathcal{M}|^{2}=\left(\mathcal{M}^{(2)}\right)^{2}+2\mathcal{M}^{(2)}\mathrm{Re}\left(\mathcal{M}^{(4)}\right)+|\mathcal{M}^{(4)}|^{2}\,,$
(55)

where $n$ inside the bracket in the superscript stands for the chiral order of
the amplitude. The well established power counting in is analyzed in most of
the cases for the amplitude itself, instead of the decay width that takes the
amplitude squared as input. By applying the chiral power counting to the
amplitude squared, the $|\mathcal{M}^{(4)}|^{2}$ term in Eq. (55) belongs to
the $O(p^{6})$ part, which would have the same order as
$\mathcal{M}^{(2)}\mathrm{Re}\left(\mathcal{M}^{(6)}\right)$. The latter term
requires the calculation of the next-to-next-to-leading order amplitude, which
is beyond the scope of the present study. However, by applying the chiral
power counting rule to the amplitude itself, one should take all the three
terms in Eq. (55), which guarantees the positivity of the amplitude squared by
construction. In this work, we will explore the consequences of the two
scenarios by taking the chiral power counting rule both to the amplitude
itself and the amplitude squared. The difference from the two scenarios can be
regarded as a measure of the higher order effects. We designate the decay rate
by keeping the first two terms in the amplitude squared in Eq. (55) as
$\Gamma_{\mathrm{NLO}}$, and denote $\Gamma_{\mathrm{NLO^{\prime}}}$ when
taking all the three terms in the previous equation. The Dalitz decay
distributions from the $\mathrm{NLO^{\prime}}$ scenario are shown in Fig. 2,
where we take three values of $m_{a}$: $m_{a}=0~{}$MeV, $m_{a}=0.5m_{\pi}$,
$m_{a}=m_{\pi}$, for illustrations. The first lesson we can learn is that the
larger the ALP mass the smoother of the Dalitz distributions result. The
$\eta\to\pi\pi a$ amplitude up to NLO becomes roughly flat in the
$m_{\pi\pi}^{2}$-$m_{a\pi}^{2}$ plane for $m_{a}\gtrsim m_{\pi}$. Another
important lesson is that almost for any given value of $m_{\pi\pi}^{2}$ the
distribution along the $m_{a\pi}^{2}$ direction is rather smooth, while in
contrast for a specific $m_{a\pi}^{2}$ value the distribution along the
$m_{\pi\pi}^{2}$ direction can be uneven. This gives hints that the $a\pi$
interaction strength is weak in the kinematical region allowed in the
$\eta\to\pi\pi a$ process and the $\pi\pi$ final-state interaction is
stronger, which validates the discussions in the previous section that focuses
on the $\pi\pi$ final-state interactions. It is also possible that in a
broader energy region there are resonant contributions, such as $a_{0}(980)$
and $a_{0}(1450)$, to the $a\pi$ channel, which are however distant from the
kinematical region of the $\eta\to\pi\pi a$ decay. This latter scenario could
also give an even distribution along the $m_{a\pi}^{2}$ direction in the
$\eta\to\pi\pi a$ decay.

Figure 3: The ratios of $|\mathcal{M}|_{\rm
NLO^{\prime}}^{2}/|\mathcal{M}|_{\rm LO}^{2}$, where $|\mathcal{M}|_{\rm
LO}^{2}$ and $|\mathcal{M}|_{\rm NLO^{\prime}}^{2}$ stand for the first and
last two terms in Eq. (55), respectively.

The ratios between the higher order contribution and the LO part of the
amplitude squared (55) in the $m_{\pi\pi}^{2}$-$m_{a\pi}^{2}$ plane are
visualized in Fig. 3, from which one can have an intuitive idea about the
convergence of the chiral expansion of the amplitude squared
$|\mathcal{M}|^{2}$. The fact that the contour lines are roughly parallel to
the $m_{a\pi}^{2}$ axis semi-quantitatively illustrates the even distributions
of the amplitude squared along the $m_{a\pi}^{2}$ direction. For $m_{a}\simeq
m_{\pi}$, the $\mathrm{NLO^{\prime}}$ correction to the amplitude squared,
with respect to the LO result, lies around at the level of $20\%\sim 40\%$ in
the whole $m_{\pi\pi}^{2}$-$m_{a\pi}^{2}$ plane. In contrast, for the QCD
axion case with $m_{a}\sim 0$, the $\mathrm{NLO^{\prime}}$ correction tends to
rapidly increase when $m_{\pi\pi}\gtrsim 0.48$ GeV. In the small area for
$0.5~{}{\rm GeV}<m_{\pi\pi}<m_{\eta}$, the ratios between the
$\mathrm{NLO^{\prime}}$ correction and the LO result can be 100% and even
larger.

Figure 4: The invariant-mass spectra of the $\pi\pi$ (upper row) and $a\pi$
(lower row) systems in the $\eta\to\pi\pi a$ process. The spectra are shown by
taking three different values of $m_{a}$, as explicitly indicated in the
plots. In each panel, we show three different curves by taking the amplitudes
squared in the LO, NLO, and NLO′ cases. The lighter and darker shaded regions
correspond to the error estimations for the NLO′ results from the Set-I and
Set-II scenarios, respectively. See the text for further explanations.

By integrating over the $t$ or $s$ in the Dalitz distribution via Eq. (50),
one can obtain the $\pi\pi$ or $a\pi$ spectra in the $\eta\to\pi^{+}\pi^{-}a$
decay, and the results from the LO, NLO and $\mathrm{NLO^{\prime}}$ cases are
summarized together in Fig. 4 with different values of $m_{a}$. By integrating
over both $s$ and $t$ variables in Eq. (50), one gets the partial decay width
of $\eta\to\pi^{+}\pi^{-}a$ and the results with varying $m_{a}$ are shown in
Fig. 5. Large uncertainties represented by the shaded areas are observed in
Figs. 4 and 5. In Appendix A, it is demonstrated that by using the lowest
resonance saturation advocated in Ref. [29] the parts of the
$\mathcal{O}(p^{4})$ LECs in the $\eta\to\pi\pi a$ decay are found to be
mainly saturated by the scalar resonances at large $N_{C}$. Phenomenologically
speaking, the couplings and the masses of the scalar resonances in the large
$N_{C}$ case are still under debate [30], and this provides an alternative way
to understand the large uncertainties seen in the $\eta\to\pi\pi a$ amplitudes
at $\mathcal{O}(p^{4})$.

According to the curves in Fig. 5, the partial decay widths of the
$\eta\to\pi^{+}\pi^{-}a$ process calculated from the amplitudes squared from
the LO, NLO and $\rm{NLO^{\prime}}$ cases look roughly similar when varying
$m_{a}$ from 0 to $m_{\eta}-2m_{\pi}$. In fact, the differences resulting by
taking the three types of amplitudes squared are obviously smaller than the
large uncertainties caused by the error bars of the $\mathcal{O}(p^{4})$ LECs.
This also tells us that the large corrections from the NLO amplitude in the
small kinematical regions with $m_{\pi\pi}>0.5$ GeV in Fig. 3 do not sizably
affect the $\eta\to\pi^{+}\pi^{-}a$ decay width.

Figure 5: The partial decay widths (left axis) and branching ratios (right
axis) at different values of $m_{a}$ for the $\eta\to\pi^{+}\pi^{-}a$ process.
The curves labeled as LO, NLO and NLO′ are obtained by taking the perturbative
amplitudes, while the curves labeled as Uni-LO and Uni correspond to the
results by taking the unitarized amplitudes in Eqs. (58) and (39),
respectively. See the text for the details of different notations.

###  5.2 Results by including the $\pi\pi$ final-state interaction

In this part, we analyze the results after including the $\pi\pi$ final-state
interactions via the unitarization procedure explained in Sec. 4. The
amplitude squared $|\mathcal{M}^{\rm CPW}|^{2}$ will be evaluated by taking
the combined partial-wave amplitude $\mathcal{M}_{\eta;\pi\pi a}^{I=0,{\rm
CPW}}$ given in Eq. (49). In this case, the energy dependence of the $a\pi$
amplitude is totally encoded in the $D$-wave part, i.e. the second term of Eq.
(49).

By construction, the phases of the unitarized $S$-wave decay amplitude in Eq.
(39) equal to the $\pi\pi$ scattering phase shifts given by Eq. (48). The
comparison between the $\pi\pi$ phases from the unitarized amplitudes and the
results from the rigorous Roy equation analyses [31] is given in Fig. 6. The
remarkable agreements of the two results justify the use of the neat formalism
in Eq. (39) to unitarize the $\eta\to\pi\pi a$ decay amplitude.

Figure 6:  Comparison between the $\pi\pi$ phases in the unitarized
$\eta\to\pi\pi a$ amplitude (39) and the $\pi\pi$ scattering phase shifts from
Roy equation analyses [31].

After the unitarization procedure, the $f_{0}(500)$ or $\sigma$ resonance is
generated in the $\eta\to\pi\pi a$ decay amplitude and the $\pi\pi\to\pi\pi$
scattering amplitude. The resonance pole position can be calculated by
extrapolating the amplitudes into the second Riemann sheet (RS). The
expressions in Eqs. (39) and (48) denote the amplitudes on the first/physical
RS. The corresponding amplitudes on the second RS can be obtained by replacing
the $G_{\pi\pi}(s)$ function with its expression on the second RS
$G_{\pi\pi}^{\rm II}(s)$, whose explicit formula is

$G_{\pi\pi}^{\rm II}(s)=G_{\pi\pi}(s)-2i\rho_{\pi\pi}(s)\,.$ (56)

To replace the $G_{\pi\pi}(s)$ with its counter part $G_{\pi\pi}^{\rm II}(s)$
in the unitarized amplitudes in Eqs. (39) and (48), one can calculate the pole
position and its residues in the complex energy plane. The $\sigma$ pole
position is found to be $457-i251$ MeV, which is close to the rigorous
dispersive analysis [32]. Around the $\sigma$ resonance pole, the unitarized
$\eta\to\pi\pi a$ amplitude on the second RS can be written as

$\mathcal{M}_{\eta;\pi\pi a}^{00,{\rm Uni,II}}(s)\big{|}_{s\to
s_{\sigma}}\sim-\frac{g_{\sigma\pi\pi}g_{\sigma a\eta}}{s-s_{\sigma}}\,,$ (57)

where the product of the couplings $g_{\sigma\pi\pi}g_{\sigma a\eta}$
corresponds to the residue. Furthermore, the $g_{\sigma\pi\pi}$ coupling can
be determined from the unitarized $\pi\pi\to\pi\pi$ scattering amplitude on
the second RS. Therefore, we can predict the axion-$\eta$-$\sigma$ coupling
$|g_{\sigma a\eta}|$ and the results are shown in Fig. 7. It is clear that the
value of $|g_{\sigma a\eta}|$ increases when increasing the ALP masses
$m_{a}$.

Figure 7:  The axion-$\eta$-$\sigma$ coupling $|g_{\sigma a\pi}|$ multiplied
by $f_{a}$ as a function of the ALP mass $m_{a}$.

Apart from the unitarized formula in Eq. (39), we also explore another
scenario by only including the LO $\eta\to\pi\pi a$ amplitude in the numerator
of the former equation, i.e.,

$\mathcal{M}_{\eta;\pi\pi a}^{00,{\rm Uni-
LO}}(s)=\frac{\mathcal{M}_{\eta;\pi\pi
a}^{00,(2)}(s)}{1-G_{\pi\pi}(s)T_{\pi\pi\to\pi\pi}^{00,(2)}(s)}\,,$ (58)

which resembles the setup employed in Ref. [20] that also uses the LO
$\eta\to\pi\pi a$ amplitude in the construction of final results. The
difference is that the Omnès function is used to take care of the final-state
$\pi\pi$ interaction in the former reference. The good reproduction of the
$\pi\pi$ scattering phase shifts by the unitarization factor
$1/[1-G_{\pi\pi}(s)T_{\pi\pi\to\pi\pi}^{00,(2)}(s)]$, as shown in Fig. 6,
indicates that the $\pi\pi$ final-state interactions are properly included in
our unitarized decay amplitudes. Furthermore, when including the higher order
effects in the $\eta\to\pi\pi a$ amplitudes, such as those in Eq. (44), one
can order by order match the unitarized decay amplitudes to the perturbative
results, which is a merit of our unitarized method.

The partial decay widths calculated by taking the two types of unitarized
amplitudes in Eqs. (39) and (58), together with the perturbative results, are
shown in Fig. 5. It is interesting to point out that the conclusion of Ref.
[20] is confirmed in our study: when only taking the LO $\eta\to\pi\pi a$
amplitude as input to proceed the unitarization, the unitarized decay
amplitude will lead to a clear enhancement around two times compared to the LO
result. However, according to the curves in Fig. 5, our study shows that when
including the NLO $\eta\to\pi\pi a$ amplitude in the unitarization procedure
the resulting partial decay widths are not necessarily enhanced, compared to
the LO result. In fact, when taking the central value of the LECs in Ref.
[21], the partial decay widths from the unitarized amplitudes are even
decreased. Nevertheless, we have explicitly verified that the uncertainties
from the unitarized amplitudes that are caused by the LECs are also quite
large, and the error bands are qualitatively similar to the shaded areas shown
in Fig. 5 from the perturbative analyses.

##  6 Summary and Conclusions

In this work, we have performed the complete calculation of the $\eta\to\pi\pi
a$ amplitude up to one-loop level within the $SU(3)$ axion chiral perturbation
theory in the isospin limit. A bare axion mass term is kept in the
calculation, which allows us to explore the phenomenological consequences by
varying the ALP mass from zero to $m_{\eta}-2m_{\pi}$. As a byproduct, we have
also worked out the ALP-$\eta$ mixing and the ALP mass up to the next-to-
leading order. It is verified that the $\eta\to\pi\pi a$ amplitude can be
renormalized with the conventional renormalization conditions of the $SU(3)$
low energy energy constants. The unitarized $\eta\to\pi\pi a$ decay amplitude
is further constructed to implement the $\pi\pi$ final-state interaction.

We conduct the phenomenological discussions of the $\eta\to\pi\pi a$ decay, by
taking the recent determinations of the $SU(3)$ $\mathcal{O}(p^{4})$ low
energy constants from Ref. [21]. We give the separate analyses by taking the
perturbative amplitudes and the unitarized ones with $\pi\pi$ final-state
interaction. The Dalitz analyses from the perturbative amplitudes reveal that
the distribution along the $a\pi$ energy is rather smooth, while the
distribution along the $\pi\pi$ energy is uneven. For small values of $m_{a}$,
the next-to-leading order correction with respect to the leading-order result
can be larger than 60% in the small region with $m_{\pi\pi}\gtrsim 0.5$ GeV.
The perturbative analyses of the Dalitz distributions indicate that the
$\pi\pi$ final-state interactions could play some relevant roles in the
$\eta\to\pi\pi a$ decay, while the $a\pi$ interaction could be treated
perturbatively at least in the kinematical region of the $\eta\to\pi\pi a$
process. The unitarized $\eta\to\pi\pi a$ decay amplitudes are constructed
under these two requirements, i.e., by nonperturbatively including the
$\pi\pi$ strong final-state interaction and perturbatively including the weak
$a\pi$ interaction. The $\pi\pi$ phases in the unitarized $\eta\to\pi\pi a$
decay amplitude are consistent with the $\pi\pi$ scattering phase shifts from
the rigorous Roy equation study. The $\eta\to\pi\pi a$ decay widths are
predicted as a function of $m_{a}$ in different scenarios by including the
leading-order, the next-to-leading order perturbative amplitudes, and also the
two different unitarized amplitudes. We find that the unitarized decay
amplitude by only taking the leading-order chiral amplitude as input indeed
can obviously enhance the $\eta\to\pi\pi a$ decay width, as observed in Ref.
[20]. However, the more sophisticated unitarized amplitude by including the
next-to-leading order perturbative expressions does not always enhance the
decay width, depending on the input values of the low energy constants. The
uncertainties in the $\eta\to\pi\pi a$ decay caused by nowadays determinations
of the $SU(3)$ low energy constants turn out to be quite large. To further pin
down the error bars of such low energy constants will be very helpful to give
more precise predictions to the $\eta\to\pi\pi a$ decay process.

## Acknowledgements

We would like to thank Sergi Gonzàlez-Solís for communications. This work is
funded in part by the National Natural Science Foundation of China (NSFC)
under Grants Nos. 12150013, 11975090, and the Science Foundation of Hebei
Normal University with contract No. L2023B09.

##  Appendix A Full $\mathcal{O}(p^{4})$ amplitude for the $\eta\to\pi\pi a$
decay

The calculation of the amplitudes is carried out with the help from the
computational tools of FeynCalc [33, 34, 35] and Package-X [36]. The explicit
expressions of the relevant two-point 1PI amplitudes to this work are

$\displaystyle\Sigma_{\pi\pi}^{(4)}(p^{2})=$
$\displaystyle-\frac{12L_{4}p^{2}\left(m_{\eta}^{2}+m_{\pi}^{2}\right)}{F_{\pi}^{2}}-\frac{8L_{5}m_{\pi}^{2}p^{2}}{F_{\pi}^{2}}+\frac{24L_{6}m_{\pi}^{2}\left(m_{\eta}^{2}+m_{\pi}^{2}\right)}{F_{\pi}^{2}}+\frac{16L_{8}m_{\pi}^{4}}{F_{\pi}^{2}}$
(59)
$\displaystyle+\frac{m_{\pi}^{2}A_{0}(m_{\eta})}{96\pi^{2}F_{\pi}^{2}}+\frac{\left(3m_{\eta}^{2}-4m_{K}^{2}+5m_{\pi}^{2}-4p^{2}\right)A_{0}(m_{K})}{192\pi^{2}F_{\pi}^{2}}+\frac{\left(m_{\pi}^{2}-4p^{2}\right)A_{0}(m_{\pi})}{96\pi^{2}F_{\pi}^{2}}\,,$
$\displaystyle\Sigma_{\eta\eta}^{(4)}(p^{2})=$
$\displaystyle-\frac{12L_{4}p^{2}\left(m_{\eta}^{2}+m_{\pi}^{2}\right)}{F_{\pi}^{2}}-\frac{8L_{5}m_{\eta}^{2}p^{2}}{F_{\pi}^{2}}+\frac{24L_{6}m_{\eta}^{2}\left(m_{\eta}^{2}+m_{\pi}^{2}\right)}{F_{\pi}^{2}}$
(60)
$\displaystyle+\frac{24L_{7}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)^{2}}{F_{\pi}^{2}}+\frac{8L_{8}\left(3m_{\eta}^{4}-2m_{\eta}^{2}m_{\pi}^{2}+m_{\pi}^{4}\right)}{F_{\pi}^{2}}+\frac{m_{\pi}^{2}A_{0}(m_{\pi})}{32\pi^{2}F_{\pi}^{2}}$
$\displaystyle-\frac{\left(-9m_{\eta}^{2}+12m_{K}^{2}+m_{\pi}^{2}+12p^{2}\right)A_{0}(m_{K})}{192\pi^{2}F_{\pi}^{2}}+\frac{\left(4m_{\eta}^{2}-m_{\pi}^{2}\right)A_{0}(m_{\eta})}{96\pi^{2}F_{\pi}^{2}}\,,$

$\displaystyle\Sigma_{a\eta}^{(4)}(p^{2})=$
$\displaystyle-\frac{2\sqrt{3}L_{4}p^{2}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)\left(m_{\eta}^{2}+m_{\pi}^{2}\right)}{F_{\pi}f_{a}m_{\eta}^{2}}-\frac{4L_{5}p^{2}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)}{\sqrt{3}F_{\pi}f_{a}}$
(61)
$\displaystyle-\frac{4m_{\pi}^{2}(3L_{7}+L_{8})\left(3m_{\eta}^{4}-4m_{\eta}^{2}m_{\pi}^{2}+m_{\pi}^{4}\right)}{\sqrt{3}F_{\pi}f_{a}m_{\eta}^{2}}-\frac{m_{\pi}^{2}\left(m_{\pi}^{2}-3m_{\eta}^{2}\right)A_{0}(m_{\pi})}{64\sqrt{3}\pi^{2}F_{\pi}f_{a}m_{\eta}^{2}}$
$\displaystyle+\frac{\left[-6p^{2}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)-3m_{\eta}^{2}m_{\pi}^{2}+m_{\pi}^{4}\right]A_{0}(m_{K})}{96\sqrt{3}\pi^{2}F_{\pi}f_{a}m_{\eta}^{2}}+\frac{m_{\pi}^{2}\left(m_{\pi}^{2}-3m_{\eta}^{2}\right)A_{0}(m_{\eta})}{192\sqrt{3}\pi^{2}F_{\pi}f_{a}m_{\eta}^{2}}$
$\displaystyle+\frac{m_{a}^{2}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)}{\sqrt{3}m_{\eta}^{2}\left(m_{a}^{2}-m_{\eta}^{2}\right)}\Bigg{\\{}\frac{6L_{4}p^{2}\left(m_{\eta}^{2}+m_{\pi}^{2}\right)}{F_{\pi}f_{a}}+\frac{4L_{5}m_{\eta}^{2}p^{2}}{F_{\pi}f_{a}}-\frac{12L_{6}m_{\eta}^{2}\left(m_{\eta}^{2}+m_{\pi}^{2}\right)}{F_{\pi}f_{a}}$
$\displaystyle-\frac{12L_{7}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)^{2}}{F_{\pi}f_{a}}-\frac{4L_{8}\left(3m_{\eta}^{4}-2m_{\eta}^{2}m_{\pi}^{2}+m_{\pi}^{4}\right)}{F_{\pi}f_{a}}-\frac{m_{\pi}^{2}A_{0}(m_{\pi})}{64\pi^{2}F_{\pi}f_{a}}$
$\displaystyle+\frac{\left(-9m_{\eta}^{2}+12m_{K}^{2}+m_{\pi}^{2}+12p^{2}\right)A_{0}(m_{K})}{384\pi^{2}F_{\pi}f_{a}}+\frac{\left(m_{\pi}^{2}-4m_{\eta}^{2}\right)A_{0}(m_{\eta})}{192\pi^{2}F_{\pi}f_{a}}\Bigg{\\}}\,,$
$\displaystyle\Sigma_{aa}^{(4)}(p^{2})=$
$\displaystyle\frac{2L_{6}m_{\pi}^{2}\left(3m_{\eta}^{2}-m_{\pi}^{2}\right)\left(m_{\eta}^{2}+m_{\pi}^{2}\right)}{f_{a}^{2}m_{\eta}^{2}}+\frac{2(3L_{7}+L_{8})m_{\pi}^{4}\left(3m_{\eta}^{2}-m_{\pi}^{2}\right)^{2}}{3f_{a}^{2}m_{\eta}^{4}}$
(62)
$\displaystyle+\frac{m_{\pi}^{2}\left(3m_{\eta}^{2}-m_{\pi}^{2}\right)^{2}A_{0}(m_{\pi})}{384\pi^{2}f_{a}^{2}m_{\eta}^{4}}+\frac{m_{\pi}^{2}\left(3m_{\eta}^{2}-m_{\pi}^{2}\right)\left(3m_{\eta}^{2}+m_{\pi}^{2}\right)A_{0}(m_{K})}{576\pi^{2}f_{a}^{2}m_{\eta}^{4}}$
$\displaystyle+\frac{m_{\pi}^{2}\left(3m_{\eta}^{2}-m_{\pi}^{2}\right)\left(m_{\eta}^{2}+m_{\pi}^{2}\right)A_{0}(m_{\eta})}{1152\pi^{2}f_{a}^{2}m_{\eta}^{4}}$
$\displaystyle+\frac{m_{a}^{2}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)}{2\sqrt{3}m_{\eta}^{2}\left(m_{a}^{2}-m_{\eta}^{2}\right)}\Bigg{\\{}\frac{4\sqrt{3}L_{4}p^{2}\left(m_{\eta}^{4}-m_{\pi}^{4}\right)}{f_{a}^{2}m_{\eta}^{2}}+\frac{8L_{5}p^{2}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)}{\sqrt{3}f_{a}^{2}}$
$\displaystyle+\frac{8(3L_{7}+L_{8})m_{\pi}^{2}\left(3m_{\eta}^{4}-4m_{\eta}^{2}m_{\pi}^{2}+m_{\pi}^{4}\right)}{\sqrt{3}f_{a}^{2}m_{\eta}^{2}}+\frac{\left(m_{\pi}^{4}-3m_{\eta}^{2}m_{\pi}^{2}\right)A_{0}(m_{\pi})}{32\sqrt{3}\pi^{2}f_{a}^{2}m_{\eta}^{2}}$
$\displaystyle+\frac{\left[3m_{\eta}^{2}\left(m_{\pi}^{2}+2p^{2}\right)-m_{\pi}^{2}\left(m_{\pi}^{2}+6p^{2}\right)\right]A_{0}(m_{K})}{48\sqrt{3}\pi^{2}f_{a}^{2}m_{\eta}^{2}}+\frac{m_{\pi}^{2}\left(3m_{\eta}^{2}-m_{\pi}^{2}\right)A_{0}(m_{\eta})}{96\sqrt{3}\pi^{2}f_{a}^{2}m_{\eta}^{2}}\Bigg{\\}}$
$\displaystyle+\frac{m_{a}^{4}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)^{2}}{12m_{\eta}^{4}\left(m_{a}^{2}-m_{\eta}^{2}\right)^{2}}\Bigg{\\{}\frac{24L_{6}m_{\eta}^{2}\left(m_{\eta}^{2}+m_{\pi}^{2}\right)}{f_{a}^{2}}-\frac{12L_{4}p^{2}\left(m_{\eta}^{2}+m_{\pi}^{2}\right)}{f_{a}^{2}}-\frac{8L_{5}m_{\eta}^{2}p^{2}}{f_{a}^{2}}$
$\displaystyle+\frac{24L_{7}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)^{2}}{f_{a}^{2}}+\frac{8L_{8}\left(3m_{\eta}^{4}-2m_{\eta}^{2}m_{\pi}^{2}+m_{\pi}^{4}\right)}{f_{a}^{2}}+\frac{m_{\pi}^{2}A_{0}(m_{\pi})}{32\pi^{2}f_{a}^{2}}$
$\displaystyle-\frac{\left(m_{\pi}^{2}+3p^{2}\right)A_{0}(m_{K})}{48\pi^{2}f_{a}^{2}}-\frac{\left(m_{\pi}^{2}-4m_{\eta}^{2}\right)A_{0}(m_{\eta})}{96\pi^{2}f_{a}^{2}}\Bigg{\\}}\,,$

where

$A_{0}(m)=-(\mu^{2})^{\frac{d-4}{2}}m^{2}\left(R+\log\frac{m^{2}}{\mu^{2}}\right)\,,\quad
R=\frac{2}{d-4}-\Gamma^{\prime}(1)-\log(4\pi)-1\,,$ (63)

with $d\to 4$. The UV divergence in the loop corrections is renormalized in
the physical quantities with the bare LECs [23]

$L_{i}=L_{i}^{r}(\mu)+\Gamma_{i}\frac{(\mu^{2})^{\frac{d-4}{2}}}{2(4\pi)^{2}}R\,,$
(64)

with $\Gamma_{1}=3/32$, $\Gamma_{2}=3/16$, $\Gamma_{3}=0$, $\Gamma_{4}=1/8$,
$\Gamma_{5}=3/8$, $\Gamma_{6}=11/144$, $\Gamma_{7}=0$, $\Gamma_{8}=5/48$.

The full perturbative $\chi$PT amplitude of the $\eta\to\pi\pi a$ process is
given by

$\mathcal{M}_{\eta;\pi\pi a}^{(2+4)}(s,t,u)=\mathcal{M}_{\eta;\pi\pi
a}^{(2)}(s,t,u)+\mathcal{M}_{\eta;\pi\pi a}^{(4)}(s,t,u),$ (65)

where the LO part $\mathcal{M}_{\eta;\pi\pi a}^{(2)}(s,t,u)$ is given in Eq.
(31) and the NLO part will be decomposed into four pieces: the $s$-channel
loop diagram (Fig. 1(d)) contribution $\mathcal{M}_{\eta;\pi\pi a}^{s}$, the
$t$\- and $u$-channel loop diagrams (Figs. 1(e) and 1(f)) contributions
$\mathcal{M}_{\eta;\pi\pi a}^{t,u}$, and the $\mathcal{M}_{\eta;\pi\pi
a}^{\mathrm{c+tad}}$ term that includes the contributions from the
$\mathcal{O}(p^{4})$ contact diagram (Fig. 1(b)), tadpole diagram (Fig. 1(c)),
the wave function renormalizations of $\pi$ and $\eta$, and the NLO ALP-$\eta$
mixing, as well as the correction of Eq. (32). the The explicit expression for
the NLO part is

$\mathcal{M}_{\eta;\pi\pi a}^{(4)}=\mathcal{M}_{\eta;\pi\pi
a}^{s}+\mathcal{M}_{\eta;\pi\pi a}^{t}+\mathcal{M}_{\eta;\pi\pi
a}^{u}+\mathcal{M}_{\eta;\pi\pi a}^{\mathrm{c}+\mathrm{tad}}\,,$ (66)

where

$\displaystyle\mathcal{M}_{\eta;\pi\pi
a}^{s}=\frac{m_{\eta}^{2}\left[9s\left(3s-m_{a}^{2}\right)-48m_{\pi}^{4}+81m_{\pi}^{2}s\right]}{2304\sqrt{3}\pi^{2}f_{a}F_{\pi}^{3}m_{\eta}^{2}}-\frac{27m_{\eta}^{4}s-16m_{\pi}^{6}+18m_{\pi}^{4}s}{2304\sqrt{3}\pi^{2}f_{a}F_{\pi}^{3}m_{\eta}^{2}}$
$\displaystyle+\frac{9m_{\pi}^{2}s\left(m_{a}^{2}-3s\right)}{2304\sqrt{3}\pi^{2}f_{a}F_{\pi}^{3}m_{\eta}^{2}}+\frac{\left(3m_{\eta}^{2}-m_{\pi}^{2}\right)\left(7m_{\pi}^{2}-6s\right)\mu_{\pi}}{18\sqrt{3}f_{a}F_{\pi}m_{\eta}^{2}}-\frac{m_{\pi}^{4}\left(m_{\pi}^{2}-3m_{\eta}^{2}\right)\mu_{\eta}}{18\sqrt{3}f_{a}F_{\pi}m_{\eta}^{4}}$
$\displaystyle+\frac{\left(2m_{K}^{2}-3s\right)\left[3m_{\eta}^{2}\left(m_{\pi}^{2}+3s\right)-3m_{a}^{2}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)-9m_{\eta}^{4}+2m_{\pi}^{4}-9m_{\pi}^{2}s\right]\mu_{K}}{72\sqrt{3}f_{a}F_{\pi}m_{\eta}^{2}m_{K}^{2}}$
$\displaystyle+\frac{m_{\pi}^{2}\left(m_{\pi}^{2}-3m_{\eta}^{2}\right)\left(m_{\pi}^{2}-2s\right)\mathcal{F}_{B_{0}}(s,m_{\pi},m_{\pi})}{192\sqrt{3}\pi^{2}f_{a}F_{\pi}^{3}m_{\eta}^{2}}+\frac{m_{\pi}^{4}\left(m_{\pi}^{2}-3m_{\eta}^{2}\right)\mathcal{F}_{B_{0}}(s,m_{\eta},m_{\eta})}{576\sqrt{3}\pi^{2}f_{a}F_{\pi}^{3}m_{\eta}^{2}}$
$\displaystyle-\frac{s\left[3m_{a}^{2}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)+9m_{\eta}^{4}-3m_{\eta}^{2}\left(m_{\pi}^{2}+3s\right)-2m_{\pi}^{4}+9m_{\pi}^{2}s\right]\mathcal{F}_{B_{0}}(s,m_{K},m_{K})}{768\sqrt{3}\pi^{2}f_{a}F_{\pi}^{3}m_{\eta}^{2}}$
$\displaystyle+\frac{m_{a}^{2}}{m_{a}^{2}-m_{\eta}^{2}}\Bigg{\\{}\frac{\left(m_{\eta}^{2}-m_{\pi}^{2}\right)}{4608\sqrt{3}\pi^{2}f_{a}F_{\pi}^{3}m_{\eta}^{2}}\Big{[}18s\left(m_{a}^{2}+2m_{K}^{2}-3s\right)-m_{\eta}^{2}\left(32m_{\pi}^{2}+9s\right)+32m_{\pi}^{4}$
$\displaystyle-45m_{\pi}^{2}s\Big{]}-\frac{\left(m_{\eta}^{2}-m_{\pi}^{2}\right)\left(7m_{\pi}^{2}-6s\right)\mu_{\pi}}{18\sqrt{3}f_{a}F_{\pi}m_{\eta}^{2}}+\frac{\left(4m_{\eta}^{4}m_{\pi}^{2}-5m_{\eta}^{2}m_{\pi}^{4}+m_{\pi}^{6}\right)\mu_{\eta}}{18\sqrt{3}f_{a}F_{\pi}m_{\eta}^{4}}$
$\displaystyle+\frac{\left(m_{\eta}^{2}-m_{\pi}^{2}\right)\mu_{K}}{144\sqrt{3}f_{a}F_{\pi}m_{\eta}^{2}m_{K}^{2}}\Big{[}3s\left(3m_{\eta}^{2}-m_{\pi}^{2}+18s\right)+4m_{K}^{2}\left(3m_{\eta}^{2}+2m_{\pi}^{2}-24s\right)$
$\displaystyle+6m_{a}^{2}\left(2m_{K}^{2}-3s\right)\Big{]}+\frac{m_{\pi}^{2}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)\left(m_{\pi}^{2}-2s\right)\mathcal{F}_{B_{0}}(s,m_{\pi},m_{\pi})}{192\sqrt{3}\pi^{2}f_{a}F_{\pi}^{3}m_{\eta}^{2}}$
$\displaystyle-\frac{m_{\pi}^{2}\left(4m_{\eta}^{4}-5m_{\eta}^{2}m_{\pi}^{2}+m_{\pi}^{4}\right)\mathcal{F}_{B_{0}}(s,m_{\eta},m_{\eta})}{576\sqrt{3}\pi^{2}f_{a}F_{\pi}^{3}m_{\eta}^{2}}-\frac{s\left(m_{\eta}^{2}-m_{\pi}^{2}\right)\mathcal{F}_{B_{0}}(s,m_{K},m_{K})}{1536\sqrt{3}\pi^{2}f_{a}F_{\pi}^{3}m_{\eta}^{2}}$
$\displaystyle\times\big{(}18s-6m_{a}^{2}+3m_{\eta}^{2}-12m_{K}^{2}-m_{\pi}^{2}\big{)}\Bigg{\\}}\,,$
(67)

$\displaystyle\mathcal{M}_{\eta;\pi\pi
a}^{t}=\frac{3m_{a}^{2}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)\left(3m_{\eta}^{2}+m_{\pi}^{2}-3t\right)+3m_{\eta}^{2}\left(2m_{\pi}^{4}+15m_{\pi}^{2}t+9t^{2}\right)}{2304\sqrt{3}\pi^{2}f_{a}F_{\pi}^{3}m_{\eta}^{2}}$
$\displaystyle-\frac{27m_{\eta}^{4}\left(m_{\pi}^{2}+t\right)}{2304\sqrt{3}\pi^{2}f_{a}F_{\pi}^{3}m_{\eta}^{2}}-\frac{11m_{\pi}^{6}-18m_{\pi}^{4}t+27m_{\pi}^{2}t^{2}}{2304\sqrt{3}\pi^{2}f_{a}F_{\pi}^{3}m_{\eta}^{2}}+\frac{m_{\pi}^{2}\left(m_{\pi}^{2}-3m_{\eta}^{2}\right)\mu_{\pi}}{18\sqrt{3}f_{a}F_{\pi}m_{\eta}^{2}}$
$\displaystyle+\frac{m_{\pi}^{4}\left(m_{\pi}^{2}-3m_{\eta}^{2}\right)\mu_{\eta}}{18\sqrt{3}f_{a}F_{\pi}m_{\eta}^{4}}+\frac{\left(3m_{\eta}^{2}+2m_{K}^{2}+m_{\pi}^{2}-3t\right)\mu_{K}}{24\sqrt{3}f_{a}F_{\pi}m_{\eta}^{2}m_{K}^{2}}\Big{[}m_{a}^{2}\left(m_{\pi}^{2}-m_{\eta}^{2}\right)+m_{\pi}^{4}$
$\displaystyle+3m_{\eta}^{2}\left(m_{\pi}^{2}+t\right)-3m_{\pi}^{2}t\Big{]}-\frac{m_{\pi}^{4}\left(3m_{\eta}^{4}-4m_{\eta}^{2}m_{\pi}^{2}+m_{\pi}^{4}\right)}{576\sqrt{3}\pi^{2}f_{a}F_{\pi}^{3}m_{\eta}^{2}t}\log\frac{m_{\eta}^{2}}{m_{\pi}^{2}}-\frac{m_{\pi}^{4}\left(m_{\pi}^{2}-3m_{\eta}^{2}\right)}{288\sqrt{3}\pi^{2}f_{a}F_{\pi}^{3}m_{\eta}^{2}}$
$\displaystyle\times\mathcal{F}_{B_{0}}(t,m_{\eta},m_{\pi})-\frac{\left(3m_{\eta}^{2}+m_{\pi}^{2}-3t\right)\mathcal{F}_{B_{0}}(t,m_{K},m_{K})}{768\sqrt{3}\pi^{2}f_{a}F_{\pi}^{3}m_{\eta}^{2}}\Big{[}m_{a}^{2}\left(m_{\pi}^{2}-m_{\eta}^{2}\right)+m_{\pi}^{4}$
$\displaystyle+3m_{\eta}^{2}\left(m_{\pi}^{2}+t\right)-3m_{\pi}^{2}t\Big{]}+\frac{m_{a}^{2}}{m_{a}^{2}-m_{\eta}^{2}}\Bigg{\\{}\frac{\left(m_{\pi}^{2}-m_{\eta}^{2}\right)}{2304\sqrt{3}\pi^{2}f_{a}F_{\pi}^{3}m_{\eta}^{2}}\Big{[}3m_{a}^{2}\left(3m_{\eta}^{2}+m_{\pi}^{2}-3t\right)$
$\displaystyle+18m_{\eta}^{4}+15m_{\eta}^{2}\left(m_{\pi}^{2}-3t\right)+11m_{\pi}^{4}-18m_{\pi}^{2}t+27t^{2}\Big{]}+\frac{m_{\pi}^{2}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)\mu_{\pi}}{18\sqrt{3}f_{a}F_{\pi}m_{\eta}^{2}}$
$\displaystyle+\frac{m_{\pi}^{4}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)\mu_{\eta}}{18\sqrt{3}f_{a}F_{\pi}m_{\eta}^{4}}+\frac{\left(m_{\eta}^{2}-m_{\pi}^{2}\right)\mu_{K}}{24\sqrt{3}f_{a}F_{\pi}m_{\eta}^{2}m_{K}^{2}}\Big{[}m_{a}^{2}\left(3m_{\eta}^{2}+2m_{K}^{2}+m_{\pi}^{2}-3t\right)+6m_{\eta}^{4}+5m_{\eta}^{2}m_{\pi}^{2}$
$\displaystyle-15m_{\eta}^{2}t+2m_{K}^{2}\left(5m_{\eta}^{2}+2m_{\pi}^{2}-5t\right)+m_{\pi}^{4}-6m_{\pi}^{2}t+9t^{2}\Big{]}+\frac{m_{\pi}^{4}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)^{2}}{576\sqrt{3}\pi^{2}f_{a}F_{\pi}^{3}m_{\eta}^{2}t}\log\frac{m_{\eta}^{2}}{m_{\pi}^{2}}$
$\displaystyle+\frac{m_{\pi}^{4}\left(m_{\pi}^{2}-m_{\eta}^{2}\right)\mathcal{F}_{B_{0}}(t,m_{\eta},m_{\pi})}{288\sqrt{3}\pi^{2}f_{a}F_{\pi}^{3}m_{\eta}^{2}}-\frac{\left(m_{\eta}^{2}-m_{\pi}^{2}\right)\left(3m_{\eta}^{2}+m_{\pi}^{2}-3t\right)\mathcal{F}_{B_{0}}(t,m_{K},m_{K})}{768\sqrt{3}\pi^{2}f_{a}F_{\pi}^{3}m_{\eta}^{2}}$
$\displaystyle\times\left(m_{a}^{2}+2m_{\eta}^{2}+m_{\pi}^{2}-3t\right)\Bigg{\\}}\,,$
(68)

$\mathcal{M}_{\eta;\pi\pi a}^{u}=\mathcal{M}_{\eta;\pi\pi a}^{t}\big{|}_{t\to
u}\,,$ (69)

$\displaystyle\mathcal{M}_{\eta;\pi\pi
a}^{\mathrm{c}+\mathrm{tad}}=-\frac{4L_{1}^{r}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)\left(s-2m_{\pi}^{2}\right)\left(m_{a}^{2}+m_{\eta}^{2}-s\right)}{\sqrt{3}f_{a}F_{\pi}^{3}m_{\eta}^{2}}-\frac{2L_{2}^{r}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)}{\sqrt{3}f_{a}F_{\pi}^{3}m_{\eta}^{2}}\Big{[}m_{a}^{4}+m_{\eta}^{4}$
$\displaystyle-s\left(m_{a}^{2}+m_{\eta}^{2}\right)-2m_{\pi}^{4}+2m_{\pi}^{2}(t+u)-t^{2}-u^{2}\Big{]}+\frac{2L_{3}^{r}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)}{3\sqrt{3}f_{a}F_{\pi}^{3}m_{\eta}^{2}}\Big{(}s^{2}+t^{2}+u^{2}-m_{a}^{4}$
$\displaystyle-
m_{\eta}^{4}-2m_{\pi}^{4}\Big{)}+\frac{4L_{4}^{r}m_{\pi}^{2}\left(m_{\pi}^{2}-m_{\eta}^{2}\right)\left(m_{a}^{2}+m_{\eta}^{2}-s\right)}{\sqrt{3}f_{a}F_{\pi}^{3}m_{\eta}^{2}}+\frac{2L_{5}^{r}m_{\pi}^{2}}{3\sqrt{3}f_{a}F_{\pi}^{3}m_{\eta}^{2}}\Big{[}m_{a}^{2}\left(m_{\eta}^{2}+m_{\pi}^{2}\right)$
$\displaystyle-6m_{\eta}^{4}+m_{\eta}^{2}m_{\pi}^{2}-m_{\pi}^{4}\Big{]}-\frac{4L_{7}^{r}\left(6m_{\eta}^{6}m_{\pi}^{2}-13m_{\eta}^{4}m_{\pi}^{4}+2m_{\eta}^{2}m_{\pi}^{6}+m_{\pi}^{8}\right)}{\sqrt{3}f_{a}F_{\pi}^{3}m_{\eta}^{4}}-\frac{4L_{8}^{r}m_{\pi}^{4}}{3\sqrt{3}f_{a}F_{\pi}^{3}m_{\eta}^{4}}$
$\displaystyle\times\big{(}2m_{\eta}^{2}m_{\pi}^{2}+m_{\pi}^{4}-15m_{\eta}^{4}\big{)}+\frac{\left(13m_{\eta}^{2}m_{\pi}^{4}+3m_{\pi}^{6}-30m_{\eta}^{4}m_{\pi}^{2}\right)\mu_{\pi}}{18\sqrt{3}f_{a}F_{\pi}m_{\eta}^{4}}-\frac{\left(m_{\pi}^{6}-3m_{\eta}^{2}m_{\pi}^{4}\right)\mu_{\eta}}{18\sqrt{3}f_{a}F_{\pi}m_{\eta}^{4}}$
$\displaystyle+\frac{\left[18m_{a}^{2}\left(m_{\eta}^{4}-m_{\eta}^{2}m_{\pi}^{2}\right)+5\left(3m_{\eta}^{4}m_{\pi}^{2}-7m_{\eta}^{2}m_{\pi}^{4}-2m_{\pi}^{6}\right)\right]\mu_{K}}{90\sqrt{3}f_{a}F_{\pi}m_{\eta}^{4}}$
$\displaystyle+\frac{m_{a}^{2}}{m_{a}^{2}-m_{\eta}^{2}}\Bigg{\\{}\frac{4L_{1}^{r}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)\left(s-2m_{\pi}^{2}\right)\left(m_{a}^{2}+m_{\eta}^{2}-s\right)}{\sqrt{3}f_{a}F_{\pi}^{3}m_{\eta}^{2}}+\frac{2L_{2}^{r}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)}{\sqrt{3}f_{a}F_{\pi}^{3}m_{\eta}^{2}}\Big{[}m_{a}^{4}+m_{\eta}^{4}$
$\displaystyle-s\left(m_{a}^{2}+m_{\eta}^{2}\right)-2m_{\pi}^{4}+2m_{\pi}^{2}(t+u)-t^{2}-u^{2}\Big{]}-\frac{2L_{3}^{r}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)}{3\sqrt{3}f_{a}F_{\pi}^{3}m_{\eta}^{2}}\Big{(}s^{2}+t^{2}+u^{2}-m_{a}^{4}$
$\displaystyle-
m_{\eta}^{4}-2m_{\pi}^{4}\Big{)}+\frac{4L_{4}^{r}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)\left[m_{a}^{2}m_{\pi}^{2}+m_{\eta}^{2}\left(3m_{\pi}^{2}-s\right)-m_{\pi}^{2}s\right]}{\sqrt{3}f_{a}F_{\pi}^{3}m_{\eta}^{2}}+\frac{2L_{5}^{r}m_{\pi}^{2}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)}{3\sqrt{3}f_{a}F_{\pi}^{3}m_{\eta}^{2}}$
$\displaystyle\times\left(m_{a}^{2}+2m_{\eta}^{2}-m_{\pi}^{2}\right)+\frac{16L_{6}^{r}\left(m_{\pi}^{4}-m_{\eta}^{2}m_{\pi}^{2}\right)}{\sqrt{3}f_{a}F_{\pi}^{3}}+\frac{4L_{7}^{r}\left(4m_{\eta}^{2}+m_{\pi}^{2}\right)\left(m_{\pi}^{3}-m_{\eta}^{2}m_{\pi}\right)^{2}}{\sqrt{3}f_{a}F_{\pi}^{3}m_{\eta}^{4}}$
$\displaystyle+\frac{4L_{8}^{r}m_{\pi}^{4}\left(2m_{\eta}^{2}m_{\pi}^{2}+m_{\pi}^{4}-3m_{\eta}^{4}\right)}{3\sqrt{3}f_{a}F_{\pi}^{3}m_{\eta}^{4}}+\frac{\left(10m_{\eta}^{4}m_{\pi}^{2}-13m_{\eta}^{2}m_{\pi}^{4}-3m_{\pi}^{6}\right)\mu_{\pi}}{18\sqrt{3}f_{a}F_{\pi}m_{\eta}^{4}}+\frac{\left(m_{\pi}^{6}-3m_{\eta}^{2}m_{\pi}^{4}\right)}{18\sqrt{3}f_{a}F_{\pi}m_{\eta}^{4}}$
$\displaystyle\times\mu_{\eta}+\frac{\left[7m_{\eta}^{4}m_{\pi}^{2}+50m_{\eta}^{2}m_{\pi}^{4}+10m_{\pi}^{6}-3m_{a}^{2}\left(m_{\eta}^{4}-m_{\eta}^{2}m_{\pi}^{2}\right)-27m_{\eta}^{6}\right]\mu_{K}}{90\sqrt{3}f_{a}F_{\pi}m_{\eta}^{4}}\Bigg{\\}}$
$\displaystyle+\frac{2L_{4}^{r}m_{a}^{2}m_{\pi}^{2}\left(m_{\pi}^{4}-m_{\eta}^{4}\right)}{\sqrt{3}f_{a}F_{\pi}^{3}m_{\eta}^{2}\left(m_{a}^{2}-m_{\eta}^{2}\right)}+\frac{4L_{5}^{r}m_{a}^{2}m_{\pi}^{2}\left(m_{\pi}^{2}-m_{\eta}^{2}\right)}{3\sqrt{3}f_{a}F_{\pi}^{3}\left(m_{a}^{2}-m_{\eta}^{2}\right)}+\frac{4m_{\pi}^{4}(3L_{7}^{r}+L_{8}^{r})}{3\sqrt{3}f_{a}F_{\pi}^{3}m_{\eta}^{2}\left(m_{\eta}^{2}-m_{a}^{2}\right)}\Big{(}3m_{\eta}^{4}+m_{\pi}^{4}$
$\displaystyle-4m_{\eta}^{2}m_{\pi}^{2}\Big{)}-\frac{m_{\pi}^{4}\left(m_{\pi}^{2}-3m_{\eta}^{2}\right)\mu_{\pi}}{6\sqrt{3}f_{a}F_{\pi}m_{\eta}^{2}\left(m_{\eta}^{2}-m_{a}^{2}\right)}+\frac{m_{\pi}^{4}\left(m_{\pi}^{2}-3m_{\eta}^{2}\right)\mu_{\eta}}{18\sqrt{3}f_{a}F_{\pi}m_{\eta}^{2}\left(m_{\eta}^{2}-m_{a}^{2}\right)}+\Big{[}m_{\pi}^{4}-3m_{\eta}^{2}m_{\pi}^{2}$
$\displaystyle-6m_{a}^{2}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)\Big{]}\frac{m_{\pi}^{2}\mu_{K}}{9\sqrt{3}f_{a}F_{\pi}m_{\eta}^{2}\left(m_{\eta}^{2}-m_{a}^{2}\right)}+\frac{m_{a}^{2}}{m_{a}^{2}-m_{\eta}^{2}}\Bigg{\\{}\frac{2L_{4}^{r}m_{\pi}^{2}\left(m_{\eta}^{4}-m_{\pi}^{4}\right)}{\sqrt{3}f_{a}F_{\pi}^{3}m_{\eta}^{2}}$
$\displaystyle+\frac{4L_{5}^{r}m_{\pi}^{2}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)}{3\sqrt{3}f_{a}F_{\pi}^{3}}+\frac{\left(m_{\pi}^{4}-m_{\eta}^{2}m_{\pi}^{2}\right)\mu_{K}}{3\sqrt{3}f_{a}F_{\pi}m_{\eta}^{2}}\Bigg{\\}}\,,$
(70)

with

$\displaystyle\mu_{p}=\frac{m_{p}^{2}}{2(4\pi
F_{\pi})^{2}}\log\frac{m_{p}^{2}}{\mu^{2}}\,,$ (71)
$\mathcal{F}_{B_{0}}(s,m_{1},m_{2})=\frac{\sqrt{\lambda(s,m_{1}^{2},m_{2}^{2})}}{s}\log\left(\frac{m_{1}^{2}+m_{2}^{2}-s+\sqrt{\lambda(s,m_{1}^{2},m_{2}^{2})}}{2m_{1}m_{2}}\right)\,.$
(72)

In the region of $s>(m_{1}+m_{2})^{2}$, the function
$\mathcal{F}_{B_{0}}(s,m_{1},m_{2})$ develops the unitary cut, which is given
by

$\mathrm{Im}\mathcal{F}_{B_{0}}(s+i0^{+},m_{1},m_{2})=\theta\left[s-(m_{1}+m_{2})^{2}\right]\pi\frac{\sqrt{\lambda(s,m_{1}^{2},m_{2}^{2})}}{s}\,.$
(73)

By assuming the lowest resonance saturation to the $\mathcal{O}(p^{4})$ LECs
and the large $N_{C}$ relations of the resonance couplings in Ref. [29]

$\displaystyle L_{1}^{r}=\frac{G_{V}^{2}}{8M_{V}^{2}}\,,\quad
L_{2}^{r}=\frac{G_{V}^{2}}{4M_{V}^{2}}\,,\quad
L_{3}^{r}=\frac{c_{d}^{2}}{2M_{S}^{2}}-\frac{3G_{V}^{2}}{4M_{V}^{2}}\,,\quad
L_{4}^{r}=0\,,$ (74) $\displaystyle
L_{5}^{r}=\frac{c_{d}c_{m}}{M_{S}^{2}}\,,\quad L_{6}^{r}=0\,,\quad
L_{7}^{r}=0\,,\quad
L_{8}^{r}=\frac{c_{m}^{2}}{2M_{S}^{2}}-\frac{d_{m}^{2}}{2M_{P}^{2}}\,,$

the $L_{i=1,\cdots,8}^{r}$ terms in Eq. (66) become

$\displaystyle\mathcal{M}_{\eta;\pi\pi a}^{\mathrm{LECs}}=$
$\displaystyle\frac{c_{d}^{2}\left(m_{\eta}^{2}-m_{\pi}^{2}\right)\left(m_{a}^{4}+m_{\eta}^{4}+2m_{\pi}^{4}-s^{2}-t^{2}-u^{2}\right)}{3\sqrt{3}f_{a}F_{\pi}^{3}M_{S}^{2}\left(m_{a}^{2}-m_{\eta}^{2}\right)}+\frac{4c_{m}^{2}m_{\pi}^{4}\left(2m_{a}^{2}-3m_{\eta}^{2}+m_{\pi}^{2}\right)}{\sqrt{3}f_{a}F_{\pi}^{3}M_{S}^{2}\left(m_{a}^{2}-m_{\eta}^{2}\right)}$
(75)
$\displaystyle+\frac{2c_{d}c_{m}m_{\pi}^{2}\left[2m_{a}^{4}-m_{a}^{2}\left(5m_{\eta}^{2}+3m_{\pi}^{2}\right)+6m_{\eta}^{4}-m_{\eta}^{2}m_{\pi}^{2}+m_{\pi}^{4}\right]}{3\sqrt{3}f_{a}F_{\pi}^{3}M_{S}^{2}\left(m_{a}^{2}-m_{\eta}^{2}\right)}$
$\displaystyle-\frac{4d_{m}^{2}m_{\pi}^{4}\left(2m_{a}^{2}-3m_{\eta}^{2}+m_{\pi}^{2}\right)}{\sqrt{3}f_{a}F_{\pi}^{3}M_{P}^{2}\left(m_{a}^{2}-m_{\eta}^{2}\right)}\,,$

which shows that the vector resonances do not contribute to the $\eta\to\pi\pi
a$ process.

## References

  * Kim and Carosi [2010] J. E. Kim and G. Carosi, Rev. Mod. Phys. 82, 557 (2010), [Erratum: Rev.Mod.Phys. 91, 049902 (2019)], arXiv:0807.3125 [hep-ph] .
  * Graham  _et al._ [2015] P. W. Graham, I. G. Irastorza, S. K. Lamoreaux, A. Lindner, and K. A. van Bibber, Ann. Rev. Nucl. Part. Sci. 65, 485 (2015), arXiv:1602.00039 [hep-ex] .
  * Irastorza and Redondo [2018] I. G. Irastorza and J. Redondo, Prog. Part. Nucl. Phys. 102, 89 (2018), arXiv:1801.08127 [hep-ph] .
  * Di Luzio  _et al._ [2020] L. Di Luzio, M. Giannotti, E. Nardi, and L. Visinelli, Phys. Rept. 870, 1 (2020), arXiv:2003.01100 [hep-ph] .
  * Choi  _et al._ [2021] K. Choi, S. H. Im, and C. Sub Shin, Ann. Rev. Nucl. Part. Sci. 71, 225 (2021), arXiv:2012.05029 [hep-ph] .
  * Sikivie [2021] P. Sikivie, Rev. Mod. Phys. 93, 015004 (2021), arXiv:2003.02206 [hep-ph] .
  * Weinberg [1978] S. Weinberg, Phys. Rev. Lett. 40, 223 (1978).
  * Peccei and Quinn [1977a] R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977a).
  * Peccei and Quinn [1977b] R. D. Peccei and H. R. Quinn, Phys. Rev. D 16, 1791 (1977b).
  * Wilczek [1978] F. Wilczek, Phys. Rev. Lett. 40, 279 (1978).
  * Kim [1987] J. E. Kim, Phys. Rept. 150, 1 (1987).
  * Di Luzio  _et al._ [2017a] L. Di Luzio, F. Mescia, and E. Nardi, Phys. Rev. Lett. 118, 031801 (2017a), arXiv:1610.07593 [hep-ph] .
  * Di Luzio  _et al._ [2017b] L. Di Luzio, F. Mescia, and E. Nardi, Phys. Rev. D 96, 075003 (2017b), arXiv:1705.05370 [hep-ph] .
  * Gaillard  _et al._ [2018] M. K. Gaillard, M. B. Gavela, R. Houtz, P. Quilez, and R. Del Rey, Eur. Phys. J. C 78, 972 (2018), arXiv:1805.06465 [hep-ph] .
  * Hook  _et al._ [2020] A. Hook, S. Kumar, Z. Liu, and R. Sundrum, Phys. Rev. Lett. 124, 221801 (2020), arXiv:1911.12364 [hep-ph] .
  * de Giorgi  _et al._ [2024] A. de Giorgi, M. F. Zamoro, and L. Merlo, (2024), arXiv:2402.14059 [hep-ph] .
  * Gan  _et al._ [2022] L. Gan, B. Kubis, E. Passemar, and S. Tulin, Phys. Rept. 945, 1 (2022), arXiv:2007.00664 [hep-ph] .
  * Landini and Meggiolaro [2020] G. Landini and E. Meggiolaro, Eur. Phys. J. C 80, 302 (2020), arXiv:1906.03104 [hep-ph] .
  * Alves [2021] D. S. M. Alves, Phys. Rev. D 103, 055018 (2021), arXiv:2009.05578 [hep-ph] .
  * Alves and Gonzàlez-Solís [2024] D. S. M. Alves and S. Gonzàlez-Solís, (2024), arXiv:2402.02993 [hep-ph] .
  * Bijnens and Ecker [2014] J. Bijnens and G. Ecker, Ann. Rev. Nucl. Part. Sci. 64, 149 (2014), arXiv:1405.6488 [hep-ph] .
  * Georgi  _et al._ [1986] H. Georgi, D. B. Kaplan, and L. Randall, Phys. Lett. B 169, 73 (1986).
  * Gasser and Leutwyler [1985] J. Gasser and H. Leutwyler, Nucl. Phys. B 250, 465 (1985).
  * Herrera-Siklody  _et al._ [1997] P. Herrera-Siklody, J. I. Latorre, P. Pascual, and J. Taron, Nucl. Phys. B 497, 345 (1997), arXiv:hep-ph/9610549 .
  * Kaiser and Leutwyler [2000] R. Kaiser and H. Leutwyler, Eur. Phys. J. C 17, 623 (2000), arXiv:hep-ph/0007101 .
  * Lu  _et al._ [2020] Z.-Y. Lu, M.-L. Du, F.-K. Guo, U.-G. Meißner, and T. Vonk, JHEP 05, 001, arXiv:2003.01625 [hep-ph] .
  * Oller and Meissner [2001] J. A. Oller and U. G. Meissner, Phys. Lett. B 500, 263 (2001), arXiv:hep-ph/0011146 .
  * Workman  _et al._ [2022] R. L. Workman _et al._ (Particle Data Group), PTEP 2022, 083C01 (2022).
  * Ecker  _et al._ [1989] G. Ecker, J. Gasser, A. Pich, and E. de Rafael, Nucl. Phys. B 321, 311 (1989).
  * Guo [2021] Z.-H. Guo, Eur. Phys. J. ST 230, 1659 (2021), arXiv:2109.14240 [hep-ph] .
  * Garcia-Martin  _et al._ [2011] R. Garcia-Martin, R. Kaminski, J. R. Pelaez, J. Ruiz de Elvira, and F. J. Yndurain, Phys. Rev. D 83, 074004 (2011), arXiv:1102.2183 [hep-ph] .
  * Peláez  _et al._ [2021] J. R. Peláez, A. Rodas, and J. R. de Elvira, Eur. Phys. J. ST 230, 1539 (2021), arXiv:2101.06506 [hep-ph] .
  * Mertig  _et al._ [1991] R. Mertig, M. Bohm, and A. Denner, Comput. Phys. Commun. 64, 345 (1991).
  * Shtabovenko  _et al._ [2016] V. Shtabovenko, R. Mertig, and F. Orellana, Comput. Phys. Commun. 207, 432 (2016), arXiv:1601.01167 [hep-ph] .
  * Shtabovenko  _et al._ [2020] V. Shtabovenko, R. Mertig, and F. Orellana, Comput. Phys. Commun. 256, 107478 (2020), arXiv:2001.04407 [hep-ph] .
  * Patel [2015] H. H. Patel, Comput. Phys. Commun. 197, 276 (2015), arXiv:1503.01469 [hep-ph] .