Self-interacting dark matter and muon $g-2$ in a gauged U$(1)_{L_{\mu} - L_{\tau}}$ model

We construct a self-interacting dark matter model that could simultaneously explain the observed muon anomalous magnetic moment. It is based on a gauged U$(1)_{L_{\mu} - L_{\tau}}$ extension of the standard model, where we introduce a vector-like pair of fermions as the dark matter candidate and a new Higgs boson to break the symmetry. The new gauge boson has a sizable contribution to muon $(g-2)$, while being consistent with other experimental constraints. The U$(1)_{L_{\mu} - L_{\tau}}$ Higgs boson acts as a light force carrier, mediating dark matter self-interactions with a velocity-dependent cross section. It is large enough in galaxies to thermalize the inner halo and explain the diverse rotation curves and diminishes towards galaxy clusters. Since the light mediator dominantly decays to the U$(1)_{L_{\mu} - L_{\tau}}$ gauge boson and neutrinos, the astrophysical and cosmological constraints are weak. We study the thermal evolution of the model in the early Universe and derive a lower bound on the gauge boson mass.


I. INTRODUCTION
Dark matter (DM) makes up 85% of the mass in the Universe, but its nature remains largely unknown. There has been great progress in studying particle DM candidates associated with extensions of the standard model (SM) of particle physics, which can be probed in high-energy and intensity-frontier terrestrial experiments (see, e.g., Ref. [1]). Astrophysical and cosmological observations can also provide important clues to the nature of DM. In fact, a number of astrophysical observations indicate that the cold dark matter (CDM) model may break down on galactic scales (see Refs. [2,3] for a review), although it works remarkably well in explaining large-scale structure of the Universe, from Mpc to Gpc scales. For example, the galactic rotation curves of spiral galaxies exhibit a great diversity in inner shape [4], which is hard to understand in CDM. It has been shown that the diverse rotation curves can be explained naturally if DM has strong self-interactions [5,6]. The analysis has been extended to the dwarf spheroidal galaxies in the Milky Way [7] and other galactic systems [8,9]. This self-interacting dark matter (SIDM) scenario has rich implications for understanding the stellar kinematics of dwarf spheroidal galaxies and galaxy clusters and interpreting DM direct, indirect, and collider search results (see, e.g., Ref. [2]).
On the particle physics side, there is also a long-standing puzzle that the measured muon anomalous magnetic moment, (g − 2) µ , is larger than predicted in the SM at the 3σ level [10,11]. This discrepancy may indicate that there is new physics beyond the SM associated with the muon sector. In this work, we propose a model that could simultaneously explain the (g − 2) µ measurement and provide a realization of SIDM. It is based on a gauged U(1) Lµ−Lτ extension of the SM (see, e.g., Refs. [12,13]). Aside from the new gauge (Z ) and Higgs (ϕ) bosons related to the symmetry, we introduce a vector-like pair of fermions (N andN ) as the DM candidate and assume that they couple to ϕ via a Yukawa interaction. The model overcomes a number of challenges in SIDM model building [14] (see also, e.g, Refs. [15][16][17][18][19][20][21][22][23][24][25][26][27][28]). Our main observations are the following: • The presence of Z could contribute to (g − 2) µ . If Z has a mass in the range of m Z ∼ 10-100 MeV, there is a viable parameter space to explain the (g − 2) µ anomaly.
• Astrophysical observations indicate that the DM self-scattering cross section per unit mass is σ/m 1 cm 2 /g in galaxies, while diminishing to σ/m 0.1 cm 2 /g in galaxy clusters [8]. The desired velocity-dependence of σ/m can be achieved if m ϕ ∼ 10-100 MeV.
• In the early Universe, the mediator ϕ is in equilibrium with Z and the SM neutrinos (ν µ and ν τ ), and its number density is Boltzmann suppressed when the temperature is below its mass. Thus, the model avoids over-closing the Universe.
• Since Z and ϕ dominantly decay to the neutrinos, the big bang nucleosynthesis (BBN), cosmic microwave background (CMB), and indirect detection bounds are significantly weaker, compared to the models with electrically charged final states [29,30].
This paper is organized as follows. We introduce the U(1) Lµ−Lτ model and discuss relevant experimental and cosmological constraints in Sec. II. We discuss DM phenomenology of the model in Sec. III and devote Sec. IV to conclusions. In the appendix, we provide detailed discussion on the neutrino-electron scattering cross section induced by Z (A), the temperature evolution of the model in the early Universe (B), and the decay width of ϕ (C).

II. THE MODEL
We introduce a vector-like pair of fermions N andN and a U(1) Lµ−Lτ Higgs Φ in addition to the U(1) Lµ−Lτ gauge boson Z µ . The charge assignments are summarized in Table I. The renormalizable Lagrangian density can be written as The covariant derivative on U(1) Lµ−Lτ -charged particles is written as where Q is a U(1) Lµ−Lτ charge. The field strengths of Z µ and the U(1) Y gauge boson B µ are denoted by Z µν and B µν , respectively. The scalar potential of Φ takes a form of We set λ ΦH = 0 so that it does not induce the DM interaction with SM particles through the Higgs portal. Φ develops a vacuum expectation value (VEV) and can be expanded as The resultant masses of Z and ϕ are given by m Z = g v Φ and m ϕ = λ Φ /2 v Φ , respectively. Note that the VEV of Φ breaks U(1) Lµ−Lτ into a Z 2 symmetry, which stabilizes the lightest state of N andN , i.e., the DM candidate. In our model, ϕ has a mass of O(10) MeV and plays a role of the SIDM mediator.
We assume that there is no kinetic mixing between Z and B at some high-energy scale, = 0. This can be achieved by imposing a charge conjugation symmetry C Lµ−Lτ : L 2 ↔ L 3 , µ ↔ τ , N ↔N , Z → −Z , and Φ → Φ * . However, since this symmetry is broken by the Yukawa mass terms of µ and τ , the mixings of Z with the photon A and with the Z boson arise at the 1-loop level at low energy. The A-Z mixing below m µ is where e is the electric charge of the electron and m µ and m τ are the mu and tau lepton masses, respectively. To obtain the canonical gauge fields, we redefine the field as A → A + AZ Z for | AZ | 1. It induces a coupling between Z and the SM electromagnetic current: While the Z-Z mixing is given by where c W = cos θ W and s W = sin θ W with θ W being the Weinberg angle. After performing the following field shifts: and the coupling of Z to the SM neutral current, Since the coupling of Z to the neutral current is suppressed by m 2 Z /m 2 Z ∼ 10 −8 , its contribution to Z phenomenology is negligible. In the rest of this section, we discuss the observational constraints on Z and ϕ.
• Υ decays. The BABAR experiment has searched for the eē → µμZ process followed by Z → µμ at the Υ resonance [42], resulting an upper limit on g for 200 MeV m Z 10 GeV.
• ν-e scattering. The Borexino experiment measures the interaction rate of the monoenergetic 862 keV 7 Be solar neutrino [43]. It puts a strong constraint on models that predict new ν-e interactions. In this work, we calculate the ν-e reaction rate by taking into account the Z contribution (see appendix A for details). We require that the total reaction rate deviates from the SM prediction no more than 8% [43] and obtain the bound.
• White dwarf (WD) cooling. The white dwarf cooling also gives a strong limit [44]. The plasmon in the white dwarf star could decay to neutrinos through off-shell Z , and this process increases the cooling efficiency. In the SM, the effective operator where G F is the Fermi constant and C V 0.964 [45], induces the plasmon decay through the electron loop. In our model, we have a similar operator, It is remarkable that given these constraints there is still a viable parameter space for explaining (g − 2) µ , i.e., 5 MeV m Z 200 MeV and 3.0 × 10 −4 g 1.1 × 10 −3 . Since Z does not couple to the electron directly, our model is much less constrained when compared to others invoking a new gauge boson, such as the U(1) B−L model [46] and the hidden photon model with the kinetic mixing (see, e.g., Ref. [47]). In what follows, we take two model examples: (m Z , g ) = (10 MeV, 5 × 10 −4 ) and (50 MeV, 7 × 10 −4 ) as indicated in Fig. 1 (green stars).

B. Cosmological constraint
The main decay channel of Z is Z → νν, and its lifetime is For the model parameters that we are interested in, τ Z is much shorter than the BBN timescale, t ∼ 1 sec (1 MeV/T ) 2 . Thus, Z is in equilibrium with ν µ and ν τ through the decay and inverse decay, after neutrinos decouple from the SM plasma, T = T ν-dec 1.5 MeV. Its energy density is suppressed by the Boltzmann factor at T m e , m Z , and its direct contribution to the effective number of neutrino degrees of freedom N eff is negligible. Note N eff is related to the total energy density (ρ tot ) and the photon density (ρ γ ) as Even in this case, the presence of Z may modify N eff , since Z injects energy only into ν µ and ν τ and it changes the ratio of the neutrino temperature to the photon temperature. In the limit of AZ =0, we follow the analysis in Ref. [48] to estimate the effect. We assume that (γ, e), (ν e ), and (ν µ , ν τ , Z ) form independent thermal baths after the neutrino decoupling, and the comoving entropy density is conserved in each bath. By following the temperature evolution in each sector, we evaluate N eff at T m e , m Z . In Fig. 1 (right), we show the lower bound on m Z (blue), i.e., m Z 6 MeV, where we take the upper limit N eff < 3.5 [38]. In the presence of a sizable AZ as in Eq. (3), Z can decay to eē, and there is heat transfer between (ν µ , ν τ , Z ) and (γ, e) baths. In this case, we implement the heat transfer through Z → eē in the evolution of the comoving entropy density as where T , T ν , and T , denote the temperatures of (γ, e), (ν e ), and (ν µ , ν τ , Z ), respectively, and Γ Z →eē is the decay width of Z → eē, where α is the fine-structure constant and m e is the electron mass. Note that ρ Z for the decay is evaluated with T , while ρ Z for the inverse decay is evaluated with T . We numerically solve Eqs. (11)-(13) with the Hubble expansion rate, where m Pl is the reduced Planck mass. See appendix B for detailed discussion on the temperature evolution. We obtain N eff as shown in Fig. 1 (right). For N eff < 3.5 and AZ 7.2×10 −6 (orange), corresponding to g = 5×10 −4 , the lower bound is m Z 10 MeV, and stronger than the limit for AZ = 0.
In this model, ϕϕ ↔ Z Z and ϕν ↔ Z ν keep ϕ in equilibrium with Z and thus with SM particles below T ∼ 20 GeV for g = 5 × 10 −4 . To be conservative, we can assume that m ϕ > m Z so that the lifetime of ϕ is less than 1 s (see appendix C). In this case, the presence of ϕ does not change the lower bound on m Z as shown in Fig. 1 (right). We comment on the possibility of m ϕ < m Z , where ϕ decouples from the thermal bath when ϕν ↔ Z ν becomes inefficient. In this case, ϕ dominantly decays to four neutrinos long after the BBN. It can also decay to two neutrinos and two electrons with a small branching ratio ∼ e 2 2 AZ /g 2 . The latter process could inject electromagnetic energy to the plasma and lead to observational consequences. If the ϕ lifetime is longer than 10 6 s, the energy injection could distort the CMB spectrum from the blackbody radiation (see, e.g., Ref. [49]). We find that m ϕ has a lower limit, i.e., m ϕ 6 MeV, to satisfy the COBE constraint [50] for (m Z , g ) = (10 MeV, 5 × 10 −4 ), while m ϕ 20 MeV for (m Z , g ) = (50 MeV, 7 × 10 −4 ). In addition, the observation of the light-element abundances can further tighten the constraint (see, e.g., Ref. [51]). A detailed study of the cosmological constraint is beyond the scope of this work, and we will take the conservative assumption in the rest of the paper.

III. DARK MATTER PHENOMENOLOGY
In this section, we discuss DM phenomenology predicted in this model. The mass matrix of the DM candidates N andN can be diagonalized by a unitary matrix U as where 0 ≤ M 1 ≤ M 2 and Then the Lagrangian density in the 4-component spinor notation, where N i denotes a Majorana fermion, becomes We consider two extreme cases to illustrate the main predictions of our model and further simplify the Lagrangian density.
• Pseudo-Majorana DM, where |yN |v Φ > |y N |v Φ |m N |. After performing the phase rotation of N andN such that both y N and yN are positive, we have Note that DM phenomenology is unchanged for |y N |v Φ > |yN |v Φ |m N |. In this limit, it turns out that we cannot simultaneously obtain the SIDM cross section, explain the (g − 2) µ discrepancy, and realize viable cosmology. We will not focus on this case in the rest of the paper.
• Pseudo-Dirac DM, where |m N | |y N |v Φ , |yN |v Φ . We choose a basis such that m N and y * N + yN are positive and further restrict our discussion to the case where y N = yN ≡ y > 0. In this limit, the DM sector respects C Lµ−Lτ (N ↔N and Φ → Φ * ) and the parity conjugation (N → iN † andN → iN † ). The mass eigenstates are N 1 = (N −N )/( √ 2i) and N 2 = (N +N )/ √ 2, and the corresponding masses are

. The gauge and Yukawa interactions are
To explore DM phenomenology predicted in this pseudo-Dirac case, we need to specify five parameters: m Z , g , m ϕ , M 1 , and y. For the first two, we take the examples motivated by the (g − 2) µ measurement as denoted by the green stars in Fig. 1 (left). We then scan the (m ϕ , M 1 ) parameter space, while fixing y by the relic density requirement accordingly, as discussed next.

A. DM relic abundance
We determine the Yukawa coupling y such that the N 1 relic density gives rise to the observed DM abundance. We assume that M 1 ∼ M 2 100 GeV so that the freeze-out occurs after the U(1) Lµ−Lτ symmetry breaking, T < v Φ = O(10) GeV, and we can neglect the thermal corrections to ∆M = M 2 −M 1 = √ 2yv Φ , m ϕ , and m Z . 2 Since ∆M = O(1) GeV, N 1 and N 2 are in equilibrium with each other during the freeze-out and we need to take into account their co-annihilation process [52]. The annihilation channels are N 1 N 1 , N 2 N 2 → Z Z , ϕϕ and N 1 N 2 → ϕZ , and the corresponding thermally-averaged cross sections times relative velocity can be written as where x = M 1 /T . In deriving Eq. (20), we have ignored the terms depending on g , since they are too small to play a role. The effective annihilation cross section is σ ann v rel eff = σ ann v rel 11 r 2 1 + σ ann v rel 22 r 2 2 + 2 σ ann v rel 12 r 1 r 2 , where r 1 = 1/ 1 + e −∆x , r 2 = e −∆x / 1 + e −∆x , and ∆ ≡ ∆M/m N . We equate this effective annihilation cross section to the canonical one, σ ann v rel can = 3 × 10 −26 cm 3 /s, at x f = 20 and derive the y value for given (m ϕ , M 1 ). We have checked that the Sommerfeld effect is negligible for the freeze-out calculation (see, e.g., Refs. [53,54]) for the DM mass range that we are interested in.

B. DM self-scattering
The elastic scattering of N 1 is dominated by the Yukawa interaction in Eq. (19), which is represented by the following Yukawa potential for non-relativistic scattering: where the y value is fixed by the relic density consideration as discussed above. The scattering amplitude f (θ) is given by the partial wave expansion as where δ is the phase shift and P is the th order Legendre polynomial. We numerically calculate δ by solving the Schrödinger equation as in Ref. [55]. Since N 1 is a Majorana fermion, the scattering particles are indistinguishable. Then the differential cross section is given by the sum of spin-singlet and spin-triplet channel amplitudes as where the initial state is assumed to be unpolarized [56]. The standard cross section, σ = dΩ (dσ/dΩ), is not an appropriate quantity for characterizing the effect of self-scattering on the structure formation, since the scattering through a light mediator receives strong enhancement at θ = 0 and π, which do not affect the DM distribution. Instead the transfer cross section is often used in the literature [2], which is written for the scattering of identical particles as [57] σ T = 4π In Fig. 2, we show the parameter regions for σ T /M 1 = 0.1-1 cm 2 /g (light blue) and 1-10 cm 2 /g (blue) in dwarf galaxies, 0.1-1 cm 2 /g (red) in Milky Way-size galaxies, and > 0.1 cm 2 /g (green) in galaxy clusters, where we take v rel = 30, 200, and 1000 km/s, respectively. We have used the Born approximation when it is valid. We see that with a light mediator, i.e., m ϕ = O(10) MeV, the SIDM cross section can be large in galaxies, while diminishing towards galaxy clusters, as shown in both panels. However, when we impose the cosmological constraint (gray), i.e., m ϕ > m Z , most of the SIDM parameter space is excluded for the case of g = 7 × 10 −4 and m Z = 50 MeV (right). Figure 3 shows the velocity-dependence of the self-scattering cross section σ T v rel /M 1 , averaged over the Maxwell-Boltzmann distribution, for the three benchmark cases marked in Fig. 2 (left), together with the inferred values from stellar kinematics of dwarf (red) and low surface brightness (blue) galaxies, and galaxy clusters (green) [8]. We see that at least two cases (solid and dashed) have a large self-scattering cross section to address the diversity The points with errors are taken from Ref. [8], corresponding to dwarf galaxies (red), low surface brightness galaxies (blue), and galaxy clusters (green).
problem of galactic rotation curves. On cluster scales, all three cases have a cross section below ∼ 0.1 cm 2 /g, as required by the cluster constraints [8]. Interestingly, this upper limit is preferred if the DM self-interactions also explain the inferred density cores in the galaxy clusters [58]. The case with M 1 = 10 GeV may explain observations of stellar kinematics in both galaxies and galaxy clusters. As indicated in Fig. 2, all three benchmark cases are in the resonance regime, where the self-scattering cross section has a strong velocity dependence [59], i.e., σ T /M 1 ∝ 1/v 2 rel . Thus, the cross sections decrease from 1 cm 2 /g in dwarf galaxies to < 0.1 cm 2 /g in galaxy clusters. This is a generic feature of the model due to the tight constraints. Combining the muon (g − 2) and the Borexino constraints, we require m Z 10 MeV. We further demand m ϕ m Z so that ϕ decays before the onset of the BBN, and the Yukawa coupling constant y to be fixed by the DM relic abundance. After imposing all these constraints, we find that our SIDM model is in the resonant regime.

C. Direct and indirect searches
If the SIDM mediator decays through the Higgs portal, DM direct detection experiments have put a strong constraint on the mixing parameter λ ΦH [14,56,60] and the lifetime of the mediator is too long to be consistent with early Universe cosmology. In the model we consider, ϕ decays to Z to avoid over-closing the Universe. Thus, it does not need to mix with the SM Higgs boson to deplete its abundance. We set λ ΦH = 0 at the tree level. A finite value of λ ΦH can be generated through a 2-loop process, but it is small and much below the experimental sensitivity. Moreover, the spin-dependent DM-nucleon scattering is inelastic, since the relevant interaction term involves two states, as shown in Eq. (19). In this model, the predicted mass splitting is ∆M > 1 GeV (see Fig. 2) which is much larger than the nuclei recoil energy in DM direct detection. Thus, DM direct searches do not constrain this model.
For indirect detection, the annihilation processes, N 1 N 1 → Z Z and ϕϕ, are p-wave dominated, and their final states decay to the SM neutrinos. Since the upper bound on the DM annihilation cross section to the neutrinos is σ ann v rel 10 −24 cm 3 /s [61], our model easily avoids the indirect detection constraints.
Before closing this section, we comment on the scenario where y N = yN . In this case, there is a pseudo-scalar interaction, iϕN 1 γ 5 N 1 , in addition to the scalar one. This new term induces a DM-DM potential that is suppressed by a factor of m 2 ϕ /m 2 N 1 , compared to the Yukawa one. But, it is strongly singular as V (r) ∝ 1/r 3 . We expect it has a subdominant effect on DM self-interactions since in the Born approximation limit the pseudo-scalar cross section vanishes for v rel → 0. For the gauge interactions, we now need to include the DM axial current, −(g κ/2)Z µ N 1 γ 5 γ µ N 1 , with κ ∼ yv Φ /m N . It induces a DM-nucleon scattering cross section that decreases with the DM velocity. The current constraint, e.g., from the LUX experiment [62] is satisfied if |κ| < O(10 2 ) for g = 5 × 10 −4 with the momentum transfer ∼ 100 MeV [63].

IV. CONCLUSION
We have constructed a U(1) Lµ−Lτ model that could explain the (g −2) µ measurement and reconcile the discrepancies of CDM on galactic scales. In this model, the U(1) Lµ−Lτ Higgs boson acts as the light dark force carrier, mediating DM self-interactions. We thoroughly studied the constraints from the high-energy and intensity-frontier experiments and cosmological and astrophysical observations, and found a viable model parameter space. Since the mediator dominantly decays to the neutrinos, the model avoids tight constraints from DM indirect searches. To be consistent with the cosmological upper bound, N eff < 3.5, we found that the mediator ϕ and gauge boson Z masses satisfy the condition of m ϕ > m Z 10 MeV. Interestingly, with the same mass range of Z , the model could also explain the dip feature [64,65] in the IceCube neutrino spectrum [48,66,67]. Since the Higgs mixing term vanishes, our model does not lead to direct detection signals via the Higgs portal. In addition, our model naturally predicts a GeV mass splitting between two DM mass states, which forbids DM-nucleon scattering via the gauge boson-photon kinetic mixing.
Our model could be tested in future collider experiments, such as the facilities in the intensity frontier [68][69][70][71][72]. In addition, the observation of a nearby supernova may also provide important hints on the properties of Z , since the non-standard neutrino self-interaction shortens the mean free path and increases the neutrino diffusion time from the supernova core [48]. On the model building aspect, we could introduce three right-handed neutrinos and realize the see-saw mechanism to reproduce the observed neutrino masses and mixings [73]. We provide details on the constraint from the neutrino-electron scattering discussed in Sec. II A (see also, e.g., Refs. [74,75]). The electron-neutrino scattering cross sections are where g V = −1/2+2s 2 W and g A = −1/2. The upper (lower) sign corresponds to the neutrino (anti-neutrino) and α = µ, τ . The incident neutrino energy is E ν = 862 keV for the 7 Be solar neutrino. T is the electron recoil energy and its maximal value is T max = 2E 2 ν /(m e + 2E ν ). The minimal value of T is T min 270 keV for the Borexino experiment [43,74]. Then, the total reaction rate R tot is given by where t exp and ρ e are the exposure time and the electron number density of the target, respectively. We assume that the neutrino spectrum dΦ νe /dE ν is mono-energetic at E ν = 862 keV. Since we are only interested in the ratio between the total reaction rates with and without Z , the overall normalization is irrelevant for our discussion. We define the total cross section used in Eq. (A4) as σ tot (νe) = P ee σ SM (ν e e) + (1 − P ee ) α [σ SM (ν α e) + σ Z (ν α e)] , with the electron neutrino survival probability P ee 0.51 [43]. For the SM prediction, one drops the σ Z contribution. Z → eē lowers T , while raises T , when compared to the case of the SM. This is because heat from eē annihilation is partially transferred from (γ, e) to (ν µ , ν τ , Z ).

Appendix C: Decay width of ϕ
For m ϕ ≥ 2m Z , the dominant decay channel is ϕ → Z Z and the decay width is While for m Z ≤ m ϕ < 2m Z , the 3-body decay process, i.e., ϕ → Z νν, dominates and the width is given by where y min = r 2 − r 4 − 2r 2 (x + 1) + (x − 1) 2 − x + 1 2 , In Fig. 4 (right), we show the lifetime of ϕ vs its mass. For m ϕ > m Z , τ ϕ is less than 1 s.