{'fraction_non_alphanumeric': 0.03488228768642919, 'fraction_numerical': 0.01723864696870983, 'mean_word_length': 4.349157012338702, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Flavor-specific Interaction Favors Strong Neutrino Self-coupling in the Early Universe

Flavor-universal neutrino self-interaction has been shown to ease the tension between the values of the Hubble constant measured from early and late Universe data. We introduce a self-interaction structure that is flavor-specific in the three active neutrino framework. This is motivated by the stringent constraints on new secret interactions among electron and muon neutrinos from several laboratory experiments. Our study indicates the presence of a strongly interaction mode which implies a late-decoupling of the neutrinos just prior to matter radiation equality. Using the degeneracy of the coupling strength with other cosmological parameters, we explain the origin of this new mode as a result of better fit to certain features in the CMB data. We find that if only one or two of the three active neutrino flavors are interacting, then the statistical significance of the strongly-interacting neutrino mode increases substantially relative to the flavor-universal scenario. However, the central value of the coupling strength for this interaction mode does not change by any appreciable amount in the flavor-specific cases. We also briefly analyze a scenario with more than three neutrino species of which only one is self-interacting. In none of the cases, we find a large enough Hubble constant that could resolve the so-called Hubble tension.


Introduction
Neutrino remains to be the most elusive particle in the Standard Model (SM) even after sixtyfour years of its discovery. Even though in the SM, neutrinos are predicted to be massless, the observed phenomena of neutrino oscillation predicts at least one state of the neutrinos to have a mass 0.03 eV [1,2]. Neutrino oscillation requires mixing between different flavors, which implies that their mass matrix is non-diagonal in the flavor basis. Several neutrino oscillation experiments have now measured the neutrino mixing angles to few percent accuracy level. Terrestrial beta decay experiments probing the absolute masses of neutrinos, however, is yet to reach the sensitivity of that of the oscillation experiments [3][4][5]. In future, beta decay experiment KATRIN is expected to reach a sensitivity of 0.2 eV for electron neutrino mass [6].
Cosmology, on the other hand, can probe certain properties of neutrino remarkably better than any terrestrial experiment. After the advent of state-of-the-art observations of the cosmic microwave background (CMB) by the Planck collaboration, the sum of neutrino masses is now constrained to be below 0.12 eV through its dependence on the matter power spectrum [7]. CMB and Big Bang Nucleosynthesis (BBN) measurements of the effective number degrees of freedom in neutrinos in the vanilla ΛCDM cosmology agrees with the theoretical prediction of the SM. Interestingly, the constraints of neutrino properties from the cosmological data are also very sensitive to the underlying cosmological model. In other words, precision cosmological measurements are sensitive probes for beyond Standard Model (BSM) interactions of neutrinos. New physics such as a secret interaction among the neutrinos is difficult to probe in terrestrial oscillation experiments because of their very feeble interaction strength. Whereas, the dense environment in the early Universe can amplify the effects of any secret interaction in the neutrino sector due to large number density.
Several BSM interaction scenarios have been invoked in the neutrino sector to explain several experimental measurements and observations, e.g., anomalous neutrino signal in shortbaseline experiments [8,9], discrepancy between the CMB and local measurement of the Hubble parameter [7, 10] etc. One of the proposals to address these observations/tensions is selfinteracting neutrino (SINU) [11][12][13][14][15][16][17][18]. The deviation of the neutrino sector from its vanilla freestreaming nature have been shown to affect various cosmological observables. In the early universe, the dense neutrino environment gives us an opportunity to test these SINU models [19][20][21]. Several studies have been done on the effects of self-interaction among neutrinos on CMB and large scale structure [22][23][24][25][26][27][28][29][30]. Most of these studies are based on a parametrization of the self-interacting neutrino sector using c eff and c vis . However, with the availability of more precise experimental data, it is now important to model the SINU using a realistic particle physics-based approach. Previous studies have also shown that the neutrino flux from a supernova will be sensitive to new interactions among the neutrinos [31,32]. However, the strongest constraints on the parameter space of the SINU models come from K-meson decay, double beta decay, invisible width of the Z-boson, and τ -decay experiments [33][34][35][36]. Therefore, modelling a BSM SINU scenario that can address the above anomalies, while respecting the phenomenological constraints, is a challenging task [37,38].
As alluded earlier, self-interacting neutrinos have been suggested to ease a crisis of modern cosmology -the so-called Hubble tension. Since the early days of the Planck experiment, a discrepancy in the value of the Hubble parameter from CMB observations and local (low redshift) measurements have become a topic of active research [39]. This discrepancy has been exacerbated with the latest data from the Planck experiment and the local H 0 -measurement using cepheids and supernovae [7, 10,40,41]. The Planck-measured value of the Hubble constant using the temperature, polarization and lensing CMB data is H 0 = 67.36 ± 0.54 km s −1 Mpc −1 , whereas the SH0ES collaboration found H 0 = 74.03 ± 1.42 km s −1 Mpc −1 . There is now a 4.4σ tension between the two values of H 0 . However, an independent local measurement of H 0 was done in Ref. [42] which found a value H 0 = 69.8 ± 1.9 km s −1 Mpc −1 . Apart from unknown experimental systematics, this discrepancy could point towards a new physics signature beyond the ΛCDM paradigm. Several such scenarios have been envisaged to address this tension [15,16,18,. Among them, Ref. [15,16] introduced new secret interaction between Standard Model active neutrinos mediated by a new heavy scalar φ. A flavor-universal four-Fermi interaction with a coupling strength G eff was assumed, where g φ is the coupling between ν and the mediator φ, and M φ is the mass of φ. A Bayesian P18-TTTEEE+lowE+Lensing 3c + 0f 2c + 1f 1c + 2f Figure 1: The 1D posteriors for log 10 (G eff /MeV −2 ) in the three scenarios-3c + 0f , 2c + 1f , and 1c + 2f with three, two, and one self-interacting neutrinos, respectively. The two-mode feature is present in all cases. The relative significance of the larger G eff mode increases with less number of interacting neutrinos, hence, is largest in the 1c + 2f scenario. See section 4 for more details. analysis of this model with the CMB data prefers two distinct regions of G eff values: strongly interacting (SI) mode with a large value log 10 (G eff /MeV −2 ) −1.711 +0.099 −0.11 (68% confidence limit), and the moderately interacting (MI) mode with an upper bound log 10 (G eff /MeV −2 ) < −3.57 at 95% confidence level using the Planck 2015 and BAO data. Note that the SI mode coupling strength is about 10 9 times stronger than the weak interaction Fermi constant G F = 1.17 × 10 −11 MeV −2 .
However, SINU scenario meets strong constraints from the laboratory experiments [34,36]. Ref. [34] showed that a simple model of flavor-universal SINU scenario is ruled out by laboratory constraints from meson decay, τ decay, and double beta decay. Whereas, Ref. [36] demonstrated that Standard Model Effective Field Theory (SMEFT) operators, giving rise to flavor universal SINU scenario, are also highly constrained by meson decays, Z decays and electroweak precision measurements etc. These studies also derive flavor-dependent SINU constraints where only one of the neutrinos is interacting. The constraints in those cases are very stringent for ν e , whereas for ν µ and ν τ the constrains are comparatively weaker.
In this paper, we study flavor-specific neutrino self-interaction scenario using latest cosmological data. We perform a detailed Bayesian analysis of SINU with flavor-specific interactions for the first time using the latest 2018 likelihood from the Planck collaboration. As such, we allow the coupling strengths for different neutrino flavors to be different from each other. Our goal is to complement the flavor-specific SINU studies from laboratory experiments using the cosmological data. However, there is a notable distinction between the framework of SINU studies in cosmology vs laboratory searches. In laboratory experiments, the flavor of neutrino can be identified via electroweak interaction. Whereas, cosmological observables such as CMB are only sensitive to the temperature, free-streaming properties and the total mass of the neu-trinos. Therefore, in cosmology, we cannot distinguish any particular flavor of neutrino and can only study their effects collectively. All three generation of neutrinos ν e , ν µ and ν τ are on equal footing for our analysis.
In this work, we considered only three massless SM neutrinos, and therefore, fixed N eff = 3.046 (N eff ≈ 1.015 for each flavor). Massless neutrinos are a fairly good approximation as the bound on the total mass of the neutrinos are an order of smaller compared to the temperature of the plasma at the last scattering surface. The effect of massive neutrinos is reasonably well understood and can be speculated from our results. We consider three possible scenarios depending on the interactions of the three neutrino species -3-coupled (3c + 0f ), 2-coupled + 1free-streaming (2c + 1f ), and 1-coupled + 2-free-streaming (1c + 2f ) respectively. We assumed same coupling strengths for the coupled species for the first two cases. This is justified because the coupled neutrino species are completely equivalent due to the massless approximation. Setting neutrino mass to zero implies that the neutrinos remain in the same flavor eigenstate while propagating and, hence, there is no mixing between different flavor eigenstates. Note that, the phenomenological analysis performed here is also applicable to the scenarios where neutrinos have tiny masses, hence non-zero mixing, but the interaction is diagonal in the mass basis. However, in this paper, we treat neutrinos to be exactly massless and thus do not consider mixing. As a result, the flavor structure of the interaction remains the same during propagation. Also note that, the 3c + 0f case is identical to the universal flavor coupling scenario which has been studied previously. Additionally, we also analyze the scenario 1c + 2f +∆N eff , where the total N eff is allowed to vary, in appendix C.
Our main result is shown in figure 1. We find that the flavor-universal case yields a low significance for the SI mode with the Planck 2018 data. However, the SI mode significance is drastically increased once the flavor-universality of the coupling strength is relaxed, becoming maximum when only one neutrino state is self-interacting. We also find that the SI mode significance even surpasses the MI mode in certain cases. We explain the origin of the SI mode as a results of a better fit of certain features of the CMB data, compared to ΛCDM, using degeneracy of G eff with other cosmological parameters. When only one or two neutrino states are self-interacting, the resulting changes are smaller compared to the scenario when all three neutrinos are interacting and hence, can be compensated by the other correlated parameters relatively easily. This results in substantial enhancement of the significance of SI mode in 2c + 1f and 1c + 2f cases. We also find that the SI mode best-fit value of G eff has a mild dependence on the number of interacting species.
The outline of the paper is as follows. In section 2, we describe our cosmology model and explain the rationale behind it. In section 3 we detail our methodology and the experimental likelihoods used in this work. We show and interpret our results in section 4, and conclude with a discussion in section 5.

The model
We consider scalar interactions between massless ν and φ below the electroweak scale as where g ij is the coupling between φ and the neutrino flavors i and j 1 . When the temperature of the neutrino bath (T ν ) cools down below the mass of φ, i.e., T ν M φ , a four-Fermi interactions among the neutrinos, similar to Eq.(1.1), is generated as follows Here we note that for a general interaction of φ with any two flavors of neutrino, the four-Fermi interaction strength G (ijkl) eff in Eq.(2.2) has four indices for a process like ν i + ν j → ν k + ν l . Therefore, the most general scenario would involve many different couplings for different flavor combinations. In such a case, the energy and momentum of individual neutrino species will not be conserved and one will need to incorporate energy and momentum transfer between different species in the perturbation equations [17]. However, in this work we limit ourselves to a simpler scenario as described below. We consider only diagonal interactions in the flavor space with different coupling strengths: This assumption introduces only three new SINU parameters: G eff and G eff . As mentioned earlier, in this work we fix the number of neutrino flavors to three with all of them having a same temperature T ν . However, because of the complete equivalence of the interacting states in the context of CMB, we need to consider only one common coupling parameter G eff for all the interacting states for a given scenario.
We also assume M φ > 1 MeV to avoid BBN constraints on extra relativistic species around plasma temperature T ∼ 1 MeV [68,69]. The four-Fermi interaction picture is valid only below T ν = M φ when φ cannot be produced from scattering of the neutrinos. However for M φ > 1 MeV, its population is Boltzmann suppressed and hence does not affect the evolution of the density fluctuations during the period (T < 100 eV) relevant to the present analysis. The annihilation and decay of φ into neutrino increases the temperature of the neutrino bath, but this is modeldependent and we do not take this extra heating into account in this work. 2 Using dimensional analysis, the thermally averaged scattering cross section between the ith neutrino state goes as σv ∼ (G eff ) 2 T 2 ν . Therefore the interaction rate scales as Γ ν ≡ n ν σv ∼ (G eff ) 2 T 5 ν because of the T 3 ν -scaling of the neutrino number density. We absorb all other modeldependent prefactors of the interaction rate into G eff . The comoving neutrino self-interaction opacityτ ν is defined asτ where a is the scale factor of the Universe. Note that the T 5 ν -scaling of the opacity is a characteristic of the four-Fermi interaction. The neutrinos are self-interacting when the interaction rateτ ν is greater than the comoving Hubble expansion rate aH, i.e.,τ ν > aH. The interaction freezes-out whenτ ν drops below the Hubble expansion rate at a redshift z dec which is given by (2.5)  The chosen values of G eff are the best-fit values for SI and MI mode in 3c + 0f scenario with Planck temperature and polarization data. The top axis shows the angular multipole H corresponding to the modes that enter horizon at redshift z. The epochs of neutrino decoupling in ΛCDM and the matter-radiation equality are shown as gray-shaded region.
In the above equation, we have set other background cosmological parameters to their ΛCDM best-fit values and assumed the decoupling to take place in radiation domination era. We show the relative strength of the neutrino self-interaction to the Hubble expansion rate in figure 2. From figure 2, we see that the neutrinos are interacting with each other until z 4000 for G eff = 10 −1.86 MeV −2 which is much later than SM neutrino decoupling (z 10 9 ). For the MI modes the decoupling happens at z > a few × 10 5 . The neutrino decoupling for the SI mode happens very close to the matter-radiation equality which has interesting implications that will be discussed later.
Note that, because of the massless approximation, the structure of the interaction matrix does not change during propagation. However, this will not be true in case of massive neutrinos because of mixing between different flavor eigenstates. In that case, interaction of a particular flavor will generate interaction among different mass eigenstates depending on the structure of the mixing matrix. However, in our case, there is no mass mixing and the flavor eigenstates are the eigenstates of the Hamiltonian due to massless approximation.

Perturbation equations
The Boltzmann hierarchy of the perturbation equations for massless, self-interacting neutrinos in the Newtonian gauge is shown below following Refs. [15,70]. (2.6) Here α are -dependent O(1) angular coefficients that depend on the details of the neutrino interaction model, and the anisotropic stress σ ν is related to F ν,2 as σ ν = F ν,2 /2. Energy and momentum conservation dictates α 0 = α 1 = 0. The values of α for ≥ 2 is approximately of order unity. In this work, we take α = 1 for ≥ 2 which is a fairly good approximation. We implement these new perturbation equations for each neutrino species in the public code CLASS and solve them numerically [71,72]. 3 For very large value of the coupling G eff when |τ ν | aH, this system of equations may become difficult to solve as the equations become stiff, and tight-coupling approximation may be necessary. However, we checked that the default ndf15 integrator in CLASS is able to solve the equations for log 10 (G eff /MeV −2 ) ≤ −0.8 without invoking tight-coupling approximation which turns out to be sufficient for this work.
In the presence of self-interaction, neutrinos are strongly coupled at very early time which modifies the initial conditions for the perturbation variables compared to ΛCDM. For strongly coupled fluids, the initial anisotropic stress is zero. However, for large scale modes entering the horizon after neutrino decoupling this characterization is not true. For those modes, the initial anisotropic stress is nonzero as neutrinos are free-streaming when the modes enter horizon. To set unique initial conditions for all the modes of interest of CMB, we start evolving each k-mode in CLASS from a very high redshift (z ∼ 10 7 ). Because at z ∼ 10 7 , |τ ν | (aH) for all SINU modes including the MI mode, we can safely set initial anisotropic stress to zero for all k values. We also modify other initial conditions accordingly. However, we found that the resulting SINU spectrum with these modified initial conditions differ very slightly ( 0.1%) compared to SINU model with ΛCDM initial spectrum where initial anisotropic stress of neutrinos is non-zero. This is because, the anisotropic stress for a k mode, that starts with a non-zero value, vanishes very quickly due to the strong self-interaction. Since the modifications due to the initial conditions are very small compared to the precision of the Planck data and also makes the code comparatively much slower due to the additional integration time, we chose not to incorporate those for the mcmc analysis.

Changes in the CMB power spectra
In this section, we shall discuss the changes in the CMB angular power spectra due to neutrino self-interaction. The self-interaction stops the neutrinos from free-streaming before decoupling. As can be seen from Eq.(2.6), the new interaction plays the role of damping in the perturbations , and ΛCDM (black, solid). The white region in each plot, which is determined by |τ ν |/(aH) ≥ 10, is approximately the region in scale factor a upto which neutrinos are tighly-coupled for From the top left plot we see that the σ ν is suppressed in the non-free streaming region due to the self-interaction. The suppression is largest for 3c + 0f where all three neutrinos are tightly coupled, and gradually decreases for 2c + 1f and 1c + 2f , respectively. The suppression of σ ν results in the enhancement of φ and ψ which in turn enhances δ γ . The insets show the absolute changes in the corresponding cases compared to ΛCDM.
for ≥ 2. Therefore, it impedes the growth of the anisotropic stress σ while the neutrinos are strongly-coupled. We show the evolution of neutrino anisotropic stress σ ν in the top left panel of figure 3 in both SINU and ΛCDM cosmology. The important difference between them is the initial suppression of F ν,2 in the SINU scenario caused by the new interaction. The suppression is maximum for 3c + 0f where all three neutrinos are interacting, and gradually decreases for 2c + 1f and 1c + 2f where the number of interacting neutrino flavor is two and one respectively. The anisotropic stress is related to the gravitational potentials φ and ψ via the following equation, where ρ i , P i , σ i are individual energy density, pressure and anisotropic stress, respectively, for the i-th species, and ρ tot is the total energy density. In the last step of the above equation, we have ignored the small anisotropic stress of photon σ γ before recombination. Also, R ν is the fractional energy density of free streaming neutrinos which, in radiation domination, is given by In the radiation domination era R ν ∼ 0.41 in ΛCDM. Therefore, the suppression of neutrino anisotropic stress due to SINU during this period plays an important role in enhancing the gravitational potentials φ and ψ, as can be seen from figure 3. The gravitational potentials in turn affect the evolution of the photon perturbations as can be seen in the bottom-right panel of figure 3.
The changes in the CMB power spectrum in figure. 4 can be understood as a result of change in the propagation speed for perturbation of the neutrinos as explained below following Refs. [73,74]. In ΛCDM cosmology, the neutrino perturbations travel at the speed of light while the perturbations in the photon-baryon fluid propagates at the speed of sound c s 1/ √ 3. Therefore, the neutrinos tend to create perturbations ahead of the sound horizon. This creates a phase shift φ ν in the acoustic oscillations of the photon [73] , which in the radiation domination, is given by The amplitude of the oscillation is also changed by a factor of (1 + ∆ ν ) where This results in a suppression of the CMB power spectrum in ΛCDM due to free streaming neutrinos for the modes that enter the horizon before matter-radiation equality. Therefore, these changes in the photon acoustic oscillations affect the CMB angular power spectra for 200. Effectively, the phase shift moves the peaks in the TT and EE power spectra towards smaller , and suppresses the power spectra [73].
This story is changed in the presence of self-interacting neutrinos. The self-interaction stops the neutrinos from free-streaming and delays the neutrino decoupling from the thermal bath until a later time (z dec ) depending on the strength of the interaction. As a result, the free-streaming neutrino fraction R ν is decreased relative to its ΛCDM value depending on the number of neutrino species which are coupled at a certain time: Decreasing R ν in Eq.(2.9) and (2.10) readily implies a phase shift of the CMB spectrum towards larger , and an enhancement of power relative to ΛCDM. The resulting changes in the CMB TT and EE spectra can be seen in figure 4. In the upper panels, we show D XX = ( where XX = T T and EE in three SINU scenarios and ΛCDM. In the lower panels, we show the fractional changes of the spectrum relative to ΛCDM. We see that there is an enhancement in the SINU spectra, and also a phase shift which shows up as wiggles in the fractional difference plot. In compliance with the explanation above, both of these effects are maximal in 3c + 0f and gradually decreases as the number of self interacting states is decreased. Taking everything into account, we can see that the overall changes in the spectrum are milder when less number of neutrinos are interacting, which allow these changes to be compensated relatively easily by changing other parameters. A higher Hubble parameter can undo the change in the neutrino induced phase shift φ ν in SINU. This can be understood in terms of the photon transfer function cos(kr * s + φ ν ) which sources the acoustic peaks in the CMB spectrum. The CMB multipole corresponding to a mode k is given by, where D * A is the angular diameter distance, and r * s is the sound horizon at recombination defined as, where H(z) is the Hubble rate, and c s (z) ≈ 1/ √ 3 is the sound speed in the photon-baryon fluid. Decrease in the phase shift φ ν from its ΛCDM value moves the spectrum towards a higher value as can be seen from Eq.(2.12), which can be compensated by an increase in θ * ≡ r * s /D * A . This shift in θ * can be accommodated by changing D * A via modifying H 0 and Ω Λ which changes the Hubble evolution at late times [18]. The increase in Hubble constant H 0 is only relevant for SI mode where the self interaction strength is significantly strong inducing a large phase shift. In other words, higher values G eff is positively correlated with H 0 . Note that, the increase of H 0 for SI mode is maximum in 3c + 0f case where all three neutrinos interact producing largest phase-shift compare to ΛCDM.

Data & Methodology
We performed a Bayesian analysis of this model using the latest version of the Markov Chain Monte Carlo (MCMC) sampler MontePython3.3 [75,76]. To analyze the MCMC chains and plot the parameter posteriors, we use the GetDist 1.1.2 software package [77]. We use MultiNest interfaced with MontePython to sample the parameter space. The operational settings for Multinest that we use are shown in table 1 [78][79][80].
To constrain the model, we used the Planck 2018 likelihoods for the temperature and polarization power spectra [81]. Here 'TT+lowE' denotes the combination of low-TT ( < 30), low-EE and high-TT plik-lite ( ≥ 30) likelihood, and 'TTTEEE+lowE' denotes the combination of low-TT, low-EE and high-TTTEEE plik-lite likelihood. In addition we also use Planck 2018 lensing likelihood. For BAO, we used the 6DF Galaxy survey, SDSS-DR7 MGS data, and the BOSS measurement of BAO scale and f σ 8 from DR12 galaxy sample [82- Table 2: Prior ranges used for all parameters except log 10 (G eff /MeV −2 ).

Parameter
Prior  [41]. We use the following combinations of likelihoods in our analysis: 'TT+lowE', 'TTTEEE+lowE', 'TTTEEE+lowE+lensing', 'TTTEEE+lowE+lensing+BAO', and 'TTTEEE +lowE+lensing+BAO+H 0 '. In table 2, we show the prior ranges used. We use a log-prior for the extra parameter G eff because the expected features in its posterior distribution span over many orders of magnitude. Also, we do not vary the relativistic degrees of freedom N eff in this analysis to disentangle the effect of only neutrino self-coupling on the CMB power spectra, fix it to N eff = 3.046. Furthermore, to analyse the MI and SI modes individually, we used separate priors for log 10 (G eff /MeV −2 ) to distinguish them. The prior ranges for the two modes for all datasets are shown in Table 3 below. The prior ranges are chosen such that the boundary between MI and SI modes is approximately at the minimum of the valley region between the two peaks in the posterior of log 10 (G eff /MeV −2 ).

Results
In this section, we present the results of our MCMC analysis. In the 3c + 0f scenario, we assume that all of three species self-interact and we use the same value of G eff for all of them. This is because all of them are massless, and hence are equivalent. We have explicitly checked that even if we let the couplings for each species vary independently, their final posteriors are identical. Therefore, we use the same value of coupling without loss of generality. Below we present the results in flavor universal and flavor-specific scenario. In figure 5 we show the posteriors for log 10 (G eff /MeV −2 ) and other relevant parameters for all three SINU scenarios and ΛCDM.

Flavor-universal scenario: Three-coupled states (3c + 0f )
This is the only scenario which has been analysed before in Refs. [16,66,86]  Our results for the posteriors for log 10 (G eff /MeV −2 ) and other relevant parameters are shown in figure 6a. We find a multimodal posterior for log 10 (G eff /MeV −2 ) in agreement with the previous analyses. The SI mode corresponds to a larger value of log 10 (G eff /MeV −2 ) = −1.92 ± 0.18, whereas the MI mode is characterized by log 10 (G eff /MeV −2 ) = −4.35 ± 0.42 for  4). These values do not change appreciably for other datasets, except for TT+lowE, the reason of which we discuss later. The most important effect of including the Planck 2018 polarization data is that it suppresses the significance of the SI mode substantially compared to TT+lowE. Inclusion of the polarization data also shifts the log 10 (G eff /MeV −2 ) posterior towards smaller value as can be seen in the marginalized posteriors in figure 6a. In the rest of this subsection, we explain different aspects of these results. First, let us try to understand the origin of the SI mode which can be explained using the degeneracy of G eff with parameters. To this end, we first show the changes in the residuals of the TT, TE, and EE power spectra in figure 7 as we change G eff and the other correlated parameters, successively, to their SI mode best-fit values for TTTEEE+lowE dataset, starting from the ΛCDM best-fit spectrum. These other parameters are H 0 , A s , τ reio , n s , and ω m , respectively. When we first incorporate G eff (red, solid), we see that the spectrum moves upwards and shifts towards larger according to our discussion in the previous section. This gives rise to the positive residuals and the oscillations. Next, we change the best-fit value of H 0 (blue, dashed) which primarily compensates for the phase shift of the spectrum. What remains after this is mostly an overall amplitude offset barring a small amount of residual phase-shift. The amplitude offset is taken care of by the best-fit values of A s and τ reio (magenta, dotted). A smaller value of A s and a larger τ reio suppress the amplitude bringing down the residuals very close to the ΛCDM values. This overall amplitude change over-compensates the modifications at very large scale which entered the horizon after neutrino decoupling, resulting in the dip in the ∼ 200 region. A smaller value of the n s (orange, dot-dashed) corrects for this and red-tilts the spectrum increasing power at low while suppressing power at large . Finally, a smaller ω m (green, solid) reduces the residuals even more, bringing it very close to the ΛCDM spectrum. A similar
behavior is shown by the EE spectrum as well which, however, has much larger error bars at high than TT, and has less constraining power. This exercise shows that the Planck CMB data allows for a larger value of G eff using its degeneracy with other parameters giving rise to the SI mode. However, this compensation mechanism works for a special range of G eff values where it is large enough to enhance all the acoustic peaks so that the effects can be compensated by other global parameters, and also the same time does not impart very large modification of the first acoustic peak and large phase shift. For very high values of G eff , neutrino remains strongly coupled till very late times, even after recombination, and behaves almost always like a perfect 68 70 H 0 (km/s/Mpc)  Figure 6: Marginalised 1D posteriors for log 10 (G eff /MeV −2 ), H 0 , A s e −2τ reio , and n s for three SINU scenarios -3c + 0f (top), 2c + 1f (middle) and 1c + 2f (bottom) for three datasets.
fluid before recombination, and this scenario is disfavored by CMB data [74]. All of these prefer G eff values for which neutrino decoupling happens slightly prior matter-radiation equality in the SI mode. The origin of the MI mode is rather easier to understand. For smaller value of G eff , the changes in the spectrum are very small compared to the ΛCDM. This results in the plateau at small value in the posterior of log 10 (G eff /MeV −2 ). Now we shall give a quantitative explanation for the significance of the SI mode relative to the MI mode. From figure 6a, clearly the significance of the former depends on the dataset Geff + H0 + As + τreio Geff + H0 + As + τreio + ns Geff + H0 + As + τreio + ns + ωm used. To better understand the origin of this variation in significance, we first define a quantitỹ where D Planck is the binned Planck data, and D BF is the power spectrum corresponding to our best-fit point, and σ is the error bar of the Planck binned data. Therefore,χ 2 carries information about the goodness-of-fit of the spectrum in different regions of , and gives an idea about which part of the spectrum prefers/penalizes the fit. We define a quantityχ 2 asχ 2 = χ 2 which gives the approximate goodness-of-fit information in the whole spectrum. Here we are using the binned Planck data for visual clarity in the plots. Inχ 2 , we ignore the bin-by-bin correlation of the CMB data. Note that, theχ 2 is not the same quantity as the χ 2 calculated using the full Planck likelihood.
In figure 8, we plot the residuals andχ 2 for TT, TE, and EE spectrum in the highregion for the SI and MI mode best-fit points in 3c + 0f for the TT+lowE (left panels) and TTTEEE+lowE (right panels) datasets 5 . We also show the corresponding plots for the best-fit point in the 'valley' (VA) between the two modes which is defined within the range  P18 : TTTEEE + lowE Figure 8: Residual plots relative to the ΛCDM best-fit for high-( = 29 − 2500 for TT and = 29 − 2000 for TE and EE) modes for SINU with TT+lowE (left) and TTTEEE+lowE(right) datasets. In the sub-panel for each residual plot, we show the distribution ofχ 2 ( see Eq.(4.1)) with for ΛCDM and SINU modes. In the legend of each plot, we show the totalχ 2 . The gray-shaded TE (middle-left)and EE (bottom-left) plots for TT+lowE signifies that high-TE and EE mode data are not included in that analysis. VA denotes the 'best-fit' point in the 'valley' between the two modes in G eff posterior. −2.6 < log 10 (G eff /MeV −2 ) < −2.4. 6 We define the VA point to help explain the origin behind the separate SI mode. Firstly, we note that the MI residuals are very small and very close to the ΛCDM points, i.e., the x-axis. This is expected as the CMB data loses any sensitivity for G eff values below the MI mode limit, and the spectrum becomes practically indistinguishable from ΛCDM. Interestingly however, the SI mode best-fit point yields a better fit to the ΛCDM residuals which is evident in the TT spectrum in the top panels of figure 8. This substantially enhances the significance of the SI mode which is reflected in theχ 2 distribution. In several -regions, e.g. around 750, 1300, 1700, and 1900, the SI modeχ 2 is below the ΛCDM or the MI mode values, which largely compensates for the other regions with relatively worse fit. These are exactly the places where the SI best-fit point fits the ΛCDM residuals. In fact, according to the simplified analysis, the SI modeχ 2 for the TT+lowE data is better than both MI mode or ΛCDM. This is, however, not true for the TE and EE spectra though. The peaks in the EE spectrum are sharper than in the TT spectrum. Therefore, the polarization spectrum is more sensitive to phase shift due to SINU than the temperature spectrum, even though the error bars are larger [74]. This is reflected in the bottom panels in figure 8, where we see that the SI modẽ χ 2 is always equal to or above the MI or ΛCDM plots, yielding a poorer fit especially to the low and intermediate-polarization data. Similar behaviour is also observed in the TE residuals (middle panels). As a result, when the polarization data is included, the significance of the SI mode is reduced.
The presence of the 'valley' between the two modes can be understood from the lowest mode that is affected by G eff . From the top axis of figure 2, we see that a typical SI mode G eff affects all modes 200 which happens to be the approximate position of the first peak in the CMB spectrum. Therefore, it is easier to compensate for the modifications by changing the other parameters, as explained earlier, yielding a reasonably good fit. On the other hand, MI mode values of G eff affects only 10 4 which are far beyond the range of the present CMB experiments. Planck measured the TT spectrum only upto = 2500, and the EE spectrum upto = 2000. Therefore, the MI mode is virtually indistinguishable from ΛCDM as far as the Planck data is concerned. For intermediate values of G eff between the two modes, the CMB spectrum is modified only in the high-part of the Planck -range. In this case the degeneracy with other parameters, which impart changes in whole spectrum, cannot be exploited to get an overall good fit. Thus, a 'valley' appears between the two modes as a result of the poor fit to the CMB data.
The polarization data (TTTEEE+lowE) also shifts of the whole posterior G eff to the left relative to TT+lowE which evident from the left panels of figure 6. This can be understood from the relative sizes of the errorbars in the TT, TE, and EE spectra. The TT spectrum has smaller errorbars at large . Whereas, the TE and EE spectra uncertainties are smaller at small (see figure 8). This essentially means that the TT data can accommodate large deviations at small-, whereas TE and EE data have more freedom at larger-. Therefore, inclusion of the polarization data penalizes any deviation at low . Because, stronger self-interaction implies a later decoupling of the neutrinos affecting relatively smaller , TTTEEE+lowE data prefers a slightly smaller value of G eff compared to TT+lowE. , and EE (bottom) relative to ΛCDM in three scenarios. We note that these curves exhibit a better fit to the ΛCDM residuals which shown in gray. circles.

Flavor-specific scenario: Two-coupled (2c + 1f ) & one-coupled states (1c + 2f )
In this section, we shall concentrate on the flavor-specific neutrino self-interaction scenarios, namely, 2c + 1f with only two self-interacting neutrino states, and 1c + 2f with only one selfinteracting state. For the results in this section, we used the same coupling strength for the two interacting states instead of having two independent couplings in the MCMC sampling. We have explicitly verified the validity of this assumption by running chains with independent couplings and same prior for the states which yielded same posterior distributions. This was rather expected as all states are massless, and are not distinguishable from one another as far as their cosmology is concerned. Hence, we chose the same coupling to derive all the results here reducing the number of free parameters without loss of generality. In passing, we also note that the one could have chosen different priors for different states inspired by the bounds from BBN and laboratory experiments which would yield distinct posteriors [34,36]. We want to point out that 1c + 2f scenario is equivalent to the flavor-specific ν e /ν µ /ν τ coupling cases discussed in Ref. [34]. The posteriors for log 10 (G eff /MeV −2 ) are shown in figure 5, and the inferred parameter values and their 68% confidence limits for the TT+lowE and TTTEEE+lowE datasets are given in table 5 for 2c + 1f . The two modes are present in both of these cases. The SI mode appears at almost at the same value of G eff as in 3c + 0f . As discussed before, this specific value of G eff is determined by the -range of the Planck data. The interplay between G eff and other degenerate parameters seeks out this value. At the same time there is a very small shift in the SI mode G eff values across three SINU scenarios. However, the most striking difference between 3c + 0f and 2c + 1f or 1c + 2f is the enhancement of the SI modes for all datasets in the latter, as can be seen from figure 1 and 5 (see also figure 14 and 15). The enhancement is more pronounced in 1c + 2f than 2c + 1f . From Using Eq. (2.9), (2.10), (2.11) and we see that the phase shift and the amplitude suppression due to neutrino free-streaming are proportional to the number of free-streaming neutrino states. With less number of self-interacting states, these changes in the CMB spectra are relatively milder compared to 3c + 0f . This is evident from figure 4. As a result, there exists more room to use the degeneracy between G eff and other parameters, especially A s and n s , to achieve a better fit. This is seen in figure 9 where we see that SI mode the residuals in the flavor-specific scenario are always smaller than the flavor universal case.
Another new feature here is that the upward rise of the valley between the two modes in the log 10 (G eff /MeV −2 ) posterior. In the previous section, we explained the origin of the valley as a result of the fact that intermediate values of G eff affect only a part of the CMB spectra observed by Planck, which cannot be compensated by varying other degenerate parameters. Hence, the intermediate values are disfavored. However in 2c + 1f or 1c + 2f , those changes are comparatively modest, and can be partly undone by varying other parameters. As a result, the valley points yield smaller χ 2 compared to 3c + 0f .

Effects on H 0 and σ 8
We know that neutrino self-coupling strength G eff has a positive correlation with H 0 through the phase-shift as explained in section 2.3. From table 4, we see that the mean value of the Hubble shifts to H 0 = 69.46 ± 0.52 km s −1 Mpc −1 for the TTTEEE+lowE dataset for the SI mode, reducing the tension with the SH0ES measurement to a modest ∼ 3σ. On the other hand, when compared with the CCHP measurement H 0 = 69.8 ± 1.9 km s −1 Mpc −1 which uses the tip of the red giant branch calibration, our value of H 0 is fully consistent with that [42]. However, the SI mode value of H 0 slightly decreases in the flavor-specific 2c + 1f and 1c + 2f scenarios because of smaller phase shift, but the significance of the mode is increased. This is evident from figure 10. In table 7, we show H 0 , Ω Λ , 100θ s , r * s , and D * A for 3c + 0f , 2c + 1f and ΛCDM. The shift in θ s due to the phase shift in SINU decreases D * A that helps increase the value of H 0 and Ω Λ . Note that, the value of r * s changes by only ∼ 1σ from ΛCDM, and plays a sub-dominant role in changing the Hubble constant.
The significance of the SI mode is enhanced when the local Hubble measurement data H 0 = 74.0 ± 1.4 km s −1 Mpc −1 from SH0ES is included [87]. At the same time, it suppresses the  Figure 10: The contours of 68%, 95%, and 99% confidence levels for log 10 (G eff /MeV −2 ) and H 0 in 1c + 2f (green), 2c + 1f (blue), and 3c + 0f (red) for the Planck TTTEEE+lowE+lensing dataset. The SI mode contour for 3c + 0f yields the largest value of H 0 because of bigger phase shift. However, the significance of the mode is much less compared to 2c + 1f and 1c + 2f .
MI mode significantly as can be seen from figure 6. This is not surprising given the fact that the SH0ES likelihood includes only Hubble data, and favors SI mode which yields a larger H 0 . However, one should be cautious about combining the Planck with the SH0ES data because of the more than ∼ 4σ discrepancy between the two. Even though we show the results with Planck+lensing+BAO+SH0ES here, we do not draw any conclusion from them. Note that, the 3c + 0f scenario yields higher H 0 but its significance is smaller. In 2c + 1f and 1c + 2f scenarios, although the significance of the SI mode is increased, the corresponding value of H 0 decreases due to smaller phase shift. Therefore, both flavor-universal and flavor-specific SINU scenarios fail to yield a large enough value of H 0 to completely solve the Hubble tension (see figure 6). The parameters inferred from the combined TTTEEE+lowE+lensing+BAO+H0 dataset are quoted in table 10 and 11. In appendix C, we show the results for the model 1c + 2f +∆N eff , i.e., when one neutrino species is self-interacting and additional free-streaming radiation with effective neutrino number ∆N eff is present. The value of H 0 increases slightly compared to the 1c + 2f scenario due to presence of the additional radiation relics which can be seen from Table 13.
Apart from H 0 , another parameter that shows a small discrepancy between CMB and low redshift data is the late-time matter clustering amplitude σ 8 . Planck 2018 results found σ 8 = 0.811 ± 0.006 using TTTEEE+lowE+lensing data [7]. However, the gravitational weak lensing observation by KiDS-1000 found σ 8 = 0.76 +0.025 −0.020 [88]. There is about ∼ 2.4σ discrepancy between the two values. The value of σ 8 in SINU cosmology is mainly determined by two competing effects. The lack of anisotropic stress in the neutrino bath boosts the perturbation modes which enter horizon before neutrino decoupling. In case of SI mode values of G eff , these modes happen to be in range of σ 8 , and thereby boosting its value. On the other hand, the lower value of best-fit A s and n s suppress σ 8 . As a result, G eff and σ 8 show a very weak correlation as can be seen in figure 13, 14, and 15. We find σ 8 = 0.834 ± 0.0067 and 0.829 ± 0.0058 in SI mode for 3c + 0f and 2c + 1f , respectively for TTTEEE+lowE (see table 4 and 5).

Implications of BAO measurements on SINU
In this subsection, we examine in details the effects of the BAO measurements on SINU cosmology. The BAO likelihood constraints the distance combinations D M (z)/r d and H(z)r d where Here, r d is the sound horizon at baryon-photon decoupling, z d is the redshift of that epoch, and and D A (z) is the angular diameter distance 7 . By the virtue of phase shift due to the strongly interacting neutrinos, the SI mode accommodates a larger Hubble constant compared to the MI mode, as discussed previously.
7 Note that, the phase shift in the CMB acoustic peaks due to interacting neutrinos will also induce phase shift in the BAO peaks [90,91] (also see [18]). For SI mode in particular, the effect of this scale dependent phase shift of BAO peaks with respect to ΛCDM cosmology can be important. The BAO scales extracted from the observed spectrum assuming fiducial ΛCDM cosmology [84], may differ from the true BAO scales in SINU cosmology due to the phase shift (see Ref. [92] for quantitative study of the neutrino phase-shift on BAO scales in different cosmologies). In this work we ignore this subtlety and use the BAO likelihood as reported in Ref. [84] −5 P18-TTTEEE+lowE+Lensing P18-TTTEEE+lowE+Lensing+BAO 3c + 0f 2c + 1f 1c + 2f Figure 12: 1-D posterior of log 10 (G eff /MeV −2 ) for 1c + 2f , 2c + 1f , 3c + 0f with TT-TEEE+lowE+lensing and TTTEEE+lowE+lensing+BAO datasets. Inclusion of the BAO data reduces the relative signification of the SI mode peak compared to the MI mode.
We show in figure 11 the evolution of H(z)r d and D M (z)/r d (with additional scaling of the scale factor a = 1/(1+z)) for the TTTEEE+lowE+lensing bestfit values in 3c + 0f and 1c + 2f scenarios along with the BOSS data for the BAO scales. It is evident from the figure that the BOSS data mildly disfavors the SI mode compared to the MI mode (which is virtually identical to ΛCDM). This can also be seen from table 8 where we compare the inferred parameters values from Planck-only and Planck+BAO data. Therefore with the inclusion of the BAO likelihood, the relative significance of the SI mode drops compared to the MI mode which can be seen from figure 12. The decrease in relative significance for SI mode is highest in 3c + 0f and lowest in 1c + 2f . This is because the increase in H 0 is the least in 1c + 2f among all three scenarios. The H(z) and D M (z) evolutions of the SI mode in 1c + 2f differ the least from the corresponding MI mode evolutions and are in better agreement with the BAO data as shown in figure 11. Thus, upon inclusion of the BAO data, the relative significance of the SI peaks is affected the least.

Mode comparison
In this subsection, we analyze the relative significance of the two modes in the SINU posterior. We compute the Bayesian evidence Z for the two modes as follows Here M is the cosmological model. The Bayesian evidence essentially is the probability of the data d given a model M . We use the Bayesian evidence provided by Multinest for the two modes. Then we compute the ratio of the probabilities Pr(M |d) of model M for given data d, which is the Bayes factor, to find the relative significance of each mode assuming equal prior.
A Bayes factor B SI > 1, therefore, means that the SI mode is more significant than the MI mode.

Discussion & Outlook
Secret self-interactions among the active neutrinos could leave an observable imprint in the CMB anisotropy power spectra. Such interaction stops the neutrinos from free-streaming much later than weak decoupling . As a result, the CMB power spectra experience a phase shift and an enhancement for the modes which entered the Hubble horizon before the decoupling of the neutrinos. We studied these effects in three cases-3c + 0f , 2c + 1f , and 1c + 2f with three, two, and one self-interacting massless neutrino species, respectively, and their cosmological implications using the latest CMB data from Planck 2018, BAO data, and local Hubble measurement data from SH0ES. This study of the flavor-specific scenario is inspired by the recent strong constraints from several laboratory experiments on the flavor universal scenario. We make the following key observations in this work.
• In all three cases and all datasets, the posterior of log 10 (G eff /MeV −2 ) has two distinct modes, namely, the SI mode characterized by a larger log 10 (G eff /MeV −2 ), and the MI mode with a smaller value. The origin of the SI mode is explained as a result of degeneracies between G eff and other ΛCDM parameters. The absence of free streaming neutrinos in SINU leads to a phase shift and enhancement of the CMB angular power spectra. These changes are compensated by the other ΛCDM parameters. Specifically, H 0 compensates for the phase shift, A s and τ reio correct the overall amplitude, and n s tilts the whole spectrum to achieve a good fit to the data. We showed that this compensation mechanism works only for those values of G eff for which neutrino decoupling happens close to matter radiation equality.
• In intermediate values in the 'valley' between the SI and MI modes are not favored by data.
The large G eff corresponding to the SI mode globally affects the whole CMB spectrum observed by Planck, which can be undone by changing other parameters. In contrast, the MI mode spectrum is virtually indistinguishable from the ΛCDM as it only affects very high-modes, due to small G eff , which are not observed by Planck. The valley in between corresponds to intermediate G eff values which modifies only the high-portion of the spectrum. These partial changes are difficult to compensate using other parameters, resulting to a worse fir to the data.
• In the 3c + 0f scenario, the significance of the SI mode is greatly diminished if CMB polarization data is included. This is a consequence of the relatively poorer fit of the model to the low-l polarization data. Even though the phase shift from the free-streaming neutrinos is the same for both temperature and polarization spectra, the peaks in the polarization spectrum are relatively sharper. As a result, the EE data is more sensitive to any change in the number of free-streaming neutrinos than the TT data [74]. The inclusion of the polarization data also shifts the whole posteriors towards smaller G eff .
• However, the SI mode is rejuvenated in the flavor-specific 2c + 1f and 1c + 2f scenarios which are in lesser conflict with laboratory experiments. We showed that even with the polarization data, the SI mode significance is comparable or sometimes even greater than that of the MI mode. In these cases, the number of self-interacting neutrinos are less than 3c + 0f , and as a result, the changes in the CMB spectra are also relatively moderate. This allows for more freedom to use other degenerate parameters to achieve a significantly better fit.
• Due to the phase shift in CMB spectrum, SINU scenario favors a larger H 0 in the SI mode. For example, in the 3c + 0f scenario, we find H 0 = 69.52 ± 0.41 km s −1 Mpc −1 for Planck temperature, polarization, and lensing data. Although in flavor-specific scenario, the significance of the SI mode increases, the value of H 0 slightly decreases due to lesser phase-shift. We do not find any strong correlation between σ 8 and G eff .
In this work, we have confined our discussion within the standard framework of only three active neutrino states. However, additional sterile neutrino states are interesting to consider for various different reasons. The presence of such an additional neutrino state(s) might change our results significantly. To begin with, the relativistic degrees of freedom N eff is expected to increase, which would in turn increase the value of H 0 [66]. As hinted by the short-baseline neutrino oscillation experiments, the mass of the additional sterile state is predicted to be around ∼ 1 eV. In appendix C, we briefly discuss the results when extra radiation (∆N eff ) is added. This is motivated by the sterile neutrino models. In this case, it is possible that only the sterile state can have large self-interaction. This scenario will evade all laboratory constraints on the active sector. However, as can be seen from table 13, we do not find a large enough value of H 0 that can reconcile its value with the low redshift measurement data. Also, in order to a more careful CMB analysis with sterile neutrinos, the mass of the sterile state needs to be taken into account which could be of the order of the recombination temperature for ∼ 1 eV sterile neutrino. The finite neutrino mass will have further implications for our analysis. In presence of massive neutrinos, the flavor-structure of the interaction matrix would certainly be non-diagonal which can have interesting signatures. Finite neutrino mass will also affect the matter power spectrum at late times. The strong neutrino self interaction will also enhance the primordial CMB B-mode at small scale which can have interesting implications for future B-mode experiments [18,93].
In this work, we have kept the helium fraction Y p fixed to its BBN value. However, there is a well-known degeneracy between Y p and the number of free-streaming neutrinos. We plan to investigate it in a future work. Note that, the present and all previous analyses on this topic have assumed a diagonal coupling matrix between different neutrino flavors. However, this is far from a realistic scenario as off-diagonal couplings will be present even in the most simple model of neutrino self-interaction. The evolution of the cosmological perturbations become much more complicated when off-diagonal couplings are present, and a more general set of evolution equations need to be used [17].
Future CMB observations like, CMB-S4, will make measurements at smaller scales and will probe even higher values than Planck. This would probe the neutrino interaction at even earlier times, and perhaps could shed light on the nature of the mediator particle φ. This work is the first step towards studying a more general and realistic flavor profile of neutrino self interaction. A more rigorous analysis including CMB, BBN, and laboratory experiment data would be worthwhile.
Government of India during the initial stage of this project. This work used the computational facility of Department of Theoretical Physics, Tata Institute of Fundamental Research.

A Parameter values for other datasets
Here we show the parameter values and their 68% confidence limits for the TTTEEE+lowE+lensing and TTTEEE+lowE+lensing+BAO+H 0 datasets.

B Posterior distributions of all parameters C Flavor dependent SINU with additional radiation species
In this section, we briefly discuss the effects of SINU cosmology with varying N eff . Flavor universal SINU with additional interacting radiation species are interesting as they can accommodate even higher value of the Hubble constant. For this analysis, we chose the scenario 1c + 2f + ∆N eff which contains additional free streaming radiation with effective neutrino number ∆N eff . One possible model interpretation of this setup is the existence of sterile neutrino state(s) among which one is self-interacting, while all other neutrinos are free-streaming. Figure 16 shows the 1-D and 2-D posteriors of all the parameters. The value of H 0 increase slightly compared to the 1c + 2f scenario due to presence of the additional radiation relics which can be seen from Table 13.