Stellar cooling anomalies and variant axion models

TL;DR

This paper revisits stellar cooling anomalies suggesting axions as an energy loss mechanism, proposes variant axion models with specific couplings, and identifies parameter regions testable by future experiments.

Contribution

It introduces variant axion models with flavor-dependent charges that explain cooling anomalies and address cosmological issues, expanding viable parameter space.

Findings

Weakening of the cooling anomaly significance to 2.4σ

Identification of axion mass range 0.45 meV to 30 meV consistent with data

Predictions for axion couplings accessible to upcoming helioscope experiments

Abstract

A number of observations of stellar systems show a mild preference for anomalously fast cooling compared with what predicted in the standard theory, which leads to a speculation that there exists an additional energy loss mechanism originated from the emission of axions in stars. We revisit the global analysis of the stellar cooling anomalies by adopting conservative assessments on several systematic uncertainties and find that the significance of the cooling hints becomes weaker but still indicates a non-vanishing axion-electron coupling at around 2.4$\,\sigma$. With the revised analysis results, we explore the possibility that such excessive energy losses are interpreted in the framework of variant axion models, which require two Higgs doublets and flavor-dependent Peccei-Quinn charge assignments. These models resolve two fundamental issues faced in the traditional KSVZ/DFSZ models by predicting a sizable axion coupling to electrons required to explain the cooling anomalies and at the same time providing a solution to the cosmological domain wall problem. We also find that a specific structure of the axion couplings to electrons and nucleons slightly relaxes the constraint from supernova 1987A and enlarges viable parameter regions compared with the DFSZ models. It is shown that good global fits to the observational data are obtained for axion mass ranges of $0.45\,\mathrm{meV} \lesssim m_a \lesssim 30\,\mathrm{meV}$, and that the predicted parameter regions can be probed in the forthcoming helioscope searches.

Figures

Tables

Table 1. Table 1: PQ charge assignments for the quark fields.

Table 2. Table 2: PQ charge assignments for the lepton fields.

Table 3. Table 3: Axion couplings to up-type ( u 𝑢 u ) and down-type ( d 𝑑 d ) quarks in the variant axion models.

(3.25)

Table 4. Table 4: Axion couplings to electrons and photons and the model-dependent coefficient for the axion-photon coupling in the variant axion models.

(3.31)

Table 5. Table 5: Best fit parameters and χ min 2 / d.o.f. subscript superscript 𝜒 2 min d.o.f. \chi^{2}_{\rm min}/\text{d.o.f.} for the interpretation of the stellar cooling anomalies in the variant axion models with general quark mixing. The column of ( V u , d ) 11 subscript subscript 𝑉 𝑢 𝑑 11 (V_{u,d})_{11} shows the value of ( V u ) 11 subscript subscript 𝑉 𝑢 11 (V_{u})_{11} for Model U, C, T and that of ( V d ) 11 subscript subscript 𝑉 𝑑 11 (V_{d})_{11} for Model D, S, B.

(4.10)

Table 6. Table 6: The coefficient ζ 𝜁 \zeta for the estimation of experimental sensitivities and minimal axion mass ranges accessible to IAXO and IAXO+ in the variant axion models.

(4.18)

Table 7. Table 7: 2  σ 𝜎 \sigma ranges in the axion mass and couplings to electrons, photons, neutrons, and protons, hinted by WD, HB, and RGB cooling anomalies and compatible with the bound from SN 1987A and the requirement of perturbativity of Yukawa interactions for different cases of the variant axion models.

(4.28)

Table 8. Table 8: Best fit parameters and χ min 2 / d.o.f. subscript superscript 𝜒 2 min d.o.f. \chi^{2}_{\rm min}/\text{d.o.f.} for the interpretation of the stellar cooling anomalies in Model U and D with the assumption of ( V u , d ) 11 = 0 subscript subscript 𝑉 𝑢 𝑑 11 0 (V_{u,d})_{11}=0 and in DFSZ models.

(4.37)

Equations

$${\mathcal{L}}_{a\gamma }=-\frac{1}{4}{g}_{a\gamma }a{F}_{\mu \nu }{\tilde{F}}^{\mu \nu },g_{a\gamma }\equiv \frac{\alpha }{2\pi f_a}{C}_{a\gamma },$$

$${\mathcal{L}}_{ae}=-i{g}_{ae}a\overline{e}{\gamma }_{5}e,g_{ae}\equiv \frac{C_{ae}{m}_e}{f_a},$$

$$R^{\mathrm{th}}=7.33Y+0.02-0.095\sqrt{21.86+21.08g_{10}}-1.61\delta {\mathcal{M}}_c-0.067{\alpha }_{26},$$

$$\delta {\mathcal{M}}_c=0.024M_{\odot }\left(\sqrt{4\pi {\alpha }_{26}+(1.23{)}^2}-1.23-0.921{\alpha }_{26}^{0.75}\right)$$

$$M_{I,\mathrm{TRGB}}^{\mathrm{th}}=-4.03-0.25\left((g_{13}^2+0.93{)}^{0.5}-0.96-0.17g_{13}^{1.5}\right)+\delta M_{I,\mathrm{TRGB}}^{\mathrm{th}}\mathrm{mag},$$

$${\mathcal{L}}_{aN}=-i{g}_{aN}a\overline{N}{\gamma }_{5}N,g_{aN}\equiv \frac{C_{aN}{m}_N}{f_a},$$

$$g_{ap}^2+g_{an}^2<3.6\times 10^{-19}.$$

$$H_k\to e^{i{X}_kϵ}{H}_k,\sigma \to e^{i{X}_{\sigma }ϵ}\sigma ,$$

$$\mathcal{L}\supset -\frac{{\alpha }_s}{8\pi }\mathcal{N}\frac{a}{v_{\mathrm{PQ}}}{G}_{\mu \nu }^a{\tilde{G}}^{a\mu \nu },$$

$$\mathcal{N}=X_1-X_2\text{or}{X}_2-X_1,$$

$$-{\mathcal{L}}_{\mathrm{Yukawa},q}$$

$$+{\Gamma }_i^d{\overline{q}}_{iL}{H}_d{d}_R+{\Gamma }_i^s{\overline{q}}_{iL}{H}_s{s}_R+{\Gamma }_i^b{\overline{q}}_{iL}{H}_b{b}_R+h.c.,$$

$$-\mathcal{L}_{\mathrm{Yukawa},\ell}=\left\{\begin{array}[]{ll}\Gamma^{\ell}_{ij}\overline{\ell}_{iL}H_{1}e_{jR}+\mathrm{h.c.}&\text{(Type I)},\vspace{1mm}\\ \Gamma^{\ell}_{ij}\overline{\ell}_{iL}H_{2}e_{jR}+\mathrm{h.c.}&\text{(Type II)},\end{array}\right.$$

$${\partial }^{\mu }{j}_{\mu }^{\mathrm{PQ}}=\mathcal{N}\frac{{\alpha }_s}{8\pi }{G}_{\mu \nu }^a{\tilde{G}}^{a\mu \nu }+\mathcal{E}\frac{\alpha }{8\pi }{F}_{\mu \nu }{\tilde{F}}^{\mu \nu },$$

$$j_{\mu }^{\mathrm{PQ}}=v_{\sigma }\partial \tilde{a}-\sum _{\psi }\sum _i\left[{\overline{\psi }}_{iL}{X}_{{\psi }_{iL}}{\gamma }^{\mu }{\psi }_{iL}+{\overline{\psi }}_{iR}{X}_{{\psi }_{iR}}{\gamma }^{\mu }{\psi }_{iR}\right],$$

$$j_{\mu }^{\mathrm{PQ}}$$

$$X_{{\psi }_{iL,R}}^{\prime }$$

$$\tan \beta \equiv \frac{v_2}{v_1}.$$

$$\mathcal{L}$$

$$-\frac{{\partial }_{\mu }a}{f_a}\sum _{\psi }\sum _i\left[{\overline{\psi }}_{iL}\frac{X_{{\psi }_{iL}}^{\prime }}{\mathcal{N}}{\gamma }^{\mu }{\psi }_{iL}+{\overline{\psi }}_{iR}\frac{X_{{\psi }_{iR}}^{\prime }}{\mathcal{N}}{\gamma }^{\mu }{\psi }_{iR}\right],$$

$$f_a\equiv \frac{v_{\mathrm{PQ}}}{\mathcal{N}}$$

$${\mathcal{L}}_{a\psi }=-\frac{{\partial }_{\mu }a}{2f_a}\sum _{\psi }\sum _{i,j}{\overline{\psi }}_i{\gamma }^{\mu }\left[{\left(C_{a\psi }^V\right)}_{ij}-{\left(C_{a\psi }^A\right)}_{ij}{\gamma }_5\right]{\psi }_j,$$

$$C_{a\psi }^V$$

$$C_{a\psi }^A$$

$$V_u$$

$$V_d$$

$$\Lambda$$

$$q^{\prime}=\left(\begin{array}[]{c}u\\ d\end{array}\right)\to e^{-i\gamma_{5}\frac{a}{2f_{a}}Q_{a}}\left(\begin{array}[]{c}u\\ d\end{array}\right),$$

$$\mathcal{L}\supset \frac{1}{2}\frac{{\partial }_{\mu }a}{f_a}{\bar{q}}^{\prime }{\gamma }^{\mu }{\gamma }_5(C_{a{q}^{\prime }}-Q_a)q^{\prime }-({\bar{q}}_L^{\prime }{M}_a{q}_R^{\prime }+h.c.),$$

$$\displaystyle M_{a}=e^{i\frac{a}{2f_{a}}Q_{a}}M_{q}e^{i\frac{a}{2f_{a}}Q_{a}},\quad M_{q}=\left(\begin{array}[]{cc}m_{u}&0\\ 0&m_{d}\end{array}\right),\quad C_{aq^{\prime}}=\left(\begin{array}[]{cc}(C_{au}^{A})_{11}&0\\ 0&(C_{ad}^{A})_{11}\end{array}\right).$$

$${\mathcal{L}}_{\mathrm{eff}}\supset \frac{1}{2}{f}_{\pi }^2\mu \mathrm{Tr}\left(U{M}_a^{†}+M_a{U}^{†}\right),$$

$$\displaystyle U=e^{i\Pi/f_{\pi}},\quad\Pi=\left(\begin{array}[]{cc}\pi^{0}&\sqrt{2}\pi^{+}\\ \sqrt{2}\pi^{-}&-\pi^{0}\end{array}\right),$$

$$m_a^2=\frac{m_u{m}_d}{(m_u+m_d{)}^2}\frac{m_{\pi }^2{f}_{\pi }^2}{f_a^2}.$$

Full text

1. 1 Introduction

2. 2 Axions and stellar energy losses

3. 2.1 Horizontal and red giant branch stars in globular clusters

4. 2.2 Tip of the red giant branch

5. 2.3 White dwarf luminosity function

6. 2.4 White dwarf variables

7. 2.5 Supernova 1987A and neutron stras

8. 2.6 Implications for axion models

9. 3 Variant axion models

10. 4 Interpretation of stellar cooling anomalies

11. 4.1 Models with general quark mixings

12. 4.2 Models with small quark mixings

13. 5 Discussion and conclusions

14. A Analysis method

15. B Higgs sector and Yukawa interactions

1 Introduction

Extreme environment realized in stars can be regarded as a good laboratory to test fundamental physics [1]. The consideration of stellar evolution particularly leads to quite strong constraints on light weakly interacting particles. A leading example of such a low mass particle is the axion [2, 3], which is a Nambu-Goldstone boson emerging from the spontaneous breaking of hypothetical global U(1) Peccei-Quinn (PQ) symmetry introduced to provide a solution to the strong CP problem of quantum chromodynamics (QCD) [4]. If the scale of the PQ symmetry breaking is sufficiently high, interactions of axions with ordinary matter are quite weak. This fact implies that axions produced in the core of a star can easily escape from the system, providing a new energy loss mechanism. Such an exotic energy loss can affect the evolution of stars, which allows us to investigate the property of axions by using the results from astrophysical observations.

Recently, it has been pointed out that observations of stars in several different evolutionary stages show systematic trends to prefer excessive energy losses over the standard cooling scenario, and that such cooling anomalies can be explained just by adding new cooling channels due to axion emission [5, 6, 7, 8, 9, 10]. Understanding such astrophysical phenomena is not straightforward, of course, and the results should be taken carefully. However, at this point we cannot discard the possibility of interpreting them as hints of the existence of axions. According to the analysis in Ref. [10], there is a more than 3 $\sigma$ preference for the axion interpretation.111The significance of the hints becomes weaker if we adopt a more conservative approach to include systematic uncertainties (see Sec. 2.6).

The hints obtained from the stellar cooling anomalies can be used as a guide for next generation axion search experiments [see e.g. Ref. [11] for reviews]. In particular, the proposed helioscope, International Axion Observatory (IAXO) [12, 13], will probe a broad parameter space with improved sensitivities and have a potential to cover a large fraction of the region motivated by the solution to the stellar cooling anomalies [14].

It should be noted that the preceding argument on the axion interpretation of the stellar cooling anomalies might remain inadequate when we consider cosmological issues of the axion models. It is known that the axion models suffer from a constraint from isocurvature fluctuations in the cosmic microwave background [15, 16], which arises if the PQ symmetry is broken before or during inflation and never restored afterwards. This constraint becomes severe unless the inflationary energy scale is sufficiently low. For instance, if the scale of the PQ symmetry breaking is of order $\sim 10^{9}\,\mathrm{GeV}$, which is preferred by the stellar cooling hints [10], the Hubble scale during inflation must be lower than $H_{\rm inf}\lesssim 10^{6}\,\mathrm{GeV}$ in order to avoid large isocurvature fluctuations [17].222The isocurvature constraint also depends on the abundance of cold axions produced by the vacuum re-alignment mechanism [18, 19, 20], and the limit $H_{\rm inf}\lesssim 10^{6}\,\mathrm{GeV}$ corresponds to the case where the fraction of the energy density of cold axions is about 1 % of the total dark matter abundance. The constraint becomes tighter if the axion fraction gets larger. A simple way to avoid this isocurvature constraint is to assume that the PQ symmetry is broken after inflation. However, in this case domain walls are created around the epoch of the QCD phase transition [21], leading to disastrous consequences if they are stable [22]. This cosmological domain wall problem poses a serious concern for the popular Dine-Fischler-Srednicki-Zhitnitsky (DFSZ) models [23, 24], which provide good fits to the stellar cooling hints [10]. On the other hand, in the Kim-Shifman-Vainshtein-Zakharov (KSVZ) models [25, 26], which are regarded as another well-motivated class of axion models, the axion interaction with electrons is too small to explain the observed data [10], though they can avoid the domain wall problem straightforwardly. This situation is alleviated in a KSVZ-like axion/majoron model [27, 28, 29], but in this case a sufficiently large radiative correction to the axion-electron coupling is required, which is likely to violate the condition of perturbativity.

In this paper, we point out that both of the above astrophysical and cosmological issues can be naturally addressed in the framework of variant axion models, which were proposed first by Peccei, Wu and Yanagida [30] and by Krauss and Wilczek [31]. These models are constructed based on two Higgs doublets and flavor-dependent PQ charge assignments for the quark fields. Although the first version of the variant axion models was excluded as it is “visible” like the original Weinberg-Wilczek model [2, 3] (i.e. the PQ scale is assumed to be as low as the electroweak scale), it is possible to make it “invisible” [32, 33] by adding a singlet complex scalar field such that the PQ scale becomes arbitrary high like all other viable axion models.333Extensions of the variant axion models were proposed recently in e.g. Refs. [34, 35, 36, 37], which are characterized by flavor-changing axion couplings. A potential problem for such models is that they generically suffer from a severe constraint coming from the non-observation of rare decays $K^{+}\to\pi^{+}+a$. We point out that in the simplest framework of the variant axion models this constraint can be straightforwardly avoided if only one of the up-type quarks has a nonzero PQ charge. The important feature of this framework is that we can arrange the model such that it leads to a unique vacuum in the low energy effective theory, which evades the domain wall problem [32]. Furthermore, it can make the axion-electron coupling naturally large as it appears in the tree level, in a similar manner to the DFSZ framework. This feature is contrasted to the axion/majoron model, where the electron coupling is induced only radiatively.

The remaining part of the paper is organized as follows. In Sec. 2, we review the observational results of cooling anomalies and previous attempts to interpret them in terms of additional cooling channels induced by axions. The global fits to the astrophysical data are revisited with slightly different assessments on the systematic uncertainties. In Sec. 3, we construct the variant axion models and derive axion couplings to the Standard Model (SM) particles. In Sec. 4, we show the results of global fits to the data of stellar cooling anomalies, specifying the preferred parameter regions. We also compare these predictions with the sensitivities of future helioscope searches such as IAXO. Section 5 is devoted to discussion and conclusions. In Appendix A, we summarize the method of the global analysis. In Appendix B, we describe the structure of the Higgs sector and Yukawa interactions in the variant axion models in order to derive typical parameter ranges compatible with requirements of perturbativity of Yukawa interactions.

2 Axions and stellar energy losses

Low mass weakly interacting particles such as axions can affect the time scale of the evolution of stars, since they can escape from the system and take away a considerable amount of energy. Since the existence of such exotic cooling channels changes various astrophysical observables, we can constrain the properties of hypothetical particles by comparing the observational results with the prediction of the standard cooling scenario [1]. Furthermore, results of some recent analyses even show a mild preference for the existence of excessive energy losses. In this section, we revisit the hints of excessive energy losses reported in previous literatures and their possible interpretations in terms of axions. We take a different approach to deal with several systematic uncertainties, shedding light on challenging aspects of the data analysis.

2.1 Horizontal and red giant branch stars in globular clusters

A strong indicator to measure the cooling rate of stars in different evolutionary stages is provided by the observation of globular clusters. Globular clusters are gravitationally bound systems composed of about $10^{6}$ stars, which show a characteristic distribution in the color-magnitude diagram according to a certain evolutionary phase. In particular, the ratio of the number of stars in the horizontal branch (HB) to that of those in the red giant branch (RGB), called the $R$-parameter, is used to estimate the relative time scale for different evolutionary stages associated with these branches. The existence of axions can affect the value of the $R$-parameter, since they can be efficiently produced via the Primakoff process $\gamma+Ze\to Ze+a$ in the HB stars but not in the RGB stars, changing the relative amount of time spent on two different branches. Hence the additional cooling due to the emission of axions in the HB stars could be observed as a decrease in the value of $R$.

Based on the above argument, we can constrain the magnitude of the axion-photon coupling $g_{a\gamma}$, which is given by

[TABLE]

where $F_{\mu\nu}$ is the electromagnetic (EM) field strength, $\tilde{F}^{\mu\nu}$ its dual, $f_{a}$ the axion decay constant, $\alpha$ the fine structure constant, and $C_{a\gamma}$ a dimensionless coefficient.

Axions can also interact with electrons through the following coupling,

[TABLE]

where $C_{ae}$ is a dimensionless coefficient and $m_{e}$ the electron mass. In this case, axions can be produced in the core of red giants through bremsstrahlung off electrons $e+Ze\to Ze+e+a$, which cools the core and causes a delay of helium ignition. The delay of helium ignition leads to a larger core mass, which affects the value of $R$-parameter. Therefore, the information on the $R$ can be used to probe not only the axion-photon coupling but also the axion-electron coupling [8].

The prediction for the effect of the axion emission on the $R$-parameter is given by [8]

[TABLE]

where $g_{10}=g_{a\gamma}\times 10^{10}\,\mathrm{GeV}$, $\alpha_{26}=10^{26}\times g_{ae}^{2}/4\pi$, $Y$ is the helium mass fraction in the globular clusters, and

[TABLE]

is the shift of the core mass due to the the axion emission during the RGB phase [38, 39].

In Ref. [40], observational value $R^{\rm obs}=1.39\pm 0.03$ was obtained and compared to a theoretical prediction to derive a bound on $g_{a\gamma}$. This analysis was revised in Ref. [8] by adding the dependence on $g_{ae}$ in Eq. (2.3). As the model depends not only on $g_{a\gamma}$ and $g_{ae}$ but also on the helium mass fraction $Y$, the results of such analyses rely on the assumption about the value of $Y$. The authors of Ref. [40] preferred to use a value of $Y=0.2535\pm 0.0036$, which corresponds to the weighted average for the helium abundance measured at different oxygen abundance (O/H) estimated in Ref. [41]. Note that this value is somewhat larger than the primordial helium abundance $Y_{p}$, which is obtained by extrapolating the observational results to $\mathrm{O/H}\to 0$. As noted in Ref. [40], such different assumptions on $Y$ lead to different bounds on the axion couplings. It was found that there is a mild preference for a non-vanishing value of $g_{a\gamma}$ when the average value $Y=0.2535\pm 0.0036$ is adopted, while such a hint is diminished for the primordial value $Y_{p}$.

The claim of Ref. [40] is that one should not use $Y_{p}$ extracted at zero metallicity but the actual helium abundance in the globular clusters, since the chemical composition there should be different from the primordial one. In this sense, the average value of $Y$ can be regarded as representative of the actual helium mass fraction in the globular clusters. In the analysis performed in this paper, we follow this suggestion and adopt $Y=0.2535\pm 0.0036$ in the global fits.

Equation (2.3) was derived in Ref. [8] based on the analysis methods developed in Refs. [40, 42]. In these analyses, the stellar evolution models were computed via a 1D hydrostatic code FUNS [43], which was developed based on the FRANEC code [44]. In Ref. [45], behavior of the intermediate mass stars in the pre-asymptotic giant branch phase computed by means of the FRANEC code was compared to other evolutionary codes, and it was found that there is a $40\textendash 50\%$ difference in the evaluation of the helium-burning lifetime between different models. Although we cannot directly apply this estimate to the present case of the FUNS code, it should be kept in mind that Eq. (2.3) could be subject to a similar amount of uncertainty.

The effect of stellar rotation is also a possible source of uncertainty in Eq. (2.3). Rotation could modify the structure of a star, which results in a change in its chemical composition. The effect of rotation was studied in Ref. [46] with a modified version of the FRANEC code. It was shown that a nonzero rotation velocity leads to a shorter helium-burning lifetime, which implies a smaller value of $R$. In the models considered in Ref. [46], the variation of the helium-burning lifetime due to the effect of rotation is about $3\%$ in low metallicity environments. In this paper, we include it as an additional systematic uncertainty of $R$ when we perform the global analysis.

2.2 Tip of the red giant branch

As mentioned in the previous subsection, the additional cooling due to the axion bremsstrahlung off electrons can lead to a delay of helium ignition in the core of red giants. As a result, the core becomes more massive and brighter at helium ignition, and such effects can be observed as a tip of the RGB (TRGB) in the color-magnitude diagram of a globular cluster. This observable can be used to constrain the axion-electron coupling $g_{ae}$.

In Refs. [38, 39], the color-magnitude diagram of the globular cluster M5 was studied, and the observational value of $M_{I,\mathrm{TRGB}}^{\rm obs}=-4.17\pm 0.13\,\mathrm{mag}$ for the TRGB brightness was obtained. This value was compared to the theoretical prediction given by [39]

[TABLE]

where $g_{13}=10^{13}\times g_{ae}$, and a shift $\delta M_{I,\mathrm{TRGB}}^{\rm th}=0.039\,\mathrm{mag}$ was introduced since the ranges of systematic errors in $M_{I,\mathrm{TRGB}}^{\rm th}$ were not symmetric. A number of possible sources of systematic uncertainties in the theoretical prediction $M_{I,\mathrm{TRGB}}^{\rm th}$ were investigated and summarized in Ref. [38]. Among them, relatively important ones are the effect of mass loss [$+(0.022\textendash 0.035)\,\mathrm{mag}$], the equation of state ($-0.0045$ to $+0.0242\,\mathrm{mag}$), the treatment of convection ($\pm 0.056\,\mathrm{mag}$), and the color transformations and bolometric corrections [$\pm(0.08+0.02g_{13})\,\mathrm{mag}$].

As in the case of the $R$-parameter described in Sec. 2.2, the stellar rotation may also affect the results of $M_{I,\mathrm{TRGB}}$. However, we can assume that such an effect is negligible in this case from the following reason. According to Ref. [46], the variation of the core mass due to the stellar rotation at the TRGB phase is just about $0.35\%$ in low metallicity environments. The analysis in Ref. [39] indicates that a change in the core mass by $3\%$ corresponds to the error of $\Delta M_{I,\mathrm{TRGB}}=\pm 0.002\,\mathrm{mag}$, which is much smaller than other dominant sources of uncertainties. Therefore, we expect that the stellar rotation would only have a minor effect in the estimation of $M_{I,\mathrm{TRGB}}$.

In Refs. [38, 39], all the systematic errors in $M_{I,\mathrm{TRGB}}^{\rm th}$ were convolved to a Gaussian statistical error, whose magnitude was estimated by adding them in quadrature. With this estimate of the theoretical error, the comparison with the observational value gives a bound $|g_{ae}|<4.3\times 10^{-13}$ (95% CL) with a small hint that the agreement between observations and theoretical predictions improves at a nonzero value of $g_{ae}$. However, the validity of this procedure of converting systematic uncertainties into statistical uncertainties is not obvious. In the analyses performed in this paper, we take a different approach to deal with the error in $M_{I,\mathrm{TRGB}}^{\rm th}$.

2.3 White dwarf luminosity function

The hint for non-vanishing value of $g_{ae}$ was also pointed out in the context of the cooling of white dwarfs (WDs). The cooling rate of WDs can be inferred from their luminosity function (WDLF), which represents the number density of WDs per luminosity interval. The dominant process of the axion emission is given by the bremsstrahlung off electrons, which hastens the cooling of WDs. This additional cooling leads to some modifications of the WDLF. The result of the analysis in Ref. [47] shows that a non-vanishing value of $|g_{ae}|\sim 1.4\times 10^{-13}$ is marginally preferred, although it excludes higher values of $g_{ae}$ and sets a bound $|g_{ae}|<2.1\times 10^{-13}$ (95% CL).

As noted in Ref. [47], there are two obstacles that prevent us from drawing credible conclusions from the analysis of the WDLF. One is the uncertainties in the stellar formation rates (SFR). In Ref. [47], theoretical WDLF curves were derived with an assumption of a constant SFR, while possible time variations of the SFR can affect the estimation of the theoretical WDLF. The other is magnitude-dependent selection effects (survey incompleteness), which could distort the observed WDLF.

In Ref. [48], it was claimed that the shape of the bright part of the WDLF is almost independent of the SFR and that its effect can be absorbed in the overall normalization. The analysis performed in Ref. [49] with the aim of the “inversion” from the WDLF to the SFR (i.e. deriving the time dependence of the SFR from the observed WDLF) appears to confirm this claim, as only a fainter part ($M_{\rm bol}\gtrsim 12$) of the WDLF is sensitive to the SFR. Since the bound on the axion-electron coupling was derived by only using the luminosity range $3\leq M_{\rm bol}\leq 12.5$ [47], we expect that the possible time variation in the SFR would have only a minor effect in the shape of the WDLF of interest.

The more important and serious issue is the uncertainties due to survey incompleteness. The data of the “observed” WDLF used in Ref. [47] were constructed in Ref. [50] by considering two data sets obtained from different surveys: One is the WDLF determined by the Sloan Digital Sky Survey (SDSS) [51, 52], and the other is that determined by the SuperCOSMOS Sky Survey (SSS) [53]. It was found that these two WDLFs do not agree with each other within their error bars. This inconsistency might be attributed to the incompleteness of the SSS-WDLF, which could be not uniform at all magnitude bins, biasing its shape [50].

Understanding the origin of the discrepancy is challenging, and further independent WDLFs are likely to be needed to deal with possible systematics. To compromise over this issue, in Ref. [50] two WDLFs are unified into an averaged WDLF, whose error bras are estimated from differences between the SDSS-WDLF and the SSS-WDLF. This procedure can be regarded as an intuitively reasonable guess for the systematic uncertainty, and we use it in the analysis performed in this paper. However, it should be noted that real systematic uncertainties can be larger than this estimate and probably highly correlated, and this situation must be reviewed once further observational data are available.

2.4 White dwarf variables

Another observable to infer the cooling rate of WDs is the period decrease of WD variables. The time scale of the change of the pulsation period of WD variables can be related to the cooling rate, and it can be used to probe the exotic cooling due to the axion emission [54]. The observations of several WD variables show that the observed rate of the period decrease is substantially faster than what predicted in the standard theory. The reported hints of the cooling excess include the observations of G117-B15A [55], R548 [56], PG1351+489 [57, 58], and L19-2 [59]. These results lead to a speculation that there exits some anomalous energy loss due to the axion-electron bremsstrahlung.

The rate of the period decrease of the pulsating WD depends on its mass and effective temperature. In the study of PG 1351+489 [57, 58], theoretical prediction for the rate of the period change was given by considering two approaches, the asteroseismological model, in which the mass and effective temperature are fixed such that it reproduces the observed pulsation periods, and the spectroscopic model, in which the mass and effective temperature are chosen to be spectroscopically inferred values. It was found that the values of the surface gravity and effective temperature found by the asteroseismological model are not consistent with those given by spectroscopic measurements. Due to this discrepancy, the bounds on the axion couplings become different according to which model is assumed [58].

The above model dependence can be regarded as a systematic uncertainty in the theoretical estimate of the rate of the period change of PG 1351+489, which we include in the analysis performed in this paper. Note that a similar caution should be given for other WD variables. In the study of L19-2 [59], the results of asteroseismological models agreed with spectroscopic inferences, and hence such a model uncertainty is not the issue in that case.

2.5 Supernova 1987A and neutron stras

Several observational results described above motivate us to study the possibility to interpret the cooling anomalies in terms of the axion emission. However, there are further caveats to be kept in mind. One is the fact that the axion couplings to nucleons can be constrained by considering the energy loss of the supernova (SN) 1987A [60, 61, 62]. The additional cooling channel due to nucleon bremsstrahlung $N+N\to N+N+a$ shortens the neutrino pulse duration from the SN, and the observed neutrino signals lead to the bounds on the axion-nucleon couplings $g_{aN}$, which are given by

[TABLE]

where $C_{aN}$ is a dimensionless coefficient, $m_{N}$ the mass of a nucleon, with $N=p,n$ representing either a proton or a neutron.

We note that there is subtlety in deriving the axion bound from SN 1987A observations. This is mainly related to the difficulty in evaluating the axion production rate in the dense nuclear medium from the first principles. It was pointed out that the axion emission can be suppressed by many-body and multiple-scattering effects [63, 64], while the degree of suppression is determined by the nucleon spin-fluctuation rate that is hard to estimate in the dense medium because of the complication of nuclear physics [1]. The results of the numerical studies [65, 66] using an axion emission rate with many-body effects estimated in Refs. [67, 68] imply the following bound [10]

[TABLE]

It should be emphasized that the above constraint should be taken as an indicative result rather than a sharp bound as the nature of the axion emission is not completely understood and SN simulations do not take account of all necessary physics. More detailed study is warranted to obtain a conclusive bound on the axion-nucleon couplings from SN observations. Recently, the SN 1987A bound has been reviewed in Ref. [69] by applying several updated nuclear physics calculations, which yields a factor $\sim 5$ weaker limit on $f_{a}$ compared with the “canonical” bound for the KSVZ models $f_{a}\gtrsim 4\times 10^{8}\,\mathrm{GeV}$ from Ref. [70]. A more recent study performed in Ref. [71] shows a similar relaxation of the SN 1987A bound, but the degree of reduction is milder than Ref. [69], amounting to a factor $\sim 2$ weaker limit on $f_{a}$ compared with the canonical bound. The discrepancy between two results may originate from the fact that the corrections from different nuclear physics effects are treated separately as multiplicative factors in Ref. [69] whereas in Ref. [71] these effects are combined in a consistent manner. However, these results are still not conclusive, as they do not take account of the feedback from the axion emission on the SN. The bound would be further modified if the effect of the axion feedback is consistently included in a full SN simulation [71].

Constraints on the axion-nucleon couplings also arise from the cooling of neutron stars (NSs). In NSs the axion emission via Cooper pair-breaking-formation (PBF) process could be important in addition to the nucleon bremsstrahlung process, if baryons in the core of the NS undergo a transition to superfluid state [72]. Simulations of NS cooling by including the axion emission via PBF process and comparison with the data from NS observations were performed in Refs. [73, 74], which gives bounds on the axion decay constant $f_{a}\gtrsim(5\textendash 10)\times 10^{7}\,\mathrm{GeV}$ for the KSVZ models and $f_{a}\gtrsim(5\textendash 15)\times 10^{7}\,\mathrm{GeV}$ for the DFSZ models. These results are obtained by taking account of the average temperature of the NS in the Cassiopeia A (Cas A) SN remnant rather than its transient behavior, and they can be regarded as conservative limits. On the other hand, some attempts focusing on the cooling rate of the Cas A NS were made in Refs. [75, 76], whose results turned out to be controversial. In Ref. [75], it was argued that a slight extra cooling due to the axion emission is required to explain the observed fast cooling of the Cas A NS, which shows a preference for a non-vanishing value of the axion-neutron coupling amounting to $f_{a}/C_{an}\sim 2.5\times 10^{9}\,\mathrm{GeV}$. Such a hint was not confirmed in the analysis performed in Ref. [76], which instead yields bounds of $f_{a}\gtrsim 5(7)\times 10^{8}\,\mathrm{GeV}$ for the KSVZ (DFSZ) models. In addition, the study of the axion emission from a different NS in HESS J1731-347 [77] also disagrees with the hint reported in Ref. [75], giving even a stronger bound $f_{a}/C_{an}\gtrsim 3.4\times 10^{9}\,\mathrm{GeV}$. However, given a limited understanding of the cooling of NSs, these results should be taken cautiously. Moreover, as for the cooling of the Cas A NS, the observational data themselves are inconclusive [78, 79, 80]. Because of these controversies, we do not include the results of observations of the NS cooling in the analyses performed in this paper.

2.6 Implications for axion models

Except for the ambiguities related to the cooling of SN and NS, various observations show a mild preference for additional cooling. Although individual anomalies are not quite significant, if combined, they lead to a clear systematic tendency of excessive energy losses. The interpretation in terms of several particle physics models including neutrino anomalous magnetic moments, minicharged particles, hidden photons, and axion-like particles (ALPs) was discussed in Ref. [8]. It turned out that axion/ALPs coupled to electrons are perfectly fit to all the anomalies while other candidates are inadequate for the explanation. The combined analysis of the hints from WD, HB, and RGB stars performed in Ref. [10] indicates non-vanishing couplings of axion/ALPs with electrons and photons.

Let us revisit the global analysis of the cooling hints by taking account of several systematic uncertainties that were not included in the analysis of Ref. [10]. In Fig. 1, we show the results for the 1, 2, and 3 $\sigma$ hinted regions in the $g_{ae}$-$g_{a\gamma}$ plane. In this analysis we include the systematic uncertainties associated with the effect of stellar rotation in the globular clusters and the model dependencies in the prediction for the rate of the period change of the WD variable PG 1351+489 (see Secs. 2.1 and 2.4). Furthermore, we take a different approach to treat the systematic uncertainties in the predicted TRGB brightness. The details of the analysis method are summarized in Appendix A.

Since we adopt conservative assessments on the systematic uncertainties, the significance of the hints becomes weaker than Ref. [10]. In particular, because of the systematic error corresponding to the effect of stellar rotation, the hint for a non-zero axion-photon coupling becomes less clear, and any value in the range of $g_{a\gamma}\lesssim 0.07\times 10^{-10}\,\mathrm{GeV}^{-1}$ gives an equally good fit. A non-zero axion-electron coupling is favored at around 2.4 $\sigma$, and the best fit value is $g_{ae}=1.56\times 10^{-13}$ with $\chi^{2}_{\rm min}/\text{d.o.f.}=14.7/15$.

Although the anomalous excessive cooling can be explained in a general framework of axion/ALPs where there is no particular relation between the mass and couplings, there remains a question of whether typical axion models that provide a solution to the strong CP problem accommodate the hints. The axion associated with the solution to the strong CP problem acquires a mass from the QCD effects, and the mass $m_{a}$ is related to the decay constant $f_{a}$. On the other hand, the prediction for coupling coefficients $C_{a\gamma}$, $C_{ae}$, $C_{ap}$, and $C_{an}$ is highly model dependent. By taking account of such model dependencies, the possibilities to interpret the stellar cooling anomalies in terms of the KSVZ and DFSZ axion models were explored in Ref. [10].

In the KSVZ models, the SM is extended by introducing one singlet complex scalar and exotic heavy quark(s). According to representations of the exotic quarks under the SM SU(2)${}_{L}\times$ U(1)Y gauge group, the low energy effective potential for the axion field can only have a single minimum. This fact is conventionally described as $N_{\rm DW}=1$, where $N_{\rm DW}$ is an integer representing the number of degenerate vacua and called the domain wall number. In this case, even though domain walls are formed around the epoch of the QCD phase transition, they decay immediately after the formation [81]. Hence there is no cosmological problem associated with the domain walls. However, since the SM leptons do not have PQ charges in these models, the electron coupling vanishes at tree level and emerges only at the loop level [82], which is too small to explain the anomalous excessive cooling. This fact excludes the possibility of a pure KSVZ axion interpretation of the stellar cooling anomalies.444If we extend the KSVZ models by adding extra three right-handed neutrinos (sometimes called a KSVZ-like axion/majoron model) [27, 28, 29], it becomes possible to obtain good fits to the data since the axion-electron coupling acquires extra loop contributions from neutrinos that can be adjusted to explain the observed anomalies [10]. However, such large loop contributions are generically in tension with the requirement of perturbativity.

In the DFSZ models, the SM is extended by introducing one singlet complex scalar and two Higgs doublets. In this case, we can obtain a suitable magnitude of the axion-electron coupling, and there exists a parameter region compatible with the observed data of the cooling anomalies [10]. However, the DFSZ models lead to multiple degenerate minima in the low energy effective potential, which causes the formation of domain walls if the PQ symmetry is broken after inflation [21]. In order to avoid the cosmological domain wall problem, we have to introduce some extra assumptions such as the existence of an explicit symmetry breaking term in the effective potential which leads to the late-time annihilation of domain walls [21, 83, 84, 85, 86, 87].

Given the fact that the pure KSVZ models are incompatible with the observations and that the DFSZ models suffer from the cosmological domain wall problem, it is reasonable to investigate whether other classes of models could give good fits to the data without causing any problem. In the following, we study the variant axion models [30, 31] as alternative possibilities to interpret the stellar cooling anomalies.

3 Variant axion models

We consider the PQ axion models with two Higgs doublets $H_{1}$, $H_{2}$ and one singlet scalar $\sigma$, and assume that the global U(1)PQ symmetry acts on them as

[TABLE]

where $k=1,2$ and $\epsilon$ is a real constant parameter.

In order to provide a solution to the strong CP problem, we need to assign the U(1)PQ charges for the quark fields to make the PQ symmetry anomalous for QCD. In other words, the low energy effective Lagrangian must contain the following term,

[TABLE]

where $\alpha_{s}$ is the strong coupling constant, $\mathcal{N}$ the QCD anomaly coefficient, $v_{\rm PQ}$ the PQ symmetry breaking scale, $G^{a}_{\mu\nu}$ the gluon field strength, and $\tilde{G}^{a\mu\nu}$ its dual. In the following, we focus on the cases where the QCD anomaly coefficient is given by

[TABLE]

which is a crucial element of variant axion models and different from the DFSZ model, for which $\mathcal{N}=3(X_{1}-X_{2})$ or $3(X_{2}-X_{1})$.

We note that the value of $\mathcal{N}$ can be related to the periodicity of the effective potential for the axion field, and hence it determines the vacuum structure of the theory. According to the PQ charge assignments, we can consider two possibilities [88]: If $2X_{\sigma}=X_{1}-X_{2}$, the Lagrangian contains an interaction term of the form $H_{1}^{\dagger}H_{2}\sigma^{2}$ and its hermitian conjugate, and we have $|\mathcal{N}|=2|X_{\sigma}|$. In this case, the periodicity of the axion field becomes twice of that of the QCD $\theta$ term, and the effective potential for the axion field has two degenerate minima, leading to the domain wall number $N_{\rm DW}=2$. On the other hand, if $X_{\sigma}=X_{1}-X_{2}$, the Lagrangian contains an interaction term proportional to $H_{1}^{\dagger}H_{2}\sigma$ and its hermitian conjugate. In this case, the periodicity of the axion field becomes the same as that of the QCD $\theta$ term, which implies $N_{\rm DW}=1$. In what follows we consider the latter case, since it avoids the cosmological domain wall problem. It should be emphasized that the cosmological domain wall problem cannot be solved straightforwardly in the DFSZ model, since it leads to $|\mathcal{N}|=6|X_{\sigma}|$ or $3|X_{\sigma}|$, which implies $N_{\rm DW}=6$ or $3$. We also note that the physical consequences depend only on the combination $X_{1}-X_{2}$, and that we can assign $X_{\sigma}=1$, $X_{1}=0$, and $X_{2}=-1$ without loss of generality.

The Yukawa interactions for the quark fields read

[TABLE]

where $H_{q}$ ($q=u,c,t,d,s,b$) take either $H_{1}$ or $H_{2}$, $\Gamma^{q}_{i}$ are Yukawa couplings with $i=1,2,3$ being a generation index, and $\tilde{H}_{k}=i\sigma_{2}H_{k}^{*}$. We consider the models where only one quark flavor has a nonzero PQ charge and couples to $H_{2}$ while the others couple to $H_{1}$ such that $|\mathcal{N}|=1$ in Eq. (3.3) is satisfied. There are six possibilities (Model U, C, T, D, S, B), and the corresponding charge assignments are shown in Table 1. For the lepton sector, we consider two types of interactions,555Although we consider flavor-blind PQ charge assignments for the lepton sector for simplicity, it is also possible to assign different PQ charges for each lepton flavor. In particular, the model can be embedded in SU(5) grand unified theory, where only one generation of $\overline{\bf 5}_{i}$ multiplet has a nonzero PQ charge. In that case, only one generation of $d_{iR}$ and $\ell_{iL}$ have nonzero PQ charges and couple to $H_{2}$.

[TABLE]

where $\Gamma^{\ell}_{ij}$ is the lepton Yukawa matrix with $i,j=1,2,3$ being generation indices. These terms are consistent with the charge assignments for the lepton fields shown in Table 2.

Let us derive axion couplings to ordinary matter in the variant axion models. The key ingredient is the PQ current $j^{\rm PQ}_{\mu}$, which is associated with the QCD and EM anomalies,

[TABLE]

where $\mathcal{E}$ is the EM anomaly coefficient. Before the spontaneous breaking of the electroweak symmetry, the PQ current is given by

[TABLE]

where $v_{\sigma}^{2}=2\langle|\sigma|^{2}\rangle$ is the vacuum expectation value (VEV) of the singlet scalar, $\tilde{a}$ is its phase direction $\sigma\propto e^{i\tilde{a}/v_{\sigma}}$, the summation over $\psi$ includes up-type quarks ($u$), down-type quarks ($d$), and leptons ($\ell$), and $X_{\psi_{iL,R}}$ represent PQ charges for left-handed ($L$) and right-handed ($R$) fermions, which are specified in Tables 1 and 2.

After the spontaneous breaking of the electroweak symmetry, the angular field $\tilde{a}$ mixes with a Nambu-Goldstone boson eaten by the $Z^{0}$ boson, and we have to redefine the axion field such that the associated PQ current gives the same QCD and EM anomalies as the original one. This amounts to the shifts of the PQ charges $X_{\psi_{iL,R}}$ for fermions in the PQ current [82, 89],

[TABLE]

where $Y_{k}$ are U(1)Y charges for two Higgs doublets $H_{k}$ ($Y_{1}=Y_{2}=1/2$), $v_{k}$ represent the VEVs of their neutral components, $v=\sqrt{v_{1}^{2}+v_{2}^{2}}\simeq 246\,\mathrm{GeV}$ is the electroweak scale, and $Y_{\psi_{iL,R}}$ are U(1)Y charges for the fermions. The PQ scale $v_{\rm PQ}$ can be related to two symmetry breaking scales $v_{\sigma}$ and $v$ as $v_{\rm PQ}=\sqrt{v_{\sigma}^{2}+v^{2}\sin^{2}\beta\cos^{2}\beta}$, where

[TABLE]

Based on the PQ current given by Eq. (3.8), we can construct the effective Lagrangian at energies below the electroweak scale but above the QCD scale,

[TABLE]

where

[TABLE]

is the axion decay constant.

So far it is implicitly assumed that the fermions $\psi_{iL,R}$ are in the weak interaction basis. Now we switch to the mass basis by performing unitary transformations $\psi_{L,R}\to U_{\psi_{L,R}}\psi_{L,R}$ that diagonalize the Yukawa matrices. After the transformations, the axion-fermion couplings become

[TABLE]

where

[TABLE]

and $X^{\prime}_{\psi_{L,R}}$ are understood as diagonal matrices. If $X^{\prime}_{\psi_{L,R}}$ are not proportional to the identity matrix, the coupling coefficients $C^{V}_{a\psi}$ and $C^{A}_{a\psi}$ depend not only on the effective PQ charges $X^{\prime}_{\psi_{L,R}}$ but also on the fermion mixing matrices $U_{\psi_{L,R}}$ [89, 36]. This is the case for the variant axion models, while in the DFSZ models we can eliminate $U_{\psi_{L,R}}$ as the PQ charges are generation-independent. In general, both $C^{V}_{a\psi}$ and $C^{A}_{a\psi}$ can have off-diagonal components due to the contribution from the fermion mixing matrices, while diagonal components of $C^{V}_{a\psi}$ can be eliminated by using the equations of motion.

Applying the PQ charge assignments for the quark fields shown in Table 1 to Eqs. (3.14) and (3.15), we obtain the explicit forms for the axion-quark couplings in the variant axion models. The results are summarized in Table 3. Here we introduced the following matrices,

[TABLE]

which represent the corrections arising from the mixings $U_{u_{R}},U_{d_{R}}$ of right-handed up-type and down-type quarks. The matrices $V_{u}$ and $V_{d}$ provide off-diagonal components of the axion-quark couplings, while diagonal components of $V_{u}$ and $V_{d}$ are positive. They satisfy the constraints $\operatorname{Tr}(V_{u})=\operatorname{Tr}(V_{d})=2$. Note that these matrices are different from the Cabibbo-Kobayashi-Maskawa mixing matrix $V_{\rm CKM}=U^{\dagger}_{u_{L}}U_{d_{L}}$, and they should be regarded as new parameters of the theory. The generation-dependence in the PQ charge assignments is now hidden in the matrices $V_{u,d}$, and in what follows we group it into two different classes of models, Model U, C, T and Model D, S, B, as shown in Table 3.

Since we consider flavor-blind PQ charge assignments for the lepton fields, there is no correction from mixing matrices to the axion-lepton couplings. Hence the coupling matrix $(C^{A}_{ae})_{ij}$ becomes proportional to the identity matrix, and the coupling strength is solely given by $\tan\beta$. Now it is straightforward to extract the coefficient for the axion-electron coupling $C_{ae}$ in Eq. (2.2) by using Eq. (3.15) and the PQ charge assignments specified in Table 2. In our setup there are four possibilities according to different PQ charge assignments for the quark and lepton fields, which we summarize in Table 4.666The difference in the sign of $C_{ae}$ between Model U, C, T and Model D, S, B originates from the sign of $\mathcal{N}$ appearing in Eq. (3.15). Note that this model-dependence does not affect the sign of $g_{ae}$, since it is proportional to $C_{ae}/f_{a}$ and $f_{a}$ contains another factor of $1/\mathcal{N}$ [see Eq. (3.12)]. Hence the electron coupling is solely determined by the Type of the PQ charge assignments for the lepton fields.

Below the QCD scale, the field $a$ mixes with light mesons such as neutral pions, and we have to find a state orthogonal to these mesons in order to identify the physical axion state at low energies. The mixing with neutral pions can be eliminated straightforwardly by considering the low energy effective theory emerging from the QCD with two lightest quarks $q^{\prime}=(u,d)$. First, we perform the following chiral transformations of up and down quarks,

[TABLE]

where $Q_{a}$ is a matrix acting on $q^{\prime}=(u,d)$. After this redefinition, the coefficient of the $G^{a}_{\mu\nu}\tilde{G}^{a\mu\nu}$ term in the effective Lagrangian is shifted by $\mathrm{Tr}(Q_{a})(a/f_{a})(\alpha_{s}/8\pi)$, and we can eliminate this coupling to gluons by taking $\mathrm{Tr}(Q_{a})=1$. Then, the terms including the lightest two quarks read

[TABLE]

where

[TABLE]

The leading order chiral Lagrangian that describes physics below the QCD scale contains the following term,

[TABLE]

where

[TABLE]

$f_{\pi}\simeq 92\,\mathrm{MeV}$ is the pion decay constant, and $\mu$ is a parameter that can be related to the pion mass. Expanding Eq. (3.39) to quadratic order in the fields, we find the axion-pion mixing, which can be eliminated by choosing $Q_{a}=M_{q}^{-1}/\mathrm{Tr}(M_{q}^{-1})$ [89]. After eliminating the axion-pion mixing and fixing the parameter $\mu$ in terms of the pion mass $m_{\pi}$, we can also extract the mass of the physical axion,

[TABLE]

The above formula for the axion mass is derived based on the tree-level axion-pion mixing in the effective theory with two quark flavors. More precise computation including the extra contributions from the strange quark was performed based on the next to leading order (NLO) [90] and next-to-next to leading order (NNLO) [91] chiral perturbation theory. Here we quote the NNLO result obtained in Ref. [91],

[TABLE]

This result is consistent with that of the direct calculation of the topological susceptibility in lattice QCD performed in Ref. [92].

The axion-photon coupling [Eq. (2.1)] consists of the model-dependent contribution given by the term proportional to $\mathcal{E}/\mathcal{N}$ in Eq. (3.11) and the model-independent contribution arising from the coupling to QCD. The tree level contribution to the latter can be extracted from a shift of the $F_{\mu\nu}\tilde{F}^{\mu\nu}$ term due to the rotation of two light quarks (3.32). A more precise result including the NLO corrections [90] reads

[TABLE]

In Table 4, we summarize the values of the model-dependent coefficient $\mathcal{E}/\mathcal{N}$ as well as the values of $C_{a\gamma}$ in the variant axion models.

The axion-nucleon couplings [Eq. (2.6)] can be derived by matching the low energy effective Lagrangian involving nucleons to the ultraviolet (UV) Lagrangian involving quarks. Similarly to the axion-photon coupling, the axion-nucleon couplings can be written in terms of the model-independent contributions from the mixing with mesons and the model-dependent contributions arising from the couplings to quarks in the UV Lagrangian (3.13). Here we adopt the results obtained in Ref. [90],

[TABLE]

Substituting the diagonal components of $(C^{A}_{au})_{ij}$ and $(C^{A}_{ad})_{ij}$ shown in Table 3 to the above equations, we obtain

[TABLE]

where $\delta_{N}$ is a correction term containing $(V_{u,d})_{22}$ and $(V_{u,d})_{33}$. This correction term is small compared to the leading terms as it originates from contributions of heavy quarks.777 Changing the value of $(V_{u,d})_{22}$ from [math] to $2$ with imposing the constraint $(V_{u,d})_{22}+(V_{u,d})_{33}=2$, we obtain $\delta_{N}\lesssim 0.024$ for Model U, C, T and $\delta_{N}\lesssim 0.076$ for Model D, S, B. Note that the constant terms in $C_{ap}$ and $C_{an}$ have uncertainty comparable to these upper limits on $\delta_{N}$. Hereafter we simply ignore this term and marginalize over the parameters $(V_{u,d})_{22}$ and $(V_{u,d})_{33}$. After this simplification, the nucleon couplings can be described by two parameters $\sin\beta$ and $(V_{u,d})_{11}$, where $(V_{u,d})_{11}$ take some values within the range of [math] to $2$.

In Fig. 2, we show the value of the combined quantity $\sqrt{C_{ap}^{2}+C_{an}^{2}}$, which is relevant to the bound from SN 1987A observations [Eq. (2.7)]. Since this quantity depends on the new parameter $(V_{u,d})_{11}$ in addition to $\tan\beta$, it shows a non-trivial structure in the two dimensional parameter space. In particular, the nucleon couplings become much smaller than unity around $\tan\beta\sim\mathcal{O}(1)$ and $(V_{u,d,})_{11}\sim 0$. To be more precise, $\sqrt{C_{ap}^{2}+C_{an}^{2}}$ takes a minimum value of 0.043 at $\tan\beta\simeq 0.70$ and $(V_{u})_{11}\simeq 0$ for Model U, C, T, and 0.042 at $\tan\beta\simeq 1.5$ and $(V_{d})_{11}\simeq 0$ for Model D, S, B. The fact that there exists a parameter region where the nucleon couplings get suppressed will play a role in the interpretation of observational results described in the next section.888It is possible to arrange the model such that the nucleon couplings are further suppressed, $C_{ap}\approx C_{an}\approx 0$ [93]. Indeed, this is the case for Model D, S, B if we take account of the small correction $\delta_{N}$ in addition to the dependence on $\tan\beta$ and $(V_{d})_{11}$. Note that, however, such a cancellation requires some tuning of the additional parameters $(V_{d})_{22}$ and $(V_{d})_{33}$. In this work, we do not consider such a tuned case and focus on the general consequences following from non-vanishing nucleon couplings.

Finally, it might be worthwhile to compare the axion couplings in the variant axion models derived above with those in the DFSZ axion models [23, 24]. The DFSZ models are built based on the generation-independent PQ charge assignments for the SM quarks and leptons, and we can consider two possibilities according to the following Yukawa interaction terms,

[TABLE]

These models lead to the following photon and electron couplings at low energies [82, 90]

[TABLE]

Comparing with Table 4, we see that $|C_{ae}|$ in Type I (Type II) of the variant axion models is a factor three larger than that of DFSZ I (DFSZ II). Furthermore, the nucleon couplings in the DFSZ models are given by [90]

[TABLE]

In contrast with the variant axion models, where the value of $\sqrt{C_{ap}^{2}+C_{an}^{2}}$ becomes as small as $\sim 0.04$ up to the choice of parameters, in the DFSZ models the suppression is milder, $\sqrt{C_{ap}^{2}+C_{an}^{2}}\gtrsim 0.24$. These differences lead to different consequences on the parameter regions hinted by the stellar cooling observations, as we discuss in the following section.

4 Interpretation of stellar cooling anomalies

Based on the axion couplings to ordinary matter derived in the previous section, we now investigate the parameter region where the variant axion models provide good fits for the stellar cooling hints. See Appendix A for the details of the analysis method.

In addition to the stellar cooling hints, we take account of the constraint from SN 1987A. In order to see how the SN 1987A bound affects the hinted parameter space, we first perform the global fits without including it, and add it afterwards. Following the procedure in Ref. [10], we include the SN 1987A bound by taking it as a $1\sigma$ hint of $g_{ap}^{2}+g_{an}^{2}=0$ with the error specified by the right-hand side of Eq. (2.7). On the other hand, we do not include the bounds or hints obtained by the observation of NSs, since there remain several controversies on the interpretation of observational results as mentioned in Sec. 2.

In the variant axion models, the axion-nucleon couplings depend on the new parameters $(V_{u,d})_{11}$ [see Eqs. (3.50) and (3.53)], and their values affect the constraint from SN 1987A. These parameters arise from the mixings of right-handed up-type and down-type quarks. In principle they can take any values within the range of 0 to 2. On the other hand, it is also possible to guess their values if we make some assumption about the UV-completion of the models. In Sec. 4.1 we first consider the case of general quark mixings where $(V_{u,d})_{11}$ are treated as extra free parameters. After that, in Sec. 4.2 we consider more specific cases where the parameters $(V_{u,d})_{11}$ are fixed to some well-motivated values.

4.1 Models with general quark mixings

If we do not consider the SN 1987A bound, the observational results can be explained in terms of the axion-photon coupling $g_{a\gamma}$ and axion-electron coupling $g_{ae}$. In the variant axion models these couplings are specified by two parameters $f_{a}$ (or $m_{a}$) and $\tan\beta$, in a similar manner to the DFSZ models. Figure 3 shows 1, 2, 3, 4 $\sigma$ hinted regions in the parameter space of $m_{a}$ and $\tan\beta$ obtained based on the data of WD, HB, and RGB cooling anomalies. Here we also show the projected sensitivities of IAXO and its upgrades (IAXO+) [10, 14]. Furthermore, the best fit parameter values are shown in blue dots, and they are also summarized in Table 5. Note that there exist upper and/or lower limits on the value of $\tan\beta$ due to the requirement of perturbativity of Yukawa interactions (see Appendix B), and they are shown as gray shaded regions.

Since most of the data prefer the electron coupling of $g_{ae}\sim 1.6\times 10^{-13}$, the overall shape of 1, 2, 3, 4 $\sigma$ contours in Fig. 3 can be understood by considering the parameter dependence of $g_{ae}$. For Type I (Type II) of variant axion models the electron coupling is given by $|g_{ae}|=m_{e}\sin^{2}\beta/f_{a}$ ($m_{e}\cos^{2}\beta/f_{a}$), and the contours become vertical at large (small) $\tan\beta$ since the coefficient becomes constant, $|C_{ae}|\simeq 1$ in that region. On the other hand, at small (large) $\tan\beta$ in Type I (Type II) models the coefficient $|C_{ae}|$ becomes suppressed and smaller values of $f_{a}$ are preferred in order to maintain the desired value of $g_{ae}$. This behavior ends when the value of $f_{a}$ becomes too small to make the photon coupling $g_{a\gamma}$ compatible with the constraint from the $R$-parameter. If we fix the type of the electron coupling and do not include the bound from SN 1987A, the difference between Model U, C, T and Model D, S, B is minor and only affects the behavior at higher $m_{a}$, which essentially comes from a difference in the model-dependent coefficient $\mathcal{E}/\mathcal{N}$ for the axion-photon coupling (see Table 4).

We note that the parameter range accessible to future helioscope experiments is also model-dependent. The sensitivity of IAXO for a generic axion model is given by [10]

[TABLE]

where $g_{a\gamma}^{\rm min}$ is the minimal value of $g_{a\gamma}$ accessible to IAXO, $d_{1}$ is the data of the minimal value of $g_{a\gamma}$ accessible to IAXO if solar axions are produced solely through the processes involving $g_{a\gamma}$, and $d_{2}$ is the data of the minimal value of $\sqrt{g_{a\gamma}g_{ae}}$ accessible to IAXO if solar axions are produced solely through the processes involving $g_{ae}$. $\zeta$ is a model-dependent dimensionless factor defined as

[TABLE]

with $g_{a\gamma}$ specified in units of $\mathrm{GeV}^{-1}$. We summarize the values of $\zeta$ in the variant axion models in Table 6.

If $\zeta$ is sufficiently small, the production of solar axions through the interaction with electrons becomes irrelevant. This case corresponds to the limit of $\tan\beta\ll 1$ ($\tan\beta\gg 1$) in Type I (Type II) variant axion models. In this case, the sensitivity is simply given by $g_{a\gamma}^{\rm min}\simeq d_{1}$, from which we can define minimal mass ranges accessible to IAXO regardless of the value of $\tan\beta$. Such mass ranges are shown in Table 6. On the other hand, if $\zeta$ becomes sufficiently large, which corresponds to the limit of $\tan\beta\gg 1$ ($\tan\beta\ll 1$) in Type I (Type II) variant axion models, the sensitivity is improved as $g_{a\gamma}^{\rm min}\simeq d_{2}/\sqrt{\zeta}$. In any case, the experimental sensitivity becomes improved if a model predicts a higher value of the photon coupling $|C_{a\gamma}|$. Since the value of $|C_{a\gamma}|$ in Type II models is larger than that in Type I models (see Table 4), a wider mass range can be probed for the former models compared to the latter ones.

Now, let us include the bound from SN 1987A. As nucleon couplings depend on new parameters $(V_{u,d})_{11}$, here we scan over three dimensional parameter space of $m_{a}$, $\tan\beta$, and $(V_{u,d})_{11}$. Since we add an extra data point for the SN constraint, the number of d.o.f. remains unchanged. The best fit parameters and the values of $\chi^{2}_{\rm min}/\text{d.o.f.}$ for the three dimensional scans are shown in Table 5. In Figs. 4 and 5, we also show the projection of hinted regions to two dimensional parameter space $(a,b)=(m_{a},\tan\beta)$, $(m_{a},(V_{u,d})_{11})$, or $((V_{u,d})_{11},\tan\beta)$ by computing $\Delta\chi^{2}(a,b)$, which is obtained by minimizing $\Delta\chi^{2}(a,b,c)=\chi^{2}(a,b,c)-\chi^{2}_{\rm min}$ over the third parameter $c$. Comparing the plots on the $m_{a}$-$\tan\beta$ plane with those for the results without including the SN data (Fig. 3), we see that every model still shows a good fit to the observed data in the lower mass ranges, while the higher mass ranges become incompatible with the SN bound.

The plots for the $m_{a}$-$(V_{u,d})_{11}$ plane in Figs. 4 and 5 show that the value of $(V_{u,d})_{11}$ does not have significant effects on the fits at lower mass ranges, where the bound from SN 1987A becomes less important. On the other hand, smaller values of $(V_{u,d})_{11}$ are preferred at higher mass ranges, since the SN bound becomes important but gets relaxed due to the reduction of the nucleon couplings for smaller values of $(V_{u,d})_{11}$ (see Fig. 2). To be more precise, in order to suppress the nucleon couplings the value of $\tan\beta$ should also be adjusted, which results in some model dependencies. In Type I models, smaller values of $\tan\beta$ are required at higher mass ranges in order to keep $g_{ae}$ close to the value hinted by the cooling anomalies, and for such smaller $\tan\beta$ values a non-vanishing value of $(V_{u,d})_{11}$ is slightly preferred over $(V_{u,d})_{11}=0$ as it leads to smaller values of the nucleon couplings in the small $\tan\beta$ region. This can be contrasted with the case of Type II models, where larger values of $\tan\beta$ are required at higher mass ranges and $(V_{u,d})_{11}=0$ is preferred such that the nucleon couplings become the smallest in the large $\tan\beta$ region.

The shape of the contours in the $(V_{u,d})_{11}$-$\tan\beta$ plane in Figs. 4 and 5 can also be understood in terms of the hinted value of $g_{ae}$ and the parameter dependence of the nucleon couplings. In Type I models, the large $\tan\beta$ region is preferred by WD, HB, and RGB cooling observations, and in the lower end of such a region the SN bound becomes relevant, which exhibits a slight preference for a non-vanishing value of $(V_{u,d})_{11}$. On the other hand, in Type II models, the small $\tan\beta$ region is preferred, and in the upper end of such a region $(V_{u,d})_{11}$ is forced to be zero in order to alleviate the SN bound by reducing the nucleon couplings.

In Fig. 6, we show the 2 $\sigma$ predicted regions in the parameter space of $m_{a}$ and $\sqrt{g_{ae}g_{a\gamma}}$ together with the projected sensitivities of IAXO and IAXO+. Numerical values of the axion mass and couplings for the corresponding ranges are also summarized in Table 7. We see that in any cases the lower end of the predicted mass ranges lie around $m_{a}\sim 0.45\,\mathrm{meV}$. This region corresponds to the limit of $\tan\beta\gg 1$ ($\tan\beta\ll 1$) in Type I (Type II) models, where the electron coupling becomes $|g_{ae}|\simeq m_{e}/f_{a}\simeq 4\times 10^{-14}\,(m_{a}/0.45\,\mathrm{meV})$ regardless of the value of $\tan\beta$, and hence the hinted value of $|g_{ae}|$ fixes the value of $m_{a}$.

From Fig. 6, we also see that the higher end of the predicted mass ranges becomes different according to the models. This difference can be understood as follows: First, we note that the contours in the $m_{a}$-$\tan\beta$ plane for Type II models (like those shown in Fig. 3) should be obtained by flipping those for Type I models vertically with respect to the line of $\tan\beta=1$ if we only take care of the structure of the electron coupling. This symmetric property along $\tan\beta$ direction is broken when we add the SN bound and consider the structure of the nucleon couplings. From left panel of Fig. 2, we see that in Model U, C, T the nucleon couplings become larger at large $\tan\beta$ region, while they become smaller at small $\tan\beta$ region. Therefore, the SN bound becomes weaker (stronger) for Type I (Type II) of Model U, C, T, where smaller (larger) values of $\tan\beta$ are favored at higher masses, which results in a less (more) stringent upper limit on the axion mass. On the other hand, in Model D, S, B the asymmetry of the nucleon couplings along $\tan\beta$ direction is less significant at $(V_{d})_{11}\simeq 0$ (see right panel of Fig. 2), and the strength of the SN bound remains almost the same for Type I and Type II models. This explains the fact that there is little difference in the upper limit of the predicted mass ranges between Type I and Type II of Model D, S, B.

In addition to the difference in the predicted mass ranges described above, Fig. 6 shows that the overall magnitude of the coupling product $g_{ae}g_{a\gamma}$ differs according to the models. Such a difference can be understood in terms of the value of the coefficient $|C_{a\gamma}|$ of the axion-photon coupling (see Table 4). For Type I models, some part of the predicted region is inaccessible to IAXO since the values of $|C_{a\gamma}|$ in these models are relatively small. On the other hand, most part of the predicted parameter region is covered by the projected sensitivity of IAXO for Type II models as they predict larger values of $|C_{a\gamma}|$, and hence we expect that these models can be tested by the forthcoming experiments.

Before closing this subsection, we note that the models can potentially be constrained by the observation of heavy meson decays [94], since they predict flavor-changing couplings with quarks given by the off-diagonal elements of the matrices $(V_{u,d})_{ij}$. Currently the strongest constraint is obtained from a search for a decay $K^{+}\to\pi^{+}+a$ [95], which leads to a bound [36]

[TABLE]

This bound is trivially satisfied in Model U, C, T, since they predict $(C^{A}_{ad})_{21}=0$ (see Table 3), while Model D, S, B can have a non-vanishing value of $|(C^{A}_{ad})_{21}|=|(V_{d})_{21}|$ and in such cases the bound could become relevant. However, at this point we should treat $(V_{d})_{21}$ as an extra free parameter in the same way as $(V_{u,d})_{11}$, and cannot deduce a definite bound from Eq. (4.29) unless we fix the pattern of quark mixings to obtain a specific value of $(V_{d})_{21}$.

4.2 Models with small quark mixings

So far we have considered the general cases of the variant axion models where the mixings of the right-handed quarks are arbitrary. Although such models already show good fits to the observational data, it is enlightening to explore a more specific scenario of the quark mixings. As a well-motivated possibility, here we consider a pattern of quark mixings inspired by the Frogatt-Nielsen mechanism [96, 97], which was introduced to provide a natural explanation to the SM quark and lepton mass hierarchy.

In the Frogatt-Nielsen scenario, the flavor structure is related to a spontaneously broken global U(1)FN symmetry.999Note that this global symmetry is different from the PQ symmetry in our setup. Some attempts to relate the U(1)FN symmetry to the PQ symmetry are found in Refs. [34, 35]. The Yukawa interactions of the SM quarks and leptons originate from non-renormalizable interactions with a singlet field $\Phi$ suppressed by a cutoff scale $M$, and the mass hierarchy of quarks and leptons is generated as a power of $\varepsilon=\langle\Phi\rangle/M$ after the U(1)FN symmetry is spontaneously broken and the singlet field acquires a VEV $\langle\Phi\rangle$. The magnitude of fermion mixings is also specified as $(U_{\psi_{L}})_{ij}\sim\varepsilon^{|Q_{\psi_{iL}}-Q_{\psi_{jL}}|}$ and $(U_{\psi_{R}})_{ij}\sim\varepsilon^{|Q_{\psi_{iR}}-Q_{\psi_{jR}}|}$ [96], where $Q_{\psi_{iL,R}}$ are U(1)FN charges of the fermions. The observed quark and lepton masses imply $\varepsilon^{2}\sim 1/300$.

Assuming that the flavor structure follows a pattern predicted by the Frogatt-Nielsen mechanism, we can guess the magnitude of the parameters $(V_{u,d})_{11}$ appearing in the axion-nucleon couplings. In particular, adopting the U(1)FN charges specified in Ref. [97], from Eqs. (3.16), (3.17), and (3.21) we obtain $(V_{u,d})_{11}\sim\mathcal{O}(\varepsilon^{2})$ for Model U and D. This fact motivates us to study the consequences of Model U and D with a negligibly small value of $(V_{u,d})_{11}$.

In Fig. 7, we show the 1, 2, 3, 4 $\sigma$ hinted regions from the stellar cooling observations including the bound from SN 1987A for Model U and D with $(V_{u,d})_{11}=0$. The values of the axion mass and couplings corresponding to the 2$\,\sigma$ regions are also summarized in Table 7. In these models $(V_{u,d})_{11}$ are no longer free parameters, and we just scan over two dimensional parameter space of $m_{a}$ and $\tan\beta$. Overall, the shapes of the contours for Model U and D with $(V_{u,d})_{11}=0$ are almost the same as those for models with general quark mixings shown in Figs. 4 and 5. The difference between general and specific cases only appear in the region with higher values of $m_{a}$ where the effect of the axion-nucleon couplings becomes relevant because of the constraint from SN 1987A. Comparing Fig. 7 with Figs. 4 and 5, we see that the contours extend to $m_{a}\lesssim 10^{-2}\,\mathrm{eV}$ in Type I of Model U and D with small quark mixings, while they reach slightly higher mass regions in Type I models with general quark mixings. This difference arises from the fact that $(V_{u,d})_{11}$ are fixed to be zero in the former cases while in the latter cases they can be adjusted to compensate the increase in the nucleon couplings due to the lower values of $f_{a}$. Such a difference is not clearly seen in the contours for Type II models, since for these models the general setup already shows a preference for $(V_{u,d})_{11}=0$ even at higher masses.

The best fit parameter values and the values of $\chi^{2}_{\rm min}/\text{d.o.f.}$ for the interpretation of the stellar cooling anomalies in the framework of the variant axion models with a Frogatt-Nielsen like flavor structure are summarized in Table 8. From Table 8 we see that the fits including WD, RGB, HB, and SN data give $\chi^{2}_{\rm min}/\mathrm{d.o.f.}=14.7/16$, which is contrasted with $\chi^{2}_{\rm min}/\mathrm{d.o.f.}=15.4/16$ for DFSZ I [10]. This improvement originates from a specific structure of the axion-nucleon and axion-electron couplings in the variant axion models. In Sec. 3, we see that the combined nucleon coupling coefficient $\sqrt{C_{ap}^{2}+C_{an}^{2}}$ can become as small as $\sim 0.04$ in the variant axion models with $(V_{u,d})_{11}\simeq 0$, while it becomes larger than $\gtrsim 0.24$ in the DFSZ models. Even if we take account of the variation due to the change of $\tan\beta$, the magnitude of the nucleon couplings remains comparable between the variant and DFSZ models as long as the value of $\tan\beta$ is not far from $\mathcal{O}(1)$. On the other hand, the coefficient of the axion-electron coupling $|C_{ae}|$ of the variant axion models is a factor three larger than that of the DFSZ models. Combining these facts, we see that the ratio $\sqrt{C_{ap}^{2}+C_{an}^{2}}/|C_{ae}|$ in the variant axion models can be smaller than that in the DFSZ models for the relevant range of $\tan\beta$. This fact slightly relaxes the bound from SN 1987A, providing better fits to the stellar cooling hints.

For the sake of comparison, in Fig. 7 we also show the results for the DFSZ models. Here we see that the difference in the coefficient of the axion-electron coupling between the variant and DFSZ axion models affects the axion mass range favored by the stellar cooling hints. The factor three difference in $|C_{ae}|$ is compensated by shifting the value of $f_{a}$ by the same amount. As a result, the variant axion models provide good fits even in the mass ranges lower than those predicted in the DFSZ models.

Finally, we note that the bound from the kaon decay [Eq. (4.29)] is trivially satisfied in Model U, since it predicts $(C^{A}_{ad})_{21}=0$, which is exactly the same as in the general cases of Model U, C, T discussed in the previous subsection. On the other hand, this bound might give rise to a tension with the interpretation of the stellar cooling anomalies in Model D with Froggatt-Nielsen like quark mixings, since we have $(V_{d})_{21}\sim\mathcal{O}(\varepsilon)$ in this case.

5 Discussion and conclusions

In this paper, we have revisited the global analysis of the stellar cooling anomalies and their axion/ALP interpretation. We have adopted a conservative approach to include possible systematic uncertainties associated with the effect of stellar rotation in the globular clusters, prediction for the TRGB brightness, and different models for the pulsating WD PG 1351+489. As a result of this procedure, the significance of the cooling hints becomes weaker, pointing to a non-vanishing axion-electron coupling of $g_{ae}\sim 1.56\times 10^{-13}$ at around 2.4 $\sigma$. Furthermore, we emphasize that there could be potentially large systematic uncertainties associated with survey incompleteness in the observed WDLF (see Sec. 2.3), which we have not included in the analysis presented in this paper. These situations should be reviewed once further observational data become available.

With the revised results of the global fits, we have considered the variant axion models as possible explanation for the anomalous cooling observed in WD, HB, and RGB stars. The models are constructed by considering six different flavor-dependent PQ charge assignments for the SM quarks and two different flavor-blind PQ charge assignments for the SM leptons, which result in four different possibilities for axion couplings to ordinary matter. For each case, coupling coefficients are derived systematically, and they are summarized in Tables 3 and 4, and Eqs. (3.50) and (3.53). By using the derived couplings, we have performed the global fits with the WD, HB, RGB data and the bound from SN 1987A by considering both the general cases for quark mixings and specific cases motivated by the Froggatt-Nielsen mechanism. Every model shows a quite good fit to the data within mass ranges around $0.45\,\mathrm{meV}\lesssim m_{a}\lesssim 30\,\mathrm{meV}$ as shown in Figs. 3, 4, 5, and 7, and Tables 5 and 8. The 2 $\sigma$ parameter regions preferred by the cooling hints and compatible with the SN bound and perturbativity requirements are shown in Fig. 6 and Table 7.

We emphasize that the models presented in this paper can resolve two fundamental issues that prevent the canonical KSVZ and DFSZ axion models from interpreting the stellar cooling anomalies straightforwardly:

The KSVZ models do not account for the large axion-electron coupling preferred by the stellar cooling hints. A KSVZ-like axion/majoron model might account for it due to extra loop contributions from neutrinos, but there is a tension with perturbativity requirements [10]. This drawback is absent in the variant axion models, since a sizable electron coupling arises naturally at tree level. 2. 2.

The DFSZ models predict a sufficiently large axion-electron coupling and hence provide good fits to the observational data, but they suffer from the cosmological domain wall problem. In the variant axion models, it is possible to avoid the domain wall problem (i.e. $N_{\rm DW}=1$) with keeping the structure of the axion-electron coupling similar to that in the DFSZ models.

In addition to the above two important points, we have also found that the models where a nonzero PQ charge is assigned to first generation quarks give a slightly better fit than the DFSZ models if we make an assumption of small mixings of right-handed quarks as suggested by the Froggatt-Nielsen mechanism. This is because the constraint from SN 1987A is slightly relaxed due to the fact that the ratio of the axion-nucleon coupling to the electron coupling becomes smaller than that in the DFSZ models in the relevant parameter regions.

Another important feature of the variant axion models is that they predict flavor-changing couplings with quarks parameterized by the matrices $(V_{u,d})_{ij}$ (see Table 3). One consequence of such flavor-changing couplings is the possibility of rare decays of heavy mesons into axions. In particular, the search for the decay process $K^{+}\to\pi^{+}+a$ could lead to a severe constraint on this class of models. We have found that such a process is absent and the constraint is trivially satisfied in the models with a PQ charge assigned for an up-type quark, while it may give rise to some tension in those with a PQ charge assigned for a down-type quark.

In addition to the above constraint from the decay of kaons, processes involving other off-diagonal elements of $(V_{u,d})_{ij}$ can be probed by various other decay channels. Present and expected future limits on the flavor-changing couplings from heavy meson decays are summarized in Ref. [36]. Furthermore, if the extra heavier states coming from two Higgs doublets (see Appendix B) have masses comparable to the electroweak scale, the models predict flavor-changing neutral-current processes that can be searched by collider experiments [98, 99, 100]. Such laboratory searches can complement the results of astrophysical observations discussed in this paper. Indeed, exploring flavor-changing interactions would be crucial to distinguish between the variant axion models and other flavor-blind models such as the DFSZ models.

The mass range $m_{a}\gtrsim 0.45\,\mathrm{meV}$ predicted by the variant axion interpretation of the stellar cooling anomalies might have intriguing implications for cosmology. In this mass range, the axion decay constant is of order $f_{a}\sim 10^{9}\,\mathrm{GeV}$, and it is preferable to assume that the PQ symmetry is broken after inflation since otherwise there is a severe constraint from isocurvature fluctuations as mentioned in Sec. 1. In such a post-inflationary PQ symmetry breaking scenario, axions can be produced by the decay of strings and domain walls, and they behave as cold dark matter in the present universe [101, 102]. Although the previous estimates of the relic axion abundance from the decay of string-wall systems for models with $N_{\rm DW}=1$ showed that such axions would comprise only a small fraction of the cold dark matter in the mass range $m_{a}\gtrsim 1\,\mathrm{meV}$ [see e.g. Refs. [103, 86, 104, 105]], it has been pointed out recently that there could be a much broader uncertainty in those estimates [106]. Within the broad uncertainty suggested in Ref. [106], there remains a possibility that axions become the main constituent of dark matter up to the mass of $m_{a}\lesssim 4.4\,\mathrm{meV}$ [14]. Indeed, preliminary results of the state of the art simulations [107, 108] show a preference for higher dark matter mass ranges close to this upper limit on $m_{a}$. Therefore, there is a possibility that the $\mathrm{meV}$ mass variant axions could account for dark matter as well as providing the explanation of the stellar cooling anomalies. Direct detection of axion dark matter in the $\mathrm{meV}$ mass range is quite challenging, but several techniques are proposed in the literature, such as a dish antenna [109], absorption in superconductors [110], and topological insulators [111].

The mass range preferred by the stellar cooling anomalies will also be probed by the proposed experiment ARIADNE [112], which is the search for long range forces mediated by axions [113]. This experiment does not rely on any assumption on astrophysics (e.g. axion emissivity from the Sun) or cosmology (e.g. local dark matter density), and hence it would play a complementary role with respect to other axion searches. However, in order to achieve a reasonable sensitivity it is necessary to assume that there exists a sizable CP-violating interaction between axions and nuclei saturating the limit on the neutron electric dipole moment. Such a large CP-violating coupling is not necessarily guaranteed in the models considered in this paper.

Crucially, the most part of the parameter regions predicated by the stellar cooling hints will be covered by the next generation helioscope IAXO. Once axions are detected in such experiments, it would be even possible to measure the axion mass and couplings to photons and electrons if the detectors have enough energy resolution [114, 115]. The predictions of Type II of the variant axion models exhibit relatively high values for the coupling product $g_{ae}g_{a\gamma}$ (see Fig. 6), and hence this class of models would potentially be distinguishable from other models through such measurements in the future helioscope experiments.

Acknowledgments

K. S. would like to thank Maurizio Giannotti and Javier Redondo for a lot of important comments. We also wish to thank the anonymous referee for his/her constructive comments and suggestions that lead to great improvements in the quality of the paper. K. S. acknowledges partial support by the Deutsche Forschungsgemeinschaft through Grant No. EXC 153 (Excellence Cluster “Universe”) and Grant No. SFB 1258 (Collaborative Research Center “Neutrinos, Dark Matter, Messengers”) as well as by the European Union through Grant No. H2020-MSCA-ITN-2015/674896 (Innovative Training Network “Elusives”). T. T. Y. was supported in part by the China Grant for Talent Scientific Start-Up Project and the JSPS Grant-in-Aid for Scientific Research No. 16H02176, and No. 17H02878 and by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan. T. T. Y. thanks to Hamamatsu Photonics.

Appendix A Analysis method

For the global fits we take a slightly different procedure than Ref. [10]. In Sec. 2, we see that there are several sources of systematic uncertainties, which might not be regarded as simple Gaussian-like fluctuations. One possible way to deal with such systematic errors is to introduce additional parameters (nuisance parameters) that are marginalized to construct a profile likelihood function [e.g. Ref. [116]].

Our general procedure to marginalize the systematic variables is as follows. Suppose that an outcome of an observation is measured as a value $M$ with an error $\sigma$. Let $T(\bm{\theta})$ denotes a theoretical prediction for this quantity, which depends on a set of parameters $\bm{\theta}=(\theta_{1},\dots,\theta_{N})$. We assume that $T(\bm{\theta})$ has some systematic uncertainty, whose typical magnitude is given by $\Delta T(\bm{\theta})$. The effect of the systematic uncertainty can be embedded as a shift of $T(\bm{\theta})-M$ by an amount $\xi\Delta T(\bm{\theta})$, where $\xi$ is an additional parameter spanning a finite interval. The values of $\bm{\theta}$ compatible with the observational data are found by evaluating the following quantity,

[TABLE]

where $\hat{\xi}(\bm{\theta})$ is a value that minimizes $\chi^{2}$ for specified values of $\bm{\theta}$. The above formula resembles the pull $\chi^{2}$ function [117, 118], but here we assume that the variable $\xi$ obeys a flat probability distribution rather than the Gaussian distribution.

One can evaluate $\hat{\xi}(\bm{\theta})$ almost trivially: It takes the boundary value $\xi_{\rm min}$ or $\xi_{\rm max}$ of the specified domain $\xi\in[\xi_{\rm min},\xi_{\rm max}]$, otherwise $\hat{\xi}(\bm{\theta})=(M-T(\bm{\theta}))/\Delta T(\bm{\theta})$, for which $\chi^{2}=0$. The appearance of the flat direction $\xi=(M-T(\bm{\theta}))/\Delta T(\bm{\theta})$ just implies that there is a degeneracy between the systematic effect and the effect of theoretical parameters, and that the values of the parameters $\bm{\theta}$ cannot be determined unless the corresponding systematic uncertainty is resolved.

We apply the above procedure to the analysis of the astrophysical data. For the $R$-parameter, the observed value is $R^{\rm obs}=1.39\pm 0.03$, and the theoretical model is given by Eq. (2.3). We use the value $Y=0.2535\pm 0.0036$ as representative of the mass fraction in low metallicity environments [41, 40], and the corresponding error $\sigma_{Y}=0.026$ is added to the observational error $\sigma_{R}=0.03$. Furthermore, we introduce a parameter $\xi_{\rm rot}$ corresponding to the systematic uncertainty due to the stellar rotation. $\xi_{\rm rot}$ is allowed to vary in the range $[-1,0]$: It is negative since the effect of rotation just shortens the helium-burning lifetime [46]. The magnitude of this uncertainty is estimated as $3\,\%$ of the observed value, $\Delta R=0.042$.

For the TRGB of the globular cluster M5, the observed value is $M_{I,\mathrm{TRGB}}^{\rm obs}=-4.17\pm 0.13\,\mathrm{mag}$. We use the theoretical model of Eq. (2.5), but here we drop the last term $\delta M_{I,\mathrm{TRGB}}^{\rm th}$, since it was introduced in order to correct the asymmetric ranges of systematic errors and convert them to a Gaussian error in Ref. [38, 39], which we do not follow in our analysis. Instead, we introduce a parameter $\xi_{\rm BC}$, which is defined in the range $[-1,1]$, to represent the systematic uncertainties. We use the error estimate for bolometric corrections as the magnitude of the uncertainty, $\Delta M_{I,\mathrm{TRGB}}=(0.08+0.02g_{13})\,\mathrm{mag}$, since it gives the largest range among the possible sources of systematic uncertainties enumerated in Refs. [38, 39].

For the WDLF, we follow the same procedure as in Ref. [10]. We take 11 binned data points in the luminosity range $7<M_{\rm bol}<12.25$ from Refs. [47, 50]. For the theoretical model, we allow the normalization factor to vary as a function of $g_{ae}$ and marginalize it such that it leads to the smallest $\chi^{2}$.

For the period change of the WD variables, we adopt the observational data of R548 [$\dot{\Pi}=(3.3\pm 1.1)\times 10^{-15}\,\mathrm{s/s}$] [56], PG 1351+489 [$\dot{\Pi}=(2.0\pm 0.9)\times 10^{-13}\,\mathrm{s/s}$] [57, 58], and two pulsation modes of L19-2 [$\dot{\Pi}=(3.0\pm 0.6)\times 10^{-15}\,\mathrm{s/s}$ for both modes] [59]. Following the suggestion in Ref. [10], we exclude the data of G117-B15A [55] and add the theoretical uncertainties for R548 [$\sigma_{\dot{\Pi}_{\rm th}}=0.09\times 10^{-15}\,\mathrm{s/s}$] and for L19-2 [$\sigma_{\dot{\Pi}_{\rm th}}=0.85\times 10^{-15}\,\mathrm{s/s}$ and $\sigma_{\dot{\Pi}_{\rm th}}=1.45\times 10^{-15}\,\mathrm{s/s}$] in quadrature. For PG 1351+489, instead of adding the $\sigma_{\dot{\Pi}_{\rm th}}^{2}$ error, we introduce a parameter $\xi_{\rm model}$ to take account of the systematic uncertainty due to different models. We extract 7 different model curves from Figs. 6 and 7 of Ref. [58] and estimate the “theoretical” prediction $\dot{\Pi}(g_{ae})$ from the median between the largest and smallest value among the 7 curves. Then the magnitude of the systematic uncertainty $\Delta\dot{\Pi}$ is identified by using half of the difference between the largest and smallest value among those curves. We marginalize it by allowing $\xi_{\rm model}$ to vary in the range $[-1,1]$.

In summary, the $\chi^{2}$ function is

[TABLE]

where $M_{I,\mathrm{TRGB}}^{0}(g_{ae})=M_{I,\mathrm{TRGB}}^{\rm th}(g_{ae})-\delta M_{I,\mathrm{TRGB}}^{\rm th}$. When we include the SN 1987A bound, we add

[TABLE]

Appendix B Higgs sector and Yukawa interactions

In this appendix, we investigate the structure of the Higgs sector and Yukawa interactions in the variant axion models for the purpose of obtaining typical ranges of the parameter $\tan\beta$ compatible with perturbativity of the Yukawa interactions. The most general renormalizable Higgs potential at energies below the PQ scale but above the electroweak scale reads [98]

[TABLE]

where the $m_{12}^{2}$ terms originate from the interaction terms with the singlet scaler (i.e. terms proportional to $H_{1}^{\dagger}H_{2}\sigma$ and its hermitian conjugate) in the UV-complete theory (see Sec. 3), and we can make the parameter $m_{12}^{2}$ real and positive by the PQ symmetry transformation. All the other parameters $m_{11}^{2}$, $m_{22}^{2}$, and $\lambda_{1,2,3,4}$ are real.

After the electroweak symmetry breaking, two Higgs doublet fields can be decomposed into their VEVs $v_{k}$ and component fields, $H_{k}=(H_{k}^{+},(v_{k}+h_{k}+iA_{k})/\sqrt{2})$. These fields can be related to the SM Higgs field $H^{\rm SM}$ and the orthogonal field $H^{\prime}$,

[TABLE]

where

[TABLE]

Among eight field degrees of freedom, $G^{\pm}$ and $G^{0}$ correspond to the Nambu-Goldstone bosons which get eaten by $W^{\pm}$ and $Z^{0}$ bosons. The charged scalar $H^{\pm}$ and neutral pseudoscalar $A^{0}$ are mass eigenstates, while neutral scalars $h^{\rm SM}$ and $h^{\prime}$ are in general not mass eigenstates. We define the rotation angle $\alpha$ which relates these neutral scalars to the lighter ($h$) and heavier ($H$) mass eigenstates,

[TABLE]

Note that all the heavier states $H^{\pm}$, $A^{0}$ and $H$ decouple in the limit $m_{12}^{2}\to\infty$, and these states become irrelevant to low energy physics, in which the lightest eigenstate $h$ is identified as the SM Higgs boson. However, even in such a case, it is possible to restrict parameters of the models by considering perturbativity of Yukawa interactions.

Yukawa interactions of neutral Higgs bosons can be parameterized as [119]

[TABLE]

where the fermions $\psi_{i}$ are taken to be mass eigenstates. We recall that the Yukawa interactions are different according to the PQ charge assignments for the quark and lepton fields [see Eqs. (3.4) and (3.5)]. The coupling coefficients read

[TABLE]

for the up-type quarks, where $V_{u}$ is defined in Eq. (3.16),

[TABLE]

for the down-type quarks, where $V_{d}$ is defined in Eq. (3.17), and

[TABLE]

for leptons.

In order to guarantee perturbativity of the Yukawa interactions, we require

[TABLE]

where

[TABLE]

for $\phi=h,H,A^{0}$. Requiring that the above constraints are satisfied for any value of $\alpha$, we obtain

[TABLE]

where the upper limit for Model U, C, T comes from bottom Yukawa interactions, that for Model D, S, B comes from top Yukawa interactions, and the lower limit for Type II comes from tau Yukawa interactions. Here we have omitted the bounds obtained from couplings that depend on $(V_{u})_{ij}$ and $(V_{d})_{ij}$, whose values are unknown.101010 Although we have specified the values of $(V_{u})_{11}$ and $(V_{d})_{11}$ in the plots shown in Sec. 4, the bounds obtained from couplings that depend on these parameters are quite weak because of the smallness of up and down quark masses, and they are irrelevant in the range of $\tan\beta$ shown in those figures. Note that, however, more stringent limits could be obtained according to their values. In the main text, we have used Eq. (B.35) as a typical range compatible with perturbativity requirements.