Stellar Evolution confronts Axion Models

Axion production from astrophysical bodies is a topic in continuous development, because of theoretical progress in the estimate of stellar emission rates and, especially, because of improved stellar observations. We carry out a comprehensive analysis of the most informative astrophysics data, revisiting the bounds on axion couplings to photons, nucleons and electrons, and reassessing the significance of various hints of anomalous stellar energy losses. We confront the performance of various theoretical constructions in accounting for these hints, while complying with the observational limits on axion couplings. We identify the most favorable models, and the regions in the mass/couplings parameter space which are preferred by the global fit. Finally, we scrutinize the discovery potential for such models at upcoming helioscopes, namely IAXO and its scaled versions.


Introduction
The most elegant solution to the strong CP problem, the so called Peccei-Quinn (PQ) mechanism [1,2], implies the existence of the axion [3,4], which also provides one of the best particle physics candidates for cold dark matter (DM) [5][6][7]. The axion arises as the pseudo Nambu-Goldstone boson of a spontaneously broken global Abelian symmetry endowed with a mixed anomaly with the color gauge group SU (3) c , and the couplings and the mass of the axion are inversely proportional to the scale at which spontaneous symmetry breaking occurs. Already shortly after the axion hypothesis was conceived it became clear that this scale had to lie much beyond the electroweak scale, given that no axion signatures could be detected in laboratory experiments, 1 thus implying that axions must be very weakly coupled to ordinary matter and very light. Such light particles could then be produced in the hot and dense plasma of star cores and then freely escape, thus providing an additional channel for energy losses from astrophysical bodies that would affect their evolution. In fact, the strongest limits on axion couplings to electrons, nucleons and electromagnetic radiation come from the requirement that stellar lifetimes and energy-loss rates should not conflict with observations [9][10][11]. Astrophysical limits have been derived from the non-observation of axion emission from our Sun, from the concordance between the prediction of stellar evolution codes and direct observations of star distributions in the color-magnitude diagram (CMD) for evolved low-mass star populations, such as red giant branch (RGB) and horizontal-branch (HB) stars in globular clusters, from limits on cooling rates of white dwarfs (WD) and neutron stars (NS), and from the duration of the neutrino burst from the collapsed core of supernova (SN) SN 1987A. These limits have been frequently revised and updated in the literature, with the coming of better observations and more detailed analyses. However, intriguingly, rather than pushing the bounds on the various couplings to smaller and smaller values, the more recent astrophysical observations have hinted at finite, though small, axion couplings. More specifically, a set of astrophysical anomalies seem to indicate a preference for non-vanishing couplings of axions to electrons and, to a lesser extent, to photons. Although individually the significance of each one of these hints is marginal (1 or 2 σ), there is concordance among all the independent observations when interpreted in terms of axions [10], which raises the overall significance at about the 3 σ level. Analyses of other new physics candidates does not show the same level of agreement among the different observations [12], making axions the most interesting candidates to explain the anomalies.
In the last few years, the advent of new observations of remarkably improved accuracy has required the revision of some astrophysical bounds on axions. In particular, the RGB bound on the axion-electron coupling has been considerably improved and the significance of the hint has been reduced [13,14]. The supernova and neutron star bounds have also been updated in the last few years [15][16][17][18][19]. In light of these recent progresses, it seems timely to revise the global analysis of the stellar bounds/hints on axions, and to provide updated regions in parameter space where the experimental search may be particularly motivated. This is the task we address in this paper.
Besides a general, model independent analysis of the stellar observables, we carry out dedicated analyses focusing on specific axion models. Going beyond the approach followed in Ref. [10], the present study is not limited to scrutinize only the canonical DFSZ models, but also considers a set of non-universal realizations of DFSZ models, wherein same-type quarks of different generations can couple to different Higgs doublets [20][21][22]. In particular, we focus on models that feature the property of nucleo-phobia, i.e. the capability to strongly suppress the axion coupling to nucleons. Given the tight constrains stemming from observations of SN and NS, which strongly bound the axion-nucleon coupling, these models perform particularly well in reproducing the global set of astrophysics data. We will constrain the parameter region and ascertain the quality of the fit for each of the axion models of our representative sample, singling out those constructions that most successfully account for the observations.
As we shall show, the parameter space preferred by stars is the meV mass region (roughly, 1 to 100 meV). Intriguingly, this parameter region is the goal of the next generation of axion helioscopes. In particular, BabyIAXO, expected to be operative in the next few years, may already be able to dig into the parameter space hinted by stars for some QCD axion models. These regions will be better accessible by the full scale IAXO helioscope. The even larger IAXO+ will be able to essentially cover the entire parameter space preferred by stars for most well motivated axion models.
The paper is organized as follows. In Sec. 2 we review the axion effective Lagrangian and we define the coupling of the axion to the SM particles. In Sec. 3 we give an updated summary of the astrophysical observables relevant to our analysis. A review of the models scrutinized in this paper is given in Sec. 4. Our results are reported in Sec. 5, with Sec. 5.1 devoted to a discussion of the constrains on the axion couplings to SM fields in the various models, and Sec. 5.2 dedicated to analyse, for each model, the discovery potential at forthcoming axion search experiments. In Sec. 6 we resume our results and draw the conclusions. A detailed calculation of the axion-photon couplings for the non-universal DFSZ models is presented in Appendix A.

Axion effective Lagrangian
In this section we focus on the most relevant axion couplings from the point of view of astrophysics and experimental sensitivities for an axion mass scale ranging in the several meV region. The axion effective Lagrangian including photons and matter fields f = p, n, e (defined at a scale Λ < m p ) can be written as where the adimensional coefficients C aγ and C af read [23,24]  In Eqs. (2.2) and (2.6) E and N correspond respectively to the electromagnetic and QCD anomaly coefficients of the PQ current The coefficients C ψ (with ψ a SM quark or lepton) can be written as where c 0 ψ , which is defined in terms of the derivative interaction terms denotes the UV axion-fermion couplings as determined from the model-dependent PQ charges, while ψ is a correction due to possible fermion mixing effects that will generally arise, after electroweak symmetry breaking, in models in which the PQ charges are generation dependent. In all the models we will consider these corrections are small in the case of the quarks of the first generation. In some cases they can become sizeable for quarks of the second and third generation (see Table 3) but since their contribution as sea quarks to C ap , C an is already suppressed by small coefficients (see Eq. (2.5)) we can expected that Eqs. (2.3)-(2.4), which are derived under the assumption of no mixing effects, will still remain valid to a good approximation. Finally, in Eq. (2.6) the first term C e = c 0 e + e corresponds to the tree level axion coupling to the electron, while the second term is the one-loop contribution originating from a triangle loop involving two photons [25,26].
In the following, we will use the following rescaled axion couplings and employ the definitions g γ10 ≡ g aγ × 10 10 GeV and g e13 ≡ g ae × 10 13 . Finally the relation between the axion mass and the axion decay constant reads [27] m a = 5.691(51) 10 12 GeV f a µeV . (2.11)

Astrophysical observables and cooling anomalies
Astrophysical considerations have played a quite significant role in the investigation of the physics of light, weakly interacting particles, such as neutrinos, dark photons, axions, axion like particles (ALPs) etc. [9]. Constraints derived from the observations of stellar populations are often much tighter than bounds from direct searches. Intriguingly, a series of anomalous astrophysical observations have led to speculations that new physics may be at play in determining the details of stellar evolution [10,12,24,[28][29][30][31]. The axion (or ALP) case is especially compelling since, contrarily to other new physics candidates, fits particularly well all astrophysical observations [10,28]. A list of updated bounds and hints from stellar evolution on the axion couplings are summarized in Table 1. The aim of this Section is to review and update the status of the impact of axion emission on stellar evolution. We begin by reviewing the most relevant axion production mechanisms in stars, and next we summarize the constraints and hints from observations of different stellar systems.

Axion production in stars
The most relevant axion production mechanisms in stellar environments are the Primakoff, Compton, and electron bremsstrahlung processes. A pedagogical introduction to these processes can be found in Ref. [9], while a set of approximate expressions and numerical results can be found in the appendix of Ref. [38]. Here we just recall a few semi-analytical formulae for axion production through the different processes that will be useful for the following discussion.
In star core plasma with a relatively low density the most relevant axion production mechanisms are the Primakoff and Compton processes. The former is induced by the axion  [13] and the second based on the NGC 4258 galaxy [14]. The three NS bounds are from independent analyses, two based on observations of CAS A [16,19] and one on observations of HESS J1731-347 [17]. We have not listed the hint from the CAS A study in Ref. [37] since it is in tension with the more recent bounds shown in the Table. In the last column we give the references to the original analyses.
coupling to photons and the latter by its coupling to electrons. The Primakoff process is the photon conversion into an axion in the electric field of electrons or ions in the plasma: Neglecting degeneracy effects and the plasma frequency (a good assumption in environments in which the Primakoff process is the dominating axion production mechanism) it is possible to provide a semi-analytical expression for the energy-loss rate per unit mass due to axion emission [39]: where T and ρ are in K and in g cm −3 . The function F depends on the Debye-Huckel screening wavenumber κ S via the variable ξ ≡ κ S /2T , and can be explicitly expressed as an integral over the photon distribution (see Eq. (4.79) in Ref. [40]). An approximate expression, which agrees with numerical results at better than 2% over the entire range of ξ, is [39] F (ξ) 1.037 ξ 2 In general, F (ξ) is O(1) for relevant stellar conditions. More general numerical recipes, valid also in a degenerate medium, can be found in the appendix of Ref. [38].
In general the Primakoff process does not play a significant role in the evolution of superdense star cores, in which case other processes dominate. The Compton process γ + e → γ + a (3.4) accounts for the production of axions from the scattering of thermal photons on electrons. The Compton axion emission rate is a steep function of the temperature: where µ e = X j Z j /A j −1 is the mean molecular weight per electron with X j the relative mass density of the j-th ion and Z j , A j its charge and mass number respectively, n e is the number density of electrons while n eff e is the effective number density of electron targets. At high densities, degeneracy effects reduce n eff e , suppressing the Compton rate. Thus, this process is particularly effective in high-temperature environments, as long as the density is still relatively low to prevent electron degeneracy (cf. Fig. 1 in [41]). At higher densities and especially when electron degeneracy conditions are reached, the most efficient axion production mechanism is the electron/ion bremsstrahlung process e + Ze → e + Ze + a . (3.6) For degenerate plasma conditions, the axion energy-loss rates per unit mass can be approximated as The mild density dependence of the degenerate rate is accounted for by the dimensionless function F B . An explicit expression for this function can be found in [42] (see also Sec. 3.5 of Ref. [9] for a pedagogical presentation). Numerically its value is of order one for the typical stellar plasma conditions in which the degenerate bremsstrahlung process dominates ρ ∼ 10 5 − 10 6 g cm −3 and T ∼ 10 7 − 10 8 K.
In a nuclear medium such that of a SN collapsed core or in NS cores, the processes discussed above are generally sub-dominant [43]. If the axion coupling to nucleons is not particularly suppressed, as is the case in most axion models, a more efficient production mechanism is axion bremsstrahlung in nucleon-nucleon collisions where N, N = n, p. If we model the nucleon-nucleon interaction with the exchange of a single pion, it is evident that the pion mass in the propagator will suppress the emission rate, unless the temperature is such that the typical momentum exchanged in the collision, which is of the order of the nucleon momentum q N ∼ (3m N T ) 1/2 , is larger than the pion mass. This demands T 10 MeV, a temperature typical of SN and NS cores. Therefore, among the various stellar objects, SN and NS provide the best environments to test the axion nucleon coupling.
Approximate expressions for the nn scattering emission rates (the pp scattering is similar) in the limit of non-degenerate (ND) and degenerate (D) nuclei are given below [9] ε ND ≈ 2.0 × 10 38 g 2 an ρ 14 T 3.5 30 erg g −1 s −1 , where T 30 = T /30MeV and ρ 14 = ρ/10 14 g cm −3 . These expression are based on a series of approximations, 2 however, they still provide a good order of magnitude estimate of the axion emission rate from a SN or NS, and in particular they put in evidence the steeper temperature dependence in the case of a degenerate medium. However, throughout this work we use the more accurate numerical results from Ref. [15] (see Sec. 3.5).

Axion Bounds from White Dwarfs
WDs represent the last stage of the evolution of low mass stars, following the exhaustion of the nuclear fuel. Therefore, during this phase the star is just cooling. WDs are characterized by a dense core of degenerate electrons, with typical density of about 10 6 g cm −3 and a core temperature dependent on the age of the star. Young WDs are hot and brighter, and cool through volume neutrino emission. At later times, the photon surface cooling dominates. The addition of exotic particles, such as axions, can have a strong impact on the evolution of the WD, accelerating the stellar cooling. Thus, testing the cooling efficiency is an indirect way to probe the existence of new physics. There are two independent ways to test the cooling efficiency of WDs. One is to observe the WD luminosity function (WDLF), which shows the WD number distribution in different luminosity bins (see Sec. 3.2.1). The other, is to measure changes in the period of WD variables (WDV), a class of WD whose luminosity periodically changes with time (see Sec. 3

.2.2).
Both methods indicate a preference for some unidentified cooling, which could well be provided by axion emission provided their coupling to electrons is of a few 10 −13 . Below we give more details on the analyses of WD cooling rates.

White Dwarf luminosity function
As we have already mentioned, the WDLF represents the number distribution of WDs as a function of their luminosity. This distribution has been used for decades as a tool to measure the WD cooling efficiency, since the number of WDs with a certain luminosity obviously depends on how efficiently the star looses energy. If axions exist, they would be produced in a WD core primarily through the bremsstrahlung process. Hence, the WDLF has provided information about the axion coupling to electrons (see Ref. [44] for a recent review). Unfortunately, observations from the Sloan Digital Sky Survey (SDSS) and the and SuperCOSMOS Sky Survey (SSD), on which the current studies are based, are not consistent within their quoted error bars, indicating that the systematic uncertainties in these observations might have been underestimated [45]. Consequently, there is no complete consensus on the exact bound on the axion-electron coupling derived from observations of the WDLF. A summary of different studies can be found in Ref. [33]. In that analysis, the authors enlarged the error bars to take into account not only their systematic uncertainties, but also the discrepancies between the SDSS and SSS observations (see discussion therein, in sec. 4.2). In the present study, we adopt the result g ae ≤ 2.1 × 10 −13 at (2 σ), derived in Ref. [33].
Furthermore, most studies of the WDLF seem to indicate an excessive cooling with respect to the standard prediction. Such cooling has been interpreted as due to an axionelectron coupling g ae (1.4 ± 0.3) × 10 −13 (at 1σ) [33]. 3

White Dwarf Variables
The WDV are a set of WDs whose luminosity changes periodically. The period P ranges from a few to several minutes, depending on the particular star. It is well known that observations of the secular change,Ṗ , of the WDV period provide information about the efficiency of the WD cooling. In fact, to a very good approximationṖ /P is directly proportional to the cooling rateṪ /T . Hence, an accurate measurement ofṖ allows to set bounds on possible sources of extra cooling (see Ref. [34] for a comprehensive review).
For over two decades, observations of the period decrease (Ṗ /P ) of particular WDVs have shown discrepancies with the expected behavior. This is clear from the data in Table 2, which shows the results of the observations of the WDVs analyzed so far [34]. The systematic tendency of the observedṖ to be larger than the expected values can be interpreted in terms of a new particle, produced in the core and efficiently carrying energy outside the star. Even though specific observations have been interpreted in various ways, for example in terms of an anomalously large neutrino magnetic moment (see, e.g., Ref. [49]), a global analysis of all the data indicates a preference for axions among other WISPy candidates [28], and identifies the coupling with electrons in the range g ae = 2.9 +0.6 −0.9 × 10 −13 (at 1σ) [10], with a 2σ bound of g ae ≤ 4.1 × 10 −13 , as reported in Table 1. The analysis in this paper is based on the recent review [34], where the most updated studies on the viable WDV are considered. The results for the axion couplings are summarized in Tab. 2.
We conclude observing that, following the approach usually employed in the literature and first outlined in Ref. [28], we will not include in our fits data relative to the WD G117-B15A: the analysis for this star shares many theoretical similarities to the one relative to R548, but the experimental results are somewhat different, with the hint on g ae stemming from the former noticeably stronger than the one inferred from the latter. Given that results relative to these two WDs are based on very similar hypothesis, we will therefore conservatively include in our fit only data pertaining to R548, along with the data from PG 1351+489 and L 19-2.

Axion Bounds from the tip of RGB stars in globular cluster
After completing the evolution through the main sequence, characterized by a H-burning core, low mass stars begin climbing the RGB, well visible in the color-magnitude diagram as a diagonal line starting at the main sequence and directed toward the colder and higher luminosity region. During the evolution in the RGB, stars are characterized by a He core 3 Such result is not free from controversy. A later study of the hot part of the WDLF [46] did not confirm this anomalous behavior. However, the hotter section of the WDLF has much larger observational errors and the axion (or ALP) production would be almost completely hidden by standard neutrino cooling in the hottest WDs. The most recent work [47] seems to confirm an excessively efficient cooling, which can be explained by the emission of axions. Alternatively, it has been recently proposed that the excess of cooling could be explained by neutrino pair synchrotron emission enhanced by extremely large magnetic fields, in excess of 10 11 Gauss, confined in the WD cores [48]. This possibility can in principle be bounded by dedicated asteroseismology surveys.
Star Table 2. Measured and expected value ofṖ for a set of WDV. Here, P is the period of the variable star andṖ its time derivative. The interpretation in terms of axion coupling to electrons is also shown. (Data from Ref. [24] except for G117-B15A, which corresponds to the updated analysis in Ref. [50]). and a burning H shell, whose ashes increase the He core mass, while the star luminosity (determined by equilibrium at the surface of the He core between thermal pressure supporting the non-degenerate envelope against the gravity pull from the core) keeps growing. The process continues until the core reaches sufficiently large temperatures and densities (T ∼ 10 8 K, ρ = 10 6 g cm −3 ) to ignite He, an event known as the He-flash. At this stage the star has reached the maximum luminosity, that is the RGB tip (RGBT), after which it shrinks and moves to the HB. If an additional core-cooling mechanism is at play, He ignition is delayed, the core would accrete a larger mass, and the star would reach higher luminosities. Therefore, measurements of the luminosity of the RGB tip allow to test the rate of cooling during the RGB phase. The method is particularly effective for constraining g ae since in red giant cores axions can be efficiently produced via electron bremsstrahlung. We denote by M I,TRGB the luminosity of the tip of the RGB in globular clusters (GC). Following Ref. [41], based on the analysis in Ref. [28,51], it is possible to derive the following analytical expression for the expected magnitude of the RGBT: that has an associated theoretical uncertainty σ 2 = 0.039 2 + (0.046 + 0.012g e13 ) 2 . This should be compared with the observational values. The latest analyses are those in Refs. [13,14]. The first, is based on the global analysis of a sample of 22 GC while the second analyzes individual clusters and galaxies. Both studies give very similar results (cf. Table 1). Here, as a reference, we consider the observational value found in the comprehensive analysis from Ref. [13]. This choice results in the bound on the axion-electron coupling g e13 ≤ 1.5 (2 σ), which is perfectly compatible with the one found in Ref. [14], but also hints to a non-vanishing value for g e13 . At any rate, we have verified that removing this additional hint from our fit does not alter significantly the results.

Helium burning stars
After He ignition, the RG core expands and the star migrates to the horizontal branch of the color magnitude diagram, characterized by a He burning non-degenerate core. The core of a HB star has a density of about ρ ∼ 10 4 g cm −3 , which is about two orders of magnitude less than that of a RGB star or a WD. In this environment, axions are efficiently produced through the Primakoff and Compton processes. The effect of the additional energy loss provided by axion production and emission is to accelerate the helium consumption in the HB core and, consequently, to reduce the lifetime of this stage. A very efficient way to probe this effect is by measuring the so called R parameter, R = N HB /N RGB , which measures the ratio between the number in the HB and in the upper portion of the RGB in GC. Historically, measurements of the R parameter have been used to derive bounds on the axion-photon coupling [35,36,52,53]. However, the axion electron coupling can also impact the R parameter. Following the results of Refs. [10,24,28], we present the expected R parameter in the following form where R 0 (Y ) is a function of the helium abundance (Y ) in the GC and F aγ , F ae are some positive-definite functions of the axion couplings. Evidently, the impact of axions on the R parameter may be induced by its couplings with photons as well as by its coupling with electrons. However, surprisingly, there are no explicit stellar numerical evaluation of the R parameter which include both couplings. An approximate expression for the functions F , based on a series of older numerical results, was presented in Refs. [10,28]. There, it was found 14) The theoretical expression in Eq. (3.11) should be compared with observational results, in order to constraint the axion couplings. Ref. [35], reported the value R = 1.39 ± 0.03 from the analysis of 39 clusters, and used the result to derive the bound g aγ ≤ 0.65 × 10 −10 GeV −1 at 2 σ [35,36], under the assumption that the axion couples only to photons. More massive He burning stars, with mass M ∼ 8−12M , can also provide some insight into the axion-photon coupling [39,54]. The (core) He burning stage of these stars is characterized by a migration towards the blue (hotter) region of the CMD and back. This journey is known as the blue loop. The existence of the loop is corroborated by many astronomical observations. In particular, this stage is essential to account for the observed Cepheid stars (see, e.g., Ref. [55]). Ref. [39], based on numerical simulations of solar metallicity stars in the 8 − 12M mass range, showed that a coupling larger than ≈ 0.8 × 10 −10 GeV −1 would cause the complete disappearance of the blue loop while a somewhat lower values of g aγ might help explain the observed deficiency of blue with respect to red supergiants, as is discussed, e.g., in Ref. [56]. The numerous uncertainties in the microphysics and in the numerical description of the blue loop stage, however, have not permitted a more quantitative assessment of this possibility [12]. Hence, these results will be ignored in our present work.

Supernovae
The observation of the neutrino signal from SN 1987A, and the recent observation of a NS associated with it [57] supported the picture of a neutrino driven SN explosion with neutrinos carrying away about 99% of the energy released in the explosion. Since new weakly coupled particles could accelerate the cooling and reduce the observed duration of the neutrino signal, SN 1987A nutrino data have been widely used to constrain models of new physics, and particularly axion properties [40,[58][59][60].
The plasma in the core of a SN, in the first few seconds after the explosion, is extremely hot (T ∼ 30 MeV) and dense (ρ ∼ 3 × 10 14 g/cm 3 ). This makes the SN an extremely interesting environment to study new physics. Axions can be produced in this environment through different processes, including the Primakoff process [61][62][63], and the electron bremsstrahlung [64]. However, for couplings allowed by other astrophysical considerations the nuclear bremsstrahlung, driven by the axion coupling to neutrons and protons, largely dominates over the latter two. 4 The most recent analysis of the axion-nucleon bremsstrahlung [15] give the bound which is the one we adopt in the present study. If the axion coupling to nucleons is large enough, axions may be trapped in the SN core [15,71]. In this case, the cooling efficiency is reduced and the axion bound weakened. Unfortunately, the exact value of the coupling for which axions are sufficiently trapped so that the observed duration of the neutrino signal is not reduced, is afflicted by several uncertainties. However, in most models in which the SM fermions carry a PQ charge the trapping regime is already excluded by other stellar bounds. Thus, in the present work we will not consider further this possibility.

Neutron stars
Observations of the cooling of NS also provide information about the axion-nucleon coupling [16-19, 72, 73]. Several NS bounds on g aN exist in the literature, however, they are not always consistent with one another. One of the most studied NS is the one in CAS A. In recent years, there has been some speculation that its anomalously rapid cooling could be a hint of axions with coupling to neutrons [37] g an 4 × 10 −10 . However, the data can also be explained assuming a neutron triplet superfluid transition occurring at the present time, t ∼ 320 years, in addition to a proton superconductivity operating at t 320 years [16]. Under these assumptions, it was possible to fit the available data well, leaving little room for additional axion cooling, corresponding to the bound [16]: A stronger bound, though only on the axion-neutron coupling, was inferred from observations of the NS in HESS J1731-347 [17]. This bound is also in good agreement with a newer analysis of the NS in CAS A carried out in Ref. [19] (cf. Table 1).
In the present work, we take the bound in Eq.(3.17) as a reference for the NS bound on the axion-nucleon coupling. It should be remarked, however, that there is no universal consensus on the NS bound. For example, Ref. [18] recently proposed a considerable less stringent result, (3.18) For this reason, we present separate analyses, with and without the NS results. 4 It was recently pointed out in Ref. [65,66] that pion processes π − + p → a + n might dominate the axion production rate, becoming even more efficient than the nuclear bremsstrahlung at some temperatures and densities. This process has been discussed for a long time in the literature [67][68][69]. Only recently, however, it was shown that pion abundance in the early stages of a SN is much larger than previously expected [70], perhaps to the point of making the pion production mechanism the dominant one [65]. We have not included these results in our analysis since the current studies, although focused on the observational consequences of the pion-induced axion production mechanism, do not derive new bounds on the axion-nucleon couplings [66].

Model-independent fit
In this subsection we describe the global fit to the relevant astrophysical observables described above, that we have performed by taking the axion couplings g ae and g aγ as uncorrelated parameters. Given the absence of any particular assumption about an underling axion model, this first analysis is denoted as model-independent. Note that this implies that g ae and g aγ are not correlated to the couplings to nucleons either, and hence it is consistent to neglect, in the model-independent fit, the data inferred from SNe and NS observations. However, this additional set of data will be taken into account in the model-dependent fits performed in Sec. 5 where, given the model, all the couplings are correlated in a well defined way. In short, in this subsection we will consider only data pertaining to WD, RGB and HB stars. We depict in Fig. 1 the 1σ, 2σ and 3σ regions allowed by the fit, showing also the isolines for the ratio C ae /C aγ . The best fit point is marked with a black dot, and corresponds to g ae 1.2×10 −13 and g aγ 0.18×10 −10 GeV −1 , with C ae /C aγ O(10 −2 ). It is interesting to notice that the present data, including the information relative to HB, can be accommodated by the introduction of a non-vanishing axion-electron coupling only. On the other hand, the SM case (g ae = g aγ = 0) is excluded by present data at the about the 3σ level.
Let us add a comment regarding a comparison with the global analysis in Ref. [10]. The results of this previous study are qualitatively similar to what is shown in our Fig. 1. The most relevant quantitative difference is the RGB contribution to the global fit. Ref. [10] was based on the RGB analysis in Ref. [51], which adopted an incorrect screening prescription for the nuclear reaction rates [74]. The correction of this numerical issue in Refs. [13,14] strengthened substantially the bound on the axion-electron coupling, and caused a shift of the entire g ae region in Fig. 1 to the left. As stated in Sec. 3.3, we employed in our fit the results obtained in Ref. [13], but we tested the stability of our findings under this choice performing also a fit where we used the results from Ref. [14]. We found that the contours of the preferred regions differ only slightly between the two choices, with the SM case always excluded at the ∼ 3σ level. The results shown in Fig. 1 therefore give a reliable representation of the region preferred by the stars.
As anticipated above, in Sec. 4 we will review a selection of motivated axion models, each of which will imply specific correlations among the various axion couplings in terms of a few model parameters. The impact of the stellar evolution data on the corresponding bounds and the effectiveness of those specific models in accommodating the hints for extra energy emission from stars are discussed in Sec. 5.

A representative sample of axion models
In this Section we introduce a set of explicit axion models which yield different axion couplings to the nucleons, the electrons and the photons. In Sec. 5.1 we will consider how well they can perform in addressing the issue of possible anomalies in stellar energy losses, accommodating the observational hints while respecting all other phenomenological bounds.
The model-independent analysis of Sec. 3.7 suggests that promising candidates among axion models should comply with a first requirement of predicting a sizeable C ae /C aγ ratio. For this reason KSVZ models [75,76], in which the C ae coupling arises radiatively from a triangle loop involving two photons, and hence is induced by C aγ , are not well suited to explain the cooling hints. Indeed, for KSVZ models one obtains where in the next-to-last relation we have selected the log-enhanced contribution and assumed E/N 1.92, and in the last relation we have chosen as reference value f a = 10 9 GeV. As it can be seen from Fig. 1, the iso-line C ae /C aγ = 10 −3 already misses the 1σ region preferred by the global fit, so that the KSVZ value of 10 −4 cannot provide a good fit to the stellar cooling hints. Hence, we will not considered the KSVZ model any further in the present study.
A better starting point for addressing the cooling hints is provided by DFSZ models [77,78] since the axion coupling to electrons arises at tree level. In particular, the nucleophobic axion models of Refs. [20][21][22]79] allow to relax astrophysical bounds from NS/SN, while keeping at the same time a sizeable coupling to electrons and photons. An additional motivation to consider the latter class of models consists in the fact that their generation dependent PQ charge assignment is engineered in such a way that the contributions of two generations to the QCD anomaly factor cancel out yielding N = N 1 + N 2 + N 3 = N 1 . This implies a DW number N DW = 2N 1 = 2 or 1 (rather than N DW = 2N = 6 or 3 as in canonical DFSZ models) [20]. The importance of the existence of DFSZ-like models with N DW = 1 consists in the fact that they are strongly preferred in post-inflationary PQ-breaking scenarios since they are free from the DW problem. In this case in fact the network of axionic strings coupled to a single DW spontaneously decays around the time when the axion acquires a mass. In this scenario, under the assumption that the shape of the instantaneous axion emission from string decays is IR dominated [80,81], cosmological considerations based on the requirement that the axion relic density will not exceed that of the DM, together with the requirement that axion interactions will be sufficiently suppressed to impede axion thermalization into dark radiation, yield a viable window for the axion mass between 0.2 and 100 meV [82]. 5 As can be seen from Table 4 in Sec. 5.1, this agrees well with the mass window favoured by the cooling hints.
In this Section we first review the basic features of the above-mentioned DFSZ-like axion models, and next we generalize further the classification of axion models with two Higgs doublets and generation dependent PQ charges by including novel constructions. In particular, we show that models featuring particularly large coupling to photons and electrons can still be nucleo-phobic, a feature that allows to optimize the fit to the cooling anomalies while bringing the corresponding axions within the discovery potential of future helioscopes such as (Baby)IAXO. In spite of the fact that we analyze a reasonably large number of different constructions, we stress that our list is far from representing the entire panorama of axion models. In particular, it leaves out models in which the axion coupling to electrons can be enhanced by specific mechanisms (see e.g. [86][87][88]). Nevertheless, it is unlikely that these alternative possibilities could yield much better fits to the hinted anomalies than our representative subset. This is because the best fit conditions are realised when g ae /g aγ ≈ 0.6×10 −2 GeV and g ap , g an ≈ 0, and indeed this parameter space point is already approached in some of the models we have considered.

DFSZ1
Leaving Yukawa coupling constants understood, the Yukawa sector contains the following operators where i, j = 1, 2, 3 are generation indices, q i , i denote the quarks and leptons SU (2) L doublets and u i , d i , e i the right-handed (RH) singlets. The corresponding axion coupling coefficients are The exact value for the lower limit is subject to systematic uncertainties in the simulations of topological defects. If instead the axion emission spectrum is not IR dominated, the leading contribution to the axion relic density comes from the misalignment mechanism, and the lower value for the axion mass drops to m a > ∼ 25 µeV [83,84]. In pre-inflationary PQ breaking scenarios N DW > 1 does not constitute a problem since all topological defects are inflated away. In that case the axion relic density is generated solely by the misalignment mechanism and, while axion masses in the meV range cannot be excluded, if DM is dominantly composed by axions then limits on iso-curvature fluctuations generated during inflation constrain m a to lie below the meV scale [85].
with c β ≡ cos β, s β ≡ sin β and tan β = H u / H d ≡ v u /v d . Requiring that the quark Yukawa couplings remain in the perturbativity domain restricts the vacuum angle β to lie withing the interval tan β ∈ [0.25, 170] (the perturbativity limits for the models considered in this work are reviewed in Sec. 4.3). Note that for both DFSZ1 and DFSZ2 the axion couples to the SM fermions in a generation-independent way so that there are no corrections to the axion couplings from intergenerational mixing effects, and hence C ψ = c 0 ψ , see Eq. (2.8).

DFSZ2
The Yukawa sector contains the following operators 4) and the corresponding axion coupling coefficients are The vacuum angle can range in the same perturbativity interval than in DFSZ1, that is tan β ∈ [0.25, 170].

Non-universal DFSZ models
We denote as "non-universal" those models that have the same scalar content than DFSZ1 and DFSZ2 (two Higgs doublets and one SM singlet Φ) but for which same-type quarks of different generations can couple to different Higgs doublets. Clearly, in this case the labels u and d for the Higgs doublets loose their meaning, and it is more convenient to employ a notation where the two scalars have the same quantum numbers, and denote them as H 1,2 ∼ (1, 2, −1/2). The vacuum angle is defined as tan β = H 2 / H 1 ≡ v 2 /v 1 . For all these models, the requirement that the PQ current is orthogonal to the hypercharge current (which implies orthogonality between the respective Goldstone bosons) fixes the PQ charges of the two Higgs doublets as X 1 = −s 2 β and X 2 = c 2 β . In the following, we review a set of models which can feature the property of being nucleo-phobic, namely for a specific choice of the vacuum angle β the coefficients C p,n can be strongly suppressed with respect to their natural DFSZ values. A non-trivial result, for which a proof is given in App. A, is that in a general non-universal DFSZ model with two Higgs doublets, the model dependent ratio E/N can only span over the following finite set of values: Although larger values of E/N can be achieved by introducing an arbitrary number of extra Higgs doublets (see e.g. Refs. [87,89,90]), here we focus only on the case of two (or three, see below) Higgs doublets. In particular, for the two Higgs doublet case we consider the nucleo-phobic models M1, M2, M3 and M4 introduced in Ref. [20], and the T (u,d) 2 models of Ref. [21]. The latter feature E/N values corresponding to the last two entries in the list above, which result in the largest coupling g aγ . For the case of three Higgs doublets we consider the nucleo-phobic 3HDM model of Ref. [22], where the suppression of the axion-nucleon couplings is strongly correlated with the value of the axion coupling to the electrons.
Given that in all these models the PQ charge assignments are generation dependent, the axion couplings to quarks and leptons will in general receive corrections from intergenerational fermion mixing. These will depend on the particular model and on specific assumptions on the mixing matrices for the left-handed (LH) and right-handed (RH) states.
As regards the mixing independent part of the couplings, they can instead be written in a model independent way in terms of the PQ charges of the SM fermions as [24] c 0 u, c, t = is the QCD anomaly factor which is not affected by corrections from fermion mixing.

Non-universal 2+1 models
The M1, M2, M3 and M4 models are characterized by a 2+1 structure of the PQ charge assignments, namely two generations replicate the same set of PQ charges. Note that as explained in Ref. [20] in this case all the entries in the up-and down-type quark Yukawa matrices are allowed and there are no texture zeros. We recall below for each model the structures of the Yukawa operators and the charge assignments for the different generations. More details can be found in Ref. [20].

M1
The Yukawa sector of the M1 model contains the following operators The PQ charges of the first generation are replicated for the second generation, which thus appears in an analogous set of operators with the generation label 1 replaced by 2. The PQ charge assignments are The anomaly coefficients and the mixing independent part of the axion couplings are with tan β ∈ [0. 25,170]. Since in all the sectors the charges of the RH fields of the three generations are the same, there are no RH mixing corrections. In the LH sector mixing effects enter because the third generation has different charges from the first two. For the quarks, if we assume that the LH mixing matrix has CKM-like numerical entries, mixing effects are small and can be neglected. Accordingly, we will simply use C u,c,t,d,s,b c 0 u,c,t,d,s,b . In the lepton sector instead there are no reasons to expect a particular suppression of e L -τ L mixing, hence C e = c 0 e + L where the expression for L ∈ [−1, 1] can be found in Ref. [20].

M2
The Yukawa sector of the M2 model contains the following operators (4.11) The PQ charge assignments are The anomaly coefficients and the mixing independent part of the axion couplings are Note that for this case the range for the vacuum angle allowed by Yukawa perturbativity is tan β ∈ [0.0024, 4.0]. As it can be seen from the charge assignments, there are no mixing effects in the LH and up-type RH quark sectors (hence C u,c,t = c 0 u,c,t ). As regards the RH down-type quarks, although there is no phenomenological information on the RH mixings, we will assume that the related effects are negligible, and we will approximate C d,s,b c 0 d,s,b while, for the leptons, we adopt the same prescription used for M1.

M3
The M3 models is defined by the following set of Yukawa operators: where, as in M1 and M2, the quark charges of the second generation replicate those of the first one, while for the leptons they replicate the charges of the third generation. Explicitly, the PQ charge assignments are The anomaly coefficients and the mixing independent part of the couplings are with tan β ∈ [0.0024, 4.0]. We can easily read off the charge assignments the following relations: C d,s,b = c 0 d,s,b , C u,c,t c 0 u,c,t (assuming RH up-type quark mixings are small) and C e = c 0 e + L .

M4
For the M4 models the Yukawa sector contains the following operators where as in M3 for the second generation the quarks replicate the charges of the first generation, while the leptons replicate the charges of the third generation. The PQ charge assignments are The anomaly coefficients and the mixing independent part of the couplings are with tan β ∈ [0.0024, 4.0]. The structure of the charge assignments in generation space implies C u,c,t = c 0 u,c,t , C d,s,b c 0 d,s,b (assuming RH down-type quark mixings are small) and C e = c 0 e + L .

1+1+1 models
In Ref. [21] two other models, denoted as T (u,d) 2 , were considered, that are characterized by a 1+1+1 structure of the PQ charges in the quark sector, namely all generations can have different PQ charges. 6 They correspond to the class of nucleo-phobic models for which the Yukawa matrices have the maximum number of texture zeros that still allows to reproduce the quark masses and CKM mixings. Below we recall the main features of these two models.
The Yukawa sector of the T (u) 2 model contains the following operators (4.20) 6 The PQ charges in the lepton sector are assumed for simplicity to be universal. However, employing a generation dependent U (1) symmetry for leptons can yield non-trivial predictions also for neutrino masses and mixings [91].
For the quarks, out of the eighteen operators allowed by the SM gauge symmetry only nine are allowed by the PQ symmetry while the remaining nine, which do not appear in Eq. (4.20), are forbidden. The PQ symmetry instead acts universally in the lepton sector so that all the Yukawa operators are permitted and there are no leptonic mixing effects. The structure of the operators in Eq. (4.20) is enforced by the following set of charges The numerical values of the non-vanishing Yukawa couplings that is required to reproduce the quark masses and CKM mixings (at the PQbreaking scale) was computed in Ref. [21] y (1,1) where (with a slight abuse of notation) v 1,2 in the denominators here stand for the dimensionless numbers H 1,2 /GeV. Diagonalization of the resulting mass matrices provides the numerical values of the LH and RH quark mixing matrices, so that the complete expressions of the axion couplings to the quarks can be readily obtained. We have (4.23) The first condition to realize nucleo-phobia is C u + C d = 1 [20] which is satisfied at the 10% level. The second condition C u − C d = 1/3 is obtained for tan β = √ 2 which is well within the perturbative Yukawa window tan β ∈ [0. 117,145].
The Yukawa sector of the T (d) 2 model contains a subset of nine quark operators, while for the leptons all Yukawa operators are allowed. They are: (4.24) The quark Yukawa structure in Eq. (4.24) is enforced by the following set of charges: The values of the Yukawa couplings required to reproduce the quark masses and CKM mixings is [21] y (1,1) where again v 1,2 stand for H 1,2 /GeV. The axion couplings, including quark-mixing corrections are: (4.28) The first nucleo-phobic condition C u + C d = 1 is satisfied at the 10% level, while the second one C u − C d = 1/3 is obtained for tan β = √ 2, well within the perturbative window which is well within the perturbative Yukawa window tan β ∈ [0.010, 8.6].

3HDM model
A three Higgs doublets model (3HDM) which, besides nucleo-phobia, also allows to accomplish electro-phobic conditions (i.e. approximate axion-electron decoupling) was introduced in Ref. [22]. In this model the leptons couple to a third Higgs doublet, H 3 , with the same charges for all generations i e jH3 . (4.29) In the quark sector the second generation replicates the PQ charges of the first generation, namely the structure is 2 + 1 like in the M 1 model.
For the scalar potential, the following couplings with the SM singlet scalar φ is assumed: These two operators fix the PQ charges of H 1,2 in terms of the charge X 3 of H 3 , as X 1 = −2 + X 3 and X 2 = 1 + X 3 . The anomaly coefficients and the mixing independent part of the couplings are c 0 e,µ,τ = and the condition of orthogonality between the PQ and the hypercharge scalar currents yields X 3 = 3(c 2 1 − 1)c 2 2 where c 1 = cos β 1 , c 2 = cos β 2 and the two vacuum angles are defined in terms of the doublet VEVs H 1,2,3 = v 1,2,3 as tan The values of the vacuum angles β 1,2 are subject to the following perturbativity constraints [22]: where y t,b,τ = √ 2m t,b,τ /v are the heavy fermions Yukawa couplings expressed in terms of the quark masses and v ≡ (v 2 1 + v 2 2 + v 2 3 ) 1/2 = 246 GeV.

Summary of axion models and discussion
The values of E/N and of the axion-fermion couplings relevant for the astrophysical analysis for all the models reviewed in this Section are collected in Table 3. In general, in models with non universal PQ charge assignments, axion interactions with matter fields will feature a certain amount of flavour violation. Searches for kaon decays of the type K → πa provide the strongest constraints on flavour-violating axion couplings (see for example Table 2 in Ref. [92]). For model M3 the charge assignments for the quark doublets and for the RH down type quarks respect universality, so this model automatically evades the corresponding limits. Model M1 features universality of the couplings for the RH quarks of the three generations, but for the LH quark doublets only the first two families have equal charges. In this case, if the off-diagonal entries V d 13 , V d 23 in the LH mixing matrix are not particularly suppressed, the predicted K → πa decay rate could easily conflict with the experimental limit. If instead the mixings are at most CKM-like (V d 13 , V d 23 1) model M1 would also evade the K → πa constraints. Because of the unitarity constraints on the mixing matrices, this also ensures that the diagonal couplings of the axion to the first generation quarks would be negligibly affected so that, to a good approximation, the conditions for axion nucleo-phobia remain preserved. Models M2 and M4 feature universality for the LH quark doublets and RH up-type quarks, however in the RH down-sector the charges of the first and second generation differ. In this case even for CKM-like RH mixings the decay rate would exceed the experimental limits by several orders of magnitude [20]. Nevertheless, as it was shown in Ref. [21], it is possible to assume some specific mass-matrix textures in such a way that, for example, s R does not mix with the other RH down-type quarks, in which case K → πa (and also B → Ka) limits are easily evaded. For example, the 1+1+1 models T u,d 2 incorporate precisely this feature (cf. Ref. [21]). Although in these models the PQ charges of all the three generations differ, it is the PQ symmetry itself that enforces suitable matrix textures such that dangerous off-diagonal mixings are absent and flavour-diagonal axion couplings to the quarks of the first generation are negligibly affected. Hence flavor violation limits can be evaded while preserving the nucleo-phobic property of the axion.
A second issue regards the range in which the vacuum angles that define the couplings of the physical axion are allowed to vary, and that should respect the perturbativity constraints on the relevant Yukawa couplings. In order to estimate the perturbativity domain we have employed the tool of perturbative unitarity on the 2 → 2 scatterings of the Yukawa theory in the presence of two or more Higgs doublets. In particular, taking into account gauge group factors (see e.g. [93][94][95]) one gets [22]: y q < 16π/3 for the Yukawa couplings of the quarks and y < 4 √ 2π for the charged leptons, where y q, are the Yukawa couplings defined in the multi-Higgs doublet theory that have to be matched with their SM counterparts in order to extract a bound on the electroweak vacuum angles. This matching presents some model dependency related to the scale at which the extra Higgs doublets are integrated out, which can vary between the PQ and the electroweak scale. For DFSZ1-2, M1-2-3-4 and 3HDM models we perform the matching at the electroweak scale (which allows to neglect top running effects in the low-energy axion couplings [96][97][98][99][100]), while for T (u,d) 2 models which feature textures in the quark mass matrices that are enforced by the PQ symmetry, we perform the matching at the PQ scale, and we neglect for simplicity top-related running effects.

Axion models in the light of cooling anomalies
This Section is devoted to the interpretation of the cooling anomalies discussed in Sec. 3 in a model-dependent way, putting under scrutiny the axion models presented in Sec. 4. Specifically, the results of the fits to the astrophysical data will be presented in Sec. 5.1, where the ability to reproduce data will be scrutinized and the models capable to provide a better accordance with observations will be singled out. Then, in Sec. 5.2, we will analyze the experimental potential to probe these models in the region relevant for stellar evolution.

Axion fits to cooling anomalies
In order to carry out our Bayesian analyses, we implemented all the models under investigation in the public HEPfit package [101], which performs a Markov Chain Monte Carlo (MCMC) analysis by means of the Bayesian Analysis Toolkit (BAT) [102]. This framework implements a Metropolis-Hastings algorithm, with the MCMC runs involving 20 chains with a total of 10 8 events per run, collected after all chains have achieved convergence in an adequate number of pre-run iterations.
Before performing the runs, we assigned a theoretical prior to each of the model parameters entering the fits. In particular, tan β is scanned logarithmically, in a similar fashion to what done in Ref. [29], in the broad range 0.01-10 3 . For the axion mass we also performed a logarithmic scan, choosing the range 0.1-10 3 meV which corresponds to the mass region of phenomenological interest. The mixing correction to the electron coupling L , which is not determined by the theory and in the 2+1 models represents an additional free parameter, is assumed to be distributed flatly in the physical range [−1, 1]. Similarly, the X 3 parameter of 3HDM is not constrained and is flatly distributed in the range [0, 1].
The Bayesian model comparison between the different models is performed evaluating for each of them an approximation of the Bayes factor, namely the Information Criterion (IC) [103]. 7 When comparing two models, the preferred one is the model displaying the smallest IC value [104]. Hence, after computing the IC value of the SM, we compute for each model M the score factor ∆IC ≡ IC SM − IC M , which is a quantity expressed in units of standard deviations. A larger value for this score thus signals a greater improvement in reproducing the data, compared to the SM case. The preferred models are thus the ones displaying the higher values of ∆IC.
The results of the global fits are reported in Table 4. Remembering that the constrains stemming from data relative to SN 1987A and to various NS are less sound than the ones extracted from WD, RGB and HB, and considering that, as we will see, they do affect in a strong manner the outcome of the fit, we report the results obtained for each model in three different analyses: in the first, less inclusive one, only data from WD, RGB and HB are included; in the second, more inclusive one, we fitted also the data pertaining to SN 1987A; and in the last, the most inclusive one, both data sets stemming from SN 1987A and NS are added to the likelihood. For each model and for each set of observables included in the fit we list the values of its parameters at the global mode, corresponding to the point where the multi-dimensional posterior p.d.f. is maximal. Only for the less inclusive fit for DFSZ1, for which the global mode of a parameter would fall outside the range of perturbative unitary discussed in Sec. 4, we imposed as an additional prior also a unitarity constraint. This procedure is however employed solely for the sake of reporting in Table 4 a value at the global mode within the unitarity limits, while all the remaining information pertaining to the model are inferred without imposing it. Finally, for each fit we give in the last column the value of the score factor ∆IC.
Starting from the canonical DFSZ models, we observe that the results are qualitatively similar to the ones found in the previous global analysis of Ref. [10]. Indeed, when considering only data from WD, RGB and HB, we obtain a better agreement with data compared to the SM, with the global mode for the axion mass found at m a ∼ 70-80 meV. On the other hand, including also data from SN 1987A and from the NS has the well understood effect of strongly reducing the allowed mass range for the axion, with the global mode now sitting at m a ∼ 10 meV. Given the fact that these additional measurements only put an upper bound on the axion couplings to nucleons, without hinting to new physics, the SM without axions is already in compliance with them. This is the reason why the overall preference for DFSZ models (or for any kind of axion model) over the SM is reduced when considering all experimental information at hand. Nevertheless, due to the presence of cooling anomalies, the preference for DFSZ models over the SM persists even in the most inclusive case, with DFSZ1 performing slightly better than DFSZ2.
Moving now to the non-universal cases, we start our discussion from the 2+1 models. Among all the models that we have considered these are the ones that display the better agreement with data. Thanks to their nucleo-phobic character, this feature is maintained  even when observations stemming from SN 1987A and the NS are included in the fit. As a consequence, the global modes of the axion mass, which in the less inclusive fit for the M2 model was found to be as high as m a ∼ 123 meV, do not drop below m a ∼ 50 meV even in the most inclusive cases. An interesting feature that can be inferred from Table 3 is that the models M1, M2 and M4 share many similarities: the global modes of the lepton mixing parameter are always found at L ∼ 0.75, while the inferred value for the M3 model is L ∼ 0.27; the global mode of tan β follows an analogous pattern, sitting at tan β ∼ 1.5 for the former models and at tan β ∼ 0.5 for the latter. Moreover, it is also evident that when adding SN and NS observables the global mode approaches the nucleo-phobic point tan β = √ 2 (for M2, M3 and M4) and tan β = 1/ √ 2 for M3. The fits for the 1+1+1 models T (u,d) 2 share some similarities with the universal DFSZ models, with the global mode for the axion mass ranging from m a ∼ 65 meV in the less inclusive fits to m a ∼ 10 meV in the more inclusive ones. However, given the nucleo-phobic nature of these models their overall agreement with data is improved compared to the universal models, hence suggesting a preference of T (u,d) 2 over DFSZ1 and DFSZ2, even if less strong than for the 2+1 models.
The last model analysed in our study is the 3HDM. In the corresponding fits the global mode for the axion mass is particularly stable, being found at m a ∼ 103 meV in the less inclusive fit and at m a ∼ 94 meV in the most inclusive one. This behaviour is due to the astro-phobic nature of the model: that is the nucleo-phobic conditions can be satisfied while simultaneously suppressing the axion-electron couplings (the global mode for the electron charge X 3 is always found at X 3 ∼ 0.04 well within the electro-phobic regime). For this reason the results of the fit are not strongly affected by increasing the number of observables included in the likelihood. Similarly to T (u,d) 2 , the 3HDM model is preferred by data over DFSZ1 and DFSZ2, but again not as strongly as the 2+1 models.
Complementary to Table 4, we also give in Fig. 2 the 2σ bounds extracted from the fits in the g ae vs. g aγ plane. For each model we show those regions for 3 different cases, according to the different set of observables included in the analysis. We recall that, contrarily to the model-independent case studied in Sec. 3.7 where g ae and g aγ are considered as independent parameters, here these two couplings depend on the same model parameters and hence are strongly correlated. This is the reason why the panels in Fig. 2 are somewhat different compared to the plot shown in Fig. 1, particularly concerning g aγ and when observables constraining the couplings to the nucleons are included in the fit. Indeed, as we have discussed above, the data relative to SN 1987A and the NS generally lower the upper bound on m a , which in turn decreases the upper bound for the coupling g aγ . This effect is more evident for the universal DFSZ models, but less accentuated for the non-universal ones. It is also worth recalling that, while in non universal 2+1 and 1+1+1 models the nucleo-phobic conditions are obtained only for a specific value of the free parameter tan β, they are always realised in the 3HDM model. This is the reason why the results of the latter model remain stable against the inclusion in the fit of the SN 1987A and of the NS data. As a closing note, it is intriguing to observe that the values of the axion mass which better fit the astrophysical data reported in Table 4, fall in the range 0.2-100 meV that is suggested, in post-inflationary scenarios, by cosmological data [82].

Discovery potential of meV-scale axion experiments
In this Section we discuss the perspectives to access experimentally the parameter regions hinted by the cooling anomalies. Fig. 2 shows that all models require a finite coupling to electrons to explain the observed stellar behavior. It is tempting, therefore, to look at experiments sensitive to the axion-electron coupling to test the preferred regions for this parameter. The most sensitive experiments of this kind are the large underground detectors XENON1T [105], LUX [106] and PandaX-II [107]. All of them can search for solar axions converted in the detector . Black lines correspond to the less inclusive global fit with only WD, RGB and HB data; red lines correspond to the more inclusive fit with data from SN 1987A added; blue lines are relative to the most inclusive case in which also data from NS are taken into account. through the axio-electric effect [108]. This strategy, however, has so far allowed to probe only relatively large axion-electron couplings, in a region in tension with RGB and WD observation (see, e.g., Ref. [41]). Even the next generation of underground detectors, such as Darwin [109], will not have a sufficient sensitivity to probe the couplings favored by stellar evolution (see, e.g., Ref. [110,111]). As we shall see, probing the axion-photon coupling offers a more efficient way to dig into the parameter region preferred by stars. Indeed, Tab. 4 shows that all our representative models give a maximal agreement with the stellar observations for axion masses in a range from a few meV to ∼ 100 meV, a mass range accessible to the next generation of axion helioscopes. As the name indicates, helioscopes [112] search for axions produced in the Sun. In the Sun core the main axion production channels are the Primakoff and the ABC (atomic transitions, bremsstrahlung and Compton) processes. The solar axion flux on earth is (see, e.g., Ref. [41]), We notice that this flux gets a similar contribution from Primakoff and ABC for values of the couplings of phenomenological interest. The exact weight of the two contributions depends on the specific axion model. There are additional contributions to the solar flux induced by the axion-nucleon coupling. These are, however, normally peaked at energies too large for standard helioscopes and will be ignored here. 8 Helioscopes exploit a strong laboratory magnetic field to convert solar axions into X-ray photons. The conversion probability is where q = q γ −q a is the momentum transfer provided by the magnetic field and L is the length of the magnet. In vacuum, q m 2 a /2ω. Coherence is ensured when qL 1. Therefore, the helioscope sensitivity is practically mass-independent up to a certain mass threshold, m th , above which it drops rapidly. The mass threshold depends on the specific helioscope. From the above expression, we find m th ≈ 10 meV L −1/2 10 , where L 10 is the magnet length in units of 10 m and we are using 3 keV as a reference solar axion energy. Whenever the coherence condition is realised, the conversion probability scales as (g aγ BL) 2 , rapidly increasing with the magnetic field and with the size of the magnet. Above m th , the sensitivity is reduced proportionally to m −2 a . The sensitivity may be regained using a buffer gas in the magnet beam pipes [115]. In this case the momentum transfer is q (m 2 a − m 2 γ )/2ω. Tuning the effective photon mass, m γ , to the axion mass allows to effectively regain coherence.
Below, we present a detailed study of the potential of the next generation of axion helioscopes to probe, in the regions of astrophysical interest, the various models that we have been discussing. In particular, we consider the proposed International Axion Observatory (IAXO) [116][117][118][119], and its scaled versions BabyIAXO and IAXO+. Of particular interest is 8  BabyIAXO, a scaled down (and significantly less expensive) version of IAXO, which will likely start operations in the mid of the current decade at the Deutsches Elektronen-Synchrotron (DESY) [120,121]. The beginning of operations for IAXO are presently not known. However, even this larger helioscope does not present particular technical challenges, besides better components [119] and may become operational in a not so distant future. Finally, IAXO+ represents a more aggressive version of the IAXO helioscope, with a larger magnet and better optics. The instrumental characteristics used in this work are extracted from Ref. [119].
The helioscope potential to detect the axion models discussed in the text is shown in the six panels of Fig. 3. In all cases, we show also the regions preferred by stellar evolution, making different assumptions for the astrophysical observables. The areas corresponding to the analysis of WDs, RGB and the R parameter (black contours) are the most solid, since the physics of these stars is relatively well known. The regions drawn by including also SN 1987A and data from various NS, for which the physics of axion emission is still not completely understood, are shown with separate curves (red and blue contours respectively).
In all models, the preferred parameter region appears as a relatively narrow mass band, spanning from a few meV for large electron couplings, to a few 10 meV in the opposite limit. The strong dependence of the mass boundaries on the axion-electron coupling should be clear from our previous discussion: as argued in Sec. 3.7 (see, in particular, Fig. 1) increasing the axion-electron coupling requires a smaller axion-photon coupling (and thus a smaller mass) to preserve the consistency with the observations. The axion-coupling with nucleons comes into the game when the SN or the SN+NS observables are also included. In the case of nucleo-phobic models their role is obviously reduced. In general, however, they play a more important role at low values of the axion-electron coupling, when the weight of WD and RGB data is decreased.
It is clear from the figures that the helioscopes may have a chance to detect the high mass end of the region favoured by stellar evolution. Not surprisingly, the models most accessible to the helioscopes are those with a large coupling to photons, particularly T (u,d) 2 (cf. Table  3). However, in general, the helioscope sensitivity increases (covering smaller masses) when the axion coupling to electrons increases, since this guarantees a larger solar flux through the ABC production processes. This feature is discernible for all the models.
A rather interesting outcome of our analysis is the potential to explore the astrophysically interesting regions of the axion parameter space for some models, and in particular for T (u,d) 2 , already with BabyIAXO. The upgraded versions, IAXO and especially IAXO+, will be clearly able to explore much larger sections of the parameter space.
Besides helioscopes, other experiments may also probe these regions of the axion parameter space. Normally, however, such experiments rely on a set of additional assumptions, for example that axions are a substantial fraction of the local DM density. In this case, preliminary studies [122,123] show that a new generation of haloscope detectors, based on Axionic Topological Antiferromagnets, might have the potential to explore the mass region of a few meV ( 10 meV), possibly down to the axion photon couplings expected in some realistic axion models (cf. Fig. 23 in Ref. [123]). Another experiment with the potential to probe the axion mass up to several meV is ARIADNE [124]. This is a long range force experiment, which relies on the axion coupling to nucleons. The experiment sensitivity, however, depends on the CP-violating axion-scalar coupling, whose phenomenologically allowed value spans several orders of magnitude [125][126][127]. Assuming the maximal experimentally allowed CP-violation beyond the SM, ARIADNE could potentially probe the DFSZ axion up to masses of several meV (cf. also Fig. 2 in Ref. [10]).  Figure 3. Contours of the 2σ regions from global fits to stellar cooling in the m a vs. C ae plane for different models, following the colour scheme defined under Fig. 2 (continuous lines). The colored regions show the helioscope sensitivities of BabyIAXO (in green), IAXO (in brown) and IAXO+ (in gray).

Summary and conclusions
In this paper we have revisited the global analysis of the stellar observations in relation to the axion couplings to the SM fields, providing an update with respect to Ref. [10] that includes the most recent observational results, and enlarging considerably the collection of axion models that are confronted with the data. The set of observations used in the present work is discussed in Sec. 3 while the statistical methodology is described in Sec. 5.1. As regards the observational information, the most significant update has been the inclusion of two recent analyses of the RGB bound on the axion-electron coupling, Refs. [13,14]. These two studies have revised the previous bound and have refined the analysis taking advantage of the new Gaia Data Release 2 [128], which has significantly improved over the previous knowledge of the distance determinations. Furthermore, we have also included new analyses of the SN [15] and NS [17] bounds on the axion-nucleon couplings.
Our study confirms a preference for non-vanishing axion couplings to electrons and photons (cf. Fig. 1). The hint is particularly strong for the axion-electron coupling, for which the best fit value lies away from zero with a significance of about ∼ 3 σ. Our best fit values for the axion-electron and axion-photon couplings are g ae 1.2×10 −13 , g aγ 1.8×10 −11 GeV −1 , which are in good agreement with previous results, although, because of the new RGB analyses, the g ae best fit point is slightly shifted to a lower value.
To investigate the theoretical implications of these hints, we have studied a set of well motivated QCD axion models whose relevant features have been summarised in Sec. 4. On top of the two benchmark DFSZ axion models, we have considered a general class of nonuniversal DFSZ-like models, featuring generation-dependent PQ charges. The logic behind the models' selection was to have in first place a sizeable C ae /C aγ ratio, a favorable condition that was highlighted by the model-independent fit in Fig. 1, and next to allow for the possibility of suppressing sufficiently the axion-nucleon couplings. This latter condition allows to circumvent to some extent the strong bounds from SN 1987A and from NS. The results of our analysis for the different models are presented in Fig. 3. The qualitative feature is the same for all models: for large electron couplings, the axion photon coupling is pushed towards relatively small values (cf . Fig 1), and correspondingly the preferred region for the axion mass is also driven to relatively small values m a ∼ 1 meV. In the case of small axionelectron coupling, on the other hand, there is a preference for a more sizeable axion-photon coupling and thus for a larger axion mass, up to 100 meV or so, depending on the model. In all cases, the inclusion of constraints on the axion-nucleon couplings from SN 1987A or from the SN and the NS pushes the mass to lower values, since these observations prefer a vanishing axion-nucleon coupling (and hence a large value of the axion decay constant f a which suppresses the axion mass). In the case of nucleo-phobic models however, axion-nucleon decoupling is assisted by a conspiracy between the values of the PQ charges of the quarks, and accordingly the mass suppression effect is considerably less important.
An important outcome of our study is that star evolution observations indicate a clear preference for an axion mass in the meV region, a range that is at least partially accessible to the next generation of axion helioscopes. The discovery potential of helioscopes is quite sensitive to the axion-photon coupling, which is largest in the T (u,d) 2 models (the complete list of axion-photon couplings for all non-universal DFSZ axion models with two Higgs doublets is given in Tab. 5). It is quite remarkable that already BabyIAXO will be able to probe some of these models, and that IAXO and IAXO+ will have a chance to cover the hinted region entirely.
It should be also stressed that a direct exploration of this region of parameter space with dedicated axion experiments will have an impact not only for axion physics but also for astrophysics. Although the stellar hints individually do not have a large significance, they all show a preference for an increased efficiency in star energy loss. Assuming no new physics is at play, this would signal a systematic problem in our understanding of stellar cooling, a possibility that would gain strength if some preferred particle physics explanations could be ruled out. Certainly, the improvements in the astrophysical instrumentation expected in the next decade or so will have a strong impact for clarifying some of these issues, and it might eventually strengthen the case for new energy loss channels. In such scenario, a new physics solution, perhaps in the form of axions, would be a most exciting result. In any case, whether this problem requires new physics, or is just a matter of understanding better the details of star evolution, or is merely an instrument calibration issue, will be clarified in the (hopefully) near future. In the meanwhile, we keep being fascinated by the synergy between astrophysics and particle physics, and by the possibility that the first evidences for the need of new particle physics could eventually come from new careful observations of the sky. the non-vanishing product of the three matrix elements it follows where in the last step we assumed only two Higgs doublets, and we have defined as H 1 the doublet that occurs more times in the product of the three mass matrix entries of the nonvanishing term of the determinant. Doing the same for the down-type quarks and leptons is straightforward. However, now we have to keep the distinct possibilities for X 1 ↔ X 2 . So in the two Higgs doublet case we get Considering together Eq. (A.1) and Eq. (A.2) we see that we have 2 × 4 different possibility for the QCD anomaly (some will give the same result). For the numerator of the second term of our E/N expression, that reads i X q i − i X u i + i X i − i X e i , one gets the same 2 × 4 possibilities, since Eq. (A.2) holds also for the leptons. Then, by combining in all the ways the possible numerators and denominators, one gets corresponding to the results displayed in Table 5.