Running effects on QCD axion phenomenology

We study the impact of renormalization group effects on QCD axion phenomenology. Focusing on the DFSZ model, we argue that the relevance of running effects for the axion couplings crucially depends on the scale where the heavier Higgs scalars are integrated out. We study the impact of these effects on astrophysical and cosmological bounds as well as on the sensitivity of helioscopes experiments such as IAXO and XENONnT, showing that they can be sizable even in the most conservative case in which the two Higgs doublets remain as light as the TeV scale. We provide simple analytical expressions that accurately fit the numerical solutions of the renormalization group equations as a function of the mass scale of the heavy scalars.


Introduction
Axions are an intrinsic prediction of the Peccei-Quinn (PQ) mechanism [1,2], which remains, after over four decades, the most appealing solution to the strong CP problem.This problem arises because Quantum Chromodynamics (QCD) predicts CP violating effects that are not observed experimentally.The PQ mechanism involves a new global chiral symmetry U(1) PQ , which is anomalous under QCD and spontaneously broken at a large energy scale f a (PQ scale).The axion is the Nambu-Goldstone boson associated with the spontanous breaking of this symmetry [3,4], and is characterised by the fact that all its interactions are inversely proportional to f a .Although the original Weinberg-Wilczek model [3,4], in which the scale of PQ breaking coincides with the electroweak (EW) symmetry breaking scale, was quickly ruled out, new viable models emerged early on, in which the PQ scale can be arbitrarily high so that all axion interactions can be sufficiently suppressed, yielding the so-called invisible axion.Two examples are particularly appealing for their simplicity: the KSVZ or hadronic axion [5,6] and the DFSZ axion [7,8].The main difference between KSVZ and DFSZ-type axions is that the former do not couple to ordinary quarks and leptons at the tree level.Though many other possible axion models have been considered in the literature (see Ref. [9] for a comprehensive overview), the two above-mentioned models are by far the most studied ones and are universally regarded as benchmark QCD axion models.In recent years continuous progress in experimental technologies has brought within reach the possibility of detecting in terrestrial experiments the invisible axions arising in these models.This has stimulated a tremendous interest in this field, with several new theoretical and phenomenological studies, as well as a wealth of new experimental proposals (see [10,11] for recent reviews).In the meanwhile, ongoing experiments have already started to probe the benchmark KSVZ/DFSZ axion models, and in the coming decades they will dig deep into the relevant parameter space region.From the theory side, this calls for the development of "precision axion physics", which will turn out to be crucial in the case of an axion discovery.Indeed, from a given determination of the low-energy axion couplings to photons and other matter fields (such as electrons and nucleons) 1 arXiv:2305.11958v2[hep-ph] 27 Dec 2023 one would like to infer the structure of the high-energy theory, that is the ultraviolet (UV) completion of the axion effective field theory (EFT).This step is highly non-trivial, since it entails a large separation of scales, from the typical low-energy scale of axion experiments up to the PQ scale, f a ≳ 10 8 GeV.Hence, axion related physical quantities, as for example the axion couplings to Standard Model (SM) fermions, are potentially affected by large radiative corrections, which can induce large deviations from the tree-level expressions.In the case of the KSVZ model, it was pointed out long ago [12,13] that although the axion coupling to electrons is zero at tree level, a non-zero electron coupling can be sourced via loop corrections by the axion-photon coupling and, more recently, it was shown that the leading correction to this coupling is generated at even higher orders via the anomalous axion coupling to gluons [14].Nowadays, the full one-loop anomalous dimensions for the d = 5 axion effective Lagrangian have been computed [15][16][17][18], while running effects have been investigated for the benchmark DFSZ/KSVZ axion models in Ref. [14], and for the so-called astrophobic axion models (which feature nonuniversal PQ charges [19,20]) in Ref. [21].
The purpose of this work is to study QCD axion phenomenology in light of renormalization group (RG) effects, focussing for definiteness on the mass window m a ∈ [meV, eV].This region of parameter space shows a remarkable complementarity among the existing bounds on the different axion couplings, namely to photons, electrons, nucleons and pions, stemming from helioscope searches, as well as from astrophysics and cosmology.It is therefore an ideal playground where to investigate the consequences of running effects for QCD axion phenomenology.In particular, we will focus on the large corrections induced by the top Yukawa coupling, which apply to a large class of axion models where the SM fermions are charged under the U(1) PQ symmetry.A paradigmatic example is the universal DFSZ model, which features two Higgs doublets and one SM singlet scalar, and whose axion parameter space at tree level depends solely on m a and tan β.However, top-Yukawa radiative corrections induce a logarithmic dependence of the effective axion couplings on the mass scale of the heavy scalar degrees of freedom of the two Higgs doublet model (2HDM) (the issue of large logarithmic corrections is well known in the 2HDM literature, see e.g.[22,23]) that can range from about 1 TeV up to the PQ scale, f a .As we shall see, these corrections are often large and may skew the parameter space region that is effectively probed by terrestrial experiments and by astrophysical/cosmological observations.The paper is organized as follows.In Sect. 2 we discuss the structure of top-Yukawa radiative corrections, and we provide approximate analytical expressions for the dependence of the axion couplings on these corrections.Sect. 3 is devoted to study the impact of running effects on QCD axion phenomenology, including the consequences for astrophysical and cosmological limits as well as for the sensitivity of future axion experiments.We conclude in Sect. 4. Details on the solutions to the RG equations are provided in Appendix B.

Running QCD axion couplings
Of central interest for axion phenomenology are the axion couplings to photons and matter fields (electrons, nucleons, as well as other hadrons relevant for axion production).They are defined by the following interaction Lagrangian: where F µν denotes the electromagnetic field strength, Fµν = 1 2 ϵ µνρσ F ρσ (with ϵ 0123 = −1) its dual, f = p, n, e runs over low-energy matter fields, and C γ, f, π, πN, N∆ are O(1) dimensionless coefficients.The axion-pion coupling in the second line of Eq. ( 1) (with f π = 92.1(8)MeV [24] the pion decay constant) is of phenomenological relevance for thermal axion production in the early Universe, while the axion contact interactions with pions and nucleons (third line) and with ∆-resonances (fourth line) are important for axion production in Supernovae (SNe).The ellipses stand for other possible axion interaction terms which will not be considered in this paper.
In the context of axion phenomenology, one usually employs the dimensional couplings g aγ = α 2π C γ / f a and g a f = C f m f / f a .In particular, C γ = E/N − 1.92 (4), where E/N is the ratio between the electromagnetic and QCD anomalies of the PQ current (for typical values in concrete axion models, see e.g.[25,26]).Note that to a very good approximation the anomalous axion-photon coupling is insensitive to running effects, with first corrections appearing at three loops [17].Moreover, mass dependent corrections to the effective axion-photon coupling are safely negligible for m a ≪ m e since they scale at most as (m a /m e ) 2 -see e.g.Ref. [27].Hence, in the following, we will only focus on radiative corrections to the axion couplings to electrons and hadrons.
Axion-hadron interactions can be expressed in terms of the model-independent axion gluon coupling (which fixes the absolute normalization in terms of f a ) and the axion couplings to quark fields, q = u, d, s, c, b, t, defined via the Lagrangian term1 In terms of the latter, the axion couplings to hadrons defined in Eq. ( 1) read (see e.g.[9,[28][29][30]) where and ∆ s = −0.035(6)(7)(at 2 GeV in the MS scheme) are the N f = 2 + 1 FLAG average [31], that is dominated by the results of [32], and z = m u (2 GeV)/m d (2 GeV) = 0.49(2) [33].
Running effects on the low-energy couplings of the axion to first generation SM fermions can be parametrized as2 (see e.g.[14]) C e (m e ) = C e ( f a ) + ∆C e , with and Ψ = u, d, e.The parameter r t Ψ (m BSM ) encodes the RG correction approximated by taking only the top-Yukawa contribution, and depends logarithmically on the parameter m BSM ≃ m H, A, H + that denotes collectively the mass scale of the heavy scalar degrees of freedom (we implicitly assume for the heavy modes of the scalar doublets the decoupling limit [34], in which all the heavy masses are approximately degenerate).The m BSM scale depends on the structure of the DFSZ scalar potential, whose details (see e.g.[35,36]) are not crucial for the calculation of the axion RG equations, and we take it to range from about 1 TeV (the approximate lower bound as set by LHC searches for new heavy scalars) up to f a .
Note that as long as the couplings are considered at a renormalization scale µ above m BSM there are no top-Yukawa running effects.This is because in this regime the axion couplings to the SM fermions correspond to the global charges of the PQ current, which is classically conserved, and thus they do not renormalize.For µ < m BSM we enter a different regime, in which Higgs doublets with different PQ charges mix to give rise to heavy scalars (which are integrated out) and to the light Higgs, that has no well-defined charge.In this effective theory there is no more a conserved PQ current, and running effects for the axion-fermion couplings can kick in.This is the reason why the largest RG effects appear when the BSM scale is taken at the largest possible scale m BSM ∼ f a .Contrary, when the 2HDM structure keeps holding all the way down to the TeV scale, running effects are much less sizeable.
In Appendix B we provide a fit to r t Ψ (m BSM ) obtained by interpolating the numerical solution to the RG equations (cf.Eqs.(B.6)-(B.7)and Tab.B.4). Taking, for instance, m BSM = f a = 10 10 GeV one finds

Analytical understanding of RG running effects
To understand the phenomenological impact of the RG corrections to the axion couplings, and to compare it with the current experimental sensitivity, it is convenient to provide some analytical approximations.To this aim, it is advantageous to introduce the iso-scalar (C 0 ) and iso-vector (C 3 ) nuclear couplings (see also [21]), defined as follows: where the right-hand sides are obtained from the expressions for C p, n given in Eqs.(3)-(4).From Eqs. ( 5)- (6) we see that all the other couplings are proportional to the iso-vector combination The RG correction to the iso-vector combination ∆C 3 ≃ 0.64 C t ( f a ) (r t u − r t d ) may be sizeable, with the exact value depending on the m BSM scale.An excellent fit to the combination r t u − r t d , for m BSM in the range 1 TeV to 10 18 GeV, is given by with x = log 10 (m BSM /GeV).This expression reproduces our numerical results with a precision better than 2% (see Appendix B).Then, in the relevant range for m BSM , we have 0 Since in universal axion models we expect C 3 ∼ C t (the exact relation depending on the model parameters), we can conclude that ∆C 3 /C 3 can be of the order of a few 10%, and even larger.For example, in the case of the DFSZ axion (to be discussed below), we find which can become quite significant at small tan β.
On the other hand, the RG correction to the iso-scalar coupling C 0 is, in general, very small.From Eq. ( 14), we see that this coupling combination gets contributions from (r t u + r t d ) and from ∆ s C s .As it was pointed out in Ref. [21], in the leading approximation in which only the contribution of the top-quark Yukawa coupling is kept, the combination (r t u + r t d ) is characterized by a strong cancellation, see for example Eqs. ( 11)-( 12), and is numerically very small ∼ 0.2%. 3Hence, eventually the leading correction to C 0 comes from the RG correction ∆ s ∆C s to the last term in Eq. ( 14).It is easy to estimate this contribution from our general results, using r t s = r t d that holds for universal models.In the end, we find that the RG corrections to C 0 are only about 3% of the corresponding corrections to C 3 and hence this combination of couplings (and the corresponding iso-scalar axion coupling to nucleons, g aN,0 = C 0 m N / f a ) is practically unaffected by RG running effects.

DFSZ axion couplings beyond tree level
The scalar sector of DFSZ models [7,8] features a SM singlet complex scalar Φ and two Higgs doublets H u,d that couple respectively to up-and down-type quarks in a generationindependent way.Under SU(3 The threefold re-phasing symmetry of the scalar sector U(1 Y by a renormalizable non-Hermitian operator that can be chosen as 4 There are two possible variants of the model, depending on whether the lepton sector couples to H d (DFSZ1) or to Hu = iσ 2 H * u (DFSZ2).For a review see Sec. 2.7.2 in Ref. [9].The Yukawa sector of the DFSZ1 model contains the following operators where a sum over generation indices i, j = 1, 2, 3 is left understood, and q i , ℓ i denote the quarks and leptons SU(2) L lefthanded (LH) doublets while u j , d j , e j the right-handed (RH) singlets.The corresponding coefficients for the axion coupling at the UV scale f a are with The domain in which tan β is allowed to vary is obtained by requiring that the DFSZ Yukawas remain perturbative up to scales of O( f a ).This corresponds to imposing perturbative unitarity on Higgs-mediated 2 → 2 SM fermion scatterings (see e.g.[37]) up to f a .The perturbative domain is evaluated by evolving the values of the gauge couplings and of the SM Yukawa couplings at m Z given in Ref. [38] up to the scale m BSM employing twoloop RG equations.For m BSM ∼ f a the SM Yukawa couplings are RG-evolved from m Z to f a , and upon matching with the DFSZ couplings . This yields the perturbative domain 3 From the more accurate numerical analysis in Appendix B we obtain that for any value of the m BSM scale |r t u + r t d |/|r t u − r t d | ≲ 0.5%. 4The first possibility yields a number of domain walls N DW = 6 while the second N DW = 3, but they remain otherwise indistinguishable from the point of view of low-energy phenomenology.
Note that the perturbative domain of tan β has a mild (logarithmic) dependence on the PQ scale f a .This is shown in Fig. 1 for the low tan β region (a similar dependence is present also for the large tan β region, where running effects are however less important).The Yukawa sector of the DFSZ2 model contains instead the following operators and the corresponding axion coupling coefficients are with tan β defined in the same perturbative domain as in DFSZ1.
Let us now proceed to discuss the impact of running effects in the DFSZ model.Approximate RG corrections to the axion couplings are collected in Tab. 1.Note that in the case of DFSZ1 the iso-vector combination C 3 , as well as the axionelectron coupling C e , receive large corrections at small tan β,

Perturbative unitarity Perturbative unitarity
Tree-level coupling

Perturbative unitarity Perturbative unitarity
Tree-level coupling that is when the tree-level coupling vanishes. 5On the other hand, the RG corrections on the iso-scalar combination C 0 remain small in the whole tan β range.This is also displayed in Fig. 2, where the dashed line corresponds to the tree-level result, while the red band encodes the range of RG corrections obtained by varying m BSM between f a = 10 9 GeV (lower border of the red region) and 1 TeV (upper border of the red region).We see that although for m BSM = 1 TeV the running couplings trace closely the tree-level couplings as long as tan β > 1, also in this case RG corrections become non-negligible at small tan β.

RG effects on QCD axion phenomenology
In the following section, we will discuss the phenomenological impacts of the RG corrections.In particular, we will see how these affect axion astrophysical and cosmological bounds, as well as the sensitivity of terrestrial experimental searches.For clarity, we will refer to the DFSZ axion models, even though our results can be applied to other axion models.
Several observables depend dominantly (or entirely) on C 3 and thus, as discussed in Sect.2, are subjected to large RGinduced modifications.These include for example the coupling to pions, which is responsible for the axion thermalization in the early Universe (via ππ ↔ πa), which controls the hot dark matter (HDM) bound discussed in Sect.3.2.More recently, it has also been shown that the pion-nucleon scattering may be responsible for a large contribution to the axion emission rates in dense media, particularly in SNe [39][40][41].Finally, the isovector coupling C 3 is entirely responsible for the nuclear reaction process p + d → 3 He + a, which is one of the most efficient and widely studied production mechanisms of axionlike particles from solar nuclear reactions [42][43][44][45][46]. On the other hand, the nucleon coupling to 57 Fe, relevant for axion production through nuclear de-excitations in the Sun [42,43,47], turns out to be less sensitive to RG corrections (see Tab. 2).
The axion-electron coupling C e , which plays a significant role in astrophysics (see Section 3.1) as well as in terrestrial experimental searches (see Section 3.3), is also subjected to large RG corrections.In fact, r t e ≈ r t d . 6Thus, using r t u + r t d ≈ 0 and Red Giants Red Giants Eq. ( 16), we find where x = log 10 (m BSM /GeV) parameterizes the new physics scale. 7Therefore the corrections to C e can also be expressed in terms of ∆C 3 .Specifically, from our previous results we see that Thus, at our level of approximation all running effects can be expressed in terms of ∆C 3 , which for DFSZ axions is ∝ cos 2 β.A general consequence of this observation is that, in the case of DFSZ axions, the running effects are mostly relevant only at small tan β, something that is apparent from our numerical analysis. 8 Before moving to the phenomenological study, it is instructive to anticipate how RG corrections will modify the usual DFSZ bands, obtained by varying the value of tan β within the perturbative unitarity limits, for the g ae , g ap , g an couplings (see for example the section on axions in the Review of Particle Physics [48]).This can be easily estimated by considering the RG corrections to the electron and to the nucleon couplings C p,n = C 0 ±C 3 given in Table 1.The results are shown in Fig. 3, where we have taken m BSM = f a ∝ 1/m a to maximize the RG effects.In each panel the bands correspond to varying tan β in the interval [0.14, 100].The first figure shows the modification of the usual g ae band in DFSZ1.For this case we obtain the most dramatic effect, that is a marked reduction of the width of the band that, after including RG corrections, shrinks down to the green hatched region.This is due to the fact that the treelevel suppression of g DFSZ1 ae in the limit sin β → 0 is cut-off at it was argued that r t u + r t d ≃ 0. Deviations from these relations are due to subleading contributions of Yukawa terms other than Y t .A detailed discussion can be found in Section 3 of Ref. [21]. 7This result should not be surprising since it holds in the limit of |r t u +r t d | ≈ 0 and, as discussed above, |r t u + r t d |/|r t u − r t d | ≲ 0.5%. 8The only exception is in cases where the cos β dependence cancels, as for ∆C e /C e in the DFSZ2 model (see the Red Giant bound on DFSZ2 axion in the right panel of Fig. 4).small tan β by the RG correction proportional to cos 2 β.Instead there is not such a dramatic effect for g DFSZ2 ae since both the tree level coupling and the RG correction are proportional to cos 2 β.In the second panel in Fig. 3 we show the RG effect on the band for g ap (that is the same in DFSZ1 and in DFSZ2).In this case we see that the shrinking of the allowed band is much less pronounced.Finally, the third panel shows the RG effects on the g an band.In this case the allowed band is sizeably widened, which this is due to a cancellation in the tan β independent part of the coupling, which enhances the overall dependence on this parameter.
The most important phenomenological consequences of the RG corrections to the axion couplings will be analyzed in the following sections.

Astrophysical constraints
In this Section, we discuss the impact of RG corrections to astrophysical observables.For reference, we will mostly focus on the DFSZ1 axion model.The analysis for DFSZ2 goes along similar lines.
Axions can be copiously produced in stars, mostly due to thermal processes (see Ref. [53] for a recent review).Here we will not consider astrophysical bounds on the axion-photon coupling [60][61][62][63], since C γ does not receive any relevant RG correction.We focus instead on the axion-electron and on the axion-nucleon couplings.
The most stringent astrophysical bound on the axion-electron coupling is derived from observations of the tip of the red giant branch (RGB) in globular clusters.The production of axions during the RGB evolution cools the core, playing a role similar to that of neutrinos, and thus delays the helium ignition.The delay leads to a larger helium core and, consequently, to a higher stellar luminosity.Thus, comparison between observations and predictions for the luminosity of the RGB tip (the brightest stars in the RGB) is an efficient way to test anomalous channels of stellar cooling.The most recent analysis has set the constraint |g ae | ≤ 1.48 × 10 −13 (at 95% C.L.) [64] Axion production/detection in 57 Fe [42,43,47] C e ∆C e = −0.78∆C 3 Axion production in stars [53] Axion detection (Xenon etc.) [56][57][58] C γ ∆C γ = 0 Axion production in stars and labs [9] Most axion detection experiments [10,11] Axion coupling for Sikivie helioscopes [9,59] Table 2: RG corrections (second column) to axion couplings listed in the first column, in the approximation of keeping only the contribution from the top Yukawa coupling Y t .Note that in this approximation all the various corrections can be expressed just in terms of ∆C 3 given in the second line, with r t u − r t d given in Eq. ( 16).
this coupling, g ae = C e m e / f a , the RGB bound translates into This relation provides an upper bound for the axion mass at any given value of tan β and x = log 10 (m BSM /GeV).The full numerical results are shown in Fig. 4, where the red-shaded bands incorporate the range fixed by the possible values of We can gain some intuition about these effects using our approximate results, shown in Table .1.In the case of the DFSZ1 model (left panel in Fig. 4) we can conveniently rewrite the RGB bound on the axion mass as follows 9 where l(x) = ln √ x − 0.52 .The first important observation is that in the limit l(x) → 0, that is ignoring the RG corrections, the RGB bound on the mass is a monotonic function of tan β and disappears in the limit of small tan β.This result is modified by the RG corrections, which in the limit tan β → 0 still provide a useful limit on the axion mass, m a ≤ 0.018 eV/l(x).From these considerations, we can conclude that the RG correction to the RGB bound becomes particularly important in the low tan β limit, a result confirmed by the full numerical result shown in Fig. 4. The most conservative value for the RGB bound corresponds to m BSM = 1 TeV, for which we obtain, in 9 Notice that in the range of m BSM we are considering here, the absolute value in Eq. ( 26) is unnecessary.our approximation, m a ≤ 0.018 eV/l(3) ≃ 8.75×10 −2 eV, which agrees well with the complete numerical result shown in the left panel of Fig. 4. In the case of the DFSZ2 model instead, there is not such a striking feature, and this is because in this case both the coupling g ae and its RG correction depend on cos β.The RGB bound for this case is shown in the right panel of Fig. 4.
Let us now move to the axion-nucleon coupling and analyze the axion bound from SN1987A. 10 This is a quite more complex problem since axion production in a SN environment, at temperatures of order ∼ 30 MeV and densities in excess of ∼ 10 14 g/cm 3 , gets contributions also from pions [29,39,40] and from ∆ baryon resonances [30].Following the notation of Ref. [41], the effective low-energy axion-nucleon interaction relevant for the axion production processes in a SN environment is given in Eq. ( 1).The second term in the first line describes the usual axion-nucleons interactions.The third line contains the axion-pion interactions.The four-particle interaction vertex in the third line accounts for the pion-nucleon contact term recently discussed in Ref. [29], while the last line accounts for the axion couplings to the ∆-resonances, whose contribution to the axion emissivity has been recently calculated in Ref. [30].The interaction Lagrangian in Eq. ( 1) can be used to compute the axion emissivity due to nucleon-nucleon (NN) bremsstrahlung, NN → NNa, as well as the Compton-like pion scattering processes, π − p → na, including also the contribution from the ∆ resonances (see Ref. [41] for an updated overview).
In general, the axion luminosity from a SN, L a , depends only  27) and (28).The coefficients are calculated at a post-bounce time of 1s (see Ref. [41]).
The first row refers to the NN bremsstrahlung contribution only [66], ignoring the pion scattering processes and the ∆ resonance contribution.The second row gives the total contribution, calculated from the results in Ref. [41].The mass parameter m is defined in Eq. ( 30).
on a particular combination of C 0 and C 3 , and thus only this combination can be constrained.The luminosity can be expressed as [41] L where ϵ 0 is a numerical factor and The numerical values of the coefficients ϵ 0 , a and b can be found in Table 3.In the table we present, in the first line, the results obtained by considering only the purely NN bremsstrahlung production, which corresponds to the first line of Eq. ( 1).The results for the total emission rate are given in the second line (we remind to the reader that up until very recently, the NN bremsstrahlung production was the only process considered for estimating SN axion emission rate.)Notice that, as evident from Table 3, the addition of the pion-induced scatterings increases the relative importance of C 3 (controlled by the coefficients b and c) and thus enhances the effects of the RG corrections.More specifically, from Eq. ( 28), and assuming ∆C 0 ≈ 0, we find As expected, the RG effects are reduced in the case of purely NN-bremsstrahlung production, due to the partial cancellation between the b and c terms. 11mposing L a ≤ L ν = 3 × 10 52 erg/s [41], we find the bound on the axion mass In the case of for DFSZ axions, this bound is shown in Fig. 4. Specializing on the DFSZ axion case, we immediately find from Tab. 1, and Note that the above expression is never larger than about 10%-15%.Thus, for the SN bound RG effects are somewhat less

Thermal axion cosmology
If axions are in thermal equilibrium until below the quarkhadron phase transition (which can occur for m a ≳ 0.1 eV) the axion thermal population will give a sizeable contribution to the effective number of extra relativistic degrees of freedom [67], ∆N eff , that is constrained by Big Bang Nucleosynthesis (BBN) [68] and cosmic microwave background (CMB) observations [69,70].The highest attainable axion mass from such cosmological constraints is also known as the Hot Dark Matter (HDM) bound.The forecast sensitivity of the planned CMB-S4 [71] and Simons Observatory (SO) [72] surveys will fully cover the mass range in which the axion decouples below or across the QCD crossover, thus a precise determination of the axion-pion thermalization rate, including running effects, would be necessary to set definite targets. 12or T ≲ T c , where T c ≃ 155 MeV is the QCD deconfinement temperature, the dominant thermalization channels is aπ ↔ ππ [13,49].It has been recently shown, however, that the standard computation of this process, that is based on chiral perturbation theory (ChPT), breaks down for T ≳ 70 MeV [50,52].Phenomenological extensions of the validity of the chiral expansion, based on unitarization methods, have been proposed in Refs.[51,52].
In the following, we will consider the unitarized thermal rate based on the Inverse Amplitude Method (IAM), recently discussed in Ref. [52], which gives the thermal scattering rate: with C π given in Eq. ( 5) and m π = 137 MeV representing the average neutral/charged pion mass.The numerical function h IAM is provided in Ref. [52] (cf.Fig. 3 of this reference) and is normalized to h IAM (m π /T c ) = 1.
We will estimate the impact of RG effects on the HDM bound relying for simplicity on the instantaneous decoupling condition Γ a (T D ) ≃ H(T D ), with Γ a (T ) the axion-pion scattering rate given in Eq. ( 33) and H(T ) = 4π 3 g ⋆ (T )/45 T 2 /m pl the Hubble rate, where m pl = 1.22 × 10 19 GeV is the Planck mass and g ⋆ (T ) the effective number of relativistic degrees of freedom. 13he axion contribution to the effective number of extra relativistic degrees of freedom is given by [67] ∆N eff ≃ 4 7 with T a /T ν the ratio of the axion to neutrino temperature at T ≪ 1 MeV (i.e.well after ν-decoupling) and g S (T D ) the number of entropy degrees of freedom at axion decoupling, that in the last relation has been normalised to the total number of SM degrees of freedom g S (T > m t ) = 106.75.We then confront Eq. ( 34) with the bound on ∆N eff from Planck's 2018 data [69,70], and from this we extract a bound on the axion mass.
Our results for the HDM bound in the DFSZ1/2 models are summarised in Fig. 5, where we show the tree-level results compared with the RG corrections included.Again we see that RG effects are especially important at small tan β.In Fig. 5 the DFSZ1 and DFSZ2 cases coincide because the (subleading) effects of scattering off leptons have been neglected.In Ref. [75] it was argued that thermalization channels involving axion scattering off leptons can become relevant in DFSZ2 at small tan β.However, since RG corrections keep the axion-pion coupling sizeable also in this regime, the effect of lepton scattering becomes accordingly less important.

Helioscope experiments
One of the most appealing result of the RG correction analysis is the implication for the next generation of experiments hunting for solar axions.The main reason is that the solar flux is strongly dependent on the axion-electron coupling and, as we have seen, this can receive large RG corrections.As a consequence, helioscope sensitivities to DFSZ axions, that have been so far estimated using tree level electron-axion couplings, have been underestimated.
Here, we focus mostly on the Sikivie type of axion helioscopes [59].This kind of experiment is designed to detect solar axions by converting them in X-ray photons using a large laboratory magnetic field.The importance of the axion-electron coupling for Sikivie's helioscope sensitivity to solar axions is expressed by the following relation [9] where g γ10 = g aγ /10 −10 GeV −1 , g e12 = g ae /10 −12 , and g γ10 is the helioscope sensitivity to g aγ (again, in units of 10 −10 GeV −1 ).Notice that g aγ is, in general, a function of the axion mass.Defining the effective coupling the above expression leads to the following heliscope sensitivity relation which, in the case of the DFSZ axion, can be readily translated into a limit on the tan β accessible to the helioscope as a function of the axion mass.Notice that, according to this expression, the DFSZ sensitivity to tan β (which enters only through C e ), should disappear for C γ ≫ 37C e , which is fulfilled for tan β ≪ 0.25.In general, if the helioscope sensitivity is good enough, there could be mass regions where the entire range of tan β is accessible.
To give an example of an application of Eq. ( 37), let us consider the case of BabyIAXO, a next-generation axion helioscope presently under construction [76].Its sensitivity at m a = 0.1 eV is expected to reach g γ10 = 0.33.Using this value and C γ = 8/3 − 1.92 for the DFSZ1 model, if RG effects are ignored, one would conclude that at this mass value BabyIAXO could be sensitive to the region tan β ≳ 0.62.The results of our complete numerical analysis for all values of the axion mass are plotted in the left panel in Fig. 6, where the dashed line contours correspond to the estimated BabyIAXO sensitivity if RG effects are ignored.The reach of the more advanced helioscope experiment IAXO [77] is also shown in the left panel in Fig. 6.In this case we see that there is a mass region for which the experiment is sensitive to all values of tan β.When RG corrections are ignored, this region extends to masses between ∼ 50 meV and ∼ 200 meV.The reason of this is that IAXO is sensitive enough to see the solar axion flux even in models in which axions are only coupled to the photon and not to the electron.
Let us now consider the effects of RG corrections on the projected sensitivities.As shown in Tab. 2, the RG corrections to the effective helioscope coupling is which is valid in the limit of ∆C e /C e ≪ 1.This condition is always verified in DFSZ2, while for DFSZ1 it holds for tan β ≫ 0.5 l(x) 1/2 (Cf.Tab. 1).Since in the case of BabyIAXO the expected sensitivity is not sufficient to detect DFSZ axions unless tan β ≳ 0.6 (see Fig. 6) which implies C γ /(37C e ) ≪ 1, we can simplify the correction to the effective coupling to which can be readily estimated using our results from the Tab 1.
Notice that this correction can be quite sizable and implies that the reach of BabyIAXO to small electron couplings (low tan β values) can be pushed down significantly, as is shown by the blue region in the left panel of Fig. 6.
The impact of RG corrections on the IAXO sensitivity to DFSZ1 axions is also shown in the left panel of Fig. 6, and corresponds to the red regions.In this case we notice an interesting effect, that is that the IAXO reach in the region of small tan β is sizeably enlarged for all values of m BSM , since the solar axion flux is necessarily larger than what predicted ignoring the RG corrections.As a result, the mass region for which IAXO is sensitive to the entire range of tan β is extended.
Finally, the correction to the axion-electron coupling has also an obvious impact on experiments which detect axions through the axio-electric effect.Such experiments include large underground detectors such as Panda-X [57], LUX [56], or XENON-nT [58], originally designed for dark matter searches.The RG modification of the axion-electron coupling extends the potential of these experiments to explore the DFSZ parameter space.Our numerical results in the case of XENON-nT are shown in the right panel of Fig. 6.A fundamental difference with respect to the previous results is that, because of the RG-induced corrections, in principle XENON-nT could be sensitive to DFSZ1 axions for any value of tan β.However, the current experimental sensitivity is insufficient to reach inside the mass region allowed by the RGB bound discussed in Sect.3.1 (see the left panel in Fig. 4).

Conclusions
In this paper, we have studied the impact of RG effects on QCD axion phenomenology, focusing on DFSZ models.We have shown that running effects on axion couplings depend crucially on the scale at which the heavy Higgs states are integrated out, and the 2HDM effectively reduces to the SM with a single light Higgs.We have discussed the implications of running axion couplings on astrophysical and cosmological bounds, as well as the sensitivity of helioscope experiments such as (Baby)IAXO.We have found that running effects are sizable even in the most conservative case in which the 2HDM structure keeps holding down to the TeV-scale, and thus they can never be neglected.We have also provided simple analytic expressions fitted to reproduce the numerical solutions of the RG equations, which can be a useful tool for studying the implications of running axion couplings.In the case of an axion discovery, running effects might prove to be crucial in order to reconstruct the axion UV completion.
In order to take into account running effects it is convenient to adopt the Georgi-Kaplan-Randall (GKR) field basis [78], where the PQ symmetry is realised non-linearly, so that under a U(1) PQ symmetry transformation all fields are invariant except for the axion field, which changes by an additive constant a → a + α f , that is where ) † H u,d and c q L , . . .are diagonal matrices in generation space.Note that in the effective field theory below f we have neglected the heavy O( f ) radial mode of Φ and we focused for simplicity on the 2HDM.In order to match an explicit axion model to the effective Lagrangian in Eq. (A.1) at the high scale µ ∼ O( f ), we perform an axion dependent field redenfinition: ψ → e −iX ψ a/ f ψ, where ψ spans over all the fields, and X ψ is the corresponding PQ charge.Due to U(1) PQ symmetry, the non-derivative part of the renormalizable Lagrangian is invariant upon this field redefinition, while the d = 5 operators in Eq. (A.1) are generated from the variation of the kinetic terms and from the chiral anomaly.The couplings are then identified as where in the second equation c ψ R,L refer to the charges of the chiral fermion fields. 14For the DFSZ1 model introduced in Sect.2.2, with e.g. the operator H u H d Φ † in the scalar potential, the universal charges X ψ can be set to where X H u = c 2 β and X H d = s 2 β .The DFSZ2 model features a similar charge assignment, with instead X e = X H u .For the anomaly coefficients in Eq. (A.Running effects induced by Yukawa couplings (and in particular by the top Yukawa which is the most relevant one) only occur below the scale of the heavy radial modes of the 2HDM, denoted as m BSM ≃ m H, A, H ± , with the heavy scalars assumed to be degenerate in the decoupling limit (see e.g.[34]).This is due to the fact that as long as the complete set of Higgs doublets appear in the effective field theory, the PQ current is conserved (up to anomalous effects) and thus the couplings, which correspond to PQ charges, do not renormalize.Once the heavy scalar components are integrated out, the sum-rule of PQ charges set by U(1) PQ invariance breaks down, and non-vanishing contributions to the running of the couplings arise (see e.g.[14]).We can now directly match Eq. (A.1) at the scale µ = O(m BSM ) with a GKR basis featuring only one SM-like Higgs doublet where In particular, by employing global U(1) Y invariance, it is convenient to cast the RG equations in a form that does not depend explicitly on c H .This can be achieved via the axion-dependent field redefinition: ψ → ψ ′ = e −ic H β ψ a/ f ψ, with β ψ = Y ψ /Y H the ratio of the corresponding hypercharges, which redefines the effective couplings as c ′ ψ = c ψ − c H β ψ (so in particular c ′ H = 0).
to the conservation of the vector current.In our case, since the models that we have considered enjoy flavour universality, the matrices of couplings c ′ f R,L are proportional to the identity, then L f simplifies to: After including matching corrections at the weak scale [17], the running for µ < m Z is given by where Θ(x) is the Heaviside theta function, while N f c and Q f denote respectively the colour number and EM charge of the fermion f .The axion-nucleon couplings, neglecting the tiny contributions of the matrix elements ∆ t,b,c of the heavy flavours, can be calculated by using Eqs.

Appendix B. Numerical fit to RG effects
Running of axion couplings is examined in detail in Refs.[14,15,17,80], where a complete set of one-loop (and partially two-loop) anomalous dimensions are derived including matching corrections at the EW scale [17].The leading contribution to the running axion couplings arises from top loop diagrams induced by the axion-top coupling C t .The RG evolved couplings at µ = 2 GeV are thus expressed to a good approximation by where Ψ = u, d, e.Note that the running occurs below the heavy Higgs scale m BSM ≃ m H,A,H + , where in the decoupling limit the heavy scalars are assumed to be approximately degenerate, and r t Ψ (m BSM ) is a function only of m BSM .Keeping only the top Yukawa and the strong gauge couplings, the running of C Ψ below µ = m BSM is governed by [14,15,17,80]  where T 3,Ψ is the weak isospin of Ψ, a Ψ = 1 for quarks and 0 for leptons, µ w = O(m Z ) is a matching scale at which weak gauge bosons, Higgs boson and top quark are integrated out, with c G defined in Eq. (A. 16).We see from Eq. (B.2) that the RG corrections to the axion couplings consist of one-loop iso-vector contribution, proportional to the weak isospin T 3,Ψ , and a two-loop level iso-scalar contribution generated from c G , 15 and can be expressed in the form r t Ψ (m BSM ) ≃ T Note that r t 3,0 are independent of Ψ to a good precision, even after including the threshold corrections at the EW scale, which turn out to be iso-vector (numerically (r t 0 ) th /(r t 3 ) th ∼ 10 −6 ).Let us now derive approximate formulae for r t 3,0 (m BSM ).To this end, we first evaluate the running effects by numerically solving the full set of the RG equations including the threshold corrections at the EW scale [17].In the calculation the two-loop running for the SM gauge and Yukawa couplings is implemented, with their input values at µ w = m Z taken from Ref. [38] GeV) are low-energy couplings evaluated at the scale µ = 2 GeV by numerically solving the RG equations from the boundary values C u, d, s ( f a ) (see below), ∆ u, d, s denote the nucleon matrix elements of the quark axialvector currents.In particular, g A ≡ ∆ u − ∆ d = 1.2754(13) from β-decays[24], ∆ u = 0.847(18)(32), ∆ d = −0.407(16)(18)

Figure 1 :
Figure 1: f a dependence of the perturbative unitarity bounds on tan β at small tan β values.

Table 1 :
RG corrections (third column) to the DFSZ1 and DFSZ2 couplings (listed respectively in the first and second column) in the approximation of keeping only the contribution from the top Yukawa coupling Y t .The corrections are given in terms of β and of l(x) = ln( √ x − 0.52), where x = log 10 (m BSM /GeV) parameterizes the new physics scale.While the corrections to the quark couplings in DFSZ1/2 are the same, for the leptons the relative corrections differ: ∆C e /C e ≃ cot 2 β (DFSZ1) and ∆C e /C e ≃ const.(DFSZ2).

Figure 2 :
Figure 2: Running axion coupling combinations (C 3 and C 0 ) in DFSZ as a function of tan β.The red band encompasses the range of the corrections for m BSM /GeV ∈ [10 3 , 10 9 ].Perturbativity bounds on tan β also depend on m BSM , with the thick (thin) grey line corresponding to m BSM = 10 9 GeV (1 TeV).

Figure 4 :Table 3 :
Figure 4: RG effects on astrophysical axion bounds from Red Giants (red bands) and SN1987A (blue bands) for the DFSZ1 (left panel) and DFSZ2 (right panel) models, compared to the tree level results (black dashed lines).The gray line corresponds to the perturbative unitarity bound on tan β for m BSM = f a .

Figure 5 :
Figure 5: HDM bound in the DFSZ1/2 models.The red region shows the effect of RG corrections, for m BSM ranging from f a (left border) to 1 TeV (right border).The gray line corresponds to the perturbative unitarity bound on tan β for m BSM = f a .

Figure 6 :
Figure 6: Experimental sensitivity to DFSZ1 axions, including RG effects.Left Panel: IAXO and BabyIAXO.Right Panel: XENON-nT.In both plots the gray line corresponds to the perturbative unitarity bound on tan β for m BSM = f a .

Table B .
. A set of numerical values for r t 3,0 (m BSM ) are tabulated in Tab.B.4.These values are accurately fitted by the following fitting functions: × 10 −4 ln 2 (x − 1.25) , (B.7) with x = log 10 (m BSM / GeV).Eq. (B.6) agrees with the numerical results within 2% accuracy in the 1 TeV ≤ m BSM ≤ 10 18 GeV range.The precision of Eq. (B.7) is better than 6%.However, since |r t 0 /r t 3 | ≲ 0.5%, this function does not affect numerically r t Ψ .4: The numerical values of r t Ψ (m BSM ) and r t 3,0 (m BSM ), which are obtained by numerically solving the full RGEs, with the threshold corrections at the EW scale included.