Revisiting Neutrino Self-Interaction Constraints from $Z$ and $\tau$ decays

Given the elusive nature of neutrinos, their self-interaction is particularly difficult to probe. Nevertheless, upper limits on the strength of such an interaction can be set by using data from terrestrial experiments. In this work we focus on additional contributions to the invisible decay width of $Z$ boson as well as the leptonic $\tau$ decay width in the presence of a neutrino coupling to a relatively light scalar. For invisible $Z$ decays we derive a complete set of constraints by considering both three-body bremsstrahlung as well as the loop correction to two-body decays. While the latter is usually regarded to give rather weak limits we find that through the interference with the Standard Model diagram it actually yields a competitive constraint. As far as leptonic decays of $\tau$ are concerned, we derive a first limit on neutrino self-interactions that is valid across the whole mass range of a light scalar mediator. Our bounds on the neutrino self-interaction are leading for $m_\phi \gtrsim 300$ MeV and interactions that prefer $\nu_\tau$. Bounds on such $\nu$-philic scalar are particularly relevant in light of the recently proposed alleviation of the Hubble tension in the presence of such couplings.


I. INTRODUCTION
Recent studies have revealed a discrepancy between local measurements of the Hubble constant [1][2][3] and those obtained by analyzing the Cosmic Microwave Background (CMB) data [4] at a 4σ level. This has sparked an ongoing controversy in cosmology and the search for potential solutions is currently ongoing. At the moment the origin of the Hubble tension is unclear; potential solutions include, for example, early dark energy [5], dark matter neutrino interactions [6], certain classes of non-Gaussian primordial fluctuations [7] or, more prosaically, underestimated systematics [8]. Most of these ideas fall clearly into the realm of cosmology and astrophysics and cannot be tested in laboratory experiments. However, it was proposed recently that strong neutrino selfinteractions (νSI) can alleviate this tension [9]. The preferred value of the interaction strength is in the ballpark of 10 7 ∼ 10 9 in units of Standard Model (SM) weak interaction strength G F . In this regime neutrino free-streaming is suppressed at high red-shift and it is not surprising that such an interaction can have remarkable consequences for the physics of the early Universe.
Large νSI present a challenge from a particle physics perspective and it is expected that terrestrial experiments can help scrutinize this option. In [10], the authors explored different options for enhanced neutrino interaction. While they found that the vector forces of the aforementioned strength are already disfavored from laboratory experiments, light (below O(10 2 ) MeV) bosons strongly coupled to neutrinos remain viable. The only surviving option which alleviates the Hubble tension is ν τ -philic light scalar. This is kind of expected since it is well known that new interactions of ν τ are generically the least constrained compared to other flavors. It is therefore timely to revisit the constraints on such interactions from particle physics processes and pay particular attention to the interactions of τ neutrino.
For a light scalar interacting with ν τ many of the most sensitive probes of new physics connected to ν e and ν µ are not sensitive and the two most relevant laboratory bounds arise from τ and Z decays. The authors of Ref. [12] were the first to estimate the bound on the neutrino coupling to light scalar by studying the former process. The reported limit only applies to a particular choice of m φ and cannot be extrapolated to the mass range of interest easily. One of our goals in this paper is to derive this limit as a function of scalar mass by using state of the art numerical tools.
In Ref. [14] the authors present a comprehensive analysis of constraints on light neutrinophilic scalars. What is very interesting for us, in light of couplings to ν τ , is the constraint arising from invisible Z decay, namely the process Z → ννφ, where φ is light scalar. In addition to this process, we will also consider Z invisible decay (Z →νν) via a triangle loop diagram. Naively such a contribution may appear subdominant since it contains two powers of scalar coupling to neutrinos already at the amplitude level. However, it interferes with the SM tree level diagram and, therefore, the leading contribution is of the same order in the new physics coupling as the φ-bremsstrahlung and should be expected to give a competitive constraint.
The paper is organized as follows. In Section II we present the main results of our investigation of the new physics contribution to invisible Z decays while relegating the details of the calculation to appendices. In Section III we discuss the procedure for obtaining limits on new physics from leptonic τ -decays. We analyze the implications of our results for the allowed interaction strength of a neutrino-philic light mediator and comment on the implication for the proposed solution of the Hubble tension in Section IV. While the motivation for our study is mostly connected to τ neutrino flavor, for completeness we also present limits for a ν e and ν µ -philic scalar as well as flavor universal coupling scenario. In Section V we summarize our results and present our conclusions.

II. Z DECAY
The new neutrino interactions to be considered in this work are parameterized by where ν α and ν β are Dirac spinors of neutrinos (ν c α is the charge conjugate of ν α ) and α and β stand for flavor indices. Furthermore, y αβ = y βα is a symmetric Yukawa matrix and φ is a scalar field. Finally, note that the left projector P L = (1 − γ 5 )/2 ensures that only left-handed neutrinos are involved in the interaction.
In the presence of this interactions two new physics processes contribute to invisible Z decays and we show the corresponding Feynman diagrams in Fig. 1. The loop contribution contains two Yukawa vertices being proportional to |y αβ | 2 while the bremsstrahlung diagram is proportional to y αβ . Therefore, at the amplitude level, the left diagram is suppressed with respect to the right one by a higher power of Yukawa coupling. However, the loop diagram can interfere with the SM invisible Z decay and this substantially enhances its contribution to the decay width. Consequently, there is no obvious hierarchy between the two processes and both contributions to the invisible width should be considered in a complete analysis.

A. Loop Contribution
Let us first consider that two neutrino species with flavors denoted by α and β (α = β) are coupled to φ and other couplings in Eq. (1) are absent. The result obtained for this simple case can be easily generalized to the most general Yukawa matrix.
In the presence of a ν c α P L ν β φ interaction with α = β, there are new physics contributions to both Z →ν α ν α and Z →ν β ν β which have identical amplitude and therefore it is enough to only consider Z →ν α ν α . The 1-loop amplitude for Z →ν α ν α reads where g Z is the gauge coupling of Z boson to neutrinos; µ (q), u(p 2 ) and v(p 1 ) denote the external legs associated to the Z boson, neutrino and antineutrino, respectively while m φ (m Z ) is the mass of φ (Z). The 1-loop diagram for this process is UV divergent and to this end we have adopted gauge invariant dimensional regularization. The typical terms appearing in such calculation are abbreviated by Here, = (4−d)/2 with d representing the number of dimensions, while γ E is the Euler-Mascheroni constant. In a complete model, the UV-divergence arising from the considered diagram is expected to be canceled by other diagrams (including counter terms). Here "complete" means not only that all operators should be of dimension 4 or lower, but also that gauge invariance has to be respected. We would like to stress that the UV cancellation is model-dependent and, consequently, the finite part can not be predicted fully without committing to a specific UV-completion. Nevertheless, the log(m φ /m Z ) term is a generic feature and is independent of the regularization scheme. This can for instance be seen by considering only the loop integral with the loop momentum running between the scales of m φ and m Z which yields a result proportional to log(m φ /m Z ). This implies that this term can be physically interpreted as the contribution of the loop momentum running in the intermediate scale and being insensitive to the UV or IR behavior of the underlying complete models. Hence, we will retain only the log(m φ /m Z ) term and ignore other contributions in Eq. (2) to study the loop effect and also when setting the constraints. This allows to estimate the impact of the loop diagram on the invisible Z decay while being as model independent as possible. We refer the interested reader to appendix A where we show the cancellation explicitly in a toy model and find the behavior detailed above.
In the SM, the tree-level amplitude for Z →ν α ν α is which leads to the decay width 1 [35] By comparing Eq. (2) and Eq. (4), we can obtain the decay width including the loop contribution which yields One can check that the final result for the case α = β turns out to be the same as Eq. (6) with β → α. Therefore, in the presence of the most general Yukawa matrix, one only needs to replace |y αβ | 2 with β |y αβ | 2 in Eq. (6). If we sum over α indices and restrict ourselves to terms proportional to the second power of Yukawa coupling or lower we get where Y is the 3 × 3 Yukawa matrix with y αβ elements. One can also see that Eq. (7) is invariant under ν → U ν, Y → U Y U † basis transformations where U is an arbitrary unitary matrix.

B. Bremsstrahlung
The bremsstrahlung process is depicted by the right diagram in Fig. 1. Again, we first consider Z → ν α φν β with α = β. The decay width in case of α = β reads where For details of the derivation including the expression with the full m φ dependence of Γ new (Z → ν α φν β ) see Appx. B. The expression for the total width of the φ bremsstrahlung with the most general Yukawa couplings is given by where the first 1/2 factor is due to double counting of α =β , and the last 1/2 factor accounts for the phase space of identical particles. In the last term of Eq. (10), Γ new (Z → ν α φν β ) takes the same expression as Eq. (8). Eq. (10) allows the formulation in a basis-independent form similar to Eq. (7) 1 The relation between neutrino coupling to Z boson and the weak interaction strength GF reads The bremsstrahlung diagram in Fig. 1 represents Z → νφ * ν process. By flipping the arrows in the diagram one obtains a similar diagram for Z → νφν with identical decay width.
Also note that there is the charge conjugate process Z → νφ * ν, which has the same decay width as Z → νφν. Therefore, upon combining all the bremsstrahlung processes we reach the expression for the total contribution to invisible Z decay Lastly, we would like to stress that we simulated this three-body decay numerically in CalcHEP [36] and found excellent agreement with our analytic results.

III. TAU DECAY
Another relevant probe of neutrino self-interactions that is particularly relevant for ν τ is the decay of τ leptons. Similar to the φ-bremsstrahlung from Z decays the basic idea here is to constrain the scalar-neutrino coupling by investigating the impact of attaching a scalar line to the final state neutrino line; this for instance turns the diagram for the standard three-body decay into a charged lepton (electron or muon) and a pair of neutrinos into a 4-body process containing an extra light scalar boson in the final state. We illustrate this process in the left panel of Fig. 2. Most τ leptons decay hadronically but with a leptonic branching ratio Br l=e,µ ≈ 34% the leptonic final states are hardly suppressed. As the leptonic channels are much cleaner we focus on them in the following. A similar process has been considered previously in the context of a model with light majorons [12]. In principle the majoron limits from the literature could be used to estimate the bounds in the model under consideration here. We believe, however, that it is worthwhile to revisit the limits since the existing constraints are only available for m φ = 1 keV. As our analysis shows, the bound cannot simply be extrapolated to higher m φ and limits derived from rescaling the results of [12] become unreliable in the mass range of interest here.
We supplement the interaction term in Eq. (1)   of interest with MadGraph5 aMC@NLO [39]. As a cross check we first calculate the partial width for τ − → l −ν l ν τ where l = µ − or e − in the SM and find good agreement with the observed values. We determine the decay rate of the process τ − → l −ν lντ φ as a function of m φ numerically and construct a fit functions to derive the limit on the Yukawa coupling. In principle, the rates for the decay into electrons and muons are different due to the different masses of final state particles. In practice, the discrepancies are expected to be rather small due to the large hierarchy of charged lepton masses. We find that the differences between electron and muon channels are within the numerical uncertainties. The obtained partial width for the Yukawa coupling y τ τ equal to 1 is shown in Fig. 3 as a function of scalar mass, m φ . This decay rate is used for obtaining the limit as will be demonstrated in Section IV.
Finally, we would like to comment on another channel that can be constrained from τ decays. In the diagram for SM process τ − → l −ν l ν τ , the two neutrino lines could be connected with a scalar similar to the loop correction to Z → νν, see the right panel of Fig. 2 for an illustrative diagram. Note, however, that in contrast to Z decays only off-diagonal components of the Yukawa lead to a contribution that can interfere with the SM amplitude. These off-diagonal couplings are already strongly constrained by meson decays and, therefore, we do not consider this process further.

IV. CONSTRAINTS OF NEUTRINO INTERACTION WITH LIGHT SCALAR
With all the necessary ingredients available we can now turn to actual observables and derive limits on the parameters of the neutrino self-interaction model under consideration here. We will first consider the impact of the measurement of invisible Z decay before turning to τ decays.  τ d e c a y FIG. 4. Constraints on νSI from Z invisible decay (blue) and τ decay (orange) shown together with other known constraints taken from Ref. [14]. For the case of ν τ -philic scalar we also show the preferred region to relax the Hubble tension [10].
Combining the results in Eqs. (7) and (12), the total Z invisible width is given by Conveniently, the experimental measurement of Z invisible width can be expressed in terms of the number of light neutrino species [40] N ν = 2.984 ± 0.008, (14) which means that the observed invisible width is about 2σ lower than the SM prediction. Since both L and F in Eq. (13) are positive, the new physics we introduce can only enhance the Z invisible width. To get our limits we set the confidence level to 3σ so that the exclusion bound can be obtained by requiring In the case of τ decays the situation is more subtle. Since a ν τ is emitted in every τ decay a correction to all decay modes is expected for y τ α = 0. Naively, one could assume that the correction of the different decay modes is very similar since a φ emitted from the ν τ is only sensitive to the total momentum of the remaining final state. Consequently, the branching ratios remain similar to the SM prediction while the total width/lifetime of the τ changes. In contrast, a coupling to ν e or ν µ only affects the partial width of the leptonic decay modes. In order to derive reliable bounds on y τ τ and y µµ we make use of the partial width Γ τ µ for the three-body decay τ − → µ − ν τνµ which can be determined by combining the measured lifetime (290.6 ± 1.0) × 10 −15 s with the observed branching ratio of (17.41 ± 0.04)% [41]. The central value for the partial decay rate reads 3.94 × 10 −13 GeV. In order to get an estimate of the relative error on the leptonic partial width we add the relative errors of the lifetime and the branching ratio in quadrature and find δΓ τ µ /Γ τ µ ≈ 0.004. Therefore, we set the 3σ exclusion limit on the couplings by requiring A similar procedure utilizing τ − → e −ν τνe φ leads to essentially identical results for y ee since δΓ τ µ /Γ τ µ ≈ δΓ τ e /Γ τ e . When there are contributions to both τ − → µ −ν τνµ φ and τ − → e −ν τνe φ, we combine both channels together to set our limit.
In Fig. 4, we present our results; constrains on the diagonal elements of Y are calculated assuming the other elements of the Yukawa matrix are zero. More specifically, for the case of nonvanishing y ee , y µµ and y τ τ , we take tr[Y Y † ] = |y ee | 2 , |y µµ | 2 and |y τ τ | 2 , respectively. For the y ee = y µµ = y τ τ figure (lower right), we take tr[Y Y † ] = |y ee | 2 + |y µµ | 2 + |y τ τ | 2 . For flavor offdiagonal elements (y αβ with α = β), one can simply interpret bounds from any of these figures as the bounds on tr[Y Y † ] and convert it to the bounds on y αβ . For the coupling of ν e , ν µ as well as flavor universal scenario (upper panels as well as lower right panel) we also superimpose limits from meson decays [14]. As can be seen these bounds are stronger than those derived in this work for m φ 300 MeV. In the cosmologically most interesting case of a ν τ -philic scalar (lower left panel) we also show the preferred region for alleviating the Hubble tension (green) as well as a constraint from BBN [10]. While the derived laboratory constraints are certainly a relevant player for excluding the parameter space in y τ τ 0.1 range, the viable region still remains in the range 0.1 y τ τ 0.01. This points towards m φ ∼ O(10) MeV.

V. SUMMARY AND CONCLUSIONS
In this work we revisited constraints on neutrino self-interactions arising from a neutrino-philic light scalar φ. The employed probes are invisible Z decays and the leptonic decay modes of the τ . For invisible Z decays we consider two contributions: one withνν in the final state where we find that 1-loop diagram interfering with the usual SM contribution yields rather significant limit; the other, complementary, contribution to the invisible width arises from bremsstrahlung where two neutrinos (or antineutrinos) appear in the final state alongside φ. Summing both contributions we derive bounds on the new interactions for the case where the light scalar interacts with all flavors individually as well as a flavor universal scenario. In addition, we derive new limit from leptonic τ decays. To the best of our knowledge these are the first results that take the dependence on the φ mass and the coupling fully into account while previous calculations in the literature only apply for a restricted set of parameters. We provide a full picture of our results in Fig. 4 and compare them to constraints from meson decays. Our results constitute the leading bound on scalars with m φ 300 MeV irrespective of the preferred flavor. In the case of ν τ self-interactions, which is a particularly relevant scenario in light of recently proposed solution to the Hubble tension, these constraints constitute the leading laboratory limit throughout the considered mass range. However, a scalar in the mass range 10 − 100 MeV remains a viable option for large neutrino self-interactions and we are not able to exclude the whole parameter space preferred by cosmology. In this appendix, we discuss a toy model which is complete and rather minimal containing only a chiral fermion ν L , a gauged U (1) with the gauge boson denoted as Z µ and a scalar boson φ. Although the toy model is not realistic, it illustrates how the UV cancellation works explicitly and, in addition, shows the potential difference between the results in an incomplete model with respect to the complete one.
The model is formulated by the following Lagrangian Here, all terms are gauge invariant with the charge assignments ν L ∼ Q ν = −1 and φ ∼ Q φ = +2, except for the gauge boson mass term m 2 Z Z µ Z µ , which can be easily generated by, e.g., introducing another scalar that has a charge of +3 and a nonzero VEV. Note that such details are irrelevant for our discussion below. The covariant derivatives can be explicitly expressed as In Fig. 5, we present the Feynman diagrams involved in our analyses. We will show explicitly that the UV divergent parts in these diagrams cancel each other, as long as the U (1) charge is conserved (2Q ν + Q φ = 0). First, we compute the 1PI diagram generated by the Yukawa interaction, which will only lead to renormalization of the wave function of ν L . It will not lead to mass renormalization as one can expect from the chiral symmetry, so ν L remains massless after the loop corrections. The self-energy generated by the top left diagram in Fig. 5 reads with Here, we used Package-X [42] to evaluate the loop integral. The UV divergence in the neutrino self-energy is canceled by wave function renormalization The wave function renormalization generates a counter term δ Z ν L i / Dν L which then can be split to two counter terms δ Z ν L i / ∂ν L and δ Z ν L gQ ν Z µ ν L . The first term, corresponding to the top right diagram in Fig. 5, cancels the UV divergence in Eq. (A4); while the second term, corresponding to the bottom right diagram in Fig. 5 cancels the UV divergences of the two triangle diagrams in Fig. 5 2 . Now, by adding the counter term iδ Z / pP L to Eq. (A3) and requiring the UV cancellation, we obtain Next, we compute the Feynman diagrams for the Z µ → ν L ν L decay. The amplitudes of the three bottom diagrams in Fig. 5 are Note that Q ν c = −Q ν because instead of ν L (g Q ν )Z µ ν L , the Z-vertex should take the charge conjugate −ν c L (g Q ν )Z µ ν c L in the left bottom diagram. After computing the loop integrals and expanding the results in m 2 φ /q 2 (q 2 = m 2 Z ), we obtain Now we can clearly see that the the UV divergent parts in the above expressions cancel out if This corresponds to Q φ = −2Q ν , which can be understood from symmetry: ν c L φν L in Eq. (A1) respects the U (1) symmetry only if Q φ = −2Q ν .
Taking Q φ = −2Q ν , we get Setting q 2 = m 2 Z since we are interested in a decay we recover the log( m φ m Z ) term found in the simplified computation discussed in the main text while the constant term differs from the one in Eq. (2) as expected.
Appendix B: Analytical Calculation of three-body invisible Z decay The amplitude for the process Z(q) → ν α (p 1 )ν β (p 2 )φ(k) reads which leads to where we used the expression for the massive vector polarization sum polarizations (q) * (q) = −g µν + q µ q ν m 2 while E 1 and E 2 are energies of particles with 4-momenta p 1 and p 2 , respectively. By employing energy conservation E 1 = m Z − E 2 − E k the square matrix element |M| 2 can be expressed only in terms of 2 energies -one of massive (E k ) and one of effectively massless (E 2 ) final state particle. This allows for a straightforward evaluation of non-trivial three-body phase space integrals.
The differential decay rate reads where f (E 2 , E k ) is Eq. (B2) with the aforementioned substitution for E 1 . The integral g(E k ) = f (E 2 , E k )dE 2 can be evaluated analytically. We obtain the following result