Cosmology With a Very Light $L_\mu - L_\tau$ Gauge Boson

TL;DR

This paper investigates how a very light $L_-L_ au$ gauge boson can influence early universe cosmology, potentially resolving the Hubble tension and explaining the muon g-2 anomaly through specific parameter space regions.

Contribution

It identifies two distinct parameter space regions where a light $Z'$ can relax the Hubble tension and account for the muon g-2 discrepancy, including new decay and production mechanisms.

Findings

A $10-20$ MeV $Z'$ can heat neutrinos and delay decoupling.

Very light ($1 eV$ to MeV) $Z'$ produced via freeze-in can increase radiation density.

The model can relax the Hubble tension with $ \u2202 N_{eff} \u2248 0.2$.

Abstract

In this paper, we explore in detail the cosmological implications of an abelian $L_\mu-L_\tau$ gauge extension of the Standard Model featuring a light and weakly coupled $Z'$. Such a scenario is motivated by the longstanding $\sim \, 4 \sigma$ discrepancy between the measured and predicted values of the muon's anomalous magnetic moment, $(g-2)_\mu$, as well as the tension between late and early time determinations of the Hubble constant. If sufficiently light, the $Z'$ population will decay to neutrinos, increasing the overall energy density of radiation and altering the expansion history of the early universe. We identify two distinct regions of parameter space in this model in which the Hubble tension can be significantly relaxed. The first of these is the previously identified region in which a $\sim \, 10-20$ MeV $Z'$ reaches equilibrium in the early universe and then decays, heating the neutrino population and delaying the process of neutrino decoupling. For a coupling of $g_{\mu-\tau} \simeq (3-8) \times 10^{-4}$, such a particle can also explain the observed $(g-2)_{\mu}$ anomaly. In the second region, the $Z'$ is very light ($m_{Z'} \sim 1\,\text{eV}$ to $\text{MeV}$) and very weakly coupled ($g_{\mu-\tau} \sim 10^{-13}$ to $10^{-9}$). In this case, the $Z'$ population is produced through freeze-in, and decays to neutrinos after neutrino decoupling. Across large regions of parameter space, we predict a contribution to the energy density of radiation that can appreciably relax the reported Hubble tension, $\Delta N_{\rm eff} \simeq 0.2$.

Equations

$$\mathcal{L}={\mathcal{L}}_{\mathrm{SM}}-\frac{1}{4}{Z^{\prime }}^{\alpha \beta }{Z}_{\alpha \beta }^{\prime }+\frac{m_{Z^{\prime }}^2}{2}{Z}_{\alpha }^{\prime }{Z^{\prime }}^{\alpha }+Z_{\alpha }^{\prime }{J}_{\mu -\tau }^{\alpha },$$

$$J_{\mu -\tau }^{\alpha }=g_{\mu -\tau }\left(\mathrm{bar}\mu {\gamma }^{\alpha }\mu +\mathrm{bar}{\nu }_{\mu }{\gamma }^{\alpha }{P}_L{\nu }_{\mu }-\mathrm{bar}\tau {\gamma }^{\alpha }\tau -\mathrm{bar}{\nu }_{\tau }{\gamma }^{\alpha }{P}_L{\nu }_{\tau }\right),$$

$${\Gamma }_{Z^{\prime }\to {\ell }^{+}{\ell }^{-}}=\frac{g_{\mu -\tau }^2{m}_{Z^{\prime }}}{12\pi }\left(1+\frac{2m_{\ell }^2}{m_{Z^{\prime }}^2}\right)\sqrt{1-\frac{4m_{\ell }^2}{m_{Z^{\prime }}^2}},{\Gamma }_{Z^{\prime }\to \mathrm{bar}{\nu }_i{\nu }_i}=\frac{g_{\mu -\tau }^2{m}_{Z^{\prime }}}{24\pi },$$

$$\Delta a_{\mu }^{Z^{\prime }}=\frac{g_{\mu -\tau }^2}{8{\pi }^2}{\int }_0^{1}dz\frac{2m_{\mu }^{2}z(1-z{)}^2}{m_{\mu }^2(1-z{)}^2+m_{Z^{\prime }}^{2}z}\simeq 1.3\times 10^{-10}{\left(\frac{g_{\mu -\tau }}{10^{-4}}\right)}^2,$$

$$\mathcal{L}\supset Z_{\alpha }^{\prime }\left(J_{\mu -\tau }^{\alpha }+ϵJ_{\mathrm{EM}}^{\alpha }\right),J_{\mathrm{EM}}^{\alpha }=e\sum _f{Q}_f\mathrm{bar}f{\gamma }^{\alpha }f,$$

$$ϵ=-\frac{e{g}_{\mu -\tau }}{2{\pi }^2}{\int }_0^{1}dxx(1-x)\log \left[\frac{m_{\tau }^2-x(1-x)q^2}{m_{\mu }^2-x(1-x)q^2}\right]\underset{m_{\mu }\gg q}{\to }-\frac{e{g}_{\mu -\tau }}{12{\pi }^2}\log \frac{m_{\tau }^2}{m_{\mu }^2}\simeq -\frac{g_{\mu -\tau }}{70}.$$

$${\Gamma }_{Z^{\prime }\to e^{+}{e}^{-}}=\frac{(ϵe{)}^2{m}_{Z^{\prime }}}{12\pi }\left(1+\frac{2m_e^2}{m_{Z^{\prime }}^2}\right)\sqrt{1-\frac{4m_e^2}{m_{Z^{\prime }}^2}},$$

$${\mathrm{Br}}_{Z^{\prime }\to e^{+}{e}^{-}}=\frac{{\Gamma }_{Z^{\prime }\to e^{+}{e}^{-}}}{{\Gamma }_{Z^{\prime }\to \mathrm{bar}{\nu }_{\mu }{\nu }_{\mu }}+{\Gamma }_{Z^{\prime }\to \mathrm{bar}{\nu }_{\tau }{\nu }_{\tau }}}\simeq {\left(\frac{eϵ}{g_{\mu -\tau }}\right)}^2\simeq 2\times 10^{-5},$$

$${\dot{n}}_{Z^{\prime }}+3H{n}_{Z^{\prime }}=⟨{\Gamma }_{Z^{\prime }}⟩\left(n_{Z^{\prime }}^{(\mathrm{eq})}-n_{Z^{\prime }}\right),$$

$$⟨{\Gamma }_{Z^{\prime }}⟩\equiv {\Gamma }_{Z^{\prime }}\frac{K_1(x)}{K_2(x)},$$

$${\rho }_R=\left[1+\frac{7}{8}{\left(\frac{4}{11}\right)}^{4/3}{N}_{\mathrm{eff}}\right]{\rho }_{\gamma },$$

$$n_{Z^{\prime }}^{(\mathrm{eq})}={\int }_0^{\infty }\frac{d^3\vec{p}}{(2\pi {)}^3}\frac{g_{Z^{\prime }}}{e^{E/T}-1},$$

$$N_{\mathrm{eff}}=\frac{8}{7}{\left(\frac{11}{4}\right)}^{4/3}\frac{{\rho }_{\nu }}{{\rho }_{\gamma }},$$

$$\frac{d{T}_{{\nu }_e}}{dt}$$

$$\frac{d{T}_{{\nu }_{\mu }}}{dt}$$

$$\frac{d{T}_{\gamma }}{dt}$$

$$\frac{\delta {\rho }_{Z^{\prime }}}{\delta t}$$

$$\frac{\delta {\rho }_{{\nu }_e}}{\delta t}$$

$$\frac{\delta {\rho }_{{\nu }_{\mu }}}{\delta t}$$

$$F(a,b)\equiv 16(a^9-b^9)+7a^4{b}^4(a-b),$$

$$\frac{d{T}_{\nu }}{dt}$$

$$\frac{d{T}_{\gamma }}{dt}$$

$$g_{\mu -\tau }\text{raise}1.29167pt>\text{lower}4.30554pt\sim {\left(\frac{1.66\sqrt{g_{\ast }}}{\alpha }\frac{m_{\mu }}{m_{\mathrm{Pl}}}\right)}^{1/2}\simeq 4\times 10^{-9},$$

$$3{\rho }_{\nu }(T_{\nu },0)$$

$$3n_{\nu }(T_{\nu },0)$$

$$\frac{T_{\mathrm{eq}}}{T_{\nu }}=1.2076,$$

$$\frac{\partial f_{Z^{\prime }}}{\partial t}-H{p}_{Z^{\prime }}\frac{\partial f_{Z^{\prime }}}{\partial p_{Z^{\prime }}}$$

$$\frac{\partial f_{\nu }}{\partial t}-H{p}_{\nu }\frac{\partial f_{\nu }}{\partial p_{\nu }}$$

$$F_{\mathrm{dec}}\left(E_{Z^{\prime }},E_{{\nu }_1},E_{{\nu }_2}\right)=f_{Z^{\prime }}(E_{Z^{\prime }})\left[1-f_{\nu }(E_{{\nu }_1})\right]\left[1-f_{\nu }(E_{{\nu }_2})\right]-f_{\nu }(E_{{\nu }_1})f_{\nu }(E_{{\nu }_2})\left[1+f_{Z^{\prime }}(E_{Z^{\prime }})\right].$$

$$\displaystyle\frac{dT_{\gamma}}{dt}=-\biggl{[}4H\rho_{\gamma}+3H(\rho_{e}+p_{e})\biggr{]}\left(\frac{\partial\rho_{\gamma}}{\partial T_{\gamma}}+\frac{\partial\rho_{e}}{\partial T_{\gamma}}\right)^{-1}.$$

$${\Gamma }_{Z^{\prime }}\ge H(T_{\nu }=m_{Z^{\prime }})⟹$$

$${\Gamma }_{Z^{\prime }}=H(T_{\nu }=m_{Z^{\prime }})/5^2⟹$$

$${\Gamma }_{Z^{\prime }}=H(T_{\nu }=m_{Z^{\prime }})/17^2⟹$$

$$\frac{d{n}_{Z^{\prime }}}{dt}+3H{n}_{Z^{\prime }}=C[f_{Z^{\prime }}],$$

$$C[f_{Z^{\prime }}]=\int \frac{d^3{\vec{p}}_1}{(2\pi {)}^{3}2E_1}\frac{d^3{\vec{p}}_2}{(2\pi {)}^{3}2E_2}\frac{d^3{\vec{p}}_{Z^{\prime }}}{(2\pi {)}^{3}2E_{Z^{\prime }}}|{\mathcal{A}}_{Z^{\prime }\to \nu \mathrm{bar}\nu }{|}^2\left[f_{\nu }({\vec{p}}_1)f_{\nu }({\vec{p}}_2)-f_{Z^{\prime }}({\vec{p}}_{Z^{\prime }})\right](2\pi {)}^4{\delta }^4(p_1+p_2-p_{Z^{\prime }}),$$

Full text

††preprint:

FERMILAB-PUB-19-001-A, LPT-Orsay-18-15, IFIC-19-02, KCL-19-01, IFT-UAM/CSIC-19-7

1 Introduction

Light and weakly coupled particles are found within many well-motivated extensions of the Standard Model (SM) Essig:2013lka ; Alexander:2016aln . Among the motivations for such states is the muon’s anomalous magnetic moment, whose measured value constitutes a $4.1\sigma$ discrepancy with respect to the predictions of the SM Tanabashi:2018oca . This anomaly has inspired many explanations involving new sub-GeV particles Pospelov:2008zw ; Altmannshofer:2014pba . Although many of these explanations have already been ruled out (including both visibly and invisibly decaying dark photons with kinetic mixing $\epsilon\sim 10^{-3}$), scenarios in which new states couple predominantly to muons remain viable.

In addition to the muon’s magnetic moment, there is also a $\sim$ 3$\sigma$ tension between the values of the Hubble constant, $H_{0}$, as determined from local measurements Riess:2018byc ; Riess:2016jrr and from the temperature anisotropies of the cosmic microwave background (CMB) Aghanim:2018eyx . Such a discrepancy could be ameliorated if the expansion rate of the universe departed from the predictions of standard $\Lambda$CDM cosmology at early times Bernal:2016gxb ; Aylor:2018drw . Particularly well motivated within this context are scenarios in which the energy density in relativistic particles exceeds that predicted by the SM, generally parameterized in terms of a non-zero contribution to the effective number of neutrino species, $\Delta N_{\rm eff}$ Weinberg:2013kea ; Shakya:2016oxf ; Berlin:2018ztp ; DEramo:2018vss ; Dessert:2018khu . Scenarios involving early dark energy Poulin:2018cxd ; Poulin:2018dzj ; Poulin:2018zxs or a component of decaying dark matter Bringmann:2018jpr have also been proposed to address this tension. Upcoming CMB measurements will be significantly more sensitive to the value of $\Delta N_{\rm eff}$ Abazajian:2016yjj , providing us with further motivation to consider this class of scenarios.

Particularly interesting with this context are models with a broken $U(1)_{L_{\mu}-L_{\tau}}$ gauge symmetry, corresponding to a new massive gauge boson that couples to muons, taus, and their corresponding neutrinos. This is one of the few anomaly free $U(1)$ gauge extensions of the SM He:1990pn ; He:1991qd , and is the only such model without tree-level couplings to first generation quarks and/or leptons.111Other choices for gauging global SM quantum numbers are $B-L,L_{i}-L_{j},B-3L_{i}$, where $B$ and $L$ are respectively baryon and lepton numbers and a subscript denotes a specific lepton flavor. For this reason, $U(1)_{L_{\mu}-L_{\tau}}$ models are relatively unconstrained and lead to qualitatively different phenomenological and cosmological consequences Altmannshofer:2014pba ; Altmannshofer:2014cfa ; Kamada:2015era ; Kamada:2018zxi ; Foldenauer:2018zrz ; Asai:2018ocx ; Carena:2018cjh ; Arcadi:2018tly ; Bauer:2018egk ; Altmannshofer:2016jzy ; Biswas:2016yan ; DiFranzo:2015qea ; Bauer:2018onh ; Ilten:2018crw ; Asai:2018ocx .

In this paper, we revisit in detail the impact of a $L_{\mu}-L_{\tau}$ gauge boson on the particle content and expansion rate of the early universe. We solve the full set of Boltzmann equations that describe the evolution of the $Z^{\prime}$ and neutrino populations, calculating the value of $\Delta N_{\rm eff}$ across the parameter space of this model. We identify two distinct regions in which the tension between early and late-time determinations of the Hubble constant can be substantially relaxed. In the first of these regions, a $\sim$ $10-20$ MeV ${Z^{\prime}}$ reaches equilibrium in the early universe and then decays, heating the neutrino population and delaying the process of neutrino decoupling. For a gauge coupling of $g_{\mu-\tau}\simeq(3-8)\times 10^{-4}$, such a ${Z^{\prime}}$ is also capable of explaining the measured value of the muon’s anomalous magnetic moment. In the second region, a very light and weakly coupled ${Z^{\prime}}$ is produced in the early universe through freeze-in, and decays to neutrinos only after neutrino decoupling. We identify a large plateau of parameter space in which the contribution to the energy density is near $\Delta N_{\rm eff}\simeq 0.2$, in good agreement with the value required to relax the tension associated with the Hubble constant.

The remainder of this paper is organized as follows. In Sec. 2 we describe the $U(1)_{L_{\mu}-L_{\tau}}$ model, including kinetic mixing between the ${Z^{\prime}}$ and the photon. In Sec. 3 we solve the full system of Boltzmann equations in order to calculate the contribution to $N_{\rm eff}$ in this model, in both the thermalized and freeze-in regimes. In Sec. 4 we review the experimental and astrophysical constraints on this model. Finally, in Sec. 5 we summarize our results and conclusions. We also include a series of appendices in which we describe the Boltzmann equations for ${Z^{\prime}}$ production, the bounds from Supernova 1987A, the contribution to ${Z^{\prime}}$-photon kinetic mixing induced by additional scalars, the ${Z^{\prime}}$ mediated corrections to the energy transfer rates responsible for delaying the process of neutrino decoupling, and a detailed description of the solutions to the Boltzmann equation in the freeze-in regime.

2 Model Description

In this study, we extend the SM to include a spontaneously broken $U(1)_{L_{\mu}-L_{\tau}}$ gauge symmetry. In the broken phase222The scalar field whose vacuum expectation value is responsible for the spontaneous breaking of $U(1)_{L_{\mu}-L_{\tau}}$ is generally expected to be much heavier than the $Z^{\prime}$ mass and other energy scales considered in this work. Details regarding this symmetry breaking are beyond the scope of this paper and are not expected to affect our results., the Lagrangian is given by

[TABLE]

where $g_{\mu-\tau}$ is the gauge coupling, $m_{Z^{\prime}}$ is the mass of the gauge boson, and $Z^{\prime}_{\alpha\beta}\equiv\partial_{\alpha}Z^{\prime}_{\beta}-\partial_{\beta}Z^{\prime}_{\alpha}$ is the field strength tensor. The $\mu-\tau$ current is

[TABLE]

where $P_{L}=\frac{1}{2}(1-\gamma_{5})$ is the left chirality projector and

[TABLE]

are the rest frame partial widths for ${Z^{\prime}}$ decays to charged leptons $\ell=\mu,\tau$ and neutrinos $\nu_{i}=\nu_{\mu,\tau}$.

The contribution of the ${Z^{\prime}}$ to the anomalous magnetic moment of the muon, to leading order in $g_{\mu-\tau}$, is given by Pospelov:2008zw

[TABLE]

where $a_{\mu}\equiv\frac{1}{2}(g-2)_{\mu}$. In the last step, we have taken the $m_{{Z^{\prime}}}\ll m_{\mu}$ limit. The measured value of this quantity is $\Delta a_{\mu}\equiv a_{\mu}({\rm obs})-a_{\mu}({\rm SM})=(28.8\pm 8.0)\times 10^{-10}$ Tanabashi:2018oca , and thus requires a gauge coupling of $g_{\mu-\tau}\simeq(3-8)\times 10^{-4}$ in order to resolve the anomaly.

At tree-level, the ${Z^{\prime}}$ in this model couples only to heavy leptons and their neutrino flavors. Muon and tau loops, however, can lead to kinetic mixing between the ${Z^{\prime}}$ and the photon, inducing an effective coupling for the ${Z^{\prime}}$ to electromagnetically charged fermions. For low energy processes we can integrate out $\mu$ and $\tau$, resulting in an off-diagonal kinetic term, $F^{\alpha\beta}Z^{\prime}_{\alpha\beta}$, between the ${Z^{\prime}}$ and the photon. Diagonalizing these fields and restoring canonical normalization induces the following ${Z^{\prime}}$ coupling to the electromagnetic current

[TABLE]

where $e$ is the electron charge, $f$ is a SM fermion with charge $Q_{f}$, and the quantity $\epsilon$ quantifies the degree of kinetic mixing. The irreducible contributions to $\epsilon$ from the loops shown in Fig. 1 are given by Kamada:2015era

[TABLE]

This calculation provides us with a benchmark value for $\epsilon$, which we will refer to throughout this paper as “natural kinetic mixing”. One should keep in mind, however, that other model dependent contributions could potentially arise, in particular if there exist heavy particles that are charged under both electromagnentism and $L_{\mu}-L_{\tau}$.

As a result of this kinetic mixing, the ${Z^{\prime}}$ can decay to $e^{+}e^{-}$ with a partial width given by

[TABLE]

leading to the following branching fraction

[TABLE]

where in the last step we have used the value of $\epsilon$ given in Eq. (6) and taken the $m_{e}\ll m_{Z^{\prime}}\ll m_{\mu}$ limit. Note that this result is independent of the gauge coupling, which cancels in the absence of additional contributions to $\epsilon$.

3 Contributions to $N_{\rm eff}$

In the early universe, the ${Z^{\prime}}$ number density, $n_{{Z^{\prime}}}$, is governed by the following Boltzmann equation

[TABLE]

where $H\equiv\dot{a}/a$ is the expansion rate of the universe and $n^{\rm(eq)}_{{Z^{\prime}}}$ is the equilibrium value of the ${Z^{\prime}}$ number density. The quantity $\Gamma_{{Z^{\prime}}}$ is the rest frame width, for which the thermally averaged value is given by the following

[TABLE]

where $K_{1,2}$ are Bessel functions of the first and second kind and $x\equiv m_{{Z^{\prime}}}/T$. Although many processes can affect $n_{{Z^{\prime}}}$, in the weakly coupled regime ($g_{\mu-\tau}\ll 1$) it suffices to consider only decays and inverse decays in the collision term. For a derivation and other details, see Appendix A.

We are interested in the effect of ${Z^{\prime}}$ decays on the total radiation density just prior to matter-radiation equality at $T_{\gamma}\simeq 0.8\,\eV$, which can be written in terms of $N_{\rm eff}$, the effective number of neutrino species:

[TABLE]

where $\rho_{\gamma}$ is the photon energy density, the factor of $7/8$ accounts for the fact that neutrinos are fermions, and the $(4/11)^{1/3}=T_{\nu}/T_{\gamma}$ in the SM. Note that the SM prediction for $N_{\rm eff}^{\rm SM}=3.045$ deSalas:2016ztq ; Mangano:2005cc is slightly larger than 3 because of the entropy transferred to the neutrinos through $e^{+}e^{-}$ annihilations, the non-instantaneous nature of neutrino decoupling, finite temperature corrections, and neutrino oscillations deSalas:2016ztq ; Mangano:2005cc ; Dolgov:1997mb ; Dicus:1982bz .

The evolution of the ${Z^{\prime}}$ population in the early universe depends on the values of its mass and coupling. Broadly speaking, we will consider two qualitatively distinct regions of parameter space:

• •

Early Universe Equilibrium: If $g_{\mu-\tau}\gtrsim 4\times 10^{-9}$, the $Z^{\prime}$ population thermalizes with the SM bath at early times and decays into neutrinos when $T\sim m_{Z^{\prime}}/3$. If these decays occur predominantly after the neutrinos and photons decouple, they contribute to the neutrino energy density and thereby increase the value of $N_{\rm eff}$. Furthermore, in the presence of non-negligible kinetic mixing with the photon, ${Z^{\prime}}$ interactions with charged particles can delay the neutrino-photon decoupling, quantitatively affecting $N_{\rm eff}$ .

• •

Freeze In (Late Equilibration): If $g_{\mu-\tau}\lesssim 4\times 10^{-9}$, the ${Z^{\prime}}$ population will not have initially been in equilibrium with the SM in the very early universe, but is instead produced through the freeze-in mechanism. For a wide range of masses and couplings, the ${Z^{\prime}}$ production rate is slower than Hubble expansion at very early times, but then becomes comparable as the Hubble rate decreases. Across this broad region of parameter space, the ${Z^{\prime}}$ population eventually thermalizes with neutrinos, but only after the latter decouple from photons, inducing a contribution of $\Delta N_{\rm eff}\simeq 0.21$ through ${Z^{\prime}}\to\mathrm{bar}\nu\nu$ decays, provided that $m_{Z^{\prime}}\gtrsim 1\,\text{eV}$ so that the $Z^{\prime}$ decays prior to CMB formation.

Each of these parameter space regions can be easily identified in Fig. 2.

3.1 Early Universe Equilibrium Regime

If $g_{\mu-\tau}$ is sufficiently large, the rate of ${Z^{\prime}}$ production will exceed that of Hubble expansion, $\langle\Gamma_{Z^{\prime}}\rangle\gg H$, keeping the ${Z^{\prime}}$ population in equilibrium with the SM plasma at early times. In this limit, the ${Z^{\prime}}$ population satisfies $n_{{Z^{\prime}}}=n_{{Z^{\prime}}}^{\rm(eq)}$, where

[TABLE]

is the equilibrium number density and $g_{Z^{\prime}}=3$ is the number of spin states. As the temperature of the universe drops below $m_{{Z^{\prime}}}$, inverse decays become kinetically forbidden and the entropy of this population is transferred to other species. If only neutrinos and photons are present, we can write the effective number of neutrino species as

[TABLE]

where $\rho_{\nu}$ generically differs from its SM value due to entropy transferred from ${Z^{\prime}}$ decays. In the equilibrium regime, these decays take place while ${Z^{\prime}}$ mediated $\nu$ scattering is in chemical equilibrium, so the decay daughter neutrinos thermalize with the existing neutrino population and increase their common temperature.

The presence of a light $U(1)_{L_{\mu}-L_{\tau}}$ $Z^{\prime}$ can substantially alter the decoupling history of the neutrinos through processes of the kind shown in Fig. 3. Of these diagrams, the most important are the decays and inverse decays of the ${Z^{\prime}}$ to neutrinos and electrons. Observations from Borexino constrain the scattering rate between neutrinos and electrons to be within 8% of the SM prediction Harnik:2012ni ; Kamada:2015era , and thus such $2\rightarrow 2$ processes cannot significantly impact the process of neutrino decoupling. The impact of $\nu_{\mu}\mathrm{bar}{\nu}_{\tau}\to Z^{\prime}Z^{\prime}$ scattering is to suppress any chemical potential that might exist in the $\nu_{\mu}-\nu_{\tau}-Z^{\prime}$ sector, an effect that will be efficient for $g_{\mu-\tau}\gtrsim 10^{-5}$ Kamada:2015era . Throughout the phenomenologically viable parameter space of this model (with N_{\text{eff}}\mathrel{\raise 1.29167pt\hbox{<\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}4), however, the chemical potential will be negligible as a result of $e^{+}e^{-}\leftrightarrow\nu_{\mu,\,\tau}\mathrm{bar}{\nu}_{\mu,\,\tau}$ scattering, which is efficient for T\mathrel{\raise 1.29167pt\hbox{>\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}3\,\text{MeV}. Finally, the process $e^{+}e^{-}\leftrightarrow Z^{\prime}\gamma$ is suppressed by a factor of $\alpha\sim 1/137$ relative to $e^{+}e^{-}\leftrightarrow Z^{\prime}$ and thus can be safely neglected.

From these considerations, we conclude that it is sufficient to calculate the evolution of the following three temperatures: $T_{\gamma}(=$ $T_{e}),\,T_{\nu_{e}}$ and $T_{\nu_{\mu}}(=$ $T_{\nu_{\tau}}$ $=$ $T_{Z^{\prime}})$. We adopt the assumption of Maxwell-Boltzmann statistics in the collision terms to derive the following set of coupled differential equations which describe the evolution of these quantities Escudero:2018mvt :

[TABLE]

where the energy transfer rates are

[TABLE]

where $G_{F}=1.16\times 10^{-5}\,\GeV^{-2}$ is the Fermi constant, $s_{W}^{2}=0.223$ is the square of the weak mixing angle Tanabashi:2018oca , $K_{2}$ is a modified Bessel function of the second kind and we have defined the function

[TABLE]

where the first and second terms arise from annihilation and scattering processes respectively; for a derivation of the ${Z^{\prime}}$ induced corrections, see Appendix D. In the above equations, the $\delta\rho_{\nu}$ terms account for SM and ${Z^{\prime}}$ mediated neutrino-electron and neutrino-neutrino interactions, whereas the $\delta\rho_{Z^{\prime}}$ term accounts for $e^{+}e^{-}\leftrightarrow Z^{\prime}$ inverse decays. Here we have safely assumed that the $Z^{\prime}$ population is strongly coupled to the $\mu-\tau$ neutrinos due to the large decay rate, $\Gamma_{Z^{\prime}\to\mathrm{bar}{\nu}_{\mu}\nu_{\mu}}\gg H$. The derivation of SM collision terms here can be found in Hannestad:1995rs ; Kawasaki:2000en ; Dolgov:2002ab ; Escudero:2018mvt and the collision terms for the ${Z^{\prime}}$ mediated processes are derived in Appendix D.

Finally, we note that the mass splittings inferred from atmospheric and solar neutrino observations correspond to oscillations that are active for temperatures $T$ $\sim 3$ MeV and $T$ $\sim 5$ MeV, respectively Hannestad:2001iy ; Dolgov:2002ab ; Dolgov:2002wy . Since neutrino oscillations imply a rapid change in the neutrino flavor, they act as to equilibrate the $\nu_{e}$ and $\nu_{\mu}-\nu_{\tau}$ distribution functions, leading to $T_{\nu_{e}}\simeq T_{\nu_{\mu}}$. To properly account for such oscillations would require the use of the density matrix formalism. Since the effect of such oscillations is to equilibrate the temperatures of the different neutrino species, however, we can simply work in terms of a common neutrino temperature, $T_{\nu}$. In this case, the evolution of the temperatures $T_{\nu}$ and $T_{\gamma}$ can be simplified as follows

[TABLE]

where $\rho_{\nu}\equiv\rho_{\nu_{e}}+2\rho_{\nu_{\mu}}$ and the energy transfer rates are $\delta\rho_{\nu}/\delta t=\delta\rho_{\nu_{e}}/\delta t+2\delta\rho_{\nu_{\mu}}/\delta t$ and $\delta\rho_{Z^{\prime}}/\delta t$, as described in Eqns. (17), (18) and (19).

Using this set of equations, we can describe the process of neutrino decoupling, including the energy injection into the SM plasma from ${Z^{\prime}}$ decays and any exchange of energy between neutrinos and the electromagnetic sector that might result from $U(1)_{L_{\mu}-L_{\tau}}$ interactions. We solve this set of equations starting from an initial condition of $T_{\gamma}=T_{\nu}=20\,\text{MeV}$, for which the SM weak interactions are very efficient ($\Gamma/H\sim 300$). In Fig. 4 we plot the evolution of the photon-to-neutrino temperature ratio for several choices of the ${Z^{\prime}}$ mass and coupling. Note that when solving these equations we explicitly check that the continuity equation, $\dot{\rho}=-3H(\rho+p)$, is fulfilled to a relative accuracy of at least $10^{-5}$. In Fig. 5 we show contours for the value of $N_{\rm eff}$ predicted across a range of parameter space in this model, focusing on the range of masses and couplings that are potentially relevant for the $(g-2)_{\mu}$ anomaly. We show these results for three choices of the kinetic mixing parameter, and both including or ignoring the impact of neutrino oscillations. Also shown are the constraints on this parameter space derived from measurements of di-muon production in neutrino nucleus scattering (CCFR) Mishra:1991bv ; Altmannshofer:2014pba , the solar neutrino scattering rate (Borexino) Harnik:2012ni ; Laha:2013xua ; Kamada:2015era (see Sec. 4). For $m_{{Z^{\prime}}}\sim 10-20$ MeV, the predicted value of $N_{\rm eff}$ falls within the range capable of relaxing the tension between late and early-time determinations of the Hubble constant. Furthermore, for $g_{\mu-\tau}\simeq(3-8)\times 10^{-4}$, such a particle can also account for the observed value of $(g-2)_{\mu}$. In Fig. 6, we extend these results to smaller values of $g_{\mu-\tau}$. Although this parameter space cannot address the measured magnetic moment of the muon, the Hubble tension can be relaxed for a ${Z^{\prime}}$ with a wide range of couplings.

Before moving on, we note that the results described in this section do not agree with those presented in Ref. Kamada:2015era . This is the case for two primary reasons. First, the authors of Ref. Kamada:2015era did not consider kinetic mixing between the photon and the ${Z^{\prime}}$ and thus neglect potentially important process of the form $e^{+}e^{-}\leftrightarrow Z^{\prime}$. Second, they adopted 1.5 MeV for the decoupling temperature of the $\nu_{\mu,\tau}$ neutrinos, instead of the correct value of $T_{\nu_{\mu}}^{\text{dec}}\simeq 3.2\,\text{MeV}$ (or $T\simeq 2.4\,\text{MeV}$ after taking into account oscillations) Escudero:2018mvt ; Hannestad:2001iy ; Dolgov:2002ab ; Dolgov:2002wy . If we set $\epsilon=0$ and $T^{\rm dec}_{\nu_{\mu},\nu_{\tau}}=1.5$ MeV, we are able to reproduce the results presented in Ref. Kamada:2015era . In a follow up study, Ref. Kamada:2018zxi considered the interactions $e^{+}e^{-}\leftrightarrow Z^{\prime}$ in the particular case in which $g_{\mu-\tau}=5\times 10^{-4}$ and used entropy conservation to track the temperature evolution of the different sectors. Our results are similar to those presented in Ref. Kamada:2018zxi for $g_{\mu-\tau}=5\times 10^{-4}$. Ref. Kamada:2018zxi , however, did not include SM neutrino-electron interactions and assumed $T_{\nu_{\mu}}^{\text{dec}}=1.5\,\text{MeV}$, whereas the present work includes these interactions and calculates the value of $T_{\nu}^{\rm dec}$.

3.2 Freeze-In (Late Equilibration) Regime

If the value of $g_{\mu-\tau}$ is very small, the ${Z^{\prime}}$ population will not reach equilibrium with the SM thermal bath in the very early universe. The condition for these particles to reach equilibrium with muons and taus at early times can be written as $n_{\mu}\langle\sigma v\rangle_{\mu\mu\to\gamma{Z^{\prime}}}>H(m_{\mu})$, corresponding to the following requirement on the gauge coupling:

[TABLE]

where we have approximated $\langle\sigma v\rangle_{\mu\mu\to\gamma{Z^{\prime}}}\sim\alpha g_{\mu-\tau}^{2}/2m_{\mu}^{2}$ in the nonrelativistic limit near $T\sim m_{\mu}$, where the production rate is maximized relative to that of Hubble expansion. As long as $m_{Z^{\prime}}\ll m_{\mu}$, the condition in Eq. (23) is independent of $m_{Z^{\prime}}$.

But even if $g_{\mu-\tau}$ is too small for the ${Z^{\prime}}$ population to be thermalized through $\mu^{+}\mu^{-}\to Z^{\prime}\gamma$ annihilations, a significant abundance of $Z^{\prime}$’s can be produced through inverse decays. In fact, inverse decays may ultimately be able to generate an equilibrium ${Z^{\prime}}$ abundance. More specifically, in the case of a light ${Z^{\prime}}$ (for which inverse decays are active when $T_{\nu}\gg m_{Z^{\prime}}$), equilibrium between the ${Z^{\prime}}$ and neutrino populations will be achieved for g_{\mu-\tau}\mathrel{\raise 1.29167pt\hbox{>\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}1.3\times 10^{-10}(m_{{Z^{\prime}}}/{\rm MeV})^{1/2}.

After neutrino decoupling, neutrino oscillations are efficient enough to approximately equilibrate the neutrino distributions Dolgov:2002ab ; Hannestad:2001iy ; Dolgov:2002wy , and thus even if the $Z^{\prime}$ only interacts with muon and tau neutrinos the rapid oscillation will render the $\nu_{e}$ distribution very similar to that of $\nu_{\mu}$ and $\nu_{\tau}$. Taking this into account, and imposing energy and number density conservation, we can calculate the neutrino and $Z^{\prime}$ distributions after equilibration. These equations explicitly read

[TABLE]

where the equilibrium condition requires $\mu_{Z^{\prime}}=2\mu_{\nu}$. Solving these equations simultaneously yields the following result:

[TABLE]

Thus the neutrino temperature increases at the same time that a chemical potential is developed for the neutrino and $Z^{\prime}$ populations. Once the $Z^{\prime}$ population is in equilibrium with the neutrinos, it will maintain its distribution due to the active decays and inverse decays until the temperature drops to $T_{\nu}\sim m_{Z^{\prime}}/3$, at which time they will decay out of equilibrium and thereby generate a net contribution to $N_{\text{eff}}$. Since this decay process is occurring out of equilibrium, we cannot use Fermi-Dirac or Bose-Einstein distributions but instead have to solve the full Boltzmann equation for the distribution functions. The relevant Boltzmann equations read as follows Kawasaki:1992kg ; Dolgov:1998st :

[TABLE]

where

[TABLE]

By using the continuity equation we find that the photon temperature simply evolves as

[TABLE]

In order to solve Eqns. (27) and (28), we bin the $Z^{\prime}$ and $\nu$ distribution functions in comoving momentum, $y=a\,p/\text{MeV}$, from $y_{\text{min}}=0.01$ to $y_{\text{max}}=20$ in 100 bins (this choice has been previously shown to produce accurate results Dolgov:1997mb ; Dolgov:1998st ). We use this set of equations for both the equilibration and freeze-in cases, since out-of-equilibrium decays result in each. For the freeze-in case, we start the integration at a temperature of $T_{\nu}=100\times m_{Z^{\prime}}$ with the initial condition $f_{Z^{\prime}}=0$ and $f_{\nu}=f^{\rm FD}(T_{\nu},\mu_{\nu}=0)$, whereas in the equilibrium case we start the integration at a temperature of $T_{\nu}=10\times m_{Z^{\prime}}$ and with the initial condition $f_{Z^{\prime}}=f^{\rm BE}(T_{\rm eq},2\,\mu_{\nu})$ and $f_{\nu}=f^{\rm FD}(T_{\rm eq},\mu_{\nu})$.

In Fig. 7 we show the evolution of the neutrino and ${Z^{\prime}}$ energy densities in two representative cases. In the case of $g_{\mu-\tau}=10^{-10}$ and $m_{Z^{\prime}}=1\,\text{keV}$, the $Z^{\prime}$ population thermalizes with the neutrinos while relativistic and only decays after $T_{\nu}\sim m_{Z^{\prime}}/3$, rendering $\Delta N_{\text{eff}}\simeq 0.21$. By contrast, for the case of $g_{\mu-\tau}=1.3\times 10^{-11}$ and $m_{Z^{\prime}}=2.2\,\text{keV}$, the $Z^{\prime}$ population reaches thermal equilibrium with the neutrinos at $T_{\nu}\sim m_{Z^{\prime}}$. We see in this case that the neutrino energy density is reduced as the $Z^{\prime}$ population grows, and then increases as the ${Z^{\prime}}$ population decays, to yield $\Delta N_{\text{eff}}\simeq 0.21$333If the effect of neutrino oscillations is neglected in this computation, we find $\Delta N_{\text{eff}}\simeq 0.18$ instead.. In Appendix E we display the time evolution of $f_{Z^{\prime}}$ and the final $f_{\nu}|_{T_{\gamma}\ll m_{Z^{\prime}}}$ for some representative values of $g_{\mu-\tau}$ and $m_{Z^{\prime}}$.

In exploring the freeze-in parameter space within this model, we find the following behavior:

[TABLE]

We would like to emphasize that throughout the parameter space in which the ${Z^{\prime}}$ population reaches thermal equilibrium with the neutrinos, this model yields a prediction of $\Delta N_{\text{eff}}\simeq 0.21$. This limit is realized across a wide range of the $m_{{Z^{\prime}}}-g_{\mu-\tau}$ parameter space (see Fig. 2), making $\Delta N_{\text{eff}}\simeq 0.21$ a common feature within the context of this model.

We appreciate that this result will strike many readers as counterintuitive. To understand why this $\Delta N_{\rm eff}\simeq 0.21$ plateau arises within the parameter space of this model, consider the following summary. If g_{\mu-\tau}\mathrel{\raise 1.29167pt\hbox{<\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}4\times 10^{-9}, the ${Z^{\prime}}$ population will not thermalize with muons and taus at early times (see Eq. 23). But so long as g_{\mu-\tau}\mathrel{\raise 1.29167pt\hbox{>\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}1.3\times 10^{-10}\,(m_{{Z^{\prime}}}/{\rm MeV})^{1/2}, a light ${Z^{\prime}}$ will ultimately obtain an equilibrium abundance through inverse decays. Thus for the a sizable range of couplings that fall in between these two limits, the ${Z^{\prime}}$ population will not reach equilibrium with muons or taus at early times, but will equilibrate with neutrinos, ultimately leading to $\Delta N_{\rm eff}\simeq 0.21$ through their decays. Note that $\Delta N_{\rm eff}\simeq 0.21$ will only result from $Z^{\prime}$’s with $m_{Z^{\prime}}\gtrsim 1\,\text{eV}$, since for lighter gauge bosons the $Z^{\prime}$ population would not have decayed prior to $T_{\rm CMB}\simeq 0.26\,\text{eV}$.

4 Constraints

In this section, we consider a number of constraints on this model, including those that apply directly to the gauge coupling, and which depend on the degree of kinetic mixing between the ${Z^{\prime}}$ and the photon.444We have explicitly checked that given the small branching fraction of the $Z^{\prime}$ to electrons ${\rm Br}_{Z^{\prime}\to e^{+}e^{-}}\simeq 2\times 10^{-5}$ (8), the energy injection from the $Z^{\prime}$ population decay during the CMB or BBN Fradette:2014sza is too small to result in meaningful constraints on the model. See e.g. Bauer:2018onh ; Ilten:2018crw for an exhaustive set of constraints for MeV-GeV $Z^{\prime}$’s.

4.1 Neutrino Tridents

Light $L_{\mu}-L_{\tau}$ gauge bosons are constrained by the results of the Columbia-Chicago-Fermilab-Rochester (CCFR) Neutrino Experiment, which are consistent with the SM predictions for $\nu N\to\nu N\mu^{+}\mu^{-}$ trident production Mishra:1991bv ; Altmannshofer:2014pba . In the presence of a light ${Z^{\prime}}$, this process receives additional corrections from the three body $\nu N\to\nu N{Z^{\prime}}$ process followed by a prompt ${Z^{\prime}}\to\mu^{+}\mu^{-}$ decay. Consistency with the CCFR trident measurement requires {g_{\mu-\tau}}\mathrel{\raise 1.29167pt\hbox{<\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}10^{-3} for $m_{Z^{\prime}}\ll\MeV$, as shown in Figs. 2 and 5.

4.2 Supernova 1987A

A light ${Z^{\prime}}$ can be produced via neutrino inverse decays $\mathrm{bar}\nu\nu\to{Z^{\prime}}$ within the dense core of supernovae. If such particles are long-lived, they may efficiently remove energy from the core and thus modify the observed $\sim$ 10 second duration of the neutrino burst from Supernova 1987A Kolb:1987qy . In Figs. 2 we show the regions of parameter space in which the ${Z^{\prime}}$ luminosity introduces an order one correction to the total energy loss ($\sim$ $10^{53}$ erg) in the first $\sim$ 10 seconds of a supernova explosion, assuming a $\sim$ 10 km radius and a $\sim$ 30 MeV core temperature. The details of this estimate are described in Appendix B.

4.3 Stellar Cooling

In this model, the ${Z^{\prime}}$ can be produced in stellar plasmas and induce anomalous cooling by carrying energy away from stars. Since this process yields an invisibly decaying particle, the stellar bounds on dark photons from Refs. An:2013yfc ; Hardy:2016kme can straightforwardly be adapted to constrain the $L_{\mu}-L_{\tau}$ model. This constraint is shown in the lower frame of Fig. 2. We also note that bounds from white dwarfs cooling apply to our model Bauer:2018onh ; Dreiner:2013tja , and that they are of comparable strength to those derived from Borexino for $Z^{\prime}$’s in the MeV mass range Kamada:2018zxi . However, such constraints have been derived assuming that the $Z^{\prime}$ is more massive than the typical temperature of a white dwarf ($T\sim 5\,\text{keV}$), and therefore it is not trivial to extrapolate them for $m_{Z^{\prime}}\lesssim 10\,{\rm keV}$, and hence are not shown in the figures.

4.4 Solar Neutrino Scattering

Nonzero kinetic mixing introduces ${Z^{\prime}}$ mediated interactions between $\nu_{\mu,\tau}$ and charged fermions, which can distort the observed solar neutrino scattering rate at Borexino Harnik:2012ni ; Laha:2013xua ; Kamada:2015era . These constraints are shown in Figs. 2 and 5 (in those frames corresponding to $\epsilon>0)$.

5 Summary and Conclusions

The reported tension between early and late time determinations of the Hubble constant has motivated us to consider models that include a new light gauge boson. In particular, in light of the stringent constraints on gauge bosons with couplings to first generation fermions, we consider a model with a broken $U(1)_{L_{\mu}-L_{\tau}}$ gauge symmetry, corresponding to a massive gauge boson with tree-level couplings to muons, taus, and their corresponding neutrino species. We have studied the impact of such a particle on the evolution of the early universe, solving the full set of Boltzmann equations and determining the contribution to the energy density in radiation, parameterized in terms of $\Delta N_{\rm eff}$.

We have identified two distinct regions of parameter space in which the tension related to the Hubble constant can be substantially relaxed (corresponding to $\Delta N_{\rm eff}\sim 0.2-0.5$):

• •

For a gauge boson with g_{\mu-\tau}\mathrel{\raise 1.29167pt\hbox{>\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}4\times 10^{-9}, this particle will reach equilibrium with the Standard Model bath in the early universe. For $m_{{Z^{\prime}}}\sim$ $10-20$ MeV, the decays of these particles will heat the neutrino population and delay the process of neutrino decoupling, resulting in $\Delta N_{\rm eff}\sim 0.2-0.5$. For a coupling of $g_{\mu-\tau}\simeq(3-8)\times 10^{-4}$, such a ${Z^{\prime}}$ could also account for the measured value of the muon’s anomalous magnetic moment. We also note that the parameter space which simultaneously ameliorates both the $H_{0}$ and $(g-2)_{\mu}$ anomalies can also be tested with existing and future accelerator searches for muon specific forces with muon beams Gninenko:2018tlp ; Kahn:2018cqs ; Berlin:2018bsc ; Chen:2018vkr , measurements of rare kaon decays Ibe:2016dir ; diego , and at Belle-II Araki:2017wyg ; Banerjee:2018mnw .

• •

For a very light ($m_{Z^{\prime}}\gtrsim 1\,\text{eV}$) gauge boson with a very small coupling (g_{\mu-\tau}\mathrel{\raise 1.29167pt\hbox{<\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}4\times 10^{-9}), the ${Z^{\prime}}$ population does not reach equilibrium with muons or taus at early times. Instead, a significant abundance of such particles can be produced through $\mathrm{bar}\nu\nu\to{Z^{\prime}}$ inverse decays. In particular, for g_{\mu-\tau}\mathrel{\raise 1.29167pt\hbox{>\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}1.3\times 10^{-10}(m_{{Z^{\prime}}}/{\rm MeV})^{1/2}, inverse decays will ultimately produce a ${Z^{\prime}}$ population that reaches equilibrium with neutrinos. When this ${Z^{\prime}}$ population decays, it produces an energy density of neutrinos corresponding to $\Delta N_{\rm eff}\simeq 0.21$. This result is found across a wide range of parameter space within this model.

These results are summarized in Fig. 2, which includes contours of constant $\Delta N_{\text{eff}}$, as well as the regions favored by measurements of the muon’s anomalous magnetic moment, $(g-2)_{\mu}$. Among other motivations, this model is particularly interesting in light of the large regions of its parameter space that can relax the reported Hubble tension ($\Delta N_{\rm eff}\sim 0.2-0.5$) and that are within the reach of next-generation CMB measurements (\Delta N_{\rm eff}\mathrel{\raise 1.29167pt\hbox{>\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}0.03 Abazajian:2016yjj ).

Acknowledgments

We thank Asher Berlin, Sam McDermott for helpful conversations. Fermilab is operated by Fermi Research Alliance, LLC, under Contract No. DE- AC02-07CH11359 with the US Department of Energy. This project has received funding/support from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreements Elusives ITN No. 674896 and InvisiblesPlus RISE No. 690575. The work of MP was supported by the Spanish Agencia Estatal de Investigación through the grants FPA2015-65929-P (MINECO/FEDER, UE), IFT Centro de Excelencia Severo Ochoa SEV-2016-0597, and Red Consolider MultiDark FPA2017-90566-REDC. ME is supported by the European Research Council under the European Union’s Horizon 2020 program (ERC Grant Agreement No 648680 DARKHORIZONS). ME and MP acknowledge the Fermilab Theory and Astro groups for their hospitality when this project was initiated.

Appendix A Boltzmann Equation for ${Z^{\prime}}$ Production

The number density evolution of the $Z^{\prime}$ in the presence of the $\nu(p_{1})+\nu(p_{2})\longleftrightarrow Z^{\prime}(p_{3})$ process is

[TABLE]

where $n_{{Z^{\prime}}}$ is the ${Z^{\prime}}$ number density, $H$ is the Hubble rate, and the collision term satisfies

[TABLE]

where $|{\cal A}_{i\to f}|^{2}$ is the squared, spin-averaged matrix element for the process in the subscript. By detailed balance, we can simplify the inverse process using

[TABLE]

so the collision terms can be simplified

[TABLE]

Using the definition of the energy dependent width

[TABLE]

where the prefactor is $1/E_{Z^{\prime}}$ instead of $1/m_{{Z^{\prime}}}$ (as in the rest frame expression), the Boltzmann equation becomes

[TABLE]

where we have defined the comoving dimensionless yield $Y_{{Z^{\prime}}}\equiv n_{{Z^{\prime}}}/s$ and

[TABLE]

which is valid in the limit of Boltzmann statistics. Note that Eq. (A.6) has the familiar freeze out form, except that there is a linear dependence on the phase space distributions (as opposed to quadratic) and the thermally averaged cross section for annihilation is replaced with a thermally averaged decay rate. Defining $z\equiv m_{{Z^{\prime}}}/T$, and $Y_{Z^{\prime}}=n_{Z^{\prime}}/s$, where $s$ is the total entropy density, we have

[TABLE]

and we can rewrite the full Boltzmann equation as

[TABLE]

which is equivalent to Eq. (9) expressed in terms of the yield, $Y_{Z^{\prime}}$. Note that this equation accounts for both production and decay processes.

Appendix B Supernova 1987A Bounds

The energy-density loss rate per unit volume from $\nu(p_{1})+\nu(p_{2})\to Z^{\prime}(p_{{Z^{\prime}}})$ production in an isotropic, homogeneous supernova is

[TABLE]

Since the neutrinos are in equilibrium in the core, but the ${Z^{\prime}}$ population is not, we neglect the reverse process. By detailed balance, we can substitute $f^{\rm(eq)}_{\nu}(\vec{p}_{1})f^{\rm(eq)}_{\nu}(\vec{p}_{2})=f^{\rm(eq)}_{{Z^{\prime}}}(\vec{p}_{3})$, which yields

[TABLE]

Since we are mainly interested in very small ${Z^{\prime}}$ couplings in the non-equilibrium regime, we can neglect the reverse process which depends on $f_{{Z^{\prime}}}$, so using the definition of the rest frame width

[TABLE]

the total production rate can be written

[TABLE]

where the integral corresponds to the equilibrium number density of ${Z^{\prime}}$s inside the supernova. To calculate the total energy carried away by the ${Z^{\prime}}$ population, we need to include the probability that a given ${Z^{\prime}}$ will survive out to a distance of $R\sim 10$ km at $T\sim 30$ MeV, so we have

[TABLE]

where $\gamma=E_{{Z^{\prime}}}/m_{{Z^{\prime}}}$ is the boost factor and the factor of 3 accounts for the number of polarization states in the convention $f^{\rm(eq)}=3\,e^{-E_{{Z^{\prime}}}/T}$ since we performed a spin average in the definition of the rest frame width. To ensure that we don’t modify the observation of Supernova 1987A by an order one amount, we demand that the total energy loss in a 10 second interval not exceed $\sim$ $10^{53}$ erg, so our criterion is

[TABLE]

where we have dropped ${Z^{\prime}}$ subscripts and multiplied by the core volume, $V=\frac{4}{3}\pi R^{3}$.

Comparing Other Processes

The ${Z^{\prime}}$ can also be produced via $\nu$ scattering with radiative emission $\nu\nu\to\nu\nu{Z^{\prime}}$ which has the same scaling in the gauge coupling as inverse decays, but doesn’t depend on $m_{Z^{\prime}}$ for $m_{{Z^{\prime}}}\ll T\sim 30$ MeV. The radiative production rate is approximately

[TABLE]

where we have taken $T\sim 30$ MeV, included a $\alpha_{\mu-\tau}$ penalty for ${Z^{\prime}}$ emission, and conservatively neglected the 3-body phase suppression for this process. Comparing this to the thermally averaged inverse decay process with the same parametric scaling,

[TABLE]

we find that the latter dominates over the full mass range we consider. Only near $m_{Z^{\prime}}\sim$ eV does radiative emission compete. Similarly, we can estimate neutrino annihilation, $\nu\nu\to{Z^{\prime}}{Z^{\prime}}$, for which the rate scales as $\sim\alpha_{\mu-\tau}^{2}T$ and yields

[TABLE]

so unless the mass is very small or the gauge coupling is large and in the thermal regime, inverse decays are the dominant form of ${Z^{\prime}}$ production. Note that in this estimate, we have taken $T=30$ MeV for the SN temperature.

Appendix C Scalar Loop Induced Kinetic Mixing

In addition to the irreducible “natural” kinetic mixing considered throughout this paper, there could also be additional contributions from hypothetical heavier particles, beyond the SM. To assess the magnitude of such contributions, we calculate in this appendix the contribution to $\gamma-{Z^{\prime}}$ kinetic mixing from an additional pair of heavy scalars with identical charges under electromagnetism and opposite charges under $L_{\mu}-L_{\tau}$.

The loop-level kinetic mixing from each scalar arises from two intermediate loop diagrams, as shown in Fig. 8. The two propagator contribution is

[TABLE]

where $m$ is a scalar mass and we combine denominators with Feynman parameters

[TABLE]

where $\ell\equiv k+xq$ and $\Delta\equiv-x(1-x)q^{2}+m^{2}$. Dropping terms linear in $\ell$ in the numerator of the integrand, we can rewrite this as

[TABLE]

which we can evaluate in dimensional regularization and extract the divergent piece

[TABLE]

where the second term has the Ward identity structure, but the first term will need to cancel against the contribution from the one propagator diagram with a four point vertex, which is simpler to evaluate

[TABLE]

which yields a divergent piece

[TABLE]

Combining the divergent $\epsilon$ pole contributions now becomes

[TABLE]

which has the requisite gauge invariant structure. In the $\overline{\rm MS}$ scheme the residual finite contribution is

[TABLE]

so performing the integrals yields the gauge invariant result

[TABLE]

To eliminate the explicit dimensionful $\log m^{2}$ dependence, we need to compare this to a reference value, so just like with the fermion case with $\mu$ and $\tau$ loops, we add a second scalar with opposite $L_{\mu}-L_{\tau}$ charge and identical SM quantum numbers to obtain a mixing with the hypercharge gauge boson

[TABLE]

where $m_{1,2}$ are the scalar masses, $x_{i}\equiv\sqrt{4m_{i}^{2}/q^{2}-1}$. We can extract the kinetic mixing from the coefficient of $\left(q^{2}g^{\mu\nu}-q^{\mu}q^{\nu}\right)$ by expanding in powers of $q^{2}$ to get

[TABLE]

where we have taken the limit $q^{2}\to 0$ and assumed $m_{1}/m_{2}=m_{\tau}/m_{\mu}$ in the log. If the scalars are integrated out above the electroweak scale, the only allowed kinetic mixing is between the ${Z^{\prime}}$ and the SM hypercharge gauge boson, $B_{\mu}\equiv\cos\theta_{W}A_{\mu}-\sin\theta_{W}Z^{0}_{\mu}$, so we would need to replace $e\to g^{\prime}$ where $g^{\prime}\equiv e\cos\theta_{W}$ is the hypercharge gauge coupling. Since the induced kinetic mixing is logarithmically sensitive to heavy particle masses, there is considerable UV dependence to the “natural” value.

Appendix D Delayed Neutrino Decoupling

In this appendix we derive the ${Z^{\prime}}$ mediated corrections to the energy transfer rates in Eqs. (17) and (19), which are responsible for delaying neutrino decoupling; the SM processes which are also present and available in Escudero:2018mvt and in references therein. The Liouville equation that evolves the phase space density $f_{i}$ of particle species $i$ in an expanding FRW universe is

[TABLE]

where $p_{i}$ is the particle’s three momentum and the collision term for $i+X\to Y$ is

[TABLE]

where $X$ and $Y$ are other initial and final state species, ${\cal A}_{i+X\to Y}$ is the matrix element for this process, $S$ is a symmetrization factor which is $1/2$ for identical particles in the initial or final states, and

[TABLE]

is the phase space factor where the $\pm$ signs correspond to Bose and Pauli statistics, respectively.

Dividing Eq. (D.1) by $(2\pi)^{3}$ and integrating over $g_{i}\,d^{3}\vec{p}_{i}$ gives the familiar Boltzmann equation for the evolution of the number density

[TABLE]

Similarly, dividing Eq. (D.1) by $(2\pi)^{3}$ and integrating over $g_{i}\,E_{i}\,d^{3}\vec{p}_{i}$ we can track the evolution of the energy density

[TABLE]

where $P_{i}$ is the pressure and we have defined an energy transfer rate $\delta\rho_{i}/\delta t$. Thus, in the context of this work, computing the neutrino decoupling temperature is tantamount to evaluating integrals over the collision term.

Inverse Decays

For the inverse decay process $e^{+}(p_{1})e^{-}(p_{2})\to{Z^{\prime}}(k)$ in Eq. (17) the energy transfer rate is

[TABLE]

where $|{\cal A}_{e^{+}e^{-}\to{Z^{\prime}}}|^{2}$ is the inverse decay matrix element, $f^{\rm(eq)}(p)$ is approximated to be the Maxwell-Boltzmann distribution, and we have neglected Pauli blocking effects. Using energy conservation, detailed balance $f^{\rm(eq)}_{e}(\vec{p}_{1})f^{\rm(eq)}_{e}(\vec{p}_{2})=f^{\rm(eq)}_{{Z^{\prime}}}(\vec{k})$, and the definition of the partial width for $\Gamma_{{Z^{\prime}}\to e^{+}e^{-}}$ we can simplify this expression

[TABLE]

where we have dropped the $k$ subscript on the integration variable. Note that this result recovers the first term on the RHS of Eq. (17); the reverse process with $T_{\gamma}\to T_{\nu_{\mu}}$ accounts for the second term, so we have

[TABLE]

which recovers the expression in Eq. (17).

Annihilation and Scattering

For the 2-2 processes in Eq. (19), we follow the formalism outlined in Refs. Kawasaki:2000en ; Dolgov:2002wy ; Escudero:2018mvt , which compute the collision terms for $\mathrm{bar}\nu\nu\leftrightarrow e^{+}e^{-}$ annihilation and $\nu e\leftrightarrow\nu e$ scattering in the SM. We are interested in using these results to compute the ${Z^{\prime}}$ mediated corrections to neutrino decoupling in our $L_{\mu}-L_{\tau}$ extension with ${Z^{\prime}}-\gamma$ kinetic mixing. To simplify our analysis, we take the massless electron limit and integrate out the ${Z^{\prime}}$ to obtain an effective Lagrangian for $\nu_{\mu,\tau}$ and electron interactions, which contains the interactions

[TABLE]

where the first term represents SM interactions with vector and axial couplings $C_{V,A}$ and the second term arises from virtual ${Z^{\prime}}$ exchange. From Kawasaki:2000en ; Dolgov:2002ab , the SM energy transfer rates for scattering and annihilation can be written as Escudero:2018mvt

[TABLE]

where we have neglected all chemical potentials and included all forwards and backwards reactions. We can use this result to compute the corresponding ${Z^{\prime}}$ induced corrections by noting that the new interactions in Eq. (D.9) have the same form as their SM electroweak counterparts with the replacements

[TABLE]

so the additional energy transfer rates in an $L_{\mu}-L_{\tau}$ extension are

[TABLE]

which recovers the $g_{\mu-\nu}$ dependent correction in Eq. (19). Note that this expression includes an additional overall factor of 2 to account for both neutrinos and antineutrinos.

Appendix E Freeze-in Solutions

In this Appendix, we describe some of the solutions to the time evolution of the $Z^{\prime}$ and neutrino distribution functions as obtained by integrating Eqns. (27) and (28) as described in Sec. 3. In the upper panels of Fig. 9 we show the evolution of $f_{Z^{\prime}}$ as a function of the photon temperature for the two scenarios considered in Fig. 7. In the upper left panel we show the case in which the $Z^{\prime}$ reaches thermal equilibrium with the neutrino population while relativistic, while the upper right panel corresponds to a scenario in which the $Z^{\prime}$ thermalizes when $T_{\nu}\sim m_{Z^{\prime}}$. The two lower panels show the evolution of the $Z^{\prime}$ distribution function for two choices of $g_{\mu-\tau}$ and $m_{Z^{\prime}}$ for which the $Z^{\prime}$ never reaches thermal equilibrium with the neutrinos, leading to $\Delta N_{\rm eff}<0.21$.

In Fig. 10 we show the neutrino distribution function after a ${Z^{\prime}}$ population that was initially generated through $\mathrm{bar}{\nu}\nu\to Z^{\prime}$ inverse decays has complete decayed away (at $T_{\nu}\ll m_{Z^{\prime}}$, as relevant for CMB observations). The $y$ axis has been normalized such that $\Delta N_{\rm eff}$ can be computed as $\Delta N_{\rm eff}=3\times\frac{120}{7\pi^{4}}\int_{0}^{\infty}\,dy\,y^{3}(f_{\nu}-f_{\nu}^{\rm FD})$, where $f_{\nu}^{\rm FD}$ corresponds to the distribution function of a free-streaming decoupled neutrino, $f_{\nu}^{\rm FD}=1/(1+e^{y_{\nu}})$. Note that the results for a different $Z^{\prime}$ mass can be obtained by rescaling the corresponding couplings in Fig. 10 by $(2.2\,\text{keV}/m_{Z^{\prime}})^{-1/2}$.

Fig. 10 illustrates why a $Z^{\prime}$ that thermalizes with the neutrinos after they have decoupled from the SM plasma renders $\Delta N_{\rm eff}\simeq 0.21$. Once the neutrinos have decoupled from the SM plasma, the only relevant process are $\mathrm{bar}{\nu}\nu\to Z^{\prime}$ and $Z^{\prime}\to\mathrm{bar}{\nu}\nu$. At very high temperatures $f_{Z^{\prime}}=0$, and a $Z^{\prime}$ population is eventually generated at the expense of neutrinos, via $\mathrm{bar}{\nu}\nu\to Z^{\prime}$ (this can also be seen from Fig. 7). After the $Z^{\prime}$ population has thermalized with the neutrinos, it decays out of equilibrium to neutrinos. However, the decay products of the $Z^{\prime}$ population have an energy $E_{\nu}\sim m_{Z^{\prime}}/2$, which is substantially different from $3\,T_{\nu}$, leading to a final neutrino population that is more energetic than that found in thermal equilibrium (as can be appreciated from Fig. 10). This results in $\Delta N_{\rm eff}\simeq 0.21$ for g_{\mu-\tau}\mathrel{\raise 1.29167pt\hbox{>\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}1.3\times 10^{-10}\,(m_{Z^{\prime}}/\text{MeV})^{1/2}, and to $\Delta N_{\rm eff}<0.21$ for smaller values of $g_{\mu-\tau}$ (which do not enable the $Z^{\prime}$ population to reach thermal equilibrium with the neutrinos).