Stellar Recipes for Axion Hunters

There are a number of observational hints from astrophysics which point to the existence of stellar energy losses beyond the ones accounted for by neutrino emission. These excessive energy losses may be explained by the existence of a new sub-keV mass pseudoscalar Nambu--Goldstone boson with tiny couplings to photons, electrons, and nucleons. An attractive possibility is to identify this particle with the axion -- the hypothetical pseudo Nambu--Goldstone boson predicted by the Peccei--Quinn solution to the strong CP problem. We explore this possibility in terms of a DFSZ-type axion and of a KSVZ-type axion/majoron, respectively. Both models allow a good global fit to the data, prefering an axion mass around 10 meV. We show that future axion experiments -- the fifth force experiment ARIADNE and the helioscope IAXO -- can attack the preferred mass range from the lower and higher end, respectively. An axion in this mass range can also be the main constituent of dark matter.


Introduction
There is a strong physics case for the existence of the axion [1,2]. It occurs in the arguably most elegant solution of the strong CP problem as a pseudo Nambu-Goldstone boson of a spontaneously broken global Abelian symmetry [3]. The axion is also regarded as one of the best candidates of dark matter in the universe [4][5][6].
The couplings and the mass of the axion are inversely proportional to the symmetry breaking scale. Non-observation of axion signatures in laboratory experiments have pushed the latter much beyond the electroweak scale, rendering the axion a very weakly interacting slim particle (WISP) [7,8]. Such particles could be produced in hot astrophysical plasmas, thus transporting energy out of stars and other astrophysical objects. In fact, the coupling strength of these particles with normal matter and radiation is bounded by the constraint that stellar lifetimes and energy-loss rates should not conflict with observations [9]. Such astrophysical limits have been derived from our Sun; from evolved low-mass stars, such as red giants and horizontal-branch stars in globular clusters, or white dwarfs; from neutron stars; and from the duration of the neutrino burst of the core-collapse supernova SN1987A.
Intriguingly, there exist also marginal hints for additional energy losses in stars at different evolutionary stages -red giants, supergiants, helium core burning stars, white dwarfs, and neutron stars [10][11][12]. Indeed, the current observations show in all cases a preference for an additional cooling -in most cases of the same order of magnitude as the standard one.
Throughout the years these results, globally known as the cooling anomaly problem, have led to the speculation that the apparent systematic tendency of stars to cool faster than predicted originates from the production of new WISPs with tiny couplings to photons, electrons, and nucleons [10][11][12]. Clearly, given the large systematics in stellar modeling and observations and the small statistical significance of the individual hints, any conclusion about new physics is certainly premature. However, with the current status of understanding, the WISP option cannot be ruled out and, in fact, would provide possibly the simplest solution to the cooling anomaly problem. Remarkably, according to the recent investigation in ref. [12], axions and more generally axion-like particles (ALPs) would provide very good fits for the combined hints.
Most importantly, recent years have witnessed a revival in the experimental effort in axion/ALP searches with an increasing interest in the development of a new generation of laboratory experiments. As a consequence, the parameter region hinted by the cooling anomalies is finally becoming accessible to terrestrial experiments, such as the light-shiningthrough-walls experiment ALPS II [13], the long-range force experiment ARIADNE [14], and the helioscope IAXO [15]. Unfortunately, the cooling anomalies mainly hint at an axion coupling to electrons, while the mentioned experiments are critically sensitive to the axionphoton and the axion-neutron coupling 1 . In a concrete axion model, these couplings are related by a few constants, and thus we can study the detectability and complementarity of the experiments and obtain additional constraints in a model-dependent basis. The aim of this paper is to select a sample of generic axion models and study whether they can accommodate all or a part of the stellar cooling anomalies, providing the best fits of their parameters to guide their experimental discovery or refutation.
The paper is organized as follows. In section 2 we summarize the observations of cooling anomalies and their interpretation in terms of model-dependent axion/ALP couplings to electrons, photon, neutrons, and protons. In section 3 we confront the required couplings with the predictions of axion models and compare those with the capabilities of upcoming experimental searches. In section 4 we show that the axion hinted at by the anomalous stellar cooling may also constitute the dark matter in the universe. We summarize our results in section 5. We leave the discussion of some technical details to the appendixes. In appendix A, we describe how we computed the χ 2 for the global fits. We present the state-of-the-art interpretation of the supernova (SN) bound on the nucleon couplings in appendix B. The sensitivities of IAXO and ARIADNE are discussed in details in appendix C.

Axion/ALP interpretation of stellar cooling anomaly observations
Axions, and more generally ALPs, are pseudo Nambu-Goldstone bosons arising from the spontaneous breaking of an Abelian global symmetry at a scale much larger than the electroweak scale. At very low energies, their interactions with photons (γ), electrons (e), protons (p), and neutrons (n) can be described by the Lagrangian where f = e, p, n; F denotes the electromagnetic field strength tensor andF its dual. The model dependence rests in the mass m a , the dimensionless coefficients C ai , and the decay constant f a , which is proportional to the symmetry breaking scale. Hints for excessive cooling have been observed in: 1) several pulsating white dwarfs (WDs), whose cooling efficiency was extracted from the rate of the period change; 2) the WD luminosity function (WDLF), which describes the WD distribution as a function of their brightness; 3) red giants branch (RGB) stars in globular clusters, in particular the luminosity of the tip of the branch; 4) horizontal branch (HB) stars in globular clusters or, more precisely, the R-parameter, that is the ratio of the number of HB over RGB stars; 5) helium burning supergiants in open clusters, more specifically the ratio B/R of blue and red supergiants; and 6) neutron stars (NS). Table 1 lists the results and their references.
The most relevant axion/ALP production mechanism for stars with a high density core, such as WDs and RGB stars, is bremsstrahlung off electrons, which induces an additional energy loss rate proportional to g 2 ae T 4 , where As shown in ref. [12], this temperature dependence is optimal to fit the WDLF and provides a reasonably good explanation for the observed excess cooling in DA and DB WD variables, whose internal temperatures differ by a factor of a few. The combined analysis of all the observed WD variables gives a fairly good fit, χ 2 min /d.o.f.= 1.1, for g ae = 2.9 × 10 −13 and favours the axion/ALP interpretation at 2 σ. Moreover, the peculiar temperature dependence of the axion bremsstrahlung rate allows to account for the excessive cooling observed in RGB stars [26], which have a considerably larger internal temperature than WDs, with a comparable axion-electron coupling. For this paper we have reexamined the WDLF (see appendix A) and combined it with the WD pulsation, and RGB stars. Our fit gives with χ 2 min /d.o.f.= 14.9/15 = 1.0, and favours the axion/ALP solution at slightly more than 3 σ.
Axion bremsstrahlung off electrons contributes also to the evolution of low mass stars [12]. Its main effect is to delay helium ignition at the end of the RGB phase while leaving essentially unchanged the HB stage, thereby reducing the expected ratio of HB and RGB stars (the R-parameter). The most recent study of a set of globular clusters outputs a strong constraint, but also a mild preference for additional cooling [28,29].
The reduction of R could also be due to other axion emission channels. Most notably, the Primakoff process, consisting in the conversion of a photon into an axion/ALP in the electric field of nuclei and electrons, is the most discussed one in relation to the hint. This process depends strongly on the environment temperature and is suppressed at high density (e.g., those characterizing the core of WDs and RGB stars) by the plasma frequency and degeneracy effects (see, e.g., [36]). However, it could efficiently accelerate the cooling of the HB stars and consequently lower the observed R-parameter.
In figure 1, we display the 1, 2, and 3 σ hinted areas for the combined fit of the WD, RGB, and HB cooling anomalies 2 in the vs. g ae plane. Note, that the quantitative analysis with the data at hand shows that it is possible to explain all the observed cooling hints, including the R-parameter anomaly, even neglecting the axion-photon coupling, though there is a preference for non-vanishing couplings with both electrons and photons. The best fit value, g ae = 1.5 × 10 −13 and g aγ = 0.14 × 10 −10 GeV −1 , is indicated in the figure with a red dot. Taken at face value, it implies |C ae /C aγ | = 0.025, which is a somewhat unusual ratio because typical axion models have . We will explore axion models of both these typical types and consider also axion-majoron models where the above unusual ratio can be obtained. Note, however, that the 1 σ contour is compatible with a zero axion-photon coupling. An additional hint comes from the anomalously fast decline of the surface temperature of the neutron star in the supernova remnant Cassiopeia A (CAS A), which can be interpreted as due to axion emission in neutron 3 P 2 (m j = 0) Cooper pair formation with a neutron coupling [34] g 2 an ∼ (1.4 ± 0.5) × 10 −19 . (2.7) The estimate of [34] does not come with an error bar, though the study claims that values outside the quoted range cannot fit the data. We shall warn that the systematics in this calculation are very hard to estimate. Indeed, the cooling of the superfluid core in the neutron star may be also explained by other means, for instance, by neutrino emission in pair formation in a multicomponent superfluid state 3 P 2 (m j = 0, ±1, ±2) [37]. Thus, at the moment, we cannot take seriously this NS hint, and we can only ask the question whether the hinted parameters can be compatible with the other hinted regions in particular models. To this end, we will confront axion models with stellar cooling hints with and without including the NS hint. For these purposes it will suffice to interpret the quoted range as a 1 σ interval. Note that the range quoted is compatible with the upper limit g an < 8 × 10 −10 from NS cooling in refs. [38,39]. This, however, is more complicated to interpret in terms of a  Figure 1. Combined analysis of the hints from WD+RGB+HB stars in the g ae − g aγ plane. Also shown are the projected sensitivities of the light-shining-through-walls experiment ALPS II [13] and the helioscope IAXO [15]. combination of couplings. In fact, it includes also axion emission in nucleon bremsstrahlung N +N → N +N +a where N can be either a proton or a neutron and most of the simulations were done with a very small neutron coupling, so that the effect is mostly due to the proton coupling.
Finally, the axion/ALP bremsstrahlung off nucleons can shorten the prediction of the neutrino pulse duration of core collapse supernovae. In fact, the neutrino observations from SN1987A lead to a bound (see appendix B) g 2 ap + g 2 an < 3.6 × 10 −19 . (2.8) We will consider this as a 1 σ hint that g 2 aN = 0 within the error 3.6 × 10 −19 . However, we warn the reader that SN 1987A constraints are based on axion emissivities not completely understood and on simulations that at the moment do not include all necessary physics and therefore have systematic uncertainties themselves. Again, we will study axion models with and without including this constraint.

Axion interpretation of stellar cooling anomalies
Let us start by reviewing the generic features of the axion. The basic building block of an invisible axion model is a global U (1) PQ symmetry, which is broken at a high scale by the vacuum expectation value σ 2 = v 2 PQ /2 of a complex Standard Model (SM) singlet scalar field σ. In this notation, the axion field appears as the phase of this complex scalar σ = (v PQ / √ 2)e ia/v PQ or as a linear combination of this and other Higgs phases. The associated Noether current J PQ µ is required to have a color anomaly and, although not required for solving the strong CP problem, it may also have an electromagnetic anomaly: where G a µν is the color field strength tensor,G aµν its dual, while N and E are anomaly coefficients. The decay constant of the associated Nambu-Goldstone boson is then given by The axion mass m a , in units of the decay constant f a , is equal to the square root of the topological susceptibility χ in QCD. Recent precision calculations of the latter by next-toleading order chiral perturbation theory [40] or lattice simulations [41] result in Moreover, the photon coupling of the axion is given by while the proton and neutron have a model-independent part and a model dependent contribution that arises from possible axion-quark couplings of the form (C aq /2)(∂ µ a/f a )ψ q γ µ γ 5 ψ q in the high-energy theory, as found in the state-of-the-art calculation [40]. Note that the couplings g ai ∝ 1/f a of the axion to i = γ, e, p, n are effectively proportional to its mass, g ai ∝ m a . That is why, in specific axion models, the coupling strength is often parameterized with the mass 3 . ALPs generically do not feature this particular relation between their couplings and their mass.

DFSZ axion
In the DFSZ axion model [42,43], the Standard Model (SM) Higgs sector is extended to contain two Higgs doublets, H u and H d , whose vacuum expectation values v u and v d give masses to up-type and down-type quarks, respectively. There are two possibilities, dubbed DFSZ I or DFSZ II, according to whether leptons couple to H d , which occurs in familiar Grand Unified Theories (GUT), or to H u . In the usual nomenclature [44], the Higgs sector of DFSZ I is a type-II Two Higgs Doublet Model (2HDM), which is easily embedded in a GUT, while the one of DFSZ II is a type-IV, sometimes also called Flipped 2HDM. Correspondingly, the Yukawa interaction terms in DFSZ I read while, in DFSZ II,  Figure 2. 1, 2, 3, 4 σ contours in the global fit of DFSZ I (left panels) and II (right panels) with WD+RGB+HB data alone (top) and including the SN 1987A constraint (middle) and, in addition, the NS CAS A data (bottom). Perturbative unitarity in the Yukawa couplings to fermions is satisfied for values of tan β between the two dashed lines. The red bullets are the best fits compatible with perturbative unitarity. Also shown are the projected sensitivities of ARIADNE [14] and IAXO [15].
Here H u = H * u , i, j = 1, 2, 3 are flavor indices and Γ ij , Y ij , G ij are complex 3 × 3 matrices. The interactions given by eq. (3.6) (DFSZ I) or eq. (3.7) (DFSZ II) are assumed to be invariant under a U(1) PQ symmetry with symmetry breaking scale v PQ . At low energies, the effective Lagrangian is then given by eq. (2.1), with [45,46] f a = v PQ 6 , (3.8) (3.10) and [40]  Here, GeV. It is theoretically constrained from both ends by the requirement of perturbative unitarity of the Yukawa couplings, 0.28 < tan β < 140 . (3.13) Here, the lower limit arises in all 2HDMs, while the upper limit is specific to the type-II and type-IV 2HDMs [47]. The DFSZ models have only two parameters, f a and tan β, that we can extract from the global fit of the WDLF, the period decrease of 4 pulsating WDs (R548, L 19-2 (113), L 19-2 (192), and PG 1351+489), the luminosity of the tip in the RGB of M5 and the R-parameter in globular clusters, which we hereafter label as HB, see appendix A for specifics. The best fit values are recorded in table 2 and the 1, 2, 3, 4 σ contours are shown in figure 2. Note that we impose the constraint on perturbative unitarity on the best fit values but not on the contours. The resulting regions can be understood as follows.

Model
Global  Let us consider first the DFSZ I case. In figure 2 (up, left) we show the 1, 2, 3, 4 σ contours for the hints involving only the coupling to electrons and photons (WD+RGB+HB stars). The main driver of the fit is the WDLF, forcing a concrete value for g ae = m e sin 2 β/3f a ∼ 1.5 × 10 −13 which shows as a vertical line for large tan β (sin 2 β ∼ 1) with a diagonal turn towards the bottom right when the small tan β is compensated by a decrease of f a . The degeneracy is broken at f a ∼ 3 × 10 7 GeV because of the g aγ dependency of the R-parameter.
A small value of f a makes g aγ large and gives too small a value for the R-parameter. The fact that the 1 σ contour shows both the vertical and diagonal branches means that we can have large tan β and small tan β DFSZ I axions fitting the data.
However, adding the bound from SN 1987A forces the axion decay constant to be large, disfavouring the small tan β solution, tan β > 0.7 at 1 σ (the best fit coincides with the unitarity constraint tan β ∼ 140), and the regions become almost vertical, see 2 (medium, left). The WD+RGB+HB set is in tension with the SN constraint and that worsens a bit the χ 2 min /d.o.f.. In this large tan β region the axion-neutron coupling required to fit the NS hint from CAS A is at odds with the SN bound (C ap is comparable but slightly larger than C an ) so the fit worsens by including it, see table 2 and figure 2 (bottom, left). However, there is a second local minimum of χ 2 with χ 2 min /d.o.f.= 24.8/17 ∼ 1.4 at tan β ∼ 0.5 and m a ∼ 20 meV. At small tan β the SN constraint becomes weaker compared to the NS hinted region because the coupling to protons becomes of the same order as the coupling to neutrons. Of course, this solution makes sense only if we trust the NS hint, and somehow disregard or diminish the importance of the SN constraint.
Let us now discuss the DFSZ II case. The favoured region by WD+RGB+HB is still dominating the fit to give g ae ∼ 1.5×10 −13 . But now g ae = m e (1−sin 2 β)/3f a , so small tan β corresponds to a vertical line at the largest f a and large tan β can be compensated by a small f a , corresponding to the up-right diagonal branch in figure 2 (up, right). The SN constraint enforces again large f a so that now the small tan β solution is preferred, giving a better fit than in DFSZ I because C ap , C an are relatively smaller at small tan β, see table 2. As in DFSZ I, the SN constraint cancels the hint for cooling in HBs due to the axion-photon coupling and thus the best fit region is the one further away from the SN constraint, so the best fit turns out to be in the unitarity constraint tan β ∼ 0.28, although the 1σ region extends to tan β < 1. This fit is incompatible with the NS hint at more than 2σ, so the fit degrades by including it, see table 2. The reason is not the tension with SN like in DFSZ I, but with the WD+RGB+HB set. In DFSZ I with large tan β we have C ae /C an = (1/3)/(0.414 − 0.16) = 1.3 while in DFSZ II at small tan β the ratio is C ae /C an = (1/3)/(−0.16) = −2. Thus, for the same electron coupling the neutron coupling is a factor ∼ 2 weaker, which adds up to the already not good NS fit of DSFZ I at large tan β.
On a pessimistic note, f a = ∞, which corresponds to axion decoupling, is still allowed at ∼ 3 − 4 σ, which is not very significant when we take into account that there might be unidentified systematics in the cooling analysis. However, at 3 σ the best fit regions indicate always a non-zero values of f a and some axion cooling at work. We have shown the best fit as red points in figure 2, but we warn once more that they have little meaning. The best guidance is provided by the 1 σ region, which in each case points to f a ∼ 10 9 GeV (corresponding to a mass m a ∼ 6 × 10 −3 eV). The hinted values of tan β are 0.7 in DFSZ I and 1 in DFSZ II. Both have a large overlap with the region where perturbative unitarity is granted.
Finally, we see from figure 2 that the 1 σ parameter region preferred by the DFSZ axion interpretation of the stellar cooling anomalies can be entirely probed by IAXO+ (see appendix C), while the projected sensitivity of ARIADNE barely touches the 1(2) σ region for DFSZ I(II).

Pure KSVZ-type models
In pure KSVZ-type axion models [48,49], the color anomaly arises from beyond the SM coloured vector-like fermions Q = (Q L , Q R ) transforming chirally under U (1) PQ and according to representations The corresponding N and E/N values for specific representations R Q [50] are given in the second and third column of table 3, respectively.  Table 3. Best fit parameters for KSVZ-type axion interpretations of the HB, RGB, and WD cooling anomalies.
The dimensionless coupling C aγ to photons has already been shown in eq. (3.4). At tree level, quarks and leptons do not interact with the KSVZ axion. Correspondingly, the leading contributions to the dimensionless couplings to protons and neutrons arise solely from the mixing with the pion and read [40]: (3.14) Moreover, the dimensionless coupling to electrons is induced at one-loop level by the photon coupling and estimated as [45] C KSVZ ae 3α 2 2π where Λ ∼ 1 GeV. Since E/N is fixed for each KSVZ-type axion model, the only free parameter to determine the couplings g 2 The preferred values of f a and m a from a global fit of WD+RGB+HB are given in table 3.
The commonly shown value of E/N = 0 gives a rather bad fit, with χ 2 min /d.o.f. 1.6. A preference for E/N 2 is clearly seen. In fact, if E/N were used as a fitting parameter we would find the best value to be E/N = 1.95. The model exploiting the representation R Q = (3, 2, + 1 6 ) has E/N = 5/3 and thus comes closest to this value, but still with a relatively bad The reason for this is that in KSVZ-type models the axion coupling to electrons emerges only at the loop level and is strongly suppressed with respect to tree level couplings, cf. eq. (3.15). Though the cooling anomalies indicate a preference for C ae about two orders of magnitude smaller than |C aγ |, for the KSVZ axion (with E/N =0) the prediction for C ae is about an order of magnitude too small with respect to the corresponding |C aγ |. Therefore, the fit prefers a value of E/N close to 1.92 (since that lowers |C aγ |, leaving C ae unaffected), and increases both the couplings to photons and electrons through lower values of f a , corresponding to larger values of m a . As a consequence, the models fitting better are those with E/N ∼ 1.95 and a large mass ∼ eV.
An axion in this mass range is, however, severely constrained if not totally excluded by several astrophysical and cosmological arguments. The reader might remember that when the axion coupling to nucleons is very large, axions can interact so strongly with the medium as to remain trapped in the SN core, reducing the axion emissivity to a relatively thin and cold "axiosphere" that cannot compete with neutrino emission and thus does not affect the duration of the neutrino pulse of SN 1987A. The excluded region is nominally taken to be 6 × 10 −10 < g ap < 6 × 10 −8 , which we adapt from [51]. On the strongly interacting side, the fewer axions emitted can interact strongly enough to be detected in Kamiokande by the nuclear reaction a 16 [52]. This would exclude the range 10 −6 < g ap < 10 −3 [9,51]. The gap 6 × 10 −8 < g ap < 10 −6 (0.78 eV < m a < 13 eV, corresponding to 0.44 × 10 6 GeV < f a < 0.73 × 10 7 GeV) can be excluded by an independent constraint from cosmology. In fact axions with these masses are produced thermally in the early universe and constitute a sizeable hot dark matter component, in analogy to massive neutrinos. Cosmological precision data provide an upper limit on the possible hot dark matter fraction which translates into an upper limit on the axion mass, m a < 0.8 eV at 95%C.L. [53][54][55], which in principle closes the gap. 4 The search for monochromatic photon lines from axion dark matter decay excludes even more strongly the region 4.5 eV < m a < 7.7 eV [56], see also [57]. Once more we warn the reader that SN 1987A constraints have uncontrolled systematic uncertainties themselves. The numbers presented here have to be taken with a grain of salt and as an invitation to review the situation once more precise modeling is available. Note finally that we have not considered the hint from the NS in CAS A, since the interaction of the KSVZ axion with the neutron is compatible with zero.
We conclude that a KSVZ-type axion can not give a plausible explanation of the cooling anomalies without severely compromising other astrophysical constraints.

KSVZ-type axion/majoron
In refs. [46,[58][59][60], the KSVZ-like particle content is extended by three right-handed SMsinglet ("sterile") neutrinos N Ri , and the Peccei-Quinn (PQ) symmetry is unified with the lepton symmetry, explaining the smallness of the masses of the active neutrinos by the seesaw mechanism [61][62][63][64]. The model then features Yukawa couplings to the left-handed SM lepton doublets L, the SM Higgs doublet H, and the PQ complex scalar σ as follows: where V is the 6×6 mixing matrix to the states n R , which form the Majorana mass eigenstates n = n R + n c R . In these models, the axion A is at the same time the majoron J -the pseudo Nambu-Goldstone boson arising from spontaneous breaking of the lepton symmetry [65]. Importantly, the one-loop induced axion-electron coupling, gets an extra contribution from the loop involving the sterile neutrinos N Ri [58,66], whichto lowest order in the seesaw limit, m D /M M 1 -is given by [67] where the dimensionless hermitian 3 × 3 matrix κ is defined as We note that all the diagonal entries of κ are real and non-negative, and thus its trace is non-negative [67], In addition, the diagonal entries are bounded from above by perturbative unitarity, so 0 ≤ κ 4π . The fit is once more dominated by the WD+RGB that force g ae ∼ 1.5 × 10 −13 , imposing a degeneracy between |trκ − 2κ ee | and f a . This results in the diagonal shape of the contours. The degeneracy is broken at small f a when g aγ becomes large and starts to ruin the R-parameter fit. In the R Q = (3, 2, + 1 6 ) case, the C aγ is smaller, so the 1 σ best fit region extends further in m a , almost reaching m a ∼1 eV. The best fits prefer sizeable diagonal matrix elements in κ, of the order of the perturbative unitarity bound |trκ − 2κ ee | 1. At the smallest values of |trκ − 2κ ee | ∼ 0.1E the photonloop contribution to C ae produces a visible feature because it can cancel completely the RH neutrino contribution.
The KSVZ A/J may also have a larger coupling to the neutron than the pure KSVZ axion, due to the loop induced contribution to the axion-quark couplings from the sterile neutrinos,       3 ); right panels for R Q = (3, 2, + 1 6 )) with WD+RGB+HB data alone (top) and including the SN 1987A constraint (middle) and, in addition, the NS CAS A data (bottom). Perturbative unitarity for the Dirac Yukawas F is satisfied for values of trκ below the dashed lines. The red bullets are the best fits compatible with perturbative unitarity. Also shown are the projected sensitivities of ARIADNE [14] and IAXO [15] (see appendix C).
This contribution leads to a non-zero axion coupling to the neutron, C an , and thus may explain the NS cooling hint. For completeness we include the corrections to the axion proton coupling to consider the SN 1987A constraint as well. We find, The inclusion of the SN 1987A constraint in the global fit has dramatic consequences (see middle panels in figure 3). In this case, in general, three A/J parameters have to be fitted: f a , |trκ − 2κ ee |, and trκ. It turns out that our minima of χ 2 happen always at κ ee = 0, so that using κ ee = 0 is equivalent to marginalising over it. Moreover, for hierarchical values of the diagonal elements of κ, axion couplings effectively depend only on trκ. In fact, for κ µµ , κ τ τ κ ee and κ ee κ µµ , κ τ τ , we have |trκ − 2κ ee | ≈ trκ. Note that the unitarity limit of trκ is 8π and 4π respectively in these limits. For simplicity we will assume |trκ−2κ ee | ∼ |trκ| in the following and display the most conservative 8π unitarity limit.
Axion masses above m a ∼ 10 meV become excluded by the SN argument and so the 1 σ region recedes towards the region inaccessible to reliable perturbative calculations. All the best fits have smaller χ 2 min /d.o.f. We show the values for the best case R Q = (3, 1, − 1 3 ) in table 5. The worst case corresponds to R Q = (3, 2, + 1 6 ), because for N = 2 the maximum electron coupling compatible with perturbative unitarity is half the value than in the N = 1 cases. We see however that the favoured regions correspond to the highest masses m a ∼ 10 − 20 meV and that there is a sizeable region of parameters within the 1 σ region, although one has to remember that the best fit is not so good.
Global  Adding the hint from the NS in CAS A actually improves the fit, see table 5, so that the preferred regions shrink a bit, cf. figure 3 (bottom panel). This is due to a remarkable prediction of A/J models. Neglecting the −0.02(3) factor in C an we have while the ratio of the neutron and electron couplings favoured by the CAS A NS and the combination WD+RGB+HB, respectively, is C an C ae = g an /m n g ae /m e = 3.8 × 10 −10 1.5 × 10 −13 m e m n = 1.4 (3.26) Therefore, A/J models that fit the combination WD+RGB+HB tend to give a good fit to NS unless 2κ ee ∼ κ µµ + κ τ τ , the only case where the ratio trκ |trκ−2κee| can differ significantly from 1.
For R Q = (3, 1, − 1 3 ) and R Q = (3, 2, + 1 6 ), the 1 σ parameter region preferred by the KSVZ A/J interpretation of the stellar cooling anomalies can be entirely probed by IAXO+, while the case R Q = (3, 2, + 1 6 ) lies outside of sensitivity of IAXO. ARIADNE has sensitivity only in the region outside of perturbative control, cf. figure 3.

Axion hints compatible with axion dark matter?
In addition to providing a possible interpretation of the stellar cooling hints, the axion can also be a good candidate of dark matter. Although the mass range obtained by global fits to the stellar cooling hints is much higher than the conventional one [4][5][6], there is a possibility that the axion becomes the main constituent of dark matter at those large masses. This possibility arises from the fact that the relic axion abundance strongly depends on the early history of the universe.
We consider two possibilities in the context of PQ symmetry breaking in inflationary cosmology. One possibility is to assume that the U (1) PQ symmetry has been broken before inflation and is never restored. In this scenario, the present observable universe is predicted to be in the same vacuum due to inflation homogenising the universe up to scales beyond the horizon today. We can easily assume that the few domain walls present in the universe are not in our local observable patch. In this case, axions are produced by the standard vacuum re-alignment mechanism [4][5][6] with the initial axion angle θ i = a /f a fixed at the inflationary epoch. The prediction of the relic axion abundance depends not only on f a , but also on θ i . In particular, the slow-roll behavior of the axion field in the limit θ i → π opens up a possibility that the axion becomes the main constituent of dark matter in the m a < meV range [68]. A potential problem with this scenario is that the axion dark matter imprints isocurvature fluctuations in the cosmic microwave background temperature anisotropies and these are enhanced in the θ i → π limit. The null observation of the isocurvature fluctuations by the Planck mission leads to a severe upper bound on the Hubble expansion rate during inflation, H inf O(100) GeV for f a = 10 10 GeV, and it becomes even tighter for θ i → π [69]. Therefore, it is unlikely that the axion whose mass fits to the stellar cooling hints becomes the main constituent of dark matter in the non-restored symmetry scenario.
Next, we consider the alternative scenario in which the U (1) PQ symmetry is restored after inflation. In this case, some class of axion models suffer from a serious cosmological domain wall problem [70]. In general, the low energy effective potential for the axion field has N degenerate minima connected by N maxima that lead to N different type of domain wall solutions. Hence, N is sometimes called domain wall number and denoted as N DW . In the DFSZ models we have N DW = 6, while in KSVZ-type axion models, the value of N DW depends on the representation R Q and the number N Q of exotic quarks introduced in the theory. For instance, N DW = N Q if we take R Q = (3, 1, − 1 3 ) or R Q = (3, 1, + 2 3 ) and assume that PQ charges of Q's have the same sign, while N DW = 2N Q if we take R Q = (3, 2, + 1 6 ) and assume that PQ charges of Q's have the same sign. If the different vacua are populated in the early universe, domain walls are created as boundaries between them. For N DW > 1, the network of domain walls is stable and its energy density redshifts slower than that of radiation [71], and hence it tends to dominate the total energy density of the universe, causing trouble with the standard cosmology [72]. There is no domain wall problem if N DW = 1, because the network is unstable in this case [73].
A domain wall number equal to one is realized for example in KSVZ-type models with N Q = 1 representation of R Q = (3, 1, − 1 3 ) or R Q = (3, 1, + 2 3 ). However, in this case, the required axion mass to explain the observed dark matter abundance is predicted to be in the range 50 µeV m a 200 µeV [60,74,75] -way too small to explain the stellar cooling anomalies 5 .
Therefore, we have to consider the alternative scenario with domain wall number larger than one, which is automatic in DFSZ-type models and can be arranged in the simplest KSVZ-type models by taking N Q > 1. In this case, the domain wall problem can be avoided if new physics breaks the degeneracy between the different vacua [78,79]. Such energy difference between vacua creates pressure that renders the domain walls unstable. Here we assume that the U (1) PQ symmetry is explicitly broken due to the following Planck-suppressed operator [80] L ⊃ gM 4 where σ is the SM singlet complex scalar field, M Pl 2.435 × 10 18 GeV is the reduced Planck mass, N is an integer 6 , and g is a complex dimensionless constant. In the low energy effective Lagrangian, the above operator induces an additional potential for the axion field which lifts N DW degenerate vacua and induces late-time annihilation of domain walls. Here, ∆ includes the phase of coupling g, andθ is the sum of the QCD θ parameter and the contribution from the quark mass phases (we have redefined the axion field to have a = 0 at the CP conserving minima). The collapse of long-lived domain walls produces an additional number of axions, that contribute to and can be the dominant component of the dark matter density [75,81]. In this case, the relic axion abundance depends not only on f a , but also on the explicit symmetry breaking parameter g and on the order of the Planck-suppressed operator N . Here we follow 5 It is argued that the theoretical uncertainty of the mass of axion dark matter in the scenario with NDW = 1 might be even larger because of an effect overlooked in the previous literature [76]. A recent alternative analysis tentatively finds even smaller mass values, by a factor of order one [77]. 6 In ref. [80], the PQ symmetry is assumed to arise as an accidental symmetry from a discrete ZN symmetry.
the analysis in refs. [75,80] and find that the hinted mass ranges obtained in the previous section can be compatible with axion dark matter if N = 9. We note that, for N ≤ 8, the axion mass is determined by V g in eq. (4.2) rather than the potential induced by the topological susceptibility in QCD unless |g| takes an extremely small value [80]. Furthermore, we find that the observed dark matter abundance cannot be explained with |g| < 1 if N ≥ 10.
In table 6, we summarise the required values of |g| that allow fitting the observed cold dark matter abundance for benchmark values of the axion mass obtained in the previous section. The benchmarks are the best fit values compatible with perturbative unitarity. For KSVZ A/J models, we consider three further choices of the domain wall number, N DW = 2, 4, and 6 with the representation R Q = (3, 1, − 1 3 ) by simply adding N DW generations in the same representation. These models are always at the border of the perturbative unitarity region, and have increasing χ 2 min /d.o.f. as we increase N DW . The tension with unitarity and the SN constraint becomes so severe at large N DW that the best fits move towards g ae < 1.5 × 10 −13 worsening the fit to WD+RGB+HB, because the likelihood decreases slower in that direction than in the small f a direction, where it faces the SN constraint.
Note that for all the models the axion decay constant takes a larger value when we include the SN bound, which increases the energy bias in our parametrisation. Therefore, a smaller value of |g| is required in order to compensate it.   Table 6. Values of the symmetry breaking parameters |g| and ∆ (for the case N = 9) that explain the present dark matter abundance for different benchmark axion models explaining the stellar cooling anomalies. The range of values corresponds to the range of uncertainties in numerical simulations of topological defects considered in [75].
The constraints on ∆ are shown in table 6. 7 The resulting tuning in the phase ∆ is relatively small compared with the one required to solve the strong CP problem. The dependence of the upper limit on ∆ on the domain wall number can be understood as follows. Basically the nEDM constraint can be avoided if f a is sufficiently small (i.e. the explicit symmetry breaking term is suppressed), except for the cases with ∆ → π/2 for N DW = 2, 6 and that with ∆ → 3π/4 for N DW = 4, where two vacua are degenerate and domain walls never collapse. In those cases the dark matter abundance blows up and a larger value of |g| is required in order to compensate this effect. This large |g| enhances the magnitude of the explicit symmetry breaking term, which violates the nEDM constraint. Therefore, the upper limit on ∆ for benchmarks with HB, RGB, WD appears around ∆ π/2 for N DW = 2, 6 and ∆ 3π/4 for N DW = 4. On the other hand, if f a is sufficiently large, the nEDM constraint is relevant and some tuning of ∆ is required. The constraint is determined from the parameter dependence is compensated by a small factor of |g|, which results in the mild N DW dependence for benchmarks with HB, RGB, WD, NS, SN. The axion features CP violating couplings proportional to θ min , like a Yukawa coupling to nucleons that will be exploited by the ARIADNE experiment to search for axion-mediated long-range forces. Of course, values of ∆ close to the limits quoted in table 6 correspond to the maximum CP violation induced by our Planck-suppressed operator and thus the largest and most optimistic effects in ARIADNE, see appendix C.

Discussion and Conclusions
We have explored quantitatively the possibility that the hints of excessive energy losses of stars in various stages of their evolution (red giants, helium burning stars, white dwarfs, neutron stars) can be explained by the axion -the pseudo Nambu-Goldstone boson predicted by the Peccei-Quinn solution of the strong CP problem. The main finding of our work is that the overall > 3 σ tension of all RG+HBS+WD data taken collectively (and presumably higher if including the NS hint) is successfully removed by taking into account the presence of axions. More specifically, good fits to the data have been obtained in two well-motivated classes of axion models: DFSZ-type models and KSVZ-type axion/majoron models 8 . These fits generically prefer an axion mass around 10 meV, if the energy loss constraint from the neutrino pulse duration of SN 1987A is taken into account, cf. figures 2 and 3, and tables 2 and 5. The predicted regions, compatible with perturbative unitarity at 2 σ, are given in table 7 and figure 4.
Both DFSZ models allow a good fit to the data but DFSZ I has larger proton couplings and thus more tension with the SN 1987A constraints. None of them gives a good fit to the NS cooling hint so we have not considered it in the 2 σ values given in table 7. The photon coupling of DFSZ II is larger, making it more suitable for detection in IAXO, cf. figure 4. The mass range in DFSZ II is narrower. Axion/majoron models also give a good fit to the data, including this time the NS cooling hint that we have considered in the values given in the table. Masses tend to be larger than in DFSZ and photon couplings as well, with exception of (3, 2, + 1 6 ) whose small value of C aγ makes it quite unfavourable to detect. The  Table 7. 2 σ ranges in masses and couplings to electrons, photons, neutrons, and protons, favored by the axion interpretation of the stellar cooling anomalies in WDs, RGBs and HBs, taking into account the constraints from SN 1987A and perturbative unitarity. In the A/J models, also the hint on anomalous cooling from the NS in CAS A is included in the fit.
reason is that we have excluded values above the perturbative limit of trκ, which cuts the large values of f a . The regions are all very close to the perturbative unitarity limit, so they must be taken with a grain of salt. For these models to give dark matter one needs to include copies of the heavy quarks. In general, this narrows the predicted ranges due to the tension with the perturbativity constraint.
Importantly, the preferred regions in the parameter space can be probed by experiments. In particular, the nominal IAXO projection will already be sufficient to probe some of the region of interest while the upgraded IAXO scenario (IAXO+) would decisively test both the DFSZ and the KSVZ A/J interpretation of the stellar energy loss anomalies, as evident from figures 2, 3, and 4. In light of this result, the experimental effort to reach the IAXO+ parameters appears highly motivated. ARIADNE, on the other hand, shows a capability to access a region of parameter space highly complementary to that of IAXO, which however barely touches the region of interest for DFSZ and does not reach it for KSVZ A/J.
The future ultimate WIMP dark matter experiment DARWIN has a projected sensitivity of g ae 10 −12 to solar axions [17] and therefore misses the required sensitivity to probe the parameter region of interest from stellar cooling by about one order of magnitude. On the other hand, neutrino detectors such as IceCube, Super-Kamiokande or a future mega-ton water cerenkov detector will probe exactly the parameter region of interest by measuring the neutrino pulse duration of the next galactic SN [92]. Another possibility to test these axions observationally could be spectral signatures of photon-axion oscillations from celestial objects [93]. Moreover, for the DFSZ models, the preferred region of tan β can be complementarily probed by experiments at the intensity and energy frontier, by searching for direct production of the additional Higgs bosons or for non-SM contributions in flavor physics, cf. e.g. refs. [94][95][96].
If one ignores the SN 1987A constraint on the axion-nucleon couplings, the best fit regions move to slightly higher masses. In these cases (first rows of figures 2 and 3) the relevant region (with the exception of mentioned R Q = (3, 2, + 1 6 ) models) will be almost entirely probed by IAXO, notably thanks to the buffer gas phase (see appendix C).
An axion in this mass range can also be the main constituent of dark matter, in particular if the Peccei-Quinn symmetry is restored after inflation. Assuming that the PQ symmetry is broken explicitely by a Planck-suppressed operator L ⊃ gM 4 Pl (σ/M Pl ) 9 , the observed dark matter abundance can be obtained for both DFSZ axions as well as KSVZ axion/majorons  with a mass in the range to explain the stellar cooling anomalies 9 , requiring relatively small tuning in g, cf. table 6. However, its direct detection will be a challenge. In this connection it is of interest that in the mass range of interest, IAXO may indeed measure the mass of the axion [98]. This information could help to design a dedicated direct detection axion dark matter experiment, for example based on the idea of a dish antenna [99], exploiting the axion-photon coupling, such as BRASS. Alternatively, direct axion dark matter experiments employing the electron coupling, e.g. searches for induced atomic transitions [100,101], searches for electron spin precession in the galactic axion dark matter wind, such as QUAX [102], and searches for axion dark matter absorption on a conduction electron, followed by emission of an athermal phonon, in a superconductor [103], may also ultimately reach the required sensitivity.
In this section we describe the way we have built our χ 2 model. For the WD luminosity function, we take the binned LF data points from [23], used in [24], and digitised the synthetic WDLFs as a function of the g ae coupling (the authors present them as function of m a with an implicit relation g ae = 0.28 × 10 −13 m a [meV]). We choose 11 points in the range 7 < M bol < 12.25. Higher luminosity points have large error bars, do not help the constraint and decrease unnecessarily the χ 2 min /d.o.f. We take the luminosity as a free-parameter, marginalising over it by finding the smallest χ 2 as a function of g ae . As in [24], we find a small positive correlation between the normalisation factor and g ae . However, it cannot be used to lower the significance of the hint below 3 σ. The normalisation is typically a few percent.
We consider 4 WD variables [19][20][21][22]. We exclude G117-B15A [18] from the fit because it is very similar to R548 (the stars are almost identical, the pulsating mode considered is the same and in both cases there is an issue with the assumption of trapping), which has much more conservative errors. For each of them we take the measured period decrease with its quoted 1 σ error and fit the provided models (again as a function of g ae ) with their reported 1 σ errors. We add each hint as a term in the χ 2 with the 1 σ errors in quadrature.
We follow the same procedure for the tip of the RGB of the globular cluster M5 from [25][26][27]. In the case of the R-parameter, whose observational value is R = 1.39 ± 0.03, we employ the theoretical model derived in [12] as a function of g aγ and g ae . We take the primordial Helium abundance Y = 0.255 ± 0.002 [104] -which propagates an additional uncertainty in R, σ Y = 0.015 -and the data from the analysis of [28,29].
The χ 2 function is then For the SN constraint and NS hint we add

B Bound on nucleon couplings from SN 1987A
The observed neutrino signal from SN 1987A has been used to set stringent constraints on the axion-nucleon coupling (see, e.g., [51] for a review). However, the bounds are based on very few observational data [105,106]. Moreover, the axion production mechanism in the SN is difficult to describe in a complete way and several simplifying (and, often, unjustified) assumptions need to be made [92]. Therefore, even the accepted bound from the SN 1987A argument [44], where g 2 aN ≡ m 2 N C 2 aN /f 2 a parameterizes the axion-nucleon interaction, has to be taken more as an indicative result than as a sharp limit on the axion-nucleon coupling [92].
In this work, we have included the SN 1987A result as an additional constraint for our fits. More specifically, as discussed at the end of sec. A, we have assumed that the axion has a zero coupling to nucleons, with a (1 σ) uncertainty corresponding to g 2 aN 3.6 × 10 −19 . Though this procedure is far from rigorous, it does quantify the effect of the SN 1897A on the hinted parameter space, and it puts it on the same level as the other stellar evolutionary considerations.
A further complication is the choice of the most appropriate parametrization of g aN in terms of g an and g ap . In fact, the axion emission rate, ε, is a function, in general, of both couplings to protons and neutrons. More specifically, one can show that for nondegenerate nucleons [92,107] ε ∝ 2 3 I nn g 2 an + I pp g 2 ap + 2 9 I np (g an + g ap ) 2 + 26 9 where the functions I ij , corresponding to I(y i , y j ) in the notation of [107], depend on the stellar plasma conditions. In the case of the KSVZ axion, the coupling to neutrons is compatible with zero (within the uncertainties). Therefore eq. (B.1) leads directly to a bound on the axion-proton coupling [44].
More general models, however, such as the DFSZ or KSVZ A/J, predict interactions to both protons and neutrons, with coupling constants depending on the value of some additional parameter, e.g. tan β in the case of the DFSZ axion model. In order to extract the dependence of the constraint (B.1) on this additional parameter and to use this result in the global fits, we studied numerically the axion rate dependence on the couplings, using the same SN models exploited in [92]. The analysis reveals that the term proportional to g 2 an + g 2 ap in Eq. (B.2) is always dominant in the region where the axion emission rate is peaked (t 1 s and R 10 km). Therefore we adopted the choice corresponding to the constraint in Eq. (2.8).

C Projected sensitivity of next generation axion experiments
The region of the axion parameter space hinted by the cooling anomalies could be largely accessible to the next generation of axion experiments. In particular, as discussed in the text, IAXO and ARIADNE show a remarkable potential for both DFSZ and hadronic A/J models. In this section we discuss our methodology to extract the projected potential of these two experiments to probe the axion parameter space.

C.1 The IAXO helioscope
The helioscope is a particularly interesting experimental setup that has already shown the capability to effectively probe axion parameters close to those required to explain the anomalous cooling of stars [108]. Solar axions are transformed into x-ray photons in a magnetic field and then observed in a low background x-ray detector.
Axions can be produced in the Sun through mechanisms involving photons (most notably, the Primakoff process) and electrons (for example, electron bremsstrahlung and Compton processes). Therefore, the axion production rate depends, in general, on both g aγ and g ae . The detection, on the other hand, relies solely on the axion-photon oscillation in a magnetic field and is proportional to g 2 aγ . IAXO is a proposal for a next generation axion helioscope [109], which will improve of the CAST experiment at CERN, the current state-of-the-art axion helioscope. IAXO would implement a large multibore superconducting magnet with extensive use of x-ray focusing and low-background detection over a large cross-sectional area [15], improving current CAST sensitivity by a factor of more than 10 4 in signal-to-noise ratio. In this paper we refer to the most updated sensitivity projections of the experiment, as recently produced by the collaboration [110]. The information consists of two data sets. The first, d 1 (m a ), gives the minimal g aγ that IAXO can detect if the interaction with electrons is suppressed. The second, d 2 (m a ), provides the minimal √ g aγ g ae that IAXO can detect if axions are produced solely through interaction with electrons. Strictly speaking, these values are the median of the 90% C.L. excluded values over an ensemble of background-only Monte Carlo simulations. In both cases, we assume g aγ to be in units of GeV −1 .
The nominal projection of the experiment following its baseline configuration [15] is labeled as "IAXO" in our plots. The collaboration has also provided a more optimistic IAXO projection that is based on a series of possible upgraded experimental parameters in both magnet and detectors [110] beyond their baseline configuration of [15]. These upgrades would allow for an additional factor ∼10 in the signal-to-noise ratio over the nominal projection. This scenario is labeled as "IAXO+" in our plots. In view of our results here, the effort of going to this upgraded scenario seems very motivated. The data sets d 1 and d 2 are shown in figure 5, where we use blue for the nominal and purple for the upgraded case.
In general, axions can be produced through interactions with both electrons and photons and so the experimental potential is generally a combination of both datasets, d 1 and d 2 . To see this, let us write the axion production rate as R = (g aγ × GeV) 2 P + g 2 ae Q , (C.1) so that P and Q have the same units, and let us indicate with (g aγ × GeV) 2 C the conversion probability in the IAXO magnet. Then, if r is the distance to the Sun, the minimal number of events expected in IAXO is for the case of axions produced only through interaction with photons (Q = 0) or only through interaction with electrons (P = 0) respectively. Hence, d 4 1 P = d 4 2 Q. Notice that this relation does not depend on the axion properties since the axion interaction has been taken out.
For a specific model, one expects g ae = ζg aγ , where we consider g aγ in units of GeV −1 , so that ζ is dimensionless. In particular, ζ 0.20 sin 2 β for the DFSZ I model, and ζ Figure 5. Projected potential of IAXO and of ARIADNE. The left and center panels show the data sets d 1 and d 2 discussed in the text. The step-wise enhancement at high masses correspond to the planned buffer gas phase [15]. The right panel shows the ARIADNE sensitivity for the conditions discussed in the text. 0.12 cos 2 β for DFSZ II. Hence, in general, R = g 2 aγ P (1 + (d 1 /d 2 ) 4 ζ 2 ) and the expected number of events is where m a is in eV and g aγ in GeV −1 . Comparing this to g min aγ in eq. (C.4) allows to study the mass interval accessible to IAXO for each tan β (which enters only through ζ), as shown in figure 2.
Notice that the accessible mass interval tends to enlarge monotonically with tan β. This is easy to understand from eq. (C.4), remembering that tan β enters only through ζ. Eventually, for extreme values of tan β, low in the case of DFSZ I and high in the case of DFSZ II, the accessible mass interval becomes constant. This feature, easily noticeable in figure 2, is evident from eq. (C.4). In fact, for ζ small enough the dependence from tan β disappears all together. It is easy to estimate when this happens. From the IAXO data sets we notice that d 2 /d 1 is always confined between 5 and 9. 10 Therefore, ζ 2 (d 2 /d 1 ) 4 1 for ζ 10 −2 , which corresponds to tan β 0.2 for DFSZ I or tan β 0.3 for DFSZ II, in good agreement with figure 2.
The arguments discussed above imply that, in the DFSZ model, there is a minimal mass range accessible to IAXO, no matter the value of tan β. These mass ranges are shown in The case of the KSVZ-type axion/majoron, discussed in sec. 3.2.2, may appear somewhat more complicated to analyze since ζ has a dependence on m 1 through the first term on the right hand side of eq. (3.18). However, this term is always negligible in the region of interest for us. Repeating the arguments discussed above, we find the IAXO sensitivity regions shown in figure 3.
Notice that also in this case there is a minimal mass region always accessible to IAXO (see table 8), no matter the values of the other parameters, except in the case of the model with R Q = (3, 2, + 1 6 ), which features a very small photon coupling, cf. eq. (3.4), lowering the sensitivity of IAXO.

C.2 ARIADNE
The most compelling challenge for the axion helioscope is pushing the sensitivity to low masses 11 . Since the signal in the helioscope is proportional to m 4 a , its potential diminishes rapidly at low mass. Thus, probing QCD axion models much below m a ∼ 10 meV is prohibited even with the next generation of axion helioscopes.
Interesting alternatives to the axion helioscope paradigm, which are expected to be better suited to probe the meV mass range, are experiments that measure axion mediated long range forces [111]. In these cases, the mass controls not only the axion couplings, but also the length of the interaction, lower masses corresponding to longer interaction lengths. Hence these experiments are expected to be particularly effective in an intermediate mass range around m a ∼ meV, corresponding to an interaction length 1/m a ∼ a few 0.1 mm. A particularly compelling example of these experiments is ARIADNE [14], which shows potential to explore sections of the DFSZ axion parameter space.
Among those novel long-range forces, axion dipole-dipole interactions are the least model dependent but very difficult to probe in the near future [112]. A more tangible possibility is to measure axion mediated monopole-dipole interactions between a nucleus and a fermion. These forces are expected to exist if axions have a CP-violating interaction 12 with nuclei [111] L CP = g aN aN N . (C.7) 11 It is also difficult to probe the high mass region, since the axion-photon oscillations loose coherence at large masses. However, using a buffer gas it is possible to probe masses up to ∼ 1 eV [108] and QCD axions with higher masses are excluded by cosmological hot dark matter bounds. Therefore, here we are mostly concerned with the lower mass challenge. 12 Here we use convention of a bar above to indicate the CP violating interaction.
The scalar coupling with nuclei, g aN , is expected to be nonzero in the SM, because of the CP-violating contribution of the weak interactions [111,113,114], while experimental limits on the nEDM constrain it from above, g aN ≤ 10 −21 (f a /10 9 GeV) −1 . The axion dipole interaction with a fermion, g af , is defined by the interaction Lagrangian (2.1). Integration by parts leads to the usual axion-fermion interaction term which defines the pseudoscalar coupling g af = C af m f /f a , where m f is the fermion mass. ARIADNE searches for forces between a rotating cylinder, made of unpolarized material, and a vessel containing hyper-polarized 3 He gas [14]. Since the 3 He magnetic moment is dominated by the neutron contribution, the experiment would be measuring the monopoledipole interaction between nucleus and neutrons, proportional to |g aN g an |.
Proceeding analogously to what we have done for IAXO, we can define the ARIADNE sensitivity AS(m a ) as the minimal value of |g aN g an | that the experiment is expected to measure. For our analysis, we have extracted the AS from [14], assuming the most optimistic value for the scalar coupling g aN = 10 −21 (f a /10 9 GeV) −1 , a total integration time of 10 6 s and a transverse relaxation time of T 2 = 1000 s (cf. figure 2 in [14]). The resulting AS is shown in the right panel of figure 5.
The axion parameter space accessible to ARIADNE is the set of m a and tan β (for the DFSZ model) or m a and tr κ (for the A/J models) satisfying log |g aN g an |/AS(m a ) ≥ 0 .

(C.9)
This region is shown in figures 2 and 3. As evident from figure 2, ARIADNE could probe particularly well the low mass end of the DFSZ axion parameter space relevant for the cooling anomalies. There is therefore an interesting complementarity between IAXO and ARIADNE. The potential for the A/J models is less impressive.
Notice that, contrarily to IAXO, ARIADNE cannot cover the whole range of tan β for the DFSZ models. In fact, the DFSZ axion interaction with neutrons vanishes for tan β 0.8, and |g aN g an | becomes inaccessibly small for tan β ∼ 0.1 − 1. The same is true for the KSVZ A/J models, in which the coupling to neutrons is suppressed at low values of tr κ, and obviously for the pure KSVZ axions, in which the coupling to neutrons is compatible with zero. Therefore, in the case of ARIADNE there is not a minimal mass range accessible for all values of the other parameters, at least for the DFSZ and KSVZ axion models.
It is interesting to point out that these problems are inherent in the choice of the 3 He gas, in which the proton contribution to the magnetic moment vanishes. The choice of 3 He has undiscussed experimental advantages, most notably a long coherence time and the possibility to be highly polarized [14]. However, as we have seen, in most realistic models the axion couples to protons more strongly than to neutrons and, more importantly, the neutron interaction vanishes for some values of the parameters in all models considered here. Therefore, it would be certainly very beneficial to develop an alternative experimental setup which does not rely exclusively on the g an coupling for detection.
Finally, proceeding as we did above, it is straightforward to show that the future projection reach of the experiment, as discussed in [14] and shown with a solid blue line in figure 2 therein, would allow to cover most of the DFSZ axion parameter space. The possibility to scale the experiment to reach such sensitivity is feasible though, most likely, not in the very near future. Remember, however, that the projected sensitivity is still based on the assumption that the CP-violating scalar coupling is the largest possible allowed by measurements of the nEDM, an assumption not well justifiable. This requirement can be removed only by measuring dipole-dipole interactions, which do not depend on the scalar coupling. This would be in principle possible, using a polarized source mass, but it would present considerable more experimental difficulties. As shown in figure 3 of [14], reaching the parameter space relevant to realistic axion models in this case is, at the moment, still prohibitive.