$a \to \pi \pi \pi$ decay at next-to-leading order in chiral perturbation theory

We discuss the construction of the two-flavour axion-pion effective Lagrangian at the next-to-leading order (NLO) in chiral perturbation theory and present, as a phenomenological application, the calculation of the decay rate of a GeV-scale axion-like particle via the channel $a \to \pi\pi\pi$. Through the NLO calculation, we assess the range of validity of the effective field theory and show that the chiral expansion breaks down just above the kinematic threshold. Alternative non-perturbative approaches are called for in order to extend the chiral description of axion-pion interactions.


Introduction
The main ingredient of the axion solution to the strong CP problem [1][2][3][4] is the axion coupling to a pseudo-scalar gluon density, which sets model-independent experimental targets for the axion mass and couplings to photons, nucleons, pions and electrons. Since the axion is much lighter than the scale of chiral symmetry breaking Λ χ 1 GeV and it has the same quantum numbers of the neutral pion, chiral perturbation theory (χPT) provides a natural framework to systematically derive axion properties. In fact, those were obtained long time ago by using leading order (LO) χPT (or equivalently current algebra) in a series of renowned papers [3,[5][6][7][8][9]. The axion chiral potential and coupling to photons at the next-to-LO (NLO) in χPT were computed in Ref. [10] (see also [11]), but it is only more recently that the program of "precision" axion physics has restarted with Ref. [12], also motivated by the booming of the axion experimental program (see e.g. [13,14]). State of the art axion mass calculations are now obtained by employing next-to-NLO (NNLO) χPT [15] or, alternatively, via lattice QCD techniques [16]. The axion-nucleon interaction Lagrangian instead has been derived in heavy baryon χPT up to NNLO [17,18]. Also CPand flavour-violating axion couplings have witnessed a resurgence of interest in the recent years, with new calculations based either on χPT or other non-perturbative approaches (see respectively Refs. [19][20][21][22] and [23][24][25]).
In this paper we focus on the axion-pion chiral Lagrangian at NLO. The latter was previously considered in Refs. [10,12] in the context of the axion potential, hence limited to non-derivative axion interactions, and more generally in Ref. [26], which included also derivative axion couplings. We here expand on the derivation of the NLO axion-pion chiral Lagrangian, by providing several details which were not presented in Ref. [26].
The most interesting application of this formalism consists in the calculation of the aπ → ππ scattering, which provides the dominant channel for axion thermalization in the early Universe [27,28], when the axion decouples from the thermal bath at temperatures below that of QCD deconfinement T c 155 MeV [29][30][31]. The highest attainable axion mass from cosmological constraints on thermally-produced axions is known as the axion hot dark matter bound. However, as shown in Ref. [26], the chiral expansion of the axionpion thermalization rate breaks down well below T c . Lacking for the moment a way to extrapolate the validity of χPT, a practical solution was given in Refs. [32,33] which proposed an interpolation of the thermalization rate starting from the high-temperature region above T c . See Refs. [34,35] for recent cosmological analyses adopting this latter approach.
Another application of the axion-pion chiral Lagrangian, which is the main subject of this paper, arises in the context of GeV-scale axion-like particles (ALPs) which dominantly decay hadronically as soon as the phase space for the channel a → πππ is kinematically open. For phenomenological studies related to this channel, see e.g. Refs. [36][37][38]. This process was computed at LO in χPT in Refs. [39,40] and the chiral expansion was claimed to be valid up to ALP masses of few GeV. However, by explicitly computing the NLO correction, we find that the effective field theory (EFT) breaks down much earlier, namely for ALP masses just above the kinematical threshold m a 3m π . Hence, in practice, χPT never yields an accurate description for the process at hand. The paper is structured as follows: in Sect. 2 we discuss the construction of the axionpion chiral Lagrangian, while the calculation of the a → πππ decay up to NLO in χPT is outlined in Sect. 3. We conclude in Sect. 4, where we advocate for possible strategies in order to extend the validity of the chiral description. Further details on the NLO calculation are provided in Apps. A-C.

Axion-pion effective field theory
The construction of the LO axion-pion Lagrangian was originally discussed in Refs. [6,9]. We first recall its basic ingredients (see also [27,[39][40][41]) in view of the extension at NLO, which was recently discussed in Ref. [26]. We here complement the latter derivation by providing several details which were omitted in Ref. [26]. In particular, we will focus on the 2-flavour formulation, which is best suited for the application to be discussed in Sect. 3. This is justified a posteriori, because the presence of strange mesons as external states is kinematically suppressed up to the energy scale at which the chiral expansion breaks down. On the other hand, the generalization to the 3-flavour case is in principle straightforward. In the following we will generically indicate both the QCD axion and the ALP as "axion", specifying when needed which case we are considering.

Axion-QCD effective Lagrangian
The 2-flavour axion effective Lagrangian in terms of quarks and gluons reads For the QCD axion m 2 a,0 = 0, while m 2 a,0 = 0 for the ALP case. 1 In the following, we will be especially interested in the case where m a,0 ∼ GeV, i.e. much larger than the pure QCD axion mass contribution. The couplings c 0 q = diag(c 0 u , c 0 d ) and g 0 aγ are model-dependent. For instance, in the case of the QCD axion, c 0 u,d = 0 and g 0 aγ = 0 in the KSVZ model [42,43], while c 0 u = 1 3 cos 2 β, c 0 d = 1 3 sin 2 β and g 0 aγ = α/(2πf a )8/3 in the DFSZ model [44,45] (with tan β = v u /v d the ratio between the vacuum expectation values of the two Higgs doublets present in the DFSZ model).
Upon an anomalous axial rotation of the quark doublet with Tr Q a = 1, the aGG term in Eq. (2.1) is shifted away, and the Lagrangian in Eq. (2.1) becomes where we have redefined the parameters as with Q EM = diag (2/3, −1/3).

Axion-pion effective Lagrangian at LO
At energies 1 GeV, the axion-QCD effective Lagrangian is replaced by the axion chiral Lagrangian, which at the LO reads (in the 2-flavour approximation, relevant for the observable studied in this paper) where f π = 92.21 MeV, χ a = 2B 0 M a (with B 0 denoting the quark condensate) and σ a (a = 1, 2, 3) the Pauli matrices. U = e iπ a σ a /f π is the pion Goldstone matrix, with The pion axial current, J a A, µ , reads at the LO (see App. A) defined in terms of the covariant derivative D µ U = ∂ µ U − ir µ U + iU µ , with r µ = r a µ σ a /2 and l µ = l a µ σ a /2 external fields which can be used to include electromagnetic or weak effects. The matching of the derivative axion term in Eq. (2.5) with the corresponding one in Eq. (2.3) has been obtained by rewriting where we used the Fierz identity σ a ij σ a kl = 2(δ il δ jk − 1 2 δ ij δ kl ). The iso-singlet current is associated to the heavy η and it can be neglected for our purposes, while the iso-triplet quark axial current is replaced with the pion axial current in Eq. (2.7).
In the following, we set Q a = M −1 q /Tr M −1 q , so that terms linear in a (including aπ 0 mass mixing) drop out from Eq. (2.5) and the only linear axion term arise from the derivative interaction with the pion axial current. Explicitly, the derivative axion coupling reads Expanding the pion axial current J a A, µ | LO = f π ∂ µ π a − 1 f π π 2 ∂ µ π a − 3 2f π π a ∂ µ π 2 + . . . , with π = π 0 π 0 + 2π + π − , the axion-pion derivative terms are given by The first operator introduces a kinetic mixing between the axion and the neutral pion, parametrized by the coefficient At the quadratic level the a-π 0 Lagrangian reads is the QCD axion mass squared at the LO. The procedure in order to diagonalize the quadratic Lagrangian in Eq. (2.12) consists of three steps: i) diagonalization of the kinetic term by an orthogonal transformation, ii) re-scaling of the fields to have a canonical kinetic term and iii) diagonalization of the mass matrix (rotated and re-scaled after steps i) and ii)). The first orthogonal rotation Therefore the re-scaling is given by (fields need to be multiplied by the inverse of W ) The action of R 1 and W on the mass matrix puts it in the form whose eigenvalues are m 2 π and m 2 a plus corrections of O( 2 ) for the pion and ALP masses, and O( 4 ) for the QCD axion mass (considering m a /m π ∼ O( ) in the QCD axion case). Denoting by R 2 the matrix that diagonalizes Eq. (2.18) as R T

2M 2
LO R 2 , one obtains that the complete rotation that needs to be applied to the fields (a, π 0 ) in order to fully diagonalize the quadratic Lagrangian in Eq. (2.12) is given by where (a phys , π 0phys ) denote fields with diagonal propagators. In the following, we drop the subscript "phys" when working in the diagonal basis. After the LO diagonalization procedure, the LO chiral Lagrangian containing the axion-pions interaction terms is given by (including the contribution due to Eq. (2.21) from the standard 4-pion Lagrangian) The QCD axion case is recovered in the m 2 a → 0 limit. Note that the correction due to the kinetic mixing in Eq. (2.21) can be safely neglected in the QCD axion case since m a m π .

Axion-pion effective Lagrangian at NLO
The axion-pion Lagrangian beyond LO requires two ingredients: the O(p 4 ) chiral Lagrangian with the axion-dressed coefficient χ a = 2B 0 M a (cf. Eq. (2.4)) and the derivative axion interaction with the NLO pion axial current. Part of the material of this Section was previously presented in Ref. [26]. The 2-flavour chiral Lagrangian at O(p 4 ) can be expressed in various equivalent bases. Here we stick to the expression given by Gasser and Leutwyler [46], which in the standard trace notation reads [47] L χ(NLO) The low-energy constants 1 , 2 , . . . , 7 are not fixed by chiral symmetry, but they need to be determined from experimental data or lattice QCD. The constants h 1 , h 2 , h 3 are coupled to pion-independent terms (see Eq. (2.26) below). The f R,L µν are the field strength tensors associated to the fields r µ and l µ appearing in the covariant derivative (see [46] for details). Since we are only interested in processes involving an even number of bosons, we neglect here the O(p 4 ) Wess-Zumino-Witten term [48,49] which features intrinsic-parityodd operators.
The NLO chiral left (right) current is obtained by differentiating the NLO Lagrangian with respect to the external field l µ (r µ ). Taking the axial combination of the chiral currents (see App. A) one obtains Being interested only in axion-pion interactions, from now on we will set to zero the field strength tensors as well as the external currents. Then the axion-pion Lagrangian up to NLO is given by the sum L χ(LO) a + L χ(NLO) a . Note that the NLO terms reintroduce a quadratic mixing of the axion field with the neutral pion. In App. B we explicitly repeat the diagonalization procedure at NLO, including as well one-loop terms from the LO chiral Lagrangian. In fact, the choice Q a = M −1 q /Tr M −1 q allows us to eliminate only some of the mass mixing terms at NLO. On the other hand, no axion-pion mixing arises from the term proportional to h 1 − h 3 in Eq. (2.24), since the latter does not depend on the pion field. This is readily seen by using the identity (2.26) The remaining axion-pion mass mixing is found to be (2.27) Considering instead derivative terms, at NLO the pion axial current gives rise to the following kinetic mixing term (2.28) Besides those tree-level mixings, the axion and the neutral pion also mix through one-loop diagrams, generated by the LO terms in Eq. (2.22).

a → πππ decay at NLO
As an application of the axion-pion chiral Lagrangian formalism at NLO we present here the calculation of the a → πππ decay rate, which shares some analogies with the case of aπ → ππ scattering discussed recently in Ref. [26]. The decay a → πππ is one of the leading hadronic channels for GeV-scale ALPs, and it has been previously computed at the LO in Refs. [39,40]. By means of the NLO correction we want to assess the convergence of the chiral expansion. Figure 1: One-loop diagrams contributing to the ALP decay a → π 0 π + π − .
The ALP decay rate in three pions is obtained by integrating the differential rate (see e.g. Sect. 48 in [50]) where there are two possible decay channels: a → π 0 π + π − and a → π 0 π 0 π 0 . In the following, we present the calculation of the ALP decay amplitudes and compare the LO to the NLO decay rate.

LO amplitude
The LO amplitudes at O(1/f a ) are obtained from the interaction terms in Eq. (2.22) and are found to be with the Mandelstam variables defined as (3.4) Note that the neutral pion channel (Eq. (3.3)) is proportional to m 2 a , since it stems entirely from a-π 0 mixing.

NLO amplitude
To compute the ALP decay into three pions at NLO we employ the Lehmann-Symanzik-Zimmermann (LSZ) formalism [51], according to which the amplitude is given by where the index α runs over the external particles, (i, j) = (+, −) or (0, 0), and Z a (Z π ) is the wave-function renormalization of the axion (pion) field defined via the residue of the 2-point Green's functions while the full 4-point Green's function is given by The first term is the amputated 4-point function, multiplied by the 2-point functions of the external legs with the axion mass set to zero. We work in a basis where the a-π 0 mixing has been diagonalized at the lowest-order, O(p 2 ), via Eqs. (2.20)-(2.21) and the remaining mixing, of O(p 4 ), is retained explicitly. Working with LO diagonal propagators, the 2-point amplitude for the a-π 0 system reads where Σ ij encodes NLO corrections including mixings. The 2-point Green's function is hence Expanding the diagonal terms around the physical masses we get (see Eq. (B.7)) with primes indicating derivatives with respect to p 2 . Then, by plugging Eq. (3.7) and (3.9) into the LSZ formula for the scattering amplitude and neglecting O(1/f 2 a ) terms, we obtain the ALP-decay amplitudes which are given by Defining the invariant mass of the two-pions systems π α -π β as (p π α + p π β ) 2 = m 2 αβ ≡ s, t, u with, respectively, (α, β) = (+, −), (0, +), (0, −), we obtain (3.13) (3.14) To carry out the renormalization procedure in dimensional regularization we shift the LECs as in Eq. (B.8) and we fix γ 1 = 1/3, γ 2 = 2/3, γ 3 = −1/2, γ 4 = 2 and γ 7 = 0, consistently with the values found in the literature for the standard chiral theory [46].

ALP decay rate: LO vs. NLO
At LO we reproduce the decay rates given in Refs. [39,40], that in our notation read  [54], m u /m d = 0.50 (2) [53], f π = 92.1(8) MeV [50] and m π = 137 MeV (corresponding to the average neutral/charged pion mass). Then the LO+NLO rates can be written as where the NLO functions g LO ij0 are obtained by numerically integrating the NLO amplitudes in Eqs. (C.1)-(C.2). Their profile is shown in the left panel of Fig. 2, for comparison with the LO counterparts. Although the expansion parameter in Eq. (3.19) is formally written as (m a /4πf π ) 2 , the actual calculation of the NLO rate shows (cf. right panel of Fig. 2) that the NLO correction becomes of the same order of the LO result already for ALP masses just above the kinematical threshold m a 3m π . This is reflected by a somewhat larger value of the NLO g-functions compared to the LO ones, as shown in Fig. 2. Thus we conclude the χPT description of the a → πππ decay rate breaks down for ALP masses much smaller than 4πf π 1.2 GeV. This earlier breakdown of χPT is also found in SM processes that are similar to the ALP decay into pions considered here, as e.g. η → πππ (see e.g. [55,56]). For instance, the NLO (NNLO) rate for η → πππ was found to be a factor ≈ 2.7 (4.5) larger than the LO one [55].

Conclusions
In this paper we have discussed the formulation of the axion-pion Lagrangian at the NLO in χPT and considered as an application of phenomenological relevance the ALP decay a → πππ, which is one of the main hadronic channels for GeV-scale ALPs. Through the inclusion of the NLO correction, we have estimated the range of applicability of the chiral expansion and found that the chiral EFT fails for ALP masses just above the kinematical threshold of 3m π (cf. right panel in Fig. 2). This result shows an earlier breakdown of the chiral EFT compared to naive expectations based on previous LO calculations, see Refs. [39,40]. We conclude that the range of applicability of the axion-pion chiral Lagrangian is rather limited for the problem at hand (similar conclusions were achieved for the case of aπ → ππ scattering in Ref. [26]) and hence alternative non-perturbative approaches (based either on dispersion relations or lattice QCD techniques) are called for in order to extend the validity of the chiral description.

Acknowledgments
We thank Guido Martinelli for many enlightening discussions on the subjects of this paper, and Ennio Salvioni for useful comments on the manuscript. We also thank Claudio Toni for spotting a mistake in the previous version of the paper. The work of L.D.L. and G.P. has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement No 860881-HIDDEN.

A Pion axial current
In this Appendix we provide the derivation of the pion axial current at the NLO. The currents associated to the Left and Right chiral rotations acting on the Goldstone matrix as U → RU L † , are easily computed promoting the global symmetries to local ones, and computing the variation δL of the Lagrangian under the given transformation. From Noether's theorem, the Left and Right currents are given by .

(A.2)
Let us consider the LO chiral Lagrangian To compute e.g. the Right current, we set Θ a L (x) = 0 and perform an infinitesimal Right transformation The variation of L χ is and therefore J µ,a R is given by With an analogous procedure one obtains The R − L combination of these two currents provides the pion axial current at LO The procedure can be repeated at the NLO by employing the shift of the O(p 4 ) chiral Lagrangian in Eq. (2.24), which yields Combining the left and right currents we obtain the axial current in Eq. (2.25).

B Axion-pion mixing at NLO
We explicitly perform here the NLO diagonalization of the axion and neutral pion propagators. The axion-neutral pion Lagrangian up to order 1/f a is given by where the subscript b stands for bare fields 3 and the interaction Lagrangian reads explicitly Note that L int contains all the terms which contribute to the two-point functions of the neutral scalar fields, i.e. LO tree-level mixings, LO terms giving the one-loop corrections and NLO terms. The latter provide the counterterms needed to reabsorb the loop divergences. We next define the renormalization conditions. Firstly, it is important to note that, since the divergences come from loops of L χ(LO) a , it is sufficient to extract the counterterms from 3 , 4 and 7 . Hence, m π and f π are the physical pion mass and decay constant at LO. Let us now denote by −iΣ ij (p 2 ) (with i, j = a, π 0 ) the 1-particle-irreducible (1PI) self-energy correction. The net effect of this correction is encoded in the effective Lagrangian where we employed the pion wave-function renormalization, π 0b → (1 + 1 2 δZ π )π 0 , defined as δZ π = ∂Σ ππ (p 2 )/∂p 2 . The one-loop self-energies Σ ij (p 2 ) can be computed from the interaction Lagrangian in Eq. (B.2). Defining with R = 2 d−4 − log(4π) + γ E − 1, and using dimensional regularization we find from which we get Therefore, we define the scale-independent parameters i and h i in such a way that the R + log(m 2 π /µ 2 ) factor is subtracted [46] i = γ i 3) we find that in order to renormalize Σ ππ and Σ aπ we need to set Thus the renormalized effective Lagrangian becomes (B.10) with We observe that 7 is not renormalized, since in the LO Lagrangian the m d − m u terms are all momentum dependent. So we are left with a non-zero off-diagonal two-point function. In order to eliminate the mixing, we can rotate the axion and the pion fields as a → a − β 1 π 0 , π 0 → π 0 − β 2 a , (B.13) yielding L eff a−π 0 → 1 2 a p 2 − m 2 a a + 1 2 π 0 p 2 −m 2 π π 0 − a β 1 (p 2 − m 2 a ) + β 2 (p 2 − m 2 π ) + p 2 3C aπf π 2f a + 4 7 (m d − m u )m u m d m 4 π f a f π (m u + m d ) 3 π 0 . (B.14) Hence, to cancel the mixing term it is sufficient to set (B.16)