Flavor-Violating Higgs Decays and Stellar Cooling Anomalies in Axion Models

We study a class of DFSZ-like models for the QCD axion that can address observed anomalies in stellar cooling. Stringent constraints from SN1987A and neutron stars are avoided by suppressed couplings to nucleons, while axion couplings to electrons and photons are sizable. All axion couplings depend on few parameters that also control the extended Higgs sector, in particular lepton flavor-violating couplings of the Standard Model-like Higgs boson $h$. This allows us to correlate axion and Higgs phenomenology, and we find that that ${\rm BR}(h \to \tau e)$ can be as large as the current experimental bound of 0.22%, while ${\rm BR} (h \to \mu \mu)$ can be larger than in the Standard Model by up to 70%. Large parts of the parameter space will be tested by the next generation of axion helioscopes such as the IAXO experiment.


Introduction and Motivation
Arguably the QCD axion is one of the best candidates for New Physics beyond the Standard Model (SM), being motivated not only by the Peccei-Quinn solution to the strong CP Problem [1][2][3][4], but also by the observed Dark Matter abundance [5][6][7]. Interestingly, the axion can also account for excessive energy losses observed in several stellar environments [8][9][10][11][12][13], which hints to new cooling channels such as a light axion with large couplings to electrons. This kind of scenarios are however constrained by other astrophysical observations that strongly constrain axion couplings to nucleons, such as the observation of the neutrino burst in SN1987A [14] or neutron star cooling [15][16][17][18][19].
Taking the stellar cooling hints seriously therefore points to a rather special structure of axion couplings, which definitely prefers the class of DFSZ-like axion models [20,21], in which the axion has large couplings to SM fermions [13]. Since all couplings are proportional to the axion mass, the required size of electron couplings puts a lower bound on the axion mass of the order of a few meV, corresponding to axion decay constants of the order of 10 9 GeV. While the standard DFSZ benchmark models discussed in Ref. [13] have some tension with perturbativity and SN1987A constraints for such heavy axions, several modifications of DFSZ models (or "axion variant models" [22,23]) have been proposed that can ease this tension as a result of suppressed couplings to nucleons [24][25][26]. Moreover, these models can feature a trivial domain wall number that elegantly avoids the cosmological domain wall problem [27].
In this article we revisit these "nucleophobic" axion models which are simple Two-Higgs-Doublet models (2HDMs) with a global Peccei-Quinn (PQ) symmetry, that are constructed analogous to the standard DFSZ axion models, but with flavor-dependent PQ charges. In general the axion coupling to electrons depends on free parameters that also control lepton flavor-violating (LFV) effects, which are mediated by LFV couplings of both the axion and the physical Higgs scalars. In contrast to earlier analyses, which have focussed on implications for axion searches, here we are interested also in the Higgs phenomenology, which is correlated to the axion phenomenology to a large extent. Thus, we have in mind a DFSZ models with a light second Higgs doublet, whose scale is subject only to present experimental constraints. This setup allows us to predict various flavor-violating Higgs decays probed at the LHC in terms of the same parameters that control the axion couplings to nucleons, electrons and photons, which are constrained by astrophysical observations and probed by dedicated axion searches such as IAXO [28]. For similar studies of the phenomenological implications of DFSZ models with light Higgs doublets see e.g. Refs. [29][30][31]. We find, in particular, that in the region explaining the stellar cooling anomaly, the branching ratio for the h → τ e decay can be as large as the current LHC upper bound, while the couplings of the Higgs to taus and muons may deviate significantly from the SM prediction.
This article is structured as follows: In Section 2 we define our basic setup and provide the expressions for axion and Higgs couplings to the SM in terms of the model parameters. In Section 3 we provide the constraints on these parameter from axion physics, and identify the region preferred by stellar cooling anomalies. In Section 4 we discuss the phenomenology of the extended Higgs sector, focussing on modifications of the SM-like Higgs couplings to leptons, in particular LFV couplings. We combine all constraints in Section 5 and study the implications for future searches at the LHC and axion helioscopes in the cooling hint region, before we conclude in Section 6.

Setup
In this section we first discuss the general effective Lagrangian for the QCD axion, before we define the special DFSZ-like UV completions that we want to consider. These models fix not only the couplings of the QCD axion, but also the couplings in the extended Higgs sector in terms of a few parameters. This structure gives rise to a correlated axion-Higgs phenomenology that we will analyze in the following sections.

Axion Effective Lagrangian
At energies much below the PQ breaking scale, the effective axion couplings to gauge fields and fermions are given by where f a is the axion decay constant, FF ≡ 1 2 µνρσ F µν F ρσ with the electromagnetic (EM) field strengths and similar in the gluon sector, E/N is the ratio of EM and color anomaly coefficients and we use the convention 0123 = −1 (for more details see Appendix A).
For the purpose of addressing the stellar cooling anomalies with axions it is helpful to have small couplings to nucleons in order to avoid the stringent constraints from SN1987A and neutron star cooling, cf. Section 3.3. From Eqs. (2.4) and (2.5) it is clear that axion couplings to nucleons can be suppressed if the UV quark couplings satisfy the approximate relations C u 1 1 + z ≈ 2/3 , C u + C d 1 . (2.6) As analyzed in detail in Ref. [34] (see also Refs. [35,36]), one can realize these conditions in the context of "nucleophobic" DFSZ models that we will now discuss in more details.

UV Lagrangian
We add two Higgs doublets h i with hypercharge Y = −1/2 and a SM singlet φ to the SM. The Lagrangian admits a U (1) PQ symmetry, with fermion charges that are consistent with a 2 + 1 flavor structure. This symmetry is broken twice: explicitly by the QCD anomaly and spontaneously by Higgs and singlet vacuum expectation values (VEVs). The QCD axion is the pseudo-Nambu-Goldstone boson of the PQ symmetry, and thus a linear combination of the CP-odd components of all PQ-charged fields with VEVs, with each coefficient given by the respective PQ charge and VEV (up to a normalization). In Appendix A we summarize the general structure of these models, which are defined by the Yukawa couplings and the Higgs-singlet interactions needed to break additional U (1) factors. These Lagrangian parameters govern the axion couplings to matter, apart from fermion flavor mixing. As discussed above, here we are only interested in models that have potentially suppressed couplings to nucleons, which requires approximately C u ≈ 2/3, C d ≈ 1/3. The axion coupling to electrons is set by a free rotation angle that also controls lepton flavor-violating axion and Higgs couplings. Moreover, all models have a trivial domain wall number (see also [26]), which evades the cosmological domain wall problem in scenarios when PQ symmetry is broken after inflation [27]. The Lagrangian is given by (apart from kinetic terms): where the first term comprises Yukawa couplings, the second term is that part of the scalar potential which only depends on the modulus of the two Higgs fields and the singlet, and the last term is needed to ensure that U (1) PQ is the only global symmetry. The Yukawa Lagrangian reads l L3 e R3hA 9 + y e 3a l L3 e RahA 10 + y e a3 l La e R3hA 11 + y e ab l La e RbhA 12 + h.c. (2.8) whereh i = iσ 2 h * i , a, b = 1, 2 runs over the first two fermion generations, and A 1...12 ∈ {1, 2} are parameters that define the Higgs field to which a given fermion bilinear structure couples to.
We consider four different structures in the quark sector, Q1-Q4, which are defined by the choice of (A 1...4 )(A 5...8 ): where we also indicate in 2+1 flavor space notation to which Higgs field the quark bilinears couple to. Together with the last term in Eq. (2.7) this choice fixes the PQ color anomaly coefficient N , the quark contribution to the electromagnetic anomaly coefficient E Q , and all couplings of the QCD axion to quarks C A,V q i q j in terms of the parameters tan β and ξ q P ij (to be defined below), which we summarize in Table 1.
The quark Yukawa Lagrangians of each model is combined with one out of the four following structures in the charged lepton sector, defined by (A 9...12 ), i.e. the Higgs to which a lepton bilinear couples in 2 + 1 flavor space  Table 1: Axion couplings in the four "nucleophobic" models Q1-Q4, see Eq. (2.9), as a function of the parameters ξ q P ij and c β ≡ cos β, s β ≡ sin β. Here E Q denotes the contribution of the quark sector to the electromagnetic anomaly coefficient E, to be added to the contribution from the charged lepton sector. In all models the domain wall number is trivial, 2N = 1. Table 2: Axion couplings in the four models EL1, E1R, E2L, E2R, see Eq. (2.10), as a function of the parameters ξ e P ij and c β ≡ cos β, s β ≡ sin β. Here E L denotes the contribution of the charged lepton sector to the electromagnetic anomaly coefficient E, to be added to the contribution from the quark sector.
This choice fixes the charged lepton contribution to the electromagnetic anomaly coefficient E L and the axion couplings to leptons C V,A i j , which we summarize in Table 2. Since each quark sector model can be combined with any charged lepton sector model, we have in total 16 different models, which we denote by e.g. "Q1E1L", which has the Higgs structure (2222)(1212)(1122), axion couplings to quarks and charged leptons as in Tables 1  and 2, and an electromagnetic anomaly coefficient E/N that is the sum of both sectors, E/N = E Q /N + E L /N = 8/3 in this example.
The couplings to quarks and leptons depend on the Higgs vacuum angle tan β and the parameters ξ f P ij , with f = u, d, e and P = L, R, which are defined by where V f P are the unitary matrices which diagonalize quark and charged lepton masses ac- . Unitarity implies the relations so apart from complex phases there are only two parameters ξ f P ii in each chiral fermion sector, which depend on the structure of quark and charged lepton masses.
By construction in all models the axion couplings to nucleons can be suppressed by choosing cos 2 β ≈ 2/3 and special values for the parameters ξ q P 11 , which are ξ u L 11 = ξ d L 11 = 0 (Q2), ξ u R 11 = ξ d R 11 = 0 (Q3), ξ d R 11 = 1 (Q1) and ξ u R 11 = 1 (Q4). The coupling to electrons is then controlled mainly by the parameter ξ e L 11 (E1L,E2L) and ξ e R 11 (E1R,E2R), and addressing the stellar cooling anomalies will generically correspond to ξ e L/R 11 = 0, 1. Therefore lepton flavorviolating axion couplings, which are controlled by ξ e L/R i =j , are a generic consequence of the cooling anomalies in these models (since two different ξ e ii are non-zero), while it is always possible to avoid quark flavor violation (i.e. having only one non-zero ξ q ii ). In the following we focus for simplicity on scenarios without quark flavor violation, obtained from choosing quark Yukawa matrices such that to very good approximation (neglecting small corrections necessary to reproduce the CKM matrix) which implies that all other ξ q P ii are vanishing (cf. Eq. (2.12)). This choice implies that models Q2 and Q3 give identical predictions for all axion couplings, while models Q1 and Q4 only give the same contribution for axion couplings to quarks in the 1st generation. In all models axion couplings to nucleons can be suppressed by taking c β ≈ 2/3 or equivalently t β ≈ 1/2 ≈ 0.7, which fixes the 1st generation quark couplings to identical values in all four models, C u ≈ 2/3, C d ≈ 1/3, while the axion couplings to heavy quarks are either close to 2/3 or 1/3. Axion couplings to leptons are controlled by tan β and the free parameters ξ e L ii [E1L,E2L] or ξ e R ii [E1R,E2R]. In the following we study the consequences of the above choices for ξ f P ii for the structure of the quark and charged lepton Yukawa sectors.

Quark Yukawa Sector
The quark Yukawa Lagrangian is given by where the structure of the couplings Y u 1,2 and Y d 1,2 depends on the model under consideration. These couplings have to be chosen appropriately, such that the diagonalization of the quark mass matrices reproduces i) the quark masses as singular values, ii) suitable left-handed rotations in order to obtain the correct CKM matrix, and iii) third rows of mixing matrices that match the parameter choice in Eq. (2.13), up to correction of small CKM angles. This gives the following parametric structure of Yukawa matrices, where λ ≈ 0.2 is of the order of the Cabibbo angle, and a numerical coefficient of order unity is understood in front of the λ, λ 2 , λ 3 entries in order to reproduce the exact values of the CKM matrix: Model Q1: (2.16) Model Q2: (2.17) Model Q3: (2.18) Model Q4:

Charged Lepton Yukawa Structure
The charged lepton Yukawa Lagrangian is defined as 20) where the structure of the couplings Y e 1 and Y e 2 depends on the model under consideration. We conveniently parametrize these couplings as follows: one matrix can be chosen to have only a 33-entry without loss of generality (for E1L and E1R Y e 1 and for E2L and E2R Y e 2 ), while the other matrix can be implicitly defined through charged lepton masses and rotations upon the relation Thus we take for models E1L and E1R and for models E2L and E2R (2.23) The six rotation angles in V EL and V ER are taken as free parameters, apart from the constraints imposed by the choice of ξ e L ii and ξ e R ii in the specific scenario.

Lagrangian in the Mass Basis
The Lagrangian in the mass basis (where the mass terms of fermions and Higgs fields are diagonal) is given by where the index k = 1, 2, 3 runs over the three neutral physical Higgs fields H 0 k = (h, H, A) while i, j are fermion flavor indices. The couplings depend on rotation angles in the scalar sector through α f k , β f k and on rotation matrices in the fermion sector through V ij , U ij , f ij . The couplings from the scalar rotations are given by and only depend on the two Higgs sector angles α and β which we treat as free parameters (obtained by a suitable scalar potential). The relation between the above Higgs mass eigenstates H 0 k and the Higgs doublet fields h 1,2 in Eq. (2.14) is spelled out in Appendix B. The couplings from the fermion rotations read and only depend on the unitary rotations V f P that diagonalize quark and lepton masses, and the Yukawa couplings The unitary rotations also fix the physical CKM and PMNS matrices V and U The explicit structure of the quark Yukawa couplings in Eqs. (2.16)-(2.19) determines the quark unitary rotations for a given model Q1-Q4, and gives for the final couplings from the fermion rotations the following matrices: Model Q1:

34)
Model Q4: where we have neglected sub-leading corrections suppressed by small mass ratios and/or by higher power of the Cabibbo angle.
Similarly the explicit structure of the lepton Yukawa couplings in Eqs. (2.22)-(2.23) gives for the final couplings from the lepton rotations the following matrices: Model E1L and E1R:

37)
Model E2L and E2R: which only depends on the third row of the rotation matrices in the left-and right-handed sectors. These can be expressed through the four real parameters ξ e L 11 , ξ e L 22 , ξ e R 11 , ξ e R 22 , up to complex phases, which we set to zero for simplicity: they do not enter h → τ and would only have a minor impact on the Higgs couplings to τ 's. Thus we finally have for models E1L and E1R and for models E2L and E2R where we remind the reader that 0 ≤ ξ e L ii , ξ e R ii ≤ 1. To summarize, we consider generalized DFSZ-like models that are defined by the Yukawa structure in Eq. (2.8). In total we have 4 different models in the quark sector (Q1-Q4) and 4 different models in the charged lepton sector (E1L, E1R, E2L, E2R). These scenarios depend on 7 parameters; the 4 angular parameters ξ e L 11 , ξ e L 22 , ξ e R 11 , ξ e R 22 controlling lepton flavor-violation, the vacuum angle β entering both axion and Higgs couplings, and two parameters that control overall decoupling: the Higgs mixing angle α for Higgs phenomenology and the axion decay constant f a for axion phenomenology. The corresponding predictions for the axion couplings to the SM are summarized in Tables 1 and 2, while the structure of neutral and charged Higgs couplings to fermions are given in Eq. (2.24), along with the predictions in Eqs. (2.29)-(2.36) and (2.39)-(2.40). In the following we study first the phenomenology of axion couplings, focussing on the possibility to address stellar cooling anomalies, before we turn to Higgs phenomenology, paying particular attention to flavor-violating decays of the SM-like Higgs.

Axion Phenomenology
In this section we provide the expressions for the low-energy axion couplings to nucleons, leptons and photons and use the present constraints in order to find the allowed regions in the parameter space. Finally we discuss the stellar cooling anomalies which hint to non-vanishing axion couplings to electrons and select a preferred region of model parameters, in particular yielding an upper bound on the axion decay constant.

Predictions for Axion Couplings
The axion has no relevant flavor-violating couplings to quarks, and the most important observables in the quark sector arise from axion couplings to nucleons constrained by SN1987A and neutron star cooling. In the different quark models the predictions for the low-energy axion couplings to protons C p and neutrons C n read [33] To derive the above expressions we have neglected Yukawa running effects (see Ref. [37][38][39][40][41][42]), which in DFSZ models are relevant only below the scale where the heavy Higgs doublet of the 2HDM decouples [42]. Since in our setup this mass scale is rather low (close to the TeV scale) we do not expect large corrections from top Yukawa running.
In the lepton sector instead there are flavor-violating couplings controlled by the parameters ξ e ii , and the most important observable in the lepton sector are µ → ea decays arising from the flavor-violating axion coupling C µe = |C A µe | 2 + |C V µe | 2 , and axion couplings to muons C µ and electrons C e , which are constrained by SN1987A and white dwarf cooling, respectively. The predictions for these couplings are Finally there is the coupling to photons, which is determined by the ratio of the electromagnetic and color anomaly coefficient that is fixed in each model, see Tables 1 and 2 The axion phenomenology thus depends on four parameters: the Higgs vacuum angle β, the relevant rotations in the charged lepton sector ξ e L 11 , ξ e L 22 (E1L, E2L) or ξ e R 11 , ξ e R 22 (E1R, E2R) and the axion decay constant f a . In the following we will simplify the notation and use the parameters ξ 11 , ξ 22 for all four models, where the chirality is left understood from the model under consideration. Thus we will distinguish only between E1 and E2 models, which represent (E1L,E1R) and (E2L,E2R) models upon the proper identification of ξ 11 , ξ 22 . The phenomenology is indeed essentially independent of the chirality of ξ ii , except bounds from LFV decays as we are going to see below.

Constraints on Axion Couplings to Photons
Axion couplings to photons are constrained mainly by the CAST experiment [43] and the evolution of horizontal branch stars in globular clusters [44], which at 95% CL require and exclude the region where The most stringent constraint arise for models with E/N = −4/3, and exclude f a ≤ 5.7 × 10 7 GeV, which is typically weaker than other constraints. However, the bound on photon couplings will be improved by helioscopes of the next generation, in particular the IAXO experiment, by about an order of magnitude [45].

Constraints on Axion Couplings to Nucleons
The couplings to neutron and protons can be bounded by the burst duration of the neutrinos observed in the SN1987A and reads [14] 0.61g 2 ap + g 2 an + 0.53g an g ap < 8.26 × 10 −19 , Since the nature of the axion emission is not completely understood and present simulations do not take into account all the relevant physics (for example the feedback from axion emission presumably modifies the bound when included in simulations [14]), these constraints should not be considered as a robust bound but rather as an indication. However, also the cooling of neutron stars provides information about the axion nucleon coupling [15][16][17][18][19] and yields upper bounds that are at least of the same order as the SN1987A bound. For example the observation of the neutron star in HESS J1731-347 [18] sets a strong bound on the axion neutron coupling as which translates to the excluded region Still, the limited understanding of the cooling of neutron stars and the lack of observational data suggest that these constraints should be taken with some grain of salt.

Constraints on Axion Couplings to Electrons
The axion coupling to electrons can be constrained by the shape of the white dwarf luminosity function, giving the 95% CL bound [47] where |g ae | = C e m e /f a , which excludes the region (E2L, E2R) (3.57)

Constraints on Axion Couplings to Muons
It has been shown recently [48,49] that SN1987A also constraints the axion couplings to muons excluding the region where (E2L, E2R) (3.60)

Constraints on LFV Axion Couplings
Finally, limits on the LFV coupling C V,A µe arise from constraints on the two-body lepton decay µ + → e + a, which depends on the chiral structure. For an isotropic decay the most stringent bound was provided by an experiment at TRIUMF, which sets the limit BR(µ + → e + a) < 2.6 × 10 −6 (at 90% CL.) [50]. If the decay has the same angular distribution as the SM, the weaker bound from the TWIST experiment [51] applies, BR(µ + → e + a) < 5.8 × 10 −5 (at 90% C.L.). The resulting bounds on the axion couplings (including a recast of the TRIUMF experiment) have been given in Ref. [52] and read for purely left-handed or right-handed couplings (at 95% C.L.) which excludes the region Note that this constraint is different for left-handed or right-handed lepton couplings, because the SM background on µ → e + invis. is purely left-handed and thus the LH models have weaker constraints. This is the only feature that allows to distinguish LH and RH models with axion physics (when ξ e L ii = ξ e R ii ).

Summary of Constraints from Axion Physics
The constraints on f a are summarized in Fig. 1 for the models Q3E1L/Q3E1R (left panel) and Q3E2L/Q3E2R (right panel) for ξ 11 = ξ 22 . The figures show the allowed regions in the tan β and ξ 11 = ξ 22 parameters space for a given f a . The constraints from white dwarf cooling (brown), neutron star cooling (orange) and SN1978A (red) allow only the region between the two solid curves (for f a = 10 8 GeV) and dashed curves (for f a = 10 9 GeV), while the constraints from µ → ea exclude the regions above the horizontal blue lines, for f a = 10 8 GeV (solid) and f a = 10 9 GeV (dashed), distinguishing between E1L/E1R and E2L/E2R models (for f a = 10 9 GeV the entire region of E1L and E2L models is allowed by µ → ea). We show only the Q3 model as a representative, because choosing a different quark model would only slightly change the bounds from supernovae and neutron stars, cf. Eqs. (3.52) and (3.54). Note that the bounds from white dwarfs are the same for the models E1L, E1R and for the models E2L, E2R, while the bounds from µ → ea are equal for the models E1L, E2L and the models E1R, E2R. Fig. 1 shows that it is not possible to have models with f a a few 10 8 GeV and ξ 11 = ξ 22 that satisfy all constraints simultaneously. This is due to the fact that the bounds from µ → ea and white dwarfs select opposite regions in ξ ii . However, when ξ 11 = ξ 22 → 0 the bound from µ → ea does not apply anymore, and it is possible to have (mild) cancellations in nucleon and electron couplings near tan β ≈ 0.7 and ξ 11 ≈ 2/3 (E1L,E1R) or ξ 11 ≈ 1/3 (E2L,E2R) such that f a 10 8 GeV is allowed. For f a 10 9 GeV the parameter space opens up, but still sizable regions of parameter space are excluded even in the limit ξ 22 → 0. While clearly all bounds can be relaxed by further increasing the axion decay constant, an upper bound on f a arises when the axion is responsible for explaining the stellar cooling anomalies.

Stellar Cooling Anomalies
Several hints for excessive cooling in stellar objects have been observed in the past years. These observations include (see e.g. Ref. [10][11][12][13]: i) the cooling efficiency of pulsating white dwarfs (WDs) extracted by the rate of period change; ii) the WD luminosity function, relating the WD distribution to their brightness; iii) the luminosity of the tip of the red giant branch (RGB) in globular clusters; iv) the ratio of the number of horizontal branch (HB) stars over RGB stars in globular clusters (R-parameter); v) the ratio of blue and red supergiants in open clusters.
The interpretation of these data in terms of BSM physics (such as millicharged particles, hidden photons and axions) was discussed in Ref. [53]. This analysis revealed that axions coupled to electrons are perfectly suited to address all anomalies, in contrast to the other candidates. A combined fit of the WD luminosity function (WDLF), WD pulsation and RGB stars is driven mainly by the WDLF and favors 2 a non-zero axion coupling to electrons [13] g ae = 1.6 +0. 29 −0.34 × 10 −13 , which translates to the region Since 0 ≤ ξ 11 ≤ 1, the fraction is always smaller than unity giving an upper bound on the axion decay constant f a 4 × 10 9 GeV at 1σ, corresponding to an axion heavier than at least 1 meV. We overlay the cooling hint region with the various constraints 3 and the IAXO sensitivity in Fig. 2, which shows the 1σ contours of the cooling hint, the reach of IAXO and IAXO+ as obtained in Ref. [13], and the bounds from supernovae and neutron stars for all the models, taking ξ 22 = 0. In particular, the solid green curves show the maximal reach of IAXO in its baseline configuration [28]. A more optimistic projection based on possible upgrades in the understanding of parameters of the magnet and the detectors is shown as a dashed line and it is labeled as IAXO+. On the other hand, the intermediate experimental stage called BabyIAXO [28] will have enough sensitivity to probe models with E/N = −4/3 (Q4E1) or 14/3 (Q1E2) up to f a ∼ 10 8 GeV, assuming that the interaction with electrons is suppressed. The presence of a sizeable interaction with electrons may increase the actual sensitivity, due to a larger axion production rate in the sun.
The cooling hint regions are the same for all the E1 or E2 models. The IAXO reach, on the other hand, depends also on the value of E/N and differs depending on both properties of the quark and lepton models, cf. Table 1-2. Fig. 2 shows that the area that satisfy both the cooling hint region and the different bounds dramatically depends on the type of the lepton model (E1 vs E2). In particular, the models E1 prefer a region where ξ 11 ∼ 1/ √ 2 (ξ 11 0.6) for f a = 10 8 (10 9 ) GeV, while the models E2 prefers a region with ξ 11 ∼ 0.3 (ξ 11 0.4) for f a = 10 8 (10 9 ) GeV. Some of these regions could be probed by IAXO or IAXO+. The cooling regions for f a ∼ 10 8 GeV can all be probed by IAXO. 4 On the contrary, IAXO will not probe  Figure 2: Projected potential of IAXO (green) and stellar cooling hint region (blue). For the indicated values of f a IAXO can probe the entire region between the solid green curves, while the dashed green contours show the region that can be tested by IAXO+ using a series of possible upgrades [55]. The blue regions are the 1σ cooling hint contours for a fixed value of f a . The SN1987A and the neutron star bound from HESS J1731-347 [18] are given by the red and yellow curves, respectively. the cooling hint region for f a = 10 9 GeV for the models Q1E1 (upper left panel), and Q2E2 and Q3E2 (middle right panel), it will partially probe the region for the models Q2E1 and Q3E1 (middle left panel) and Q4E2 (bottom right panel), and it will completely probe the models Q1E2 (upper right panel) and Q4E1 (bottom left panel). The models that will not be probed by IAXO will be challenged by its upgrade IAXO+.

Higgs Phenomenology
All models can be tested also with precision Higgs physics at the LHC, provided that the heavy Higgs bosons are not decoupled. In particular all models predict flavor-violating decays of the SM-like Higgs, h → eτ , h → µτ , h → µe, as well as modifications to the SM Higgs decays h → τ τ and h → µµ. All these observables can have large deviations from the SM prediction when the mass scale of the additional Higgs doublet is not too large. The magnitude of this deviation for these Higgs decays is then only limited by the Higgs coupling measurements at the LHC and flavor-violating charged lepton decays like µ → eγ. In the following we discuss these constraints in detail, before we combine Higgs and axion phenomenology in the next section.

Higgs coupling measurements
The effects of non-decoupled heavy Higgs bosons can be parametrized by c β−α ∝ v 2 /m 2 H , which vanishes in the decoupling limit m H → ∞. The value of c β−α affects all the couplings of the Higgs boson to SM particles. This implies that possible deviations of c β−α from zero are constrained by the Higgs coupling measurements at the LHC. In order to study the impact of those constraints on the models we perform fits using the results from the combination of the measurements of Higgs boson production and decay by the ATLAS collaboration Ref. [56]. The fits are performed in the so-called κ − λ framework where the parameters of the fit can be described in terms of measurements of the Higgs boson reduced couplings κ and λ as with (i, j) = (Z, g), (t, g), (W, Z), (γ, Z), (τ, Z), (b, Z) and κ i are the couplings of the SM-like Higgs normalized to their values in the SM, e.g. the couplings to fermions are where c h f are the couplings of the Higgs to the fermion f . In absence of exotic Higgs decays, and neglecting the Higgs decays to u, d and e, we have [56]: In the absence of new physics contributing to the effective couplings of the Higgs to gluon,  Table 3: Higgs effective couplings in the κ − λ framework from Ref. [56,57]. The root mean square (RMS) values are symmetrized in our fit procedure and are given in the table on the left (we opt for the conservative choice). The table on the right contains the correlation matrix amongst the seven free parameters defined in the text. Table 4: Approximate κ f parameters for the models Q1 -Q4.
photon and Zγ, we have the following scalings for the Higgs to gauge bosons couplings [56]: The relevant experimental results are summarized in Table 3, which can be found also in Table 12 of Ref. [56] and Fig. 30 of Ref. [57].
The relevant κ f parameters are (4.73)

Constraints from flavor-violating Higgs decays
The magnitude of the flavor-violating decays of the SM-like Higgs is controlled by the offdiagonal couplings y i j which are given by so that from Eq. (4.73) we find the following prediction for the branching ratios of the SM-like Higgs decaying to a pair of leptons: where we neglect tiny phase space effects. The total Higgs width is defined as [56] where Γ SM = 4.1 MeV [58]. Here the branching ratio B BSM denotes all the decay channels that are not present in the SM. The constraints on the decays h → i j are as follows: BR(h → eµ) < 6.1 × 10 −5 at 95% C.L. [59] , We note that both ATLAS [61] and CMS [60] have observed a slight excess in the search for h → eτ decays. We find that the weighted mean of the best-fit values reported by ATLAS and CMS is BR(h → eτ ) exp ≈ (0.09 ± 0.07)%. Even though the excess is just of the order of 1σ, it is intriguing that both experiments have seen it, and an update of these analyses in the future will be interesting for the present scenario, where large effects in this channel are possible and actually expected as we are going to show later.

Constraints from flavor-violating charged lepton decays
Flavor-violating Higgs couplings are also constrained by the flavor-violating leptonic decays µ → eγ, τ → eγ and τ → µγ [62,63]. This class of decays are induced by one-loop penguin diagrams, with internal neutral or charged Higgs bosons, and by two-loop diagrams with top, W or Z running in the loop attached to the Higgs that induce the flavour violation [64]. Neglecting contributions for the heavy neutral and charged Higgs bosons, the branching ratio for the decay i → f γ is where Γ µ is the initial state particle decay width,  [62]. The two-loop and one-loop contributions are comparable for the decay τ → eγ. On the contrary, the two-loop diagrams are the dominant contributions for the µ → eγ decay. The experimental constraints are given by [65,66] BR(τ → µγ) < 4.4 × 10 −8 at 90% C.L. , The bounds on τ → µγ and τ → eγ translate in rather weak bounds on |y τ µ | 2 + |y µτ | 2 < 1.6 × 10 −2 and |y τ e | 2 + |y eτ | 2 < 1.4 × 10 −2 , (4.81) respectively, assuming SM values for y τ τ and y tt [62]. On the other hand, the latest experimental bound on µ → eγ translates in This dramatically constrains the branching ratio for the h → µe decay to be 5 BR(h → µe) 3 × 10 −9 . Assuming instead that y eµ and y µe are zero, one can obtain a bound on If |y τ µ | ∼ |y eτ | the experimental bound on these couplings from µ → eγ is much stronger than that from τ → µγ and τ → eγ and results in BR(h → τ e) and BR(h → τ µ) below O(10 −4 ). Therefore, only one of the branching ratios BR(h → τ e) or BR(h → τ µ) can be large and close to the current LHC limits.

Summary of constraints from Higgs physics
In order to simplify the discussion we introduce the following notation: ξ ii denotes the rotation of the chirality that defines the model (i.e. ξ ii = ξ e L ii for E1L and E2L), whileξ ii denotes the rotation of opposite chirality (i.e.ξ ii = ξ e R ii for E1L and E2L). In contrast to the axion phenomenology, the Higgs phenomenology depends on the rotation angles for both chiralities, i.e. both ξ ii andξ ii . In the following we will consider only the case whereξ ii = ξ ii . This choice does not have a strong impact on the resulting phenomenology, sinceξ ii has only a subleading effect on Higgs phenomenology, as long as it is not larger than the specified one (i.e. as long asξ ii ≤ ξ ii ).
The couplings y ij of the SM-like Higgs boson depend on c β−α , t β and the rotations ξ ii ,ξ ii . Apart from the decoupling limit c β−α → 0, which we are not interested in, one can avoid strong bounds on y eµ and y µe from µ → eγ by setting ξ 11 or ξ 22 to zero, as evident from Eqs. (4.78), (4.79) and (4.74). As discussed in the axion phenomenology section, ξ 22 must be small to avoid the µ → ea constraint, while sizeable ξ 11 is preferred to explain the cooling hints. This leads to a prediction of sizeable BR(h → τ e), while BR(h → µe) and BR(h → τ e) are strongly suppressed to avoid the µ → eγ constraint, with a decay width proportional to Thus the magnitude of BR(h → τ e) is controlled not only by ξ 11 but also by c β−α , which cannot be arbitrarily large due to the constraints from the LHC Higgs coupling measurements. In Fig. 3 we show the 95% allowed region for the BR(h → τ e) (yellow) and from the Higgs coupling measurement fit (blue) as a function of c β−α and tan β, with ξ 11 = 2/3 (1/3) and ξ 22 = 0 for models E1 (E2). Note that the different lepton models do not affect the Higgs coupling measurements for the indicated choice of ξ ii . Furthermore, we show the future sensitivity of the branching ratios of 0.1%, that can be reached already at the Run 3 of the LHC, and 0.01% which is the goal of high-luminosity LHC (HL-LHC) [67]. The red band shows the region not excluded by SN1987A for axion models with f a 10 8 GeV, highlighting the most interesting region for axion phenomenology. Finally, the gray areas show regions where y t > 1, which indicates the potential loss of perturbativity. It is clear from this figure that perturbativity does not impose relevant constraints on the model parameter space consistent with Higgs couplings measurements (except for a small region for model Q3).
The result for BR(h → τ e) is the same for all the lepton models, cf. Equation (4.75). On the other hand, the bounds from the Higgs coupling measurements strongly depend on the different quark models. In particular, the allowed region for models Q2 and Q3 is much smaller than the one for models Q1 and Q4. |c β−α | for models Q2 and Q3 must be always below 0.1 while it can be about 0.2 for model Q4 or even 0.3 for model Q1 in agreement with the Higgs coupling measurement. This difference is due to the fact that, as seen from Table 4, in models Q1 and Q4 κ b and κ t are approximately equal to each other while in models Q2 and Q3 κ b and κ t are anti-correlated i.e. when one is enhanced (suppressed) the other one is suppressed (enhanced). Note that κ t controls the magnitude of the Higgs production cross-section in the gluon fusion while κ b dominates the Higgs total width. This implies that for a given value of c β−α deviations from the SM of the signal strengths for the Higgs decaying into electroweak gauge bosons (which are the best measured channels) are bigger in models Q2 and Q3 than in models Q1 and Q4 where the enhancement (suppression) of the Higgs production cross-section is partially compensated by the suppression (enhancement) of the Higgs branching ratios into gauge bosons (which is a consequence of the enhancement (suppression) of the Higgs total width).
For this reason, at present the strongest constraints on c β−α for the models Q1 and Q4 come from the bound on BR(h → τ e), and are almost independent of tan β. On the other hand, although the constraints from the Higgs coupling measurements for the models Q2 and Q3 are quite stringent, it is still possible to have BR(h → τ e) 0.1%. In the remaining part of the paper we focus on the models Q1 and Q4 since they are able to predict larger rates for flavor-violating Higgs decays while being consistent with the Higgs coupling measurements.
Finally, let us comment on the fact that as long as only Higgs phenomenology is considered, it is possible to choose ξ 11 = 0 to satisfy the µ → eγ constraint. In such a case ξ 22 can be nonzero and the contours of BR(h → τ e) in Fig. 3 would correspond to contours of BR(h → τ µ) for ξ 22 = 2/3 and ξ 11 = 0. However, as we have discussed in the previous section, the stellar cooling anomalies together with the constraints from µ → ea suggest that ξ 22 = 0 with ξ 11 free to vary.  3: Parameter space allowed by BR(h → τ e) (orange region) and Higgs signal strengths (blue region) for various quark models. The results for the Q2 model are the same as those for the Q3 model. We furthermore show the sensitivity for experiments that could prove branching ratios of 0.1% and 0.01%. The gray areas denotes regions where Higgs couplings y t > 1, indicating the potential loss of perturbativity. Also shown is the region where the SN1987A constraint is satisfied for f a = 10 8 GeV.

Interplay of Axion and Higgs Phenomenology
In this section we finally study the implications for Higgs physics in the parameter space where the axion can explain the stellar cooling hints. In this way we fix the axion decay constant f a and study the maximal possible deviations for the Higgs decays h → τ e, h → τ τ and h → µµ, which can be obtained for a suitable value of c β−α (while respecting all present constraints from precision Higgs physics).
We will show results only for the models Q1E1L and Q4E1L, since they are less constrained by the Higgs coupling measurements. A different choice of lepton model would affect only the values of κ µ , κ τ and the axion phenomenology, as discussed in the previous sections. In order to avoid the bounds from µ → ea we first fix ξ 22 = 0. As a consequence of this choice, we are left with three free parameters: ξ 11 , tan β and c β−α . For given ξ 11 and tan β we fix c β−α as the maximal value allowed by the Higgs coupling measurements, leading to a positive and a negative solution. We will show results only for the positive solution c β−α > 0 and comment on the difference with respect to the negative one.
In Fig. 4 we show contours of the maximal possible value BR(h → τ e) for the present (solid red) and future experimental reach (dashed red) consistent with the Higgs coupling measurements in the plane ξ 11 vs tan β for the models Q1E1L (left) and Q4E1L (right). These contours are overlaid with the 1σ region explaining the cooling hint (blue), and various existing and future constraints on axions for f a = 10 9 GeV. In particular, we show the bound from neutron stars (light blue) and the future reach of helioscopes with IAXO (solid green) and IAXO+ (dashed green). Notice that for f a 10 9 GeV there is no bound from SN1987A. In the region explaining the cooling hint, the branching ratio BR(h → τ e) can be as large as the current upper bound from CMS in both models. In the Q1E1L model, BR(h → τ e) can be maximal in the cooling hint region for any tan β above ∼ 0.8, while the corresponding region with maximal BR(h → τ e) in the Q4E1L model is characterised by tan β 2. A branching ratio of BR(h → τ e) ∼ 0.1 could be obtained in most of the cooling hint region, just decreasing the value of c β−α .
We furthermore show contours for the reduced Higgs couplings κ τ and κ µ . The deviations from the SM (which predicts κ τ = κ µ = 1) can be O(10)%, within reach of the HL-LHC [68]. Interestingly, in the cooling hint region κ τ and κ µ are very different from each other. In particular, the cooling hint region cannot easily accommodate a SM like value for κ τ and κ µ simultaneously. An uncertainty below 10% in the measurements of the κ parameters would be enough to probe these models. It is worth to notice that taking into account the constraints from neutron stars and the cooling hint region, the Higgs couplings to muons and taus could deviate from the SM up to 10%. This means that axion physics prefers the region tan β 1 and ξ 11 0.6, compatible with deviations in the µ and τ Yukawa couplings that may be observed already at the HL-LHC [68]. The main difference between the Q1E1L and Q4E1L models is the future reach of IAXO and IAXO+. The IAXO helioscope, in its base configuration, will be able to probe all the cooling hint region for the model Q4E1L, independently of the value of tan β or ξ 11 . On the other hand, the cooling hint region for the model Q1E1L could be partially probed by an advanced configuration of IAXO. Therefore, an interplay between axion and collider searches is needed in order to fully rule out this model.
The results for c β−α < 0 are similar to the one in Fig. 4. There are two differences that affect the results. On one hand, the change in sign affects the observables that are linear in c β−α , such as κ µ and κ τ . In particular, in Fig. 4, the contour lines for κ µ/τ = 0.9 become contour lines for κ µ/τ 1.1. On the other hand, the difference in absolute value, due to the The blue region is preferred at 1σ by stellar cooling hints, the future helioscopes IAXO (IAXO+) will probe the whole region except the one between the green (dashed green) curves, and the region between light blue curves satisfies the neutron star bound for f a = 10 9 GeV.
fact that the Higgs coupling measurements allowed region is not symmetric in c β−α (see Fig. 3), influences both the observables that are linear or quadratic in c β−α , such as the BR(h → τ e).
The results for f a = 10 8 GeV are shown in Fig. 5 for the models Q1E1L (left) and Q4E1L (right). The contours of the Higgs observables are the same as in Fig. 4 while the cooling hint region and the other constraints on the axion are modified. In this case the constraint from SN1987A (brown curve) enforces tan β to be in a small range between about 0.6 and 0.8. Although the cooling hint region is much smaller in this case it is still possible to obtain BR(h → τ e) as large as 0.22%, i.e. the current upper bound from CMS [60]. Interestingly, in the cooling hint region for f a = 10 8 GeV consistent with the constraint from SN1987A the maximal deviation of κ µ from the SM always exceeds 15%, while the deviation of κ τ can be up to 10%. Therefore, the f a = 10 8 GeV case gives very sharp prediction for the pattern of the Higgs couplings to muons and taus. The projected sensitivity of ATLAS at the HL-LHC is around 7% for κ µ and 3% for κ τ [68] which should be enough to test these models for f a = 10 8 GeV. Furthermore, this scenario will be easily probed by IAXO in its default configuration. 6 Let us emphasize that predictions for the Higgs couplings in our models substantially differ from Type-I and Type-II 2HDMs in which the Higgs couplings (normalized to the SM) to all down-type quarks and charged leptons are the same. In particular, in the Q4E1L model with c β−α > 0, κ b is smaller than κ µ (κ τ ) by about 30% (20%) in the cooling hint region consistent with the neutron star bound for f a = 10 8 GeV. Therefore, it will be easy to experimentally distinguish this model from Type-I and Type-II 2HDMs. 7 We also note that in this region BR(h → µµ) could be up to ∼ 60% larger than the SM prediction, while for f a = 10 9 GeV one can have deviations of order ∼ 70%. This is a combined effect from simultaneously enhancing κ µ and suppressing κ b . In this case a significant excess in the muon decay channel may be observed already in the Run 3 of the LHC.

Conclusions
In this article we have explored the correlation of axion and Higgs phenomenology in variant axion models with a light second Higgs doublet. We restricted to a particular class of "nucleophobic" DFSZ-like models that allow to avoid the stringent constraints from SN1987A and neutron star cooling by suppressed couplings to nucleons, while couplings to electrons can be sizable and allow to address various stellar cooling anomalies. All axion couplings are fixed in terms of three relevant parameters, the axion decay constant f a , the Higgs vacuum angle tan β and a free angular parameter ξ 11 that controls lepton flavor-violation. A compact region in this parameter space is selected by the stellar cooling hints, while imposing the astrophysical bounds on nucleon couplings and perturbativity, see Fig. 2. This in particular restricts the values of the axion decay constant to values below about 4 × 10 9 GeV, which corresponds to axion masses of the order of few meV. Large parts of this parameter space will be probed by next-generation helioscope experiments such as IAXO. Such heavy axions can still account for the observed Dark Matter abundance when produced by the decay of strings and domain walls in scenarios when PQ is broken after inflation, up to roughly m a 4 meV [55], although the abundance suffers from significant uncertainties [69]. Another production mechanism that also works for heavier axions is parametric resonance from oscillations of the Peccei-Quinn symmetry breaking field [70], which can yield the correct DM relic abundance up to axion masses of roughly 60 meV. In both scenarios it is crucial to avoid stable domain walls by having a trivial domain wall number [27], which indeed is realized in the class of DFSZ models considered here (see also Ref. [26]).
Also the Higgs sector depends on the Higgs vacuum angle β and the free angular parameter ξ 11 controlling lepton flavor-violation, in addition to the mass of the second Higgs doublet that enters Higgs couplings through the angle cos(β − α). While previous studies of similar models always decoupled the additional Higgs doublet, corresponding to the alignment limit when cos(β − α) → 0, here we analyzed in detail the resulting Higgs phenomenology when the deviation from alignment is as large as allowed by present experimental Higgs data.
In this way we can correlate axion and Higgs phenomenology, since the cooling hints essentially fix the Higgs couplings to leptons as a function of cos(β − α). Of particular relevance are the Higgs couplings to muons and tau leptons and the LFV Higgs decay h → τ e, which will be probed with upcoming LHC data for precision Higgs measurements. Maximal values of these observables can be predicted by taking the heavy Higgs doublet as light as possible consistent with present data, or equivalently maximizing cos(β − α). Our results are summarized in Figs. 4 (for f a = 10 9 GeV) and Figs. 5 (for f a = 10 8 GeV), which show that lepton flavor-violating Higgs decays h → τ e can saturate the current experimental bound BR(h → τ e) = 0.22%, while deviations from the SM prediction for BR(h → µµ) can be as large as 70% in the parameter region where the axion can explain the stellar cooling hints.
The QCD axion models that we considered in this article to address the stellar cooling anomalies might therefore be testable not only by future helioscopes but also by precision Higgs data, highlighting the interplay of dedicated axion searches with IAXO and precision tests of the SM Higgs sector at the LHC.
2} is a free parameter that selects to which Higgs field the respective fermions couple to. This gives twelve constraints, which determines all fermion charges in terms of Higgs charges Compared to the SM, the above Yukawa Lagrangian has an extra U (1) 2 h × U (1) φ global symmetry that needs to be broken to a single U (1) PQ factor by adding two couplings in the scalar sector. Since U (1) PQ = U (1) φ , we need at least one coupling of φ. On the renormalizable level we can couple h † 2 h 1 to an operator O 1 ∈ {φ, φ * , φ 2 , φ * 2 }. This constrains the charges of the Higgs fields h 2 in terms of the h 1 charge X 1 and a free parameters B that can take only discrete values with the possible values The scalar potential is assumed to induce the Higgs and singlet vevs and one is free to make a field rotation such that only one linear combination of Higgs doublets takes a vev v = 246 GeV The Goldstone boson eaten up by the Z-boson resides in this linear combination, and with it is given by The anomalous U (1) PQ current reads where we omitted the fermionic part. This current creates the axion according to which defines the axion as the linear combination with the PQ breaking scale We need to impose that the axion is orthogonal to the Goldstone eaten up by the Z, which gives the condition The rotation matrix O can be constructed simply by taking v as first row and then flipping pairwise two entries with minus signs to be orthogonal. In particular one has which is the only input in the orthogonality condition Using the axion definition in Eq. (A.11), it is easy to verify that mass terms and axion-fermion couplings arise from replacing the neutral Higgs field components by Therefore axion-fermion couplings can be removed from the Yukawa Lagrangian by the flavordiagonal fermion field redefinition (a local PQ transformation acting only on fermions) Since this transformation is anomalous, it generates axion couplings to gauge field strengths, and since it is local it modifies the fermion kinetic terms. The anomalous couplings are given by with the dual field strengthF µν = 1 2 µνρσ F ρσ , 0123 = −1 and the color and electromagnetic anomaly coefficients The kinetic terms give the following axion-fermion couplings in the flavor interaction basis L = ∂ µ a v PQ u i γ µ C q ij P L +C u ij P R u j + d i γ µ C q ij P L +C d ij P R d j + e i γ µ C l ij P L +C e ij P R e j , (A.24) , C e ij = −X A 10 δ ij + diag(0, 0, X A 10 − X A 9 ) . In this basis the axion-fermion couplings are given by Now we adopt the standard convention for the axion decay constant f a = v PQ /(2N ), and write the Lagrangian as where L = q, l and R = u, d, e. The flavor structure is controlled by the matrices (f = U, D, E; P = L, R) The absolute values of these matrices depends only on two independent real parameters in each fermion sector. Notice that the parameter B only enters the domain wall number N DW = 2N , since it is equivalent to the charge normalization of φ and thus drops out from all axion couplings which only depend on charge ratios. In flavor-universal DFSZ models this setup reduces to

B Higgs Mass Eigenstates
Starting from Eq. (2.14), we parametrise the Higgs fields h i as and change into the Higgs basis Φ i with which yields the charged-Higgs boson H ± and the pseudo scalar field A .

(B.3)
Finally, the physical Higgs fields h, H are related to the fields R 1 , R 2 of the neutral Higgs decomposition in Eq. (B.1) via the orthogonal rotation which is obtained by diagonalising the 2 × 2 mass block of the CP-even states R 1 , R 2 .