Resolving the Hubble tension in a U(1)$_{L_\mu-L_\tau}$ model with Majoron

In this paper, we explore possibilities of resolving the Hubble tension and $(g-2)_{\mu}$ anomaly simultaneously in a U(1)$_{L_\mu - L_\tau}$ model with Majoron. We only focus on a case where the Majoron $\phi$ does not exist at the beginning of the universe and is created by neutrino inverse decay $\nu\nu\to \phi$ after electron-positron annihilation. In this case, contributions of the new gauge boson $Z'$ and Majoron $\phi$ to the effective number of neutrino species $N_{\rm eff}$ can be calculated in separate periods. These contribution are labelled $N'_{\rm eff}$ for the U(1)$_{L_\mu - L_\tau}$ gauge boson and $\Delta N_{\rm eff}^\prime$ for the Majoron. The effective number $N_{\rm eff} = N'_{\rm eff} + \Delta N_{\rm eff}^\prime$ is evaluated by the evolution equations of the temperatures and the chemical potentials of light particles in each period. As a result, we found that the heavier $Z'$ mass $m_{Z^\prime}$ results in the smaller $N_{\mathrm{eff}}^\prime$ and requires the larger $\Delta N_{\mathrm{eff}}^\prime$ to resolve the Hubble tension. Therefore, compared to previous studies, the parameter region where the Hubble tension can be resolved is slightly shifted toward the larger value of $m_{Z^\prime}$.


Introduction
calculate the contribution of Z and Majoron to N eff and impose a constraint on Z and Majoron parameter space. Finally, we summarize our results in section 5.
We consider a U(1) Lµ−Lτ model which contains the global U(1) L symmetry, similarly to Ref. [27]. Such a model can have a keV Majoron as a pseudo Nambu-Goldstone boson (pNG boson) originated from the spontaneous symmetry breaking of the U(1) L . In addition, this model has a U(1) Lµ−Lτ gauge boson, which can explain the muon anomalous magnetic moment and the IceCube gap of cosmic neutrino flux if this gauge boson has O(10-100) MeV mass [28][29][30][31]. As discussed in Refs. [21,26], these particles can contribute to the expansion history of the early universe and have a possibility to resolve the Hubble tension.
In this section, we show the interactions between the electron, neutrino, U(1) Lµ−Lτ gauge boson Z , and Majoron φ, which contribute to the Hubble parameter in the early universe.
At tree level, the U(1) Lµ−Lτ gauge boson interacts only with mu and tau-type leptons.

Effective coupling with electrons
parameter space. Finally, we summarize our results in section 5.

U(1) L µ L ⌧ Model
We consider a U(1) Lµ L⌧ model which contains the global U(1) L symmetry, similarly to Ref. [25]. Such a model can have a keV Majoron as a pseudo Nambu-Goldstone boson (pNG boson) originated from the spontaneous symmetry breaking of the U(1) L . In addition, this model has a U(1) Lµ L⌧ gauge boson, which can explain the muon anomalous magnetic moment and the IceCube gap of cosmic neutrino flux if this gauge boson has O(10-100) MeV mass [26][27][28][29]. As discussed in Refs. [19,24], these particles can contribute to the expansion history of the early Universe and have a possibility to resolve the Hubble tension.
In this section, we show the interactions between the electron, neutrino, U(1) Lµ L⌧ gauge boson Z 0 , and Majoron , which contribute to the Hubble parameter in the early Universe.

The U(1) Lµ L⌧ Lagrangian
The Lagrangian related to the U(1) Lµ L⌧ gauge boson is given by where Z 0 denotes the U(1) Lµ L⌧ gauge boson with the field strength Z 0 ⇢ = @ ⇢ Z 0 @ Z 0 ⇢ , and m Z 0 and g µ ⌧ are the U(1) Lµ L⌧ gauge boson mass and gauge coupling constant, respectively. J µ ⌧ denotes the L µ L ⌧ current and is written by At tree level, the U(1) Lµ L⌧ gauge boson interacts only with mu and tau-type leptons. In this model, there can be a gauge kinetic mixing between Z 0 and the SM hypercharge gauge field B : L mix = 2 B ⇢ Z 0 ⇢ where B µ⌫ is the field strength of B. Although we assume that this kinetic mixing vanishes at some high scale for simplicity, non-zero kinetic mixing appears at one-loop level at a low energy scale. This kinetic mixing then induces an interaction between In this model, there can be a gauge kinetic mixing χ between Z and the SM hypercharge gauge field B : L mix = − χ 2 B ρσ Z ρσ where B µν is the field strength of B. Although we assume that this kinetic mixing vanishes at some high scale for simplicity, non-zero kinetic mixing appears at one-loop level at a low energy scale. This kinetic mixing then induces an interaction between Z and electrons through the mixing of Z with the SM photon γ as shown in Fig 1, and the interaction term is described as follows :

E↵ective coupling with electrons
where is calculated by where e and m are the electromagnetic charge and the mass of charged lepton . The partial decay widths of Z are given as follows: Hereafter, we assume that neutrino masses are negligible and taken to be massless. Note that the BABAR experiment excludes the U(1) Lµ−Lτ gauge boson with m Z > 2m µ as a solution of the muon anomalous magnetic moment, and thus we assumed m Z < m µ .

Majoron interactions
The spontaneous breaking of the global U(1) L symmetry gives rise to a Nambu-Goldstone boson, called the Majoron φ. If the global U(1) L symmetry is slightly broken, then the Majoron has a tiny mass : The interaction between the Majoron and neutrinos is described by where g αβ = g βα is coupling constants and ν c L,α ≡ (ν L,α ) c = Cν T L,α with the charge conjugation matrix C. As we will see later, this interaction can have a significant impact on the early universe.

Time evolution equation of temperature and chemical potential
Here we consider the thermodynamics of the early universe in the presence of new light particles, Z and the Majoron φ. In our study, we assume the following conditions : 1. As for the parameters of Z , we focus on the region of g µ−τ ∼ 10 −4 -10 −3 and m Z ∼ 10 MeV to solve the (g − 2) µ anomaly.
3. We assume that there is no primordial abundance of Majorons, and they are produced after e ± annihilation through the inverse decay process νν → φ 3 ; this assumption corresponds to looking at the parameter region satisfying Eq. (41). Boltzmann equations with simultaneous contributions from Z and φ are technically difficult to solve. We leave it for future work.
Under condition 2, the scattering and the annihilation processes of Majorons can be neglected, and only the decay and inverse decay of the Majoron φ ↔ ν α ν β ,ν ανβ are relevant to our study. Moreover, because of condition 1, Z becomes non-relativistic before e ± annihilation and decays mainly into neutrinos. On the other hand, from condition 3, Majorons are produced after e ± annihilation. Therefore, thermodynamics of Z and φ can be considered separately, before and after the temperature T γ ∼ 10 −2 MeV at which the electrons and positrons have already annihilated. In the following subsections, the evolution equations are derived for each period.

Evolution equation before e ± annihilation
We consider the evolution equations for the universe before e ± annihilation, at which photons, neutrinos, electrons, and Z exist. Following the previous studies [21,35,36], we make the following approximations in the calculation.
1. All the particles follow the thermal equilibrium distribution function.
2. In the collision terms, we use the Maxwell-Boltzmann statistics.
3. Neglect the electron mass m e in the collision terms for the weak interaction processes.
4. Neglect the chemical potentials µ i for all the particles i.

5.
The temperatures T i of a particle i in the same thermal bath are equal; T γ = T e − and T να = T Z ≡ T ν for α = e, µ, τ .
Using these approximations, we obtain the evolution equations for the temperatures of photon T γ and neutrinos T ν as follows [21]: with ρ i and P i being the energy density and pressure of particle i, respectively, and H the Hubble parameter. Here, the energy transfer rates in Eqs. (11) and (12) are given by where G F is the Fermi coupling constant, K 2 is the modified Bessel function of the second kind,

Evolution equation after e ± annihilation
We derive the evolution equations for the universe after e ± annihilation, at which photons, neutrinos, and the Majoron exist. In analogy with the previous subsection, we make the following assumptions [36].
1. All the particles follow the thermal equilibrium distribution function.
2. In the collision terms, we use the Maxwell-Boltzmann statistics.
Using these approximations, we obtain the evolution equations for temperature and chemical potential as follows 4 [36] : with n i being the number density of particle i. The number and energy density transfer rate of neutrinos are given by Since φ and ν are no longer strongly coupled to the photon in this period, their chemical potentials are no longer guaranteed to be zero. Thus, the above evolution equations for µ ν and µ φ are indispensable.

Calculation of the number and energy transfer rates
To solve the evolution equations for temperatures and chemical potentials, we need to calculate the number and the energy transfer rates. For processes φ ↔ ν α ν β and φ ↔ν ανβ , the number and the energy transfer rate of φ are described as follows [36]: Actually, in addition to the decay and inverse decay of φ, there also exist the scattering and the annihilation processes of Majoron. However, we neglect these processes because we assume |g αβ | 10 −7 as mentioned at the beginning of this section. In this case, the number transfer rate for φ is given by where Γ φ is the total decay width of φ given by where λ 2 ≡ tr(g † g). In the same way, the energy transfer rate for φ is written as The transfer rates for neutrinos, δn ν /δt and δρ ν /δt, can be obtained from the number and the energy conservation law. In the present case, the physics does not depend on a basis of neutrinos, because the neutrino masses are neglected. This is understood from the fact that Γ φ depends on g αβ only in the form tr(g † g). Thus, without loss of generality, we can assume that g αβ has only diagonal components, and the number conservation is expressed as which leads to On the other hand, the energy conservation leads to From Eq. (30), δρ ν /δt is found to be

Numerical calculation
In this section, we discuss the initial conditions and the parameters for the evolution equations of temperatures and the chemical potentials derived in the previous section and show the numerical results. The codes for calculations are partially based on NUDEC BSM [36].

Initial conditions and integration range
Before e ± annihilation We solve the system of differential equations (11) and (12) starting from T γ = T ν = 20 MeV at which all the particles are in thermal equilibrium, to T γ ∼ 10 −2 MeV where the e ± annihilation has taken place.
After e ± annihilation Let us consider solving the system of differential equations (16) This is illustrated in Figure 2 in [36]. Imposing Γ νν→φ /H < 10 −4 , we obtain the condition for T ν as follows : Here, we use the approximation K 1 (x) ∼ 1/x (for x < 1) because the situation with T ν /m φ > 1 is what we want to consider. If we set the range Γ eff ≤ 10 3 , the initial value of T ν should satisfy T ν 100m φ . Thus, we will take as the initial condition for T ν . As the initial condition for T γ /T ν , we use the numerical values after e ± annihilation (T γ 10 −2 MeV) obtained by solving equations (11) and (12). The remaining initial conditions are determined so that ρ φ /ρ ν < 10 −12 is satisfied. Since the Majoron is ultra-relativistic in T ν = 100 m φ , we can treat the Majoron as a massless particle, so ρ φ /ρ ν is written by Here, Li s (z) is Polylogarithm and a ≡ ζ(3)/ζ(4) ∼ 1.2/1.08 ∼ 1.1. Therefore, to satisfy ρ φ /ρ ν < 10 −12 , the parameters should be This means that the condition for µ φ is 8 Furthermore, since the Majoron is a boson, µ φ must satisfy from µ φ ≤ m φ . Here, the equality sign is removed because the Bose-Einstein condensation cannot occur due to the very small number density of Majoron.

Parameters
As mentioned before, we consider the case where the Majoron does not exist in the very early universe and is created after e ± annihilation (T γ 10 −2 MeV). To realize this situation, the parameters of the Majoron must satisfy the following conditions : • The Majoron production is most active after e ± annihilation.

Results
Here, we show the results of solving the evolution equations derived in the previous section. In this study, the deviation of N eff from the standard value occurs two times, namely before and after the e ± annihilation. Thus, it is convenient to write N eff as Here, N eff describes the effective number of neutrino species determined at T ν /T γ = const. soon after e ± annihilation and is defined as On the other hand, ∆N eff represents the change in the effective number of neutrino species due to the Majoron production after e ± annihilation. As we will see later, N eff and ∆N eff are not completely independent, and ∆N eff slightly depends on N eff .  Figure 2 shows the evolution of the neutrino temperature obtained by solving Eqs. (11) and (12). As can be seen from this figure, the value of N eff is slightly larger than that of the SM N SM eff 3.045 [37,38] due to the new gauge boson Z . Figure 3 shows the evolution of the neutrino energy density and the Majoron energy density for the case of N eff = 3.5. This figure shows that for Γ eff 1, the Majoron begins to be produced by νν → φ when the temperature reaches T ν m φ . After that, neutrinos and Majoron gradually reach the thermal equilibrium. This corresponds to the gently sloping area around the peak in the Figure 3. Since the net energy transfer due to φ ↔ νν is negligibly small, the evolution of the energy densities can be determined by the following Boltzmann equations: At temperature T ν m φ , the Majoron becomes non-relativistic and ρ φ becomes much larger than P φ . Consequently, the energy densities are derived as follows : where R is the scale factor. Therefore, the difference between ρ ν and ρ φ occurs as the universe expands. At temperature T ν m φ /3, Majorons start to decay into neutrinos. Since the neutrinos produced by this decay are more energetic than the existing neutrinos, the overall neutrino energy density slightly increases, resulting in a slightly larger N eff .    Figure 4 shows the evolution of the neutrino energy density for the case of N eff = 3.5, m φ = 1 keV. This figure is obtained by connecting Figure 2 and Figure 3 at T γ 10 −2 MeV smoothly. Figure 5 shows the Γ eff dependence of ∆N eff for some N eff . The parameters ∆N eff and N eff are not completely independent, and ∆N eff slightly depends on N eff . As you can see, ∆N eff becomes larger for larger N eff . The reason is as follows : A large N eff corresponds to a large number of neutrinos after e ± annihilation. For Γ eff 1, which corresponds to the case where the thermal equilibrium between the Majoron and neutrino is reached due to φ ↔ νν, this process acts to equalize the number of neutrinos and Majorons. Thus, for the larger number of neutrinos after e ± annihilation, the more neutrinos are converted to the Majorons. As a result, the neutrino energy density at T ν m φ becomes larger, yielding an increase in ∆N eff . On the other hand, for Γ eff 1, the thermal equilibrium is not achieved between ν and φ, but a small amount of Majoron is produced by νν → φ. This process occurs more often for a larger number of neutrinos after e ± annihilation. Thus, the production of Majoron increases slightly and it leads to an increase in ∆N eff . Note that the contribution of Majoron ∆N eff cannot be larger than 0.12 in the case of the SM N eff = 3.045. The dark blue region is excluded by Planck 2018 data [26]. The gray region is excluded by SN1987A [33,34], Big Bang Nucleosynthesis (BBN) [26]. The white region cannot be treated in this paper.   [39]. The brown and green regions are excluded by the Borexino and CCFR experiments, respectively [40]. 13 Using N eff and ∆N eff defined above, we can write N eff as Eq. (42). If we fix either N eff or ∆N eff , a constraint can be imposed on the other parameter by using the constraint from Planck 2018 : N eff = 3.27 ± 0.15 with 68% C.L. [11]. Although the various patterns are possible, we will only discuss the following two cases.

Summary
In this paper, we explored possibilities of resolving the Hubble tension and (g − 2) µ anomaly simultaneously in realistic U(1) Lµ−Lτ models that can explain the origin of neutrino mass. In these models, there is a new light gauge boson Z and a new light scalar, the Majoron φ. It arises from the spontaneous breaking of the global U (1) L symmetry and weakly couples to neutrinos. The parameters of Z boson are set to be 10 −3 g µ−τ 10 −4 , m Z 10 MeV, neighborhoods of region that can resolve the (g − 2) µ anomaly.
We only focused on a case where the Majoron does not exist at the beginning of the universe and is created by νν → φ after e ± annihilation. In this case, contributions of Z and φ to the effective number N eff can be calculated independently. Thus, it is convenient to write N eff as N eff = N eff + ∆N eff , a sum of the effective number after e ± annihilations N eff and its change due to the Majoron ∆N eff . The effective number N eff is evaluated by evolution equations of temperatures and the chemical potentials of light particles in each period.
For simplicity, the following two cases are discussed. First, we explored the parameter space of the Majoron in the presence of Z that realizes N eff = 3.4. In this case, the Hubble tension can be resolved (N eff 3.4 − 3.5) in the wide region of the parameter space where ∆N eff 0.1 (λ 10 −12 − 10 −14 ) holds. On the other hand, the region with ∆N eff 0.1 is excluded at more than 2σ level. In the second case, we surveyed the parameter region of Z where the Hubble tension can be resolved in the presence of Majoron that realizes ∆N eff 0.1. A choice of parameters m Z 13 − 26 MeV, g µ−τ (3.6 − 7) × 10 −4 that corresponds to N eff 3.1 − 3.4 can resolve the Hubble tension and (g − 2) µ anomaly simultaneously. On the other hand, the region with m Z 10MeV is excluded at more than 2σ level. As a result, we found that the heavier m Z results in the smaller N eff and requires the larger ∆N eff to resolve the Hubbel tension. Therefore, compared to previous studies, the parameter region where the Hubble tension can be resolved is slightly shifted toward the larger value of m Z . Note that N eff and ∆N eff are not completely independent, and ∆N eff slightly depends on N eff .
Finally, Boltzmann equations with simultaneous contributions from Z and φ are more difficult to solve. We leave it for future work.

A Derivation of the evolution equation after e ± annihilation
Here, we derive the evolution equations (16) to (19) after e ± annihilation. First of all, the evolution equations for the temperature T a and chemical potential µ a of a particle species a that follows the thermal equilibrium distribution function are given by [36] In Eqs. (47) and (48), n a , ρ a , P a are the particle number density, energy density, and pressure of a particle species a, respectively. From the approximations 3 in subsection 3.2 and T να = Tν α , µ να = µν α , Eq. (47) for the neutrino and antineutrino leads to In addition, each thermodynamic quantity for {ν α }, {ν α } is expressed by the particle number density n ν , the energy density ρ ν , and the pressure P ν for the total neutrino : By adding both sides of Eqs. (49) and (50), and summing over all flavors, we obtain the evolution equation for T ν (16). The evolution equation for µ ν (17) can also be obtained in the same way. For the Majoron evolution equation, by using eqs.(47) and (48) set to a = φ, we obtain the evolution equations for T φ (18) and µ φ (19).