Anomaly-free axion dark matter in three Higgs doublet model and its phenomenological implications

We study phenomenological implications of an axion that arises as a pseudo Nambu-Goldstone boson due to the spontaneous breaking of anomaly-free global flavor symmetry. One interesting possibility for such anomaly-free axion to explain dark matter (DM) is when it has a mass of order keV and an intermediate scale decay constant, since it can be explored through direct search experiments, X-ray observations, various stellar cooling processes, and the misalignment mechanism naturally explains the DM abundance. As a concrete renormalizable model of such axion, we consider an extended Higgs sector with global flavor symmetry, which consists of three Higgs doublet fields and three singlet Higgs fields with $U(1)_{\rm B-L}$ charges. We identify viable parameter regions that satisfy theoretical bounds on the Higgs potential and various experimental limits on this model, and evaluate the mass spectra of the axion and extra Higgs bosons. We find that even an anomaly-free axion can generally couple to photons through mixing with CP-odd Higgs, and that its strength depends on the vacuum expectation values of the Higgs doublets as well as the axion mass. As a result, the ratios of the vacuum expectation values of the Higgs doublets are tightly constrained to satisfy the X-ray constraints. We show the favored parameter region where axion DM explains the XENON1T excess. We also demonstrate that the axion-electron coupling is correlated with the extra Higgs boson masses and mixing angles for CP-even Higgs bosons. Thus, if the axion is detected in future observations, the extra Higgs boson masses and the coupling of the standard model-like Higgs boson with the weak gauge bosons are restricted. This is a good example of the synergy between searches for the axion DM and the BSM around the electroweak scale.


Introduction
The Standard Model (SM) has been scrutinized to a very high degree of accuracy by various experiments. However, it does not include a candidate for dark matter (DM), which accounts for about a quarter of the current energy density of the universe. For this reason, DM is considered to be solid evidence of physics beyond the SM.
The identity of DM is still unknown, and one plausible explanation is that it is composed of unknown particles. Various candidates for particle DM have been proposed, and one of the most promising ones is the axion [1][2][3][4]. The axion is a pseudo Nambu-Goldstone (NG) boson that appears due to the spontaneous breaking of the global U(1) symmetry. As such it can easily satisfy the stability requirement for DM because it has naturally a light mass and very weak interactions with the SM particles, if the symmetry breaking scale is sufficiently high. Because of its light mass, the axion is expected to be generally away from its potential minimum in the early universe, so it starts oscillating around the potential minimum when the Hubble parameter becomes roughly equal to the curvature of the potential. The oscillation energy is a natural explanation for the cold axion DM. This is known as the misalignment mechanism [5][6][7]. The advantage of axion DM is that it can naturally explain the stability and generation of DM. See Refs. [8][9][10][11][12][13][14] for reviews.
A number of experiments are being conducted to explore the axion DM, and various limits have been placed on its interactions with the SM particles. In particular, axions coupled to photons are often referred to as axion-like particles (ALPs). Light axions are easily produced in hot stellar interiors, which affects stellar evolution by carrying energy outside the star. For instance, the interaction between the axion and photons is restricted by horizontal branch stars [15], and the interaction between the axion and electrons is restricted by observations of the tip of red giant branch [16,17] and white dwarfs [18][19][20]. It is also known that when the axion mass exceeds O(10) eV, the observational limit on the UV and X-ray photon flux produced by axion decay becomes very tight. Assuming the anomalous coupling with photons, the decay constant of the axion of keV-scale mass must be above the GUT scale. Such axion DM with the GUT-scale decay constant is far beyond the sensitivity expected in current and near-future direct search experiments, and its contribution to the stellar cooling process is negligibly small.
The axion couplings to gauge bosons are model-dependent, and various theoretical possibilities have been discussed. 1 Among them, there is an interesting axion model that escapes the above tight limits from X-ray observations, which was proposed in Ref. [24] as the anomaly-free axion model. For simplicity let us consider a case where the axion is only coupled to leptons. Then, the axion coupling to photons is caused by one-loop diagrams in which the leptons are running in the loop. The anomalous coupling is obviously absent if the charge assignment on the leptons is such that the electromagnetic anomaly is canceled, but the effect of threshold corrections allows the axion to have a suppressed, but finite, coupling to photons. On the other hand, the interaction between the axion and the leptons is not particularly suppressed. This lack of anomalous coupling to photons significantly relaxes the severe limits from the X-ray observations, and allows the axion with a decay constant of intermediate scales to become the dominant component of DM, and to be explored in direct experiments and by cooling processes in stellar objects. Especially in the case of keV-scale mass axions, the misalignment mechanism can naturally explain the observed DM abundance if their decay constant is on the intermediate scale.
Recently, the hint of an excess in the electron recoil events observed in the XENON1T experiment that cannot be explained by the previously known background has attracted much attention [25]. The possibility that the excess event is due to a tritium contribution cannot be ruled out, but various DM candidates have been considered to explain it. Among them, the anomaly-free axion with keV mass is very interesting because it can explain the excess in the XENON1T experiment and at the same time it can explain stellar cooling anomalies [26].
With the above motivation, the purpose of this paper is to investigate a possible UV completion of the anomaly-free axion with keV-scale mass and the decay constant at an intermediate scale. In particular, we consider the three Higgs doublet model (3HDM) as such UV completion. The set-up based on the 3HDM was proposed in Refs. [24,26], but the detailed analysis of the Higgs sector and various phenomenological bounds have not been investigated so far. The 3HDM was first proposed in the context of CP violation in the Higgs sector in the pioneering works [27,28] where discrete symmetries were imposed for natural flavor conservation [29]. In the 3HDM, possible interaction terms and the number of parameters crucially depend on the imposed symmetry. Possible symmetry groups implemented in 3HDMs are surveyed in Refs. [30][31][32][33][34][35][36][37][38]. Furthermore, phenomenological studies are performed in the 3HDM with CP4 symmetry [39,40], Z 3 symmetry [41][42][43][44], Z 2 symmetries [45][46][47][48][49][50][51], Sp(6) symmetry [52], U(1)×Z 2 flavor symmetry [53], and S 3 symmetry [54]. Also, eight extra Higgs bosons are predicted in addition to the SM-like Higgs boson with the mass of 125 GeV. For this reason, the phenomenology of 3HDM is very rich, and we will focus on the symmetry and matter content that would fit with the anomaly-free axion. The anomaly-free axion has also been studied from other phenomenological aspects such as lepton-flavor violation [55][56][57] as well as inflation [58]. This paper is organized as follows. In Sec. 2 we give the set-up of the 3HDM, and we study the mass spectrum of the axion and CP-odd Higgs bosons in Sec. 3. We take into account various theoretical and experimental bounds on the model parameters and identify the viable parameter region in Sec. 4. The implications for the direct DM search experiments as well as future X-ray observations are studied in Sec. 5 and Sec. 6, respectively. In particular, we will show that the anomaly-free axion generally acquires a coupling to photons via the mixing with CP-odd Higgs bosons, and we discuss its implications for the X-ray constraints on the model parameters. The last section is devoted to conclusions.

Three Higgs doublet model with B − L Higgs bosons
We first provide the set-up of the 3HDM with B − L Higgs fields, which is essentially same as the one proposed in Ref. [24], but the matter content of the B − L Higgs sector is slightly simplified. We impose a global U (1) F flavor symmetry on leptons in such a way that its electromagnetic anomaly vanishes, and the anomalous coupling of the axion to photons is canceled if the mixing effect is negligible. In this case, the axion-photon coupling arises from the threshold corrections in the lepton loop diagrams. The effect of the mixing between the axion and the CP-odd Higgs will be discussed later in this paper. We introduce three Higgs doublet fields φ i (i = 1, 2, 3), and three Higgs singlet fields with U (1) B−L charge +2, S 0 , S 1 , S2, where the subscript of S (not φ) denotes the flavor charge. The assignment of the U (1) F flavor charge q for these Higgs fields and left-handed (right-handed) lepton fields L ( R ) ( = e, µ, τ ) are presented in Table 1. Supposing L e and e R are charged under the U (1) F symmetry, there are two possible combinations for the charge assignment of the right-handed leptons, and we call them Type-A and Type-B, respectively. Quarks are assumed to be neutral under the U (1) F symmetry.
While soft breaking terms for the U (1) F symmetry are introduced in the Higgs potential, U (1) B−L symmetry spontaneously breaks after the B − L Higgs singlet fields acquire vacuum expectation values (VEVs). In the following, we describe the Higgs potential and the Yukawa Lagrangian. Hereafter, we use shorthanded notation for trigonometric functions, s θ ≡ sin θ, c θ ≡ cos θ, and t θ ≡ tan θ as needed.
As we will see later, the charged lepton Yukawa is already diagonalized in the flavor basis. As pointed out in Ref. [24], if we further introduce three right-handed neutrinos with flavor charges (1, −1, 0), and if the VEVs of the B − L Higgs are of similar magnitude, then light neutrino masses and the large mixing angles can be realized by the seesaw mechanism [59][60][61][62]. In the following, we will focus on the correlation of the axion and heavy Higgs bosons.

Higgs potential
The Higgs potential we consider can be divided into three parts as where V 3HDM (V B−L ) denotes the potential for the three Higgs doublets φ i (B − L Higgs S j ), and V I corresponds to interaction terms between φ i and S j . The potential for the three Higgs doublets, V 3HDM , is given by where the soft breaking terms are given by 3) The first term m 2 12 (φ † 1 φ 2 ) breaks U (1) F down to its Z 6 subgroup, while the other terms In this paper, we focus on the pattern of the symmetry breaking of U (1) F → Z 6 , taking m 2 12 = 0, m 2 13 = 0 and m 2 23 = 0. While we only introduce m 2 12 , the terms similar to m 2 13 and m 2 23 are dynamically generated from the interaction terms V I as we will see shortly. Without loss of generality, one can take m 12 and λ 10 to be real by using the phase degrees of freedom of φ 2 and φ 3 . We also note that the above potential corresponds to the one in the Z 3 invariant 3HDM with the U (1) F symmetry breaking terms other than m 2 12 set to 0 in Ref. [44]. Under the U (1) F and U (1) B−L symmetries, the Higgs potential in the B − L Higgs sector and the interaction terms between the Higgs doublet fields and B − L Higgs fields are given by 2 The dimensionless parameters κ 12φ1φ3 and κ2 1φ2φ3 are taken to be real while one of them is generally complex.
In the phase after the electroweak symmetry breaking (EWSB), the component fields of the Higgs doublets can be given by where v k represents the vacuum expectation value (VEV) of the CP-even components. We parameterize the VEVs as [42] v The electroweak VEV is then obtained by On the other hand, assuming that U (1) B−L is already spontaneously broken, we parameterize the B − L Higgs fields as Our purpose is to induce the axion mass by the soft breaking terms in V3HDM, namely the Higgs sector related to the electroweak symmetry breaking. Hence, we do not consider the soft breaking term µ 12 S † 1 S2. in VB−L. In fact this is ensured by considering only the explicit breaking that preserves the Z6 subgroup.
whereã denotes the pseudo NG boson associated with the breaking of the U (1) F symmetry, Q j denotes the U (1) F charge of S j (see Table. 1), and f a denotes the decay constant for thẽ a. Note that f a is generally expressed by a linear combination of the VEVs of the B − L Higgs and φ 1,2 , but we drop the contribution of φ i by assuming the B − L breaking scale is much larger than the electroweak scale. Specifically, we assume later that v S j and f a are of O(10 10 -10 12 ) GeV. We also note that the definition of the pseudo scalarã is the one for the specific case of v S 1 = v S2 , but we use it for simplicity. Although there are the other NG boson for the U (1) B−L symmetry and a heavy CP-odd Higgs boson in the B − L sector, they are integrated out in our study. We also integrated out the three CP-even Higgs bosons ρ j after we set their tadpole conditions.
Setting the B − L Higgs fields equal to their VEVs in Eq. (2.5), one obtains terms similar to m 2 13 and m 2 Note that, although these interaction terms do not explicitly violate the U (1) F flavor symmetry, at least one of them must be nonzero for the axion to obtain nonzero mass. This is because the B − L sector where the axion resides and the EW sector where the U (1) F symmetry is explicitly broken are completely separated, otherwise. Hence we call m 2 13 and m 2 23 the soft breaking parameters in the following. Hereafter, we introduce the following rescaled soft breaking parameters, and choose them as input parameters. The replacement of Eq.(2.11) makes the expressions for the quartic couplings λ 1−9 simple as shown in Appendix A. Although the rescaled mass parameters are introduced for convenience, this parameterization may obscure the scale of the portal couplings. Before closing this section, we discuss typical size of the portal couplings. It can be estimated from the minimization conditions for Higgs doublet fields φ i (i = 1, 2, 3), which are given in Eq.(A.1)-(A.3). The mass parameters m 2 ii (ii = 11, 22, 33) should be typically at the EW scale. This follows that the order of all the portal couplings should be where the subscript X can be 12φ1φ3,21φ2φ3, or mφn (m, n = 0, 1,2). As diseased in Sec.3, the corresponding mass scale of the axion is 1 keV if f a ∼ 10 10 GeV and m 12 ∼ 100 GeV. Since all the portal couplings are small, our model is natural in the sense of 't Hooft [63] (also see the discussion of naturalness in the axion model, e.g. Ref. [64]).

Definition of physical states
In the potential obtained by integrating out heavy degrees of freedom, we have the six charged Higgs bosons w ± k , the four CP-odd Higgs bosons z k ,ã and the three CP-even Higgs bosons h k . They are related to the mass eigenstates by orthogonal transformations.
For the charged Higgs bosons, the physical states are given by where R + is a mixing matrix, G ± are NG bosons eaten by the W ± bosons, and H ± 1 and H ± 2 denote the charged Higgs bosons. The mixing matrix R + is given by a product of two rotation matrices, where they can be parameterized as and (2.16) Here we have followed the parametrization of O β given in Ref. [42]. The rotation matrix O β corresponds to the transformation from the original basis (φ 1 , φ 2 , φ 3 ) into the Higgs basis [65][66][67][68]. In other words, by the rotation matrix O β , the NG bosons G ± are identified, and the remaining two states are transformed into (H ± 1 , H ± 2 ) by O γ + . For the CP-even Higgs bosons, the physical states (H 1 , H 2 , H 3 ) are obtained by where the mixing matrix R S is given by (2.18b) We identify H 1 as the SM-like Higgs boson with the mass of 125 GeV and H 2,3 as the additional CP-even Higgs bosons.
For the CP-odd Higgs bosons, the physical sates are given by where R P is the mixing matrix, G 0 is the neutral NG boson eaten by the Z boson, a is the axion corresponding to the spontaneous breaking of U (1) F , and A 1 and A 2 are additional CP-odd Higgs bosons. The mixing matrix is expressed by Here we introduce the 3 × 3 orthogonal matrix R A . While we numerically derive R A in the following numerical calculations, the matrix R A can also be parameterized by introducing rotation matrices, The 3×3 matrices O γ i are defined by the replacement of α i → γ i (i = 1, 2, 3) in Eq. (2.18b). For later convenience, we define fields in the Higgs basis as (2.21)

Mass matrices for the Higgs sector
The mass matrix of the charged Higgs bosons H ± 1 and H ± 2 in the Higgs basis is given by where, Performing the further rotation for (B C ) 2 with O γ + yields the masses of H ± 1 and H ± 2 and the mixing angle γ + as . (2.25) The quartic couplings λ 7 , λ 8 and λ 9 are expressed in terms of these three physical param- Similarly, the masses of the CP-even Higgs bosons are represented by where the mass matrix for the CP-even Higgs bosons in the basis (h 1 , h 2 , h 3 ) T , M 2 S , is given by With Eq. (2.26), the six potential parameters, λ i (i = 1 − 6) can be expressed by the physical parameters α j and m H j (j = 1, 2, 3). Finally we express the mass matrices for the CP-odd Higgs bosons A 1 , A 2 and the axion by the mass matrix in the Higgs basis as where M P denotes the mass matrix in the basis of (G 0 , z 2 , z 3 , a ) T . The elements of B 2 P are given by, (2.28g) The matrix B 2 P are fully diagonalized by the orthogonal transformation where we have assumed the mass ordering, m a < m A 2 < m A 1 . For the CP-odd Higgs sector, the masses are chosen as the output parameters 3 4 .
Using the tadpole conditions and the mass formulae for the Higgs bosons, the original 17 parameters are replaced by the following physical parameters. (2.31) The quartic coupling constants λ i (i = 1 − 9) in the original Higgs potential can be given in terms of above the 16 parameters except for f a . The analytical formulae are given in Appendix A.  Table 3. Coefficients of the Yukawa coupling for the CP-even Higgs bosons. Table 4. Coefficients of the Yukawa coupling for the CP-odd Higgs bosons.

Yukawa Lagrangian and kinetic terms
Under the charge assignment in Table.1, the Yukawa Lagrangian for Type-A and Type-B is commonly written by where Y u and Y d correspond to 3 × 3 Yukawa coupling matrices for quarks, y e , y , and y ( , = µ or τ ) are the Yukawa coupling constants for the leptons. Due to the charge 3 If one expresses RA by Eq. (2.20c), one has the six physical parameters, γi (i = 1, 2, 3), mA 1 ,mA 2 and ma. Since the mass matrix B 2 P are determined by the five potential parameters (λ10, m12, m 13 , m 23 , fa) for the fixed EW VEVs, one cannot choose all physical parameter as inputs. 4 While we integrate out heavy Higgs bosons in the B-L sector, we mention radiative corrections to the masses for a and Ai from the B-L Higgs bosons ρm (m = 0, 1,2). Focusing on the portal interaction with κ 12φ1φ3 , we can estimate the one-loop corrections to m 2 a and m 2 where we assume κ 12φ1φ3 = v 2 /f 2 a and take fa ∼ 10 10 GeV, v1,3 ∼ v, (RP )11,13 ∼ 0.1 and (RP )44 ∼ 1. Thus, the radiative corrections to the masses from the B − L sector would not be significant. assignment of the U (1) F charge, there are no off-diagonal components of the lepton Yukawa matrix. Depending on the charge assignment, and are given by In either case, the Yukawa Lagrangian above can be expressed in terms of the mass eigenstates of the Higgs bosons as where q l denotes the effective U (1) F charge for the lepton l, and m l is the lepton mass. For instance, q e = 1 − (−2) = 3. We use Eq. (2.35) in the numerical calculations in Sec. 5 and 6. The analytical expressions for the axion coupling with charged leptons are given by Appendix E. We note that the difference between Eq. (2.35) and Eq. (2.36) comes from the breaking of U (1) F symmetry (i.e., m 2 12 ) as well as the mixing with the CP-odd Higgs. The coefficients, ξ Tables 2, 3 and 4, respectively. Hereafter, we focus on the Type-B Yukawa lagrangian. 5 From the kinetic terms of the Higgs doublet fields, one can derive the gauge-gauge-Higgs couplings as where the scaling factors κ H i V are given by (2.40)

Alignment limit
Current measurements of the coupling constants of the 125 GeV Higgs boson at the LHC Run 2 [69,70] show that the properties of the discovered Higgs boson are similar to those predicted in the SM. Theoretically, this situation can be realized in the so-called alignment limit, where the CP-even Higgs boson with the mass of 125 GeV has the same tree-level couplings as the SM. The alignment limit in two Higgs doublet models (2HDMs) was discussed in e.g., Refs. [71,72], and, in the context of 3HDM, the analytical condition for the limit was systematically derived in Ref. [42]. Symmetries for the Higgs potential that naturally lead to the alignment limit were discussed in Refs. [52,73,74].
Since we identify H 1 with the SM-like Higgs boson, the alignment limit requires κ H 1 V = 1. This can be reduced as 41) and this equation yields the condition for the alignment limit: As can be easily seen in Table. 3, ξ f H 1 becomes unity in this limit. One can also see that the mixing matrix for the CP-even Higgs bosons is expressed as when Eq. (2.42) is satisfied. This means that, similar to the charged Higgs bosons, the SM-like Higgs boson H 1 is diagonalized by the rotation O β and the remaining two CP-even Higgs states are transformed into the mass eigenstates (H 2 , H 3 ) by the rotation O α 3 . In this connection, the Yukawa couplings for (H 2 , H 3 ) have the same structure as the (H ± 1 , H ± 2 ) in the alignment limit, i.e., where the minus sign comes from the different convention in the rotation matrices (see Eqs. (2.15) and (2.18b)).

Axion-photon coupling
We here give the analytical expression for the axion photon coupling, g aγ . It is defined by the following effective Lagrangian, The axion photon-coupling g aγ can be derived from the amplitude for the axion decay into two photons. Using the Yukawa interactions of the axion, i.e., Eq. (2.35), it can be calculated as where p 1,2 denote the external momenta of the photons, and µ (p 1 ), and ν (p 2 ) are the polarization vectors. The coefficient of the tensor products in M corresponds to the axionphoton coupling g aγ , which is given in terms of the Passarino-Veltman function C 0 [75]. The approximate formula in the case of m a m is given by Using this approximate formula, we arrive at the axion-photon coupling, The first term corresponds to the anomalous coupling and the second term ξ a m 2 a /(12m 2 ) corresponds to the threshold corrections by the lepton loop diagrams [24] (see also Ref. [76]). The parameter ξ a contains the effect of the mixing between axion and CP-odd Higgs bosons and the breaking of U (1) F symmetry. We note that in the absence of these two effects, ξ a can be replaced with the charge q thorough the Eq. (2.36). If this is the case, the anomaly terms are completely canceled out and one obtains a further reduced expression of g aγ , which was used in the analysis of Refs. [24,26]. However, we would like to emphasize that in the 3HDM model for the anomaly-free axion, the axion-photon coupilng generically receives contributions from the mixing with the CP-oddd Higgs and the U (1) F breaking. We will study its implications for the X-ray constraints later in this paper. For later convenience, we formulate a relation between g aγ and g ae . From Eq. (2.48), it can be written by where ∆ is defined by Analytical expressions of the axion-lepton couplings g a are presented in the limit of m 12 = 0 in Appendix E. The parameter ∆ quantifies the deviation from the limiting case (2.49) normalized by the typical anomalous coupling. In fact, in the absence of the mixing with the CP-odd Higgs and the U (1) F symmetry breaking, i.e.,t β 1 = 1 and m 12 = 0, the anomalous coupling (the first term in ∆) vanishes, and we are left with tiny threshold corrections from heavier charged leptons, ∆ = m 2 a q µ /m 2 µ + q τ /m 2 τ /(12|q e |) 2.5 × 10 −12 at m a =1 keV. In practice, the second term in ∆ is always smaller than the threshold correction due to the electron shown by the first term in Eq. (2.50), and so, the axion-photon coupling is mainly determined by the two contributions, the threshold correction due to the electron and the (residual) anomalous coupling. One can see by using the results in Appendix E that the two contributions tend to have the opposite sign, and the axion-photon coupling can be extremely small when they are nearly canceled with each other.

Masses for axion and CP-odd Higgs bosons
In this section, we focus on the axion and CP-odd Higgs bosons and evaluate their masses and dependence on the model parameters. As discussed in the previous section, the mass eigenvalues can be calculated by diagonalizing the mass matrix for the CP-odd Higgs bosons, Eq. (2.29). However, the axion mass can be roughly estimated from the (4,4) element of the This assumption is natural since it implies that the interaction V I between the EW and B − L sectors is suppressed and does not require any extra fine tuning to obtain the EW scale. Then, the scale of the decay constant f a and the soft breaking parameter m 12 determine the mass scale of the axion. For example, if m 12 is of order the EW scale and f a is about 10 10 GeV, the axion mass is of order keV. For the anomaly-free axion, the keV-scale mass and intermediate-scale decay constant are particularly attractive from a phenomenological point of view. One reason is that, especially in the absence of the anomalous coupling to photons as Eq. (2.49), the axion is stable on cosmological time scales, making it a good candidate for DM. Through this limiting expression of g aγ , the axion mainly decays into two photons with the lifetime, which is so long that it can easily satisfy the current limit from the X-ray observations. While this is the special case where the effects of breaking of U (1) F and mixing among the axion and the CP-odd Higgs bosons vanish, we will discuss in Sec. 6 how severe the X-ray bound becomes when the general expression of g aγ , Eq. (2.48), is applied. Importantly, the misalignment mechanism can naturally produce the right amount of axion to explain DM in the case where the current limit of X-ray observations can be evaded. Also, the axion can explain the excess of electron recoil events observed in the XENON1T experiment through the axion-electron coupling g ae ∼ m e /f a ∼ 10 −14 . On the other hand, if the axion accounts for about 10% of DM, it can explain not only the XENON1T excess, but also various stellar cooling anomalies simultaneously [26]. Such interesting axion scenarios can be explored by future X-ray observations such as Theseus [77,78], Athena [79], eROSITA [80], and XRISM [81] and the direct search experiments such as LZ [82] and DARWIN [83].
On the other hand, the mass of the CP-odd Higgs bosons must be somewhat heavy in order to satisfy the limits of the direct search at collider experiments. Whether such a mass spectrum can be achieved depends on the model parameters such as the soft breaking mass m 12 in the Higgs potential. The purpose of this section is to see how the masses of the CP-odd Higgs bosons and axion depend on the model parameters and to understand their behavior intuitively. To this end, here we do not impose any theoretical and experimental constraints that we will discuss in the next section.
In Fig. 1, the masses of the CP-odd Higgs bosons are shown as a function of t β 1,2 , λ 10 and the soft breaking parameters, m 12 , m 13 , m 23 . The set of the parameters, is chosen as a bench mark point, and each model parameter is individually varied in the plot, where we identify m 12 , m 13 and m 23 as inputs not the re-scaled ones. As seen in the left top panel, while m A 1 increases as t β 1 , m A 2 becomes constant in the large t β 2 region. This behavior can be understood in terms of the diagonal elements of the mass matrix B P . The (3,3) element dominantly contributes to m A 2 , and it can be approximately expressed as the constant term (B 2 P ) 33 M 2 23 when t β 1 1. On the other hand, large t β 2 makes both A 1 and A 2 heavy, since both M 2 13 and M 2 23 are enhanced in the (2,2) and (3,3) elements of B P . As to the dependence on λ 10 , the terms with λ 10 are destructive for the other terms in (B 2 P ) 33 and (B 2 P ) 22 . The dependence on the soft breaking parameters is nontrivial, especially for the mass of A 2 . While m A 1 increases with m 2 12 , m A 2 becomes almost constant for m 12 400 GeV. To understand this behavior, let us consider the mixing matrix R A parameterized by the mixing angles γ i (i = 1, 2, 3), see (2.20c). First, note that, due to the hierarchy between f a and the EW scale, the mixing angles between the axion and A 1,2 are extremely suppressed. The remaining angle γ 1 , which is responsible for the mixing of the two CP-odd Higgs states in the basis R β (z 1 , z 2 , z 3 ), can be expressed by In terms of the mixing angle γ 1 , m A 2 can be written as As m 12 increases, only (B 2 P ) 22 increases but (B 2 P ) 23 and (B 2 P ) 33 remain the same, so the mixing angle γ 1 becomes much smaller than unity. In this case, m A 2 is approximately given by the linear combination of (B 2 P ) 33 and (B 2 P ) 23 (i.e. the last two terms in Eq. (3.4)) which is independent of m 12 . Conversely, from the middle bottom panel in the lower part of Fig. 1, one can see that the m A 2 increases with m 13 . Note however that this is due to the result of setting m 13 ∼ m 23 . When either m 13 m 23 ∼ v or m 23 m 13 ∼ v, m A 2 is bounded above, which is approximately determined by m 23 for the former case and m 13 for the later case. For the former case, one can check that terms including m 13 in the right-handed side of Eq. (3.4) are indeed canceled out when m 13 is enough large to be able to neglect other terms in the mass matrix B P .
In Fig. 2, the axion mass is shown in the plane of (t β 1 ,t β 2 ), (m 12 ,m 13 ) and (m 12 ,f a ) from left to right. We set m 2 13 = m 2 23 and m 2 12 = m 2 13 = m 2 23 in the middle and right panels, respectively. In each panel we vary only the parameters corresponding to the horizontal and vertical axis, while the other parameters are fixed as Eq. (3.2).
In the left panel, one see that negative correlation between m a and t β 1,2 . This behavior can be understood from the fact that m 12 φ † 1 φ 2 is only the soft breaking term introduced in this model and thus gives the axion mass. When we take t β 2 1, the effect of U (1) F breaking is suppressed by c 2 β 2 . On the other hand, it shrinks by c β 1 in case of t β 1 1. Thereby, compared with the dependence on t β 2 , m a slightly decreases with t β 1 .
From the middle panel, one can see the characteristic soft mass dependence of m a . Namely, the axion mass increases only when all of the soft breaking masses increase simultaneously. This is because for the axion to have mass, in addition to the explicit breaking of the U(1) flavor symmetry m 12 , either m 13 or m 23 must be nonzero to transmit the breaking in the EW sector to the B − L sector where the axion lives. This behaviors can be understood by supposing the case of t β 1 = 1 and m 2 13 = m 2 23 . In this specific parameter choice, the CP-odd Higgs field in the Higgs basis z 3 corresponds to the mass eigenstate and mixing happens only in the remaining fields z 2 and a . The approximate formula for the mass of axion can be obtained as where we have assumed that X 2 f 2 a (X = m 12 , m 13 or v) and neglected terms of the order O(X 6 /f 6 a ). The soft breaking parameters m 13 and m 23 have little effect on the axion mass in the case of m 2 is similar, and m 2 12 has little effect on the axion mass. One can see that the axion mass increases when all the soft breaking masses increase and become larger than v. In addition, the behavior of m a in the right panel obeys the last line of (3.5).
To summarize this section, we have shown the following properties for the mass scales of A 1,2 and a; • Typical scale of the mass of A 1 and A 2 is mostly determined by the soft breaking masses.
Thus, the CP-odd Higgs bosons become heavy and get decoupled by taking two of the three breaking parameters are sufficiently larger than EW VEV.
• The axion mass increases when all of the soft breaking masses, m 2 12 , m 2 13 , m 2 23 , increase simultaneously. On the other hand, the axion mass is suppressed when m 2 12 or both m 2 13 and m 2 23 are small.
• If the soft breaking masses are of order the EW scale, the axion mass is of the order of keV for f a ∼ O(10 10 -10 11 ) GeV.

Theoretical and experimental bounds on 3HDM
Here we discuss theoretical and experimental constraints on the model parameters. For the former, we take into account conditions for the potential bounded from below (BFB) and perturbative unitarity, and perturbativity on the running coupling constants. For the latter, we consider constraints from the EW oblique parameters, B meson decays, and B meson mixing. We will impose these constraints on the model parameters in the numerical calculations in Sec. 5 and 6.

The potential bounded from below and perturbative unitarity
In order to obtain the stable minimum after the EWSB, the Higgs potential should be bounded from below in any direction of the Higgs fields. While the B − L Higgs fields are involved in the original potential (2.1), the radial modes of the B − L Higgs fields are integrated out in our analysis. Hence, we focus on the conditions for the 3HDM potential (2.2) to be bounded from below. In the pioneering work of Ref. [84], the BFB conditions were derived for a potential that involves two Higgs doublets and one Higgs singlet. More recently, the BFB conditions for the 3HDM with Z 3 symmetry were derived in Ref. [43].
As mentioned in Sec. 2, the potential (2.2) can be obtained by setting the U (1) F symmetry breaking terms except for m 2 12 to zero in their Z 3 invariant potential. Therefore, we can simply read off the BFB conditions from the results of Ref. [43] as and The partial wave unitarity bound for the elastic 2 → 2 scattering processes in the high energy limit restricts scalar couplings in the Higgs potential. In Ref. [85], the unitarity bounds were applied to the SM to derive an upper limit on the mass of the Higgs boson. A similar argument can be made for the extended Higgs models. In Ref. [86], the tree-level unitarity bounds were derived in the framework of 3HDM with Z 3 symmetry. Following the results of Ref. [86], we can obtain eigenvalues Λ i (i = 1 − 21) for the partial waves of the S-matrix amplitudes for 2 → 2 scattering processes with the replacement r 1−3 = λ 1−3 , r 4−9 = 1 2 λ 4−9 , c 11 = 1 2 λ 10 , (4.5) where r i (i = 1 − 9) and c 11 denote the quartic couplings in 3HDM with the Z 3 symmetry defined in Ref. [86]. The criterion that the partial wave amplitudes satisfy the unitarity is given by The Landau pole appears at Λ < 10 10 GeV when we set t β 2 = 0.9 or smaller, and there is no Landau pole at Λ < 10 12 GeV when t β 2 = 1 or larger. Since we consider the decay constant f a of order 10 10 -10 12 GeV in the following numerical analysis, we take the lower bound on t β 2 as t β 2 1. (4.8)

Electroweak S and T parameters
Here we discuss the limits of the electroweak precision measurements for heavy particles. First, the electroweak ρ parameter does not deviate from unity at the tree level in the multiple Higgs doublet models. However, all the additional Higgs bosons contribute to it at 1-loop level. The loop corrections to the ρ parameter are described by electroweak oblique parameters.
The electroweak oblique parameters, which parameterize new physics effects for the gauge boson self-energies, were first proposed in Ref. [87]. Their analytical expressions in the multi-Higgs doublet models were calculated in Refs. [88,89]. With the definitions of S and T parameters given in Ref. [90], their analytical expressions in the 3HDM with axion can be derived from the new scalar boson loop contributions to the gauge boson two-point functions. We give the general formulae for the S and T parameters in Appendix C, which are used in the numerical calculations in Sec. 5 and 6. Since the general formulae for the S and T parameters are somewhat lengthy, we here describe the analytical expressions in the alignment limit. New physics contributions to the S and T parameters, ∆S = S 3HDM −S SM and ∆T = T 3HDM − T SM , are written in terms of the Passarino-Veltman functions [75] as where These are obtained from the general expressions given in Eqs. (C. 16) and (C.17), taking (α 1 , α 2 )=(β 1 , β 2 ), and assuming that the mixing matrix for the CP-odd Higgs bosons R A is parameterized by (2.20c) with γ 2 = γ 3 =0. The B 5 function is evaluated by using LoopTools [91]. We use the experimental values for the S and T parameters given in Ref. [92]. In Fig. 4, to illustrate the parameter space favored by the constraints from the S and T parameters, we show the χ 2 values for the S and T parameters with the correlation coefficient +0.92 [92] as a function of the mass of H 2 with α 3 = π/4 (red curve), π/2.1, (blue curve) π/2 (green curve). The other parameters are taken as The large mass difference between m H 2 and the masses of other heavy Higgs bosons tends to be in conflict with the experimental results for each value of α 3 . Setting α 3 = π/4 (π/2.1) yields the bounds for the mass of H 2 , 550 (400) m H 2 900 GeV (1700 GeV). Although = m H 3 are satisfied in Eq. (4.12), the constraint becomes much tighter in case that the mass degeneracy among these additional Higgs states is assumed. Note that the contributions from the axion to the S and T parameters are negligible due to the tiny mixing angles; i.e., γ 2 , γ 3 0.

Flavor constraints
The 3HDM model parameter space is limited by measurements of B meson rare decays and B meson mixing. A particularly strong constraint is given by B → X s γ, which is altered from the SM prediction by the additional contributions of the loop diagram of the charged Higgs boson. The Heavy Flavor Averaging Group (HFAG) gives the the experimental value for BR(B → X s γ) as [93] BR(B → X s γ) exp = (3.32 ± 0.15) × 10 −4 (4.14) with the cut off for the photon energy E γ > 1. 6 GeV. The precise evaluation of the SM prediction with QCD corrections have been performed at NLO [94][95][96][97][98][99][100][101][102] and at NNLO [103,104]. Effects of the charged Higgs boson loop contributions for B → X s γ with NLO QCD [105][106][107][108][109] and NNLO QCD [104,110,111] have been investigated in the 2HDMs (Also, see the evaluation in Type II 3HDM, in Ref. [44]). We evaluate BR(B → X s γ) for our model by the linearized formula given in Ref. [104]: where the first term contains theoretical uncertainties. The second and third terms, ∆C 7,8 , denote additional new physics contributions to the Wilson coefficients C 7 and C 8 6 at the scale of the EW theory µ 0 , which is taken to be µ 0 = 160 GeV. Making use of the explicit formula of ∆C 7,8 in Ref. [110], we include the contributions of H ± 1 and H ± 2 to Eq. (4.15). We also take into account the constraint from the mass difference ∆M s for the mass eigenstates in the B s -B s system. The experimental value is taken from the result of HFAG [93], (∆M s ) exp = 17.757 ± 0.021 ps −1 . (4.16) The analytical formula is given by [112,113]  f Bs = (228.8 ± 0.7 ± 1.9) MeV , B Bs = 1.327 ± 0.016 ± 0.030 , where the second (third) values stands for statistical and systematic uncertainties (systematic theoretical uncertainties). Contributions from the charged Higgs bosons to B → X s γ and ∆M s are controlled by their masses m H ± 1,2 and the quark Yukawa couplings ξ H ± 1,2 . Thus, these measurements We take t β2 =1.2 (red), 1.5 (blue), 2 (green), and the mixing angles for the charged Higgs bosons are fixed as γ + = π/4 (left) and γ + = π/3 (right).
give the lower bounds on the masses of charged Higgs bosons for fixed t β 2 and γ + as shown in Fig 5, where we choose t β 2 = 1.2 (red line), 1.5 (blue line), 2 (green line), and γ + = π/4 (left panel), π/3 (right panel). One can clearly see that the constraint from ∆M s is more severer than BR(B → X s γ). At (t β 2 , γ + ) = (1.5, π/4) and m H ± is inversely proportional to the t β 2 , the lower bounds for the masses are relaxed by taking larger t β 2 . On the other hand, if γ + is close to γ + = π/2 (0), loop contributions of H ± 2 (H ± 1 ) almost decouple in both of B → X s γ and B s -B s mixing. We emphasize that the Yukawa structure of the quark sector for our model is the same as Type-I and Type-X 3HDMs, since, at the LO, only the up-and down-type quark Yukawa couplings are relevant for B → X s γ. Hence, the flavor constraints for the Type-I and Type-X 3HDMs are similar to our results.

Phenomenological implications of the XENON1T excess
The axion mass is determined by the ratio of the soft breaking masses m ij . This implies that for a given decay constant, the axion mass can be related to the mass spectrum of the heavy Higgs bosons. To see if this picture is correct, we study the implications of the XENON1T excess in the electron recoil event [25] for the 3HDM with the B −L Higgs bosons. We will examine the parameter region which explains the XENON1T excess and satisfies all the constraints presented in the previous section.
The favored range of the axion mass that can explain the XENON1T excess at 95 % confidence level is given by [115]

The viable parameter space
We have performed numerical calculations to find a parameter region that would explain the XENON1T excess while satisfying the experimental limits described in the previous section. To this end we assume degenerate masses for the additional Higgs bosons to satisfy the constraints from the S and T parameters, and we take the alignment limit (2.42), α 1 = β 1 and α 2 = β 2 . We also take t β 1 =1 in order to avoid the constraint from the current X-ray observations. The detailed discussions are presented in Sec. 6. In the numerical analysis we vary the remaining input parameters in the following ranges, where the lower limits of t β 2 are chosen so as to satisfy the constraints from the RG running of the top Yukawa coupling, and the measurements of the B meson decay B → X s γ and B s -B s mixing. With the scan range Eq. (5.4), we obtain 2.3 × 10 −14 g ae 4 × 10 −14 , which is within the range indicated in Ref. [115]. Furthermore, the range of axion-photon coupling is 4.5 × 10 −20 GeV −1 g aγ 1.5 × 10 −18 GeV −1 . We have checked that there are no parameter points excluded by the current limit of X-ray observations in this case.
In Fig. 6, we show the axion mass as a function of the relevant input parameters. The region between the black dotted lines indicates the mass suggested by the XENON1T excess. The negative correlation between m a and t β 2 seen in Fig. 2 is also confirmed in Fig. 6. One of the intriguing observations here is that the XENON1T excess restricts the range of t β 2 , i.e., t β 2 3. The range of the soft breaking mass is also limited as 180 GeV m 12 300 GeV. Note that the upper bounds of these parameters are correlated with the masses of the additional Higgs bosons as will be seen shortly. ij . Note also that the lower bound on the mass of additional Higgs bosons, m Φ 510 GeV, is given when the axion mass is in the range of Eq. (5.1). This lower bound comes from the combination of the XENON1T excess and the constraint from the B s -B s mixing. If the mass of additional Higgs boson is less than 1 TeV, t β 2 should not be large to keep m a ∼ O(1) keV (see the top middle panel in Fig. 6). In contrast, to evade the constraint from the B s -B s mixing, t β 2 2 is required for m Φ 450 GeV. As a result, if the masses of the additional Higgs bosons were lower, the XENON1T excess would conflict with the measurement of the B meson mixing. From the right panel of Fig. 7, one can see that the value of m Φ is not so different from M 12 . This is mostly caused by the constraints from the RG evolution of the scalar quartic couplings. A large hierarchy among the mass of the additional Higgs bosons and the rescaled soft breaking parameters enlarges the scalar quartic couplings. If one set m Φ = M 12 = M 13 = M 23 in the alignment limit, the analytical expressions of the quartic couplings can be reduced as Namely, in this case, these couplings do not depend on the masses of the additional Higgs bosons and the soft breaking masses, and so, the constraints from the RG evolution can be evaded. We also note that a further requirement, λ 10 = 0 leads to λ 7 = λ 8 = λ 9 = 0 and λ 1 = λ 2 = λ 3 = 2λ 4 = 2λ 5 = 2λ 6 . In this limit, the 3HDM potential V 3HDM except for the soft breaking terms V soft has an Sp(6) symmetry [52], which is explicitly broken by the terms with the coefficients M 2 12 or κ mφn v 2 Sm (m = 1,2, n = 1, 2, 3). In fact, the allowed region for λ 10 is −0.1 λ 10 0.17, which contains this limit. Hence, the Sp(6) symmetry would be the desired (approximate) symmetry for the 3HDM part of the Higgs potential to satisfy all the theoretical and experimental constraints.
Another consequence for heavy Higgs bosons from the XENON1T excess is that there is a correlation between the mass m Φ and the ratio of the VEVs t β 2 as shown in Fig. 8. For m Φ = 800 GeV, the allowed range of t β 2 is 1.3 t β 2 2.4. The range is enlarged for heavier additional Higgs bosons. The Yukawa couplings for the additional Higgs bosons are controlled by t β 1 and t β 2 for fixed mixing angles. Thus, if extra Higgs bosons are found in collider experiments, the decay properties of the extra Higgs bosons may allow us to test whether they are consistent with the anomaly-free axion that explains the XENON1T excess.
We also comment that the parameter regions favored by the XENON1T excess may be explored by the future measurements of the B meson mixing. In Ref. [116], projected sensitivity for new physics effect to the B meson mixing is studied, considering planned LHCb Upgrade II [117] and a possible upgrade of Belle II [118] as well as FCC-ee as a tera-Z factory. We find that the future 95% sensitivity for B s mixing by LHCb 300 fb −1 and Belle II 250 ab −1 can probe m Φ 940 (700) GeV for (t β 2 , γ + ) = (1.75, π/2) ((2, π/2)).

Predictions for the SM-like Higgs boson decays
Here let us illustrate the extent to which the decays of the SM-like Higgs boson H 1 are deviated from the SM prediction in the allowed parameter space satisfying all the theoretical and experimental constraints. As we consider the alignment limit, the decay of the SM-like Higgs boson into weak gauge bosons and the fermions are the same as the SM predictions at the tree level. On the other hand, the decays of H 1 → γγ can deviate from the SM one through the charged Higgs boson loop diagrams. Let us express the deviation from the SM prediction for H 1 → γγ in terms of the modifier defined by where Γ(H 1 → γγ) SM denotes the decay rate in the SM. The precision measurements of the Higgs boson coupling will be performed in the future collier experiments. At the HL-LHC (ILC 250 GeV), the sensitivity of the coupling modifier for H 1 → γγ can reach 1.6% (1.4%) [119]. Furthermore, in Ref. [119] the combined sensitivity of FCC-ee, FCC-eh, FCChh is estimated as 0.31%. For the evaluation of the decay rate in the 3HDM, we use the following analytical formula, For the explicit forms of the loop functions F 1 , F 1/2 and F 0 , we refer the reader to Ref. [120,121]. In the alignment limit, the scalar couplings for the charged Higgs boson are given by We show the numerical results for ∆κ γ in Fig. 9, where the different color of the points corresponds to the values of m a and t β 2 in the left and right panels, respectively. Since we take the alignment limit, the deviations from the SM predictions purely come from the contributions from the charged Higgs boson loop diagrams. Remarkably, contributions of the charged Higgs bosons do not decouple in case of m H ± 1 600 GeV. This is due to the constraint from ∆M s . As seen in the right panel, relatively high t β 2 is required for lighter H ± w1 by the constraint. We have checked that ∆κ γ 0 can be realized even in the region m H ± 1 600 GeV if the constraint from ∆M s is switched off. The deviation ∆κ γ is maximized at m 2 H ± 1 = 200 GeV and exceeds −2%. In this case, the corresponding axion mass is less than 2.1 keV. If the axion mass lies in the range favored by the XENON1T excess, the deviation shrinks, i.e., −0.22% ∆κ γ 0%, which would be difficult to detect at the HL-LHC and ILC 250 GeV. Equivalently, if the deviation of H 1 γγ is found within the range ∆κ γ −0.25%, the axion should be lighter than 2.1 keV.
We have also calculated the deviation in the self-coupling of the SM-like Higgs boson λ H 1 H 1 H 1 using the effective potential method. In contrast to the case of H 1 → γγ, all additional Higgs bosons can contribute at the 1-loop level. We find that the magnitude of the deviation in λ H 1 H 1 H 1 is not comparable with the projected sensitivity in the future collider experiments with the order of 10% [119,[122][123][124]. The main reason for this result is the strong constraint from the RG evolution of the scalar quartic couplings. The nondecoupling effects of the additional Higgs bosons are highly suppressed.

Exploration of extra Higgs bosons from X-ray observations
We have discussed in the previous section implications for the extra Higgs bosons and predictions in the SM-like Higgs boson decays in the scenario where the axion has properties suggested by the XENON1T excess. An interesting aspect of such axion from a cosmological point of view is that it can naturally explain the observed DM abundance, and it can also be probed by the future X-ray observatories such as Theseus [77,78], Athena [79], eROSITA [80], and XRISM [81]. While in the previous section the axion decay constant f a is fixed to reproduce the value of g ae suggested by the XENON1T excess, we here vary f a and investigate the connection between the axion coupling with electron g ae and the mass of the additional Higgs bosons. We also demonstrate the parameter regions allowed by the current X-ray observations or probed by future X-ray observations. We then present the expected mass spectrum of the extra Higgs bosons if the axion is detected in the future X-ray observatories. In addition, one may be interested in the effect on the alignment parameters for the CP-even Higgs bosons α 1 and α 2 when considering the testability of the axion. To see how much these parameters can deviate from the alignment limit depending on the axion mass and the axion-electron coupling, we numerically evaluate the possible size of deviations in the SM-like Higgs boson with the weak gauge bosons κ V .

Scale of the masses of extra Higgs bosons
Bounds from X-ray observations We first show the allowed parameter regions to satisfy the constraint from X-ray observations by scanning the model parameters. We take all the dimentionful parameters degenerate, i.e., m Φ = M 12 = M 13 = M 23 , (6.1) Figure 10. The parameter region excluded by the X-ray observations is indicated by the blue shaded region in the plane of (g aγ , m a ). Future sensitivities are also shown by the colored dashed lines. The size of t β1 is shown by the color difference of the points. taking into account the constraints presented in Sec. 4. In Fig. 10, we show the current X-ray constraints on the axion-photon coupling by the blue shaded regions [76,[125][126][127][128][129][130][131][132][133][134][135][136][137][138]. In particular, the bound from XMM-Newton [139] gives the strongest constraint in 5 keV < m a < 16 keV. Also shown are the future sensitivities of the Theseus XGIS-X (blue dashed line), Theseus SXI ( red dashed line), Athena (pink dashed line), and eROSITA (brown dashed line), which are taken from the projection limits in Refs [140,141]. As can be seen, most of the parameter points with t β 1 > 1.1 are excluded by the current X-ray bounds. Furthermore, the region with t β 1 > 1.001 can be probed by the future observations. This result can be understood from the analytical expression for g aγ of Eq. (2.50). While the first term is suppressed by m 2 a /m 2 e , the second term is controlled by the quantity ∆, which involves the effect of the mixing among the axion and the CP-odd Higgs bosons. In the case of t β 1 = 1, ∆ is highly suppressed as mentioned in Sec. 2, and thus the constraint from the X-ray can be evaded. On the other hand, if ∆ is O(1), the corresponding g aγ becomes too large to satisfy the current X-ray constraints. One can see that g aγ can be extremely small when t β 1 is close to unity. This is because, as mentioned before, g aγ mainly receives two contributions with an opposite sign, and a cancellation could take place. This should be contrasted to the axion-electron coupling, whose magnitude is always determined by the ratio of the electron mass and decay constant. Note that the apparent tight constraint on t β 1 can be ameliorated by considering larger values of t β 2 (see Appendix E). In the following numerical calculations, we impose the current X-ray bounds in addition to the theoretical and experimental constraints discussed in Sec. 4.

Cosmological abundance of axion
We next discuss the production mechanism of the axion and its abundance. The axion can be produced through the misalignment mechanism and/or thermal production in the early Universe. Thermally produced axion with O(1) keV is regarded as warm DM. Hence, it suffers from bounds on galactic-scale structure formation, i.e., Lyman-α forest observations [142][143][144], if it saturates all components of DM. To safely evade the constraint, we assume that the axion is primarily produced by the misalignment mechanism [5][6][7].
When the Compton length of axion is larger than the Hubble scale, by the Hubble friction the axion is fixed at a certain field value θ i (= a i /f a ), which is called the initial misalignment angle. As the Universe cools down, the Hubble scale becomes comparable with the axion mass at a certain point. Then, at the axion starts to oscillate around the potential minimum. Shortly thereafter the axion abundance gets fixed. The temperature at the onset of the oscillation, T osc can be estimated from Eq. Then the density of the axion can be estimated as with the critical energy density ρ 0 = 3m 2 pl H 2 0 /(8π). The number density at T osc reads where F (θ i ) is the anharmonicity factor F (θ i ) = ln(e/(1−θ 2 i /π 2 )) 7/6 [145,146], with e being Napier's constant, by which the effect of quartic couplings in the axion potential is included. The factor F (θ i ) affects the number density when the θ i is not small. The cosmological abundance of the axion substantially depends on the initial condition for the axion fields, i.e., θ i , the mass m a and the decay constant f a . Thus, the correct DM abundance Ω DM h 2 0.12 can be explained by taking an appropriate value of θ i for fixed m a and f a . In particular, the initial angle is of order unity for m a = O(1) keV and f a = O(10 10 ) GeV. The numerical results for θ i to explain the observed DM abundance are shown in Fig. 11 in the plane of the axion mass and the axion-electron coupling g ae . Here the decay constant f a is converted to g ae through Eq. (E.24), and t β 1 is chosen as t β 1 = 1 for the thin solid black line. Other parameters are set as g * = g s (T osc ) = 61.75, the Hubble parameter H 0 = (9.778Gyr) −1 h and the plank mass m pl = 1.22 × 10 19 GeV. We use g s (T 0 ) = 3.91, T 0 2.35 × 10 −4 eV for the present photon temperature and effective number of entropic degree of freedom. One can see that the initial misalignment angle θ i satisfying the observed density of DM is O(1) when g ae = O(10 −14 ). For smaller g ae , θ i must be smaller than unity, which requires a mild fine-tuning. The abundance scales as Ω a ∝ m 1/2 a g 2 ae , so that it is not sensitive to m a compared with g ae .

Correlation between axion-electron coupling and extra Higgs boson masses
For the evaluation of the mass of the extra Higgs bosons, we scan the parameters in the range given in Eq.(6.2), setting all the mass parameters equal, i.e., Eq. (6.1), and taking the alignment limit. In Fig. 11, the allowed parameter points satisfying all the constraints are also shown by the green points in the plane (m a , g ae ). The intensity of the color corresponds to the range of m Φ , i.e., m Φ < 1 TeV (lightest green), m Φ < 3 TeV (lighter green), m Φ < 5 TeV (moss green), and m Φ < 10 TeV (dark green). We also show the future sensitivity of the direct DM search experiments, DARWIN [83] and LZ [82], by the dashed red and blue lines, respectively. As can be seen from the figure, in the parameter region of 10keV m a 20keV and 2 × 10 −14 g ae 1 × 10 −13 , there are few allowed parameter points. This is because many parameter points in this region are excluded by the current limits of the X-ray observations. The maximum value of m a is determined for a fixed g ae and m Φ , e.g. we obtain m a 15 keV for g ae = 5 × 10 −15 and m Φ = 10 TeV. Since m a is inversely proportional to the decay constant, the maximum value of m a decreases as g ae becomes smaller.
Interestingly, there is a correlation between the mass of the extra Higgs boson m Φ and the axion-electron coupling g ae through the axion mass. Thus, one can obtain the information of the mass spectrum for the heavy Higgs bosons from the axion searches if it is detected or some anomalies are indicated in the direct searches and/or future X-ray observations. For example, unidentified X-ray line at around 3.5 keV was reported from observations of the galaxy clusters [147,148] and galaxies [148,149] (also see the recent review in Ref. [125]), which may be originated from the decay of axion into photons [150,151]. As one can see from Fig. 11, if such a hint of axion is confirmed around e.g., (m a , g ae ) ∼ (7 keV, 7 × 10 −14 ) or (2.5 keV, 1 × 10 −14 ), the favored parameter region should be m Φ 3 or 1 TeV. In other words, the lower bound on the mass of additional Higgs bosons can be derived from the direct searches of the axion and/or the X-ray observatories in future.
In order to reveal which parameter points are probed by the future X-ray observations, we show in Fig.12 the parameter points that satisfy the current X-ray bounds and other theoretical and experimental limits and that are within the sensitivity reach of future X-ray observations (this corresponds to the parameter points located below the blue solid lines and above dashed lines in Fig. 10.). Intensity of color denotes the magnitudes of |∆| and t β 1 in the left and right panels, respectively. As can be seen from the left panel, if m a 5keV, Comparing the left and right panels, one can see that the dependence on t β 1 is similar to that of |∆|. This is because the quantity ∆ is basically controlled by t β 1 and M 12 and t β 1 ∼ 1 means that the effect of mixing among the axion and the CP-odd Higgs bosons are small. Comparing Figs. 11 and 12, one can see that the distribution of points is almost the same, with only fewer points satisfying condition |∆| > 10 −6 in Fig. 11. Thus, one can conclude that most of the parameter points are surveyed by the future X-ray observations. This can be understood by noting that small g aγ is realized only when the two contributions are nearly canceled with each other.
Before closing this section, we show in Fig. 13 another example of the correlation between the mass of the additional Higgs bosons m Φ and the axion mass. All the parameter points correspond to the ones within the reach of the future X-ray observations. Here, the intensity of the color represents different values of the decay constant f a . Thus, if we fix f a and m a , the minimal value of m Φ is determined. The minimum value of m Φ increases for heavier axion mass since it requires larger M 12 for a fixed f a .

Mixing angles
We also investigate the impact on the alignment parameters for the CP-even Higgs bosons, α 1 and α 2 in case that the axion is detected (or indicated) in the future X-ray observatories.
To demonstrate this, we numerically evaluate the scaling factor for the weak gauge boson The scan range for the other parameters is taken from Eq. (6.2), assuming that all dimensionful parameters are degenerate. The numerical results for the correlation between the axion-electron coupling and the scaling factor κ V are shown in Fig. 14, where the darkness of color corresponds to the range of axion mass; (dark blue): 10 keV < m a , (blue): 5 < m a < 10 keV, and (light blue): 1 < m a < 5 keV. We note that many parameter points with m a > 10keV are excluded by the constraint from the current X-ray observations. One sees that requiring axion mass to be larger than 10 keV makes almost alignment limit, i.e., |κ V − 1| 4 × 10 −5 . This consequence comes from the fact that the maximum size of the soft breaking parameters M ( )2 ij is limited by the constraint from the RG evolution of the scalar couplings as α 1 and α 2 deviate from the alignment limit. Also, smaller g ae makes the constraint tighter, so that the possible deviation of |κ V − 1| becomes small. We also numerically checked that the deviation over 1% is difficult even if the axion mass is within the range of m a 1 keV. This is mainly due to the lower bound of m Φ > 200 GeV in the parameter scan range Eq. (6.2). Under the assumption the perturbativity for the running scalar couplings is satisfied, larger |k V − 1| requires lighter extra Higgs bosons.
The axion with m a 5 keV can be probed by the future-X ray observatories such as eROSITA, Athena, and Theseus XGIS-X and direct searches such as LZ, and DARWIN. If there are some indications in the future observatories, one can set the upper bound for κ V − 1; e.g., κ V − 1 1 × 10 −2 % (2 × 10 −3 %) for g ae 1 × 10 −13 (2 × 10 −14 ). The predicted deviations is not so large compare with the usual extended Higgs models (see, e.g., Refs. [152,153]). Hence, if one finds the deviations in hV V coupling over 1% in the future collider experiments, such as HL-LHC, ILC, CPEPC, and FCC, 3HDM with B-L Higgs bosons can be ruled out.
In this way, axion searches by the future X-ray observatories and the direct detection potentially have the impact on probing the extra Higgs bosons and the mixing parameters for the CP-even Higgs bosons, and eventually can narrow down the structure of the Higgs sector.

Conclusions
We have investigated phenomenological implications of the axion DM based on the model with anomaly-free global flavor symmetry, which was originally proposed in Ref. [24]. To build a concrete renormalizable model that includes the anomaly-free axion, we have considered the three Higgs doublet model with three B − L Higgs fields, in which a global U (1) F flavor symmetry is imposed. In particular, we have focused on the axion with the mass of order O(1) keV. Such an axion DM scenario is promising because it can explain the reported excess in the electron recoil events of the XENON1T experiment, and because it can also be probed by the future X-ray observatories such as eROSITA, Athena, Theseus, and XRISM.
We have revealed that in this concrete model the axion-photon coupling involves the breaking of U (1) F and the mixing between the axion and the CP-odd Higgs bosons as seen in Eq.(2.48). Thus, even for the anomaly-free axion, the anomalous coupling (the first term of (2.48)) is not completely canceled out unless the effects of U (1) F breaking and the mixing are absent. This leads to the anomaly-free axion being more severely constrained by the X-ray observations depending on the model parameters as seen in Fig. 10. To put it another way, future X-ray observation experiments will be able to detect anomaly-free axion more easily.
In order to investigate the possibility of narrowing the range of model parameters based on the nature of the axion (mass and coupling), we have surveyed the mass spectrum of the axion and the extra Higgs bosons in the allowed parameter space satisfying the theoretical constraints and the experimental constraints given in Sec. 4. We have found that there are correlations among the ratio of the VEVs t β 2 and the mass of the extra Higgs bosons m Φ , provided that the axion has properties indicated by the XENON1T excess. As a result, the upper bound on t β 2 is given depending on the scale of m Φ .
We have also discussed the impact of the axion searches in the future X-ray observatories and direct detection on the extended Higgs sector. In particular, we have clarified the correlation among the axion coupling with electron g ae , the axion mass, and the extra Higgs boson mass. We have revealed that the lower bound on the mass of the extra Higgs bosons can be obtained if the axion with a mass of order keV is indeed detected (or indicated) in the future X-ray observatories and direct detection. In addition, we have demonstrated that the axion searches by the X-ray observatories and direct detection can restrict the deviations of the SM-like Higgs boson couplings with the weak gauge bosons from the SM predictions. Thus, the 3HDM with the B − L Higgs fields can be probed by the synergy of the axion and extra Higgs boson searches.

C Electroweak oblique parameters
We derive here the formulae for the electroweak S and T parameters in the 3HDM with B-L Higgs bosons. To this end, we first evaluate scalar boson loop contributions to the transverse part of the gauge boson two-point functions. For the weak gauge bosons, they are written in terms of the Passarino-Veltman functions [75] as where R z P is 4 × 3 matrix for the CP -odd scalars A i ,i.e.,  It can be extracted from the original 4 × 4 orthogonal matrix R P , removing the 4th column vector, and satisfies (R z p ) T R z p = I 3×3 while R z p (R z p ) T = I 4×4 . The parameter m A j denotes the mass for the CP-odd scalar bosons, i.e.,( is the mass for the charged scalar bosons, ). We note that the last two terms in Eqs.(C.1) and (C.2) vanish in the limit of γ 2,3 → 0. Thus, hereafter, we just drop them.
To extract new physics contributions to the oblique parameters, we need to subtract the SM contributions from Π S W W,ZZ . It is written by For the photon-photon two two-point functions and the photon and Z boson mixing two-point function, the new physics contributions just stem from the charged Higgs boon loops as ) . (C.8) We use the following definition of the electroweak oblique parameters [90], where the functions Π AB T is defined by T (q 2 ) , (C.12) Π γγ = e 2 Π QQ T (q 2 ) , (C.13) Π Zγ = eg Z Π 3Q T (q 2 ) − s 2 W Π QQ T (q 2 ) , (C.14) Π ZZ = g 2 Z Π 33 T (q 2 ) − 2s 2 W Π 3Q T (q 2 ) + s 4 W Π QQ T (q 2 ) . (C. 15) In the expressions of Eqs. (C.9)-(C.11), real part for the Π AB T (q 2 ) are taken and the relevant pinch terms are added. These yield the concrete expressions for the oblique parameters to the new physics contributions in the 3HDM, The Wilson coefficient is separated into three parts, where C W W V LL corresponds to the contribution from the box diagram with two virtual W bosons, which is the same with SM contributions [158,159]. The second (third) term comes from the box diagrams with virtual W boson and charged Higgs bosons (two virtual charged Higgs bosons). The analytical expressions for each coefficient are given by where x H ± i and x t denote the mass fraction ) and x t =m 2 t (µ 0 )/m 2 W withm t (µ 0 ) being the running mass at the scale µ 0 . The terms with x b = m 2 b /m 2 W are neglected in these expressions. We note that the limit reproduces the result of 2HDMs [112,113,160,161].

E Lepton couplings of the anomaly-free axion
In this Appendix, we derive the lepton couplings for the axion, following Ref. [162]. The scalar fields are defined by where the indices run i = 1-3, m = 1,2 . The U (1) F transformation for these scalar fields and fermion fields forms with ξ = [0, 2π]. The Nether current for the U (1) F symmetry is given by From these equations, we can write the U (1) F current Here, we have defined the axion field Sm . The physical axion field should be defined in such a way that it does not mix with the NG boson G 0 , which is absorbed by the longitudinal component of Z 0 [163,164]. The NG boson G 0 is defined by with v 2 = i (2Y i v i ) 2 = (246GeV) 2 and Y i being the hypercharge for Φ i . Hence, the physical axion field a is defined by where we have used that hypercharge for the singlet fields is zero. This modifies the U (1) F Eq. (E.7) current as [162][163][164] Using the obtained U (1) F current, the axion-lepton interaction is written by where This can be rewritten by the following form where we define f a = v a . The coefficients C V al and C A al are defined by with the unitarity matrices U l L/R for the fermion fields l L/R .
We calculate the coefficients C V al and C A al in the case of Type-B: In the end, we get We have defined the unitary matrices V l L , V l R as The diagonal component of C V al vanish by the equation of motion for fermion fields. For the model presented in sec. 2, the unitary matrices U ψ L/R are identity matrix since the off-diagonal components of the lepton Yukawa matrix are zero. One then obtains the axion-lepton couplings, (E. 26) In the limit of m 12 → 0, these agree with the numerical results of (2.35) which is obtained by diagonalizing the mass matrix M 2 P . When m 12 = 0, the actual axion-electron coupling becomes slightly smaller than the above estimate. The approximate expression for the axion-lepton coupling (2.36) is obtained in the limit of c β 2 → 0 or t β 1 → 1 since the flavor charge of the electron is q e = 1 − (−2) = 3.