Phantom fluid cosmology: impact of a phantom hidden sector on cosmological observables

Phantom scalar theories are widely considered in cosmology, but rarely at the quantum level, where they give rise to negative-energy ghost particles. These cause decay of the vacuum into gravitons and photons, violating observational gamma-ray limits unless the ghosts are effective degrees of freedom with a cutoff $\Lambda$ at the few-MeV scale. We update the constraints on this scale, finding that $\Lambda \lesssim 19$ MeV. We further explore the possible coupling of ghosts to a light, possibly massless, hidden sector particle, such as a sterile neutrino. Vacuum decays can then cause the dark matter density of the universe to grow at late times. The combined phantom plus dark matter fluid has an effective equation of state $w<-1$, and functions as a new source of dark energy. We derive constraints from cosmological observables on the rate of vacuum decay into such a phantom fluid. We find a mild preference for the ghost model over the standard cosmological one, and a modest amelioration of the Hubble and $S_8$ tensions.

Ghosts, or phantoms, are scalar fields with a wrongsign kinetic term.They have been widely studied in cosmology, as a source of matter with equation of state w < −1, as was suggested in early determinations of the dark energy equation of state [1][2][3].If such an object exists in nature, it must have not only classical behavior, but also quantum mechanical: negative energy particles which are the quanta of the ghost field. 1ven if such exotic particles have no direct couplings to the standard model (SM), they cannot avoid being gravitationally coupled.As a result, through virtual graviton exchange, the vacuum is unstable to decay into pairs of ghosts plus positive-energy particles [4], and the rate for such processes is unbounded unless there is an ultraviolet cutoff Λ on the phase space of the ghosts: in other words, any consistent theory of ghosts must be a low-energy approximation that breaks down at momenta larger than Λ.Non-observation of excess gamma rays from vacuum decay leads to a bound that was roughly estimated as Λ ≲ 3 MeV [5].
As pointed out in Ref. [5], the low-energy effective theory of ghosts must be Lorentz-violating, since Λ is a bound on the magnitude of the 3-momenta of the ghosts, which implies a preferred reference frame.Explicit, UVcomplete models that give such a momentum cutoff were constructed in Refs.[6,7].In the present work, we will not be concerned with the microphysics underlying such ghost models, but rather with the late-time cosmological implications of the low-energy effective theory.
Another possible signature of phantom theories could be decay of the vacuum into ghosts plus dark matter (DM).If both types of particles had the same dispersion relation, apart from the overall sign, there would be no net cosmological effect, at the homogeneous level, due to the cancellation of stress-energies of the two components.Here we will consider the case of massless ghosts, which can have a nonvanishing effect, since the negative energy density of the ghosts redshifts faster than that of the created dark matter.The result is a net creation of DM at late times.We will show that for purely gravitationally coupled ghosts, this process is too slow to give detectable results.But it is consistent to couple the ghosts and the DM more strongly, for example by exchange of a heavy gauge boson, allowing for significant late-time DM creation.
In the present work, we examine the effect of such latetime DM creation on the cosmic microwave background (CMB) and structure formation.Even if the contributions of the two fluids cancel exactly at the homogeneous level, in the case of massless ghosts and massless sterile neutrinos, their perturbations do not cancel.This leads to constraints on the rate of DM plus phantom fluid creation, and consequently on the energy scale of new physics linking phantom particles to dark matter.
We start by revisiting the upper bound on Λ from diffuse gamma rays in Section II, and make a fully accurate determination.In Section III we describe the particle physics model for a phantom sector coupled to dark matter, and solve for the cosmological component densities created by decay of the vacuum.There we also derive preliminary constraints based on type Ia supernovae.In Section IV, the consequences for the CMB and large scale structure are studied using a full Monte Carlo search of the parameter space, resulting in an upper bound on the vacuum decay rate and corresponding new-physics energy scales.We investigate the potential of the model for providing a new source of dark energy, and for addressing various cosmological tensions.Conclusions are given in Section V.

II. CONSTRAINT ON Λ
In Ref. [5], a rough estimate was made of the rate for vacuum decay into ghosts and photons.Here we compute the rate more quantitatively.The computation is identical to that for gravitational scattering of ghosts into photons, ϕϕ → γγ, except that the initial state parti-cles are treated as final state particles, and integrated over with momenta restricted by |⃗ p i | ≤ Λ.The matrix element has been computed in Refs.[8,9], and is given by where the Mandelstam invariants are defined in the usual way, and we can pretend that the ghost energies are positive since the momentum-conserving delta function is ) where E 1,2 are ghost energies and E 3,4 are photon energies.Here m P = 2.43 × 10 18 GeV is the reduced Planck mass.
Carrying out the phase space integrals as described in appendix A, we first compute the instantaneous spectrum of gamma rays dΓ/dE.The result is shown as the black solid curve in Fig. 1.However the observed flux is a line-of-sight integral, similar to that for decaying dark matter [10], which can be expressed as an integral over the redshift of the corresponding emission distance, where the Hubble rate is H(z) ∼ = H 0 Ω m (1 + z) 3 + Ω Λ , ignoring the small contribution from radiation.Here Ω m ≃ 0.32 and Ω Λ ≃ 0.68 are the total matter and dark energy abundance, respectively.The result is shown as the red dashed curve in Fig. 1, in units of Λ 7 /(4πH 0 m 4 P ).This must be compared to the observed diffuse electromagnetic spectrum, which in the region of interest was measured by the COMP-TEL experiment [11] to be approximately dJ/dE ∼ = 5.3 × 10 −3 (E/E 0 ) −2. 4 MeV −1 cm −2 s −1 sr −1 with E 0 = 1 MeV.By demanding that the predicted flux in Eq. ( 2) not exceed the observed one, we derive a 95% Confidence Level (C.L.) limit of Λ ≲ 18.8 MeV . ( which is somewhat weaker than the rough estimate originally derived in Ref. [5].Here we assumed the ghost particles are real scalars; the bound in Eq. ( 3) should be reduced by a factor (1/2) 1/9.4 ≈ 0.93 for complex scalar ghosts, giving Λ ≲ 17.5 MeV.
We have also estimated the constraint arising from vacuum decay to ghosts and e + e − pairs; using the matrix elements computed in Ref. [9] (see Eq. ( 4) below), we find a decay rate per unit volume Γ d = Λ 8 /(368640 π 5 m 4 P ) in the approximation m e ≪ Λ.The present density of e ± is of order n e = Γ d /H 0 , giving the rate Γ s ∼ n e σ for a CMB FIG. 1. Solid: instantaneous photon spectrum (Λ× decay rate per unit volume and energy) from vacuum → ϕϕγγ (real scalar ghosts plus photons) decay, with units of Λ 8 /m4 P .Dashed: observed flux accounting for cosmological redshift, with units of Λ 7 /(4πH0m 4 P ), as given by Eq. ( 2).

III. GHOST/DARK MATTER MODELS
We next consider the decay of the vacuum into massless ghosts ϕ and a massive hidden sector particles ν s , which for definiteness is considered to be a sterile Dirac neutrino.(The main results are not expected to be sensitive to different possible choices, such as scalar dark matter.)Such a process can increase the energy density of the Universe at late times, due to the faster redshifting of the ghost versus the massive particle contribution.

III.1. Gravitational coupling
The matrix element for the gravitational scattering between Dirac fermions with mass m νs and scalar particles is [9] where the Mandelstam invariants s, t, u have their usual definitions, after changing p 1,2 → −p 1,2 .The vacuum decay rate, per unit volume, is given by the phase-space integral Γ = (2π) 4 4 i=1 where we have taken p 1,2 → −p 1,2 , so that all energies are positive.The ghost momenta p 1 and p 2 should be integrated up to the cutoff Λ.In appendix A we carry out the numerical integration of the vacuum decay rate per unit volume.By numerically fitting, the result can be approximated as3 assuming the ghosts ϕ are real scalar particles, and that m νs ≤ Λ by energy conservation.
One can estimate the expected energy density in ν s today following Ref.[5], where n νs (t) is the number density of sterile neutrinos and t 0 ∼ 1/H 0 is the age of the Universe.Taking the maximum value Λ ∼ 19 MeV and dividing by the critical density today ρ crit , we find that Ω νs = ρ νs /ρ crit ≲ 10 −9 , which is far below the sensitivity of observational constraints.It follows that the vacuum decay rate must be ∼ 9 orders of magnitude larger than that arising from purely gravitational interactions to have an observable cosmological impact.

III.2. Scalar or vector coupling
Since the gravitational interaction of ν s is too weak to have appreciable cosmological effects, we introduce a heavy mediator coupling ϕ to ν s , which can be either a scalar Φ or a vector Z ′ .Integrating out the mediator results in an effective interaction of the form Λ where M s and M v are the new-physics mass scales, to be constrained. 4For simplicity, we take the ghost to be a complex scalar in both cases, and ν s to be a Dirac fermion.This choice insures there is no gauge anomaly for the vector mediator case.Fig. 2 shows the Feynman diagrams of the vacuum decay processes.Labeling the ghost momenta as p 1,2 and the ν s momenta as p 3,4 , we obtain the following new-physics matrix elements for the vacuum decay where again the Mandelstam invariants s, t, u have their usual definitions, after negating the 4-momenta of the ghosts so that kinematically they look like initial-state particles with positive energy.Integrating Eq. ( 5) numerically we find that the decay rate for both mediators can be approximately fit by the formulas 5 (see appendix A).

III.3. Phantom sector spectra
Once the massless ghosts ϕ and the massive sterile neutrinos ν s are produced from the vacuum decay, their cosmological evolution is described by the Boltzmann equation in an expanding universe, where p i = |⃗ p i |, t is the physical time, H is the Hubble parameter, f 0 i (t, p i ) is the (dimensionless) phase-space distribution function and S i is the collision or source term for the species i = ϕ, ν s , both at the homogeneous background level.The latter quantity is defined by [13] 6 5 See footnote 3. 6 The collision term for both species is positive since it describes the probability of collision.For massless ghosts, which carry negative energies, one has generally that E ϕ = −p ϕ ≤ 0 with p ϕ ≥ 0, and hence f 0 ϕ ≥ 0, n ϕ ≥ 0, ρ ϕ ≤ 0, P ϕ ≤ 0 and w ϕ ≥ 0, from Eqs. (19) and (21).The negative pressure indicates that ghosts anti-gravitate.
where we used the same labelling as in Eq. ( 5) (i.e., ghost energies are replaced by their absolute values), and If the occupation number of the produced particles from the vacuum decay is small, namely f 0 i ≪ 1 (we will validate this assumption a posteriori ), then F 0 ≃ 1 and the collision term in Eq. ( 12) becomes independent of the phase-space distributions.This allows us to integrate Eq. ( 11) directly, obtaining [14,15] where t in , or interchangeably the scale factor a in , is the initial time of integration, and the factor a(t)/a(t ′ ) accounts for the redshifting of the momenta between t ′ and t.We take a in ≃ 10 −8 to be well deep inside the radiation era, when the vacuum decay is completely negligible; our results are insensitive to the precise value of a in .We will compute the collision terms first.The source term for the ghosts is relatively easy to evaluate, since the integrals over final-state particles have exactly the same form as for the cross section ϕϕ → ν s ν s of positive-energy massless ϕ particles.One can as usual perform this calculation in a Lorentz-invariant manner to obtain the cross section as a function of the Mandelstam variable s = (p µ 1 + p µ 2 ) 2 = 2p 1 p 2 (1 − cos θ), where θ is the angle between ⃗ p 1 and ⃗ p 2 , with The neutrino spectrum is slightly more complicated to calculate because the phase space of the ghosts is restricted by the momentum cutoff Λ; hence one cannot write a completely analogous expression involving the cross section for ν s ν s → ϕϕ.Instead, where c ij = pi • pj , and The integrals in Eqs. ( 15) and ( 17) can be evaluated numerically and they correspond to the instantaneous spectrum dΓ n /dp i once multiplied by the factor 4π p 2 i /(2π) 3 .We can use the computed collision terms to derive the phase-space distributions at the present time, which are given by Eq. ( 14) with t = t 0 or equivalently a = a 0 = 1 and where Here, Ω m , Ω r and Ω Λ are the relative abundances of total matter (baryonic plus cold dark matter), radiation (photons plus massless SM neutrinos) and dark energy described by a cosmological constant, respectively, and the contribution from the ghost plus ν s densities is Ω g (a) = (ρ ϕ (a) + ρ νs (a))/ρ crit with the critical density defined as usual, ρ crit = 3H 2 0 /(8πG).In practice, one can neglect Ω r since the ghost fluid is produced mainly at late times.In the following, we will assume for simplicity a flat universe, which is also the expectation after a period of early-time acceleration dubbed as inflation [16].In this way, we have Ω Λ + Ω m + Ω g = 1, where we have defined Ω g ≡ Ω g (a 0 ) for simplicity.
The distribution functions in ( 14) are normalized such that the particle number densities are given by7 The shapes of the spectra at the present epoch for the vector and scalar mediator models are shown in Fig. 3.
It is straightforward to show that the redshifted spectra f 0 i (t 0 , p i ) computed via Eq.( 14) are identical to those derived through Eq. ( 2).
It can be convenient to have an explicit formula for f 0 νs in the case where m νs = 0 and vector mediator, where x = p/Λ and A = 7.7 × 10 −4 (Λ/M i ) 4 .Although this is a numerical fit, the coefficient π in the exponent turns out to be accurate to three digits. 8

III.4. Phantom fluid evolution
Treating the ghosts and sterile neutrinos as cosmological fluids, one can compute their energy density, pressure and equation of state (EOS) at the background level as where f 0 i (t, p i ) is given by Eq. ( 14), and E ϕ ≤ 0. Fig. 4 shows the scale-factor evolution of the energy densities of ϕ and ν s for different values of the sterile neutrino mass m νs , in addition to their sum.As expected, in the case 8 For other choices of mν s we find similar fits by replacing x 2 e −πx 2 → (2/3)x 1.59 e −4 x 2.28 for mν s /Λ = 0.5 and → 0.349 x FIG. 4. Scale-factor evolution of the absolute value of the energy densities for ghosts ϕ (which have ρ ϕ < 0) and sterile neutrinos νs (Left), and of their sum (Right) for the vector mediator case.The units of ρi(a) are Λ 9 /(128π 6 M 4 j H0) (j = v, s).Solid lines are the result of numerically solving Eq. ( 21).The color shaded band represents the uncertainty on the back-reaction from Ωg(a) in the Friedmann equation (18) for the values of Γρ allowed by Type Ia supernova (see Section III.5).The green dotted lines are given by the analytic formulas in Eq. ( 31), derived from the Boltzmann equations in (22) when back-reaction is neglected.For both plots, we used Ωm = 0.32 and ΩΛ is determined by the flatness condition ΩΛ = 1 − Ωm − Ωg.The numerical results are shown with solid lines, while green dotted curves correspond to the numerical-fit functions (24).Right: Scale-factor evolution of wν s for different choices of mν s .The solid curves are for the vector mediator model, whereas the dashed curves are for the scalar model.The band thicknesses and Ωm and ΩΛ are as in Fig. 4. of massless ν s there is no net contribution of the two new species to the total energy density of the Universe, while for m νs > 0, there is a clear difference between |ρ i (a)| for the two components.
The energy densities shown in Fig. 4 can be derived in a simpler way, directly from the Boltzmann equation (11), after weighting by the energy E i and integrating over the particle three-momentum.This gives where Γ ρ is the rate of change in the particle energy density due to the vacuum decay, computed in appendix B (see Eq. (B12)) which gives approximately in terms of Γ n from Eq. (10).In Eq. ( 22), the equation of state for sterile neutrinos w νs depends upon m νs , and it could a priori also be a function of time.w νs ranges between zero for nonrelativistic particles when m νs ∼ Λ to 1/3 for relativistic ones, when m νs ∼ 0. This is confirmed by Fig. 5 (left), which shows the present-day w νs computed via Eq.( 21) as a function of m νs .We find that a good fitting formula is The ν s equation of state turns out to be timeindependent, to a good approximation, as displayed in Fig. 5 (right).It can be shown that the timeindependence is exact whenever the scale factor has a simple power law behavior, a(t) = t p , because in that case the integrals in Eq. ( 21) factorize to the form ρ i = F (t)ρ i (p) and P i = F (t) Pi (p) such that the timedependence cancels in their ratio.Since most of the ghost fluid production takes place while the universe is approximately matter dominated, this analytic behavior is realized by the full numerical solution.
Taking massless ghosts and massive neutrinos rather than vice-versa is an arbitrary model-building choice.The system of equations for the alternative choice is trivially related, by exchanging the roles of Ω ϕ and Ω νs in Eqs.(25), and relabeling w νs → w ϕ .The first choice results in Ω g > 0 while the second gives Ω g < 0. We will consider the general case in the following, which results in a sign choice ± for the solutions given below.By solving the system (25) numerically, one finds that the net production Ω g (y) is approximately linear in ϵ ρ , for ϵ ρ ≲ O (10).For larger ϵ ρ , linearity breaks down because the back-reaction of the ghost plus ν s fluid appears in Ĥ and reduces the production; however such large values are ruled out by observations of late-time cosmology, namely supernovae, as we will show.Hence we focus on the linear regime.Useful approximate fits to the numerical solutions are given by Ω g (y) ∼ = ±ϵ ρ f (w νs ) (1 + 0.109 ln y) factor (1 + 0.109 ln y) 3.09 should be interpreted as zero at small y where it would become negative.The choice of sign corresponds to taking massless ghosts and massive ν s (+) or vice versa (−).For small m νs , one can read from Eq. ( 24) that (1/3 − w νs ) ∼ = C m (m νs /Λ) 2 with C m ∼ = 4/3 for the scalar mediator and C m ∼ = 1 for the vector mediator.

III.5. Constraints from type Ia Supernovae
Because the net fluid of ghosts plus sterile neutrinos is produced at late times, a sensitive probe of its effect on cosmology is the use of supernovae as standard candles [17,18].Below, in Section IV.2, we will do a full analysis including other cosmological data sets.This preliminary study serves to define the parametric region of interest for the later Monte Carlo exploration.
The luminosity distance d L (z) of type Ia supernovae can probe deviations from the standard expansion history of the universe.It is given by where H(z) is given by Eq. ( 18) with a = (1 + z) −1 .It is related to the distance modulus µ(z) = 5 log 10 (d L (z)/10 pc), which is experimentally determined from the observed and modeled absolute luminosities of the supernovae.Given a covariance matrix C ij , the loglikelihood can be computed in the standard way, where δ⃗ µ is the vector of residuals between theory and observation.One minimizes the χ 2 with respect to the cosmological parameters to find the best-fit values and confidence intervals.This is known in the literature as frequentist profile likelihood, which is complementary to the Bayesian Markov Chain Monte Carlo (MCMC) since the latter identifies large regions in the parameter space that exhibit a good fit to the data.In the present case, the parameters over which to minimize the χ 2 are the Hubble rate H 0 , the present contribution to the energy density from the ghosts and sterile neutrinos Ω g , and that from the total matter contribution, Ω m .There may also be a nuisance parameter, the correction factor of the absolute magnitude of the supernovae M , which is degenerate with H 0 for uncalibrated supernovae.We assume flatness, so that the dark energy contribution is determined as Ω Λ = 1 − Ω m − Ω g .This procedure was followed for three different (but not all mutually exclusive) data sets: Union 2.1 [19], Pantheon [20], and Pantheon+SH0ES [21], comprising 580, 1048 and 1657 supernovae, respectively. 9In an alternative analysis of Pantheon+, the smaller sample of 1590 9 The Pantheon+ sample contains 1701 spectroscopically con-SNe is used, omitting the nearby SNe whose absolute distances are calibrated using Cepheid variables; this leaves H 0 and M undetermined.The determination of the other parameters gives the same result in both cases.
The results of minimizing the χ 2 for the different data firmed type Ia supernovae up to redshift z ≃ 2.3, but only 1657 of them have z > 0.01.We followed Ref. [21] in omitting supernovae with z < 0.01 from the analysis because of their large peculiar velocities, which preclude an accurate determination of the cosmological parameters.
sets are shown in Fig. 6 (left).There is no preference for Ω g ̸ = 0, but large magnitudes are allowed at 95% C.L., down to −0.75 for the case of massive ghosts and up to 0.45-0.9for massless ghosts, depending upon the data set.Fig. 6 (right) displays the preferred values of Ω Λ and Ω m as a function of Ω g , along with the dependence of the Hubble parameter H 0 on Ω g .Two values are shown, corresponding to the Union 2.1 and Pantheon+ samples (since Pantheon does not determine H 0 ).Both are linearly increasing functions of Ω g , but are offset from each other by about 4 km/s/Mpc.The Union 2.1 fit gives a minimum value of H 0 ∼ = 69.3km/s/Mpc at Ω g ∼ = −0.75, which is still in tension with the Planck determination 67.4 ± 0.5 [22] at nearly 4σ.In Section IV.2 we will show that this preliminary result agrees with a more careful analysis where all of the standard ΛCDM parameters are allowed to vary.The 95% C.L. bounds on Ω g in Fig. 6 restrict the allowed values of the dimensionless parameter ϵ ρ (defined below Eq. ( 25)) or its dimensionful counterpart Γ ρ .Fig. 7 (left) shows these limits for massless ghosts and massive sterile neutrinos as a function of the ν s equation-of-state parameter w νs .As already anticipated, values of ϵ ρ larger than O( 10) are excluded by supernova data, except when w νs ∼ 1/3.However, for such values of w νs the magnitude and shape of the energy densities of the new species become irrelevant since their contributions cancel, giving a negligible effect in the Friedmann equation (18).Therefore, the use of the approximate fits in Eq. ( 26), which assume linearity in ϵ ρ , is justified even for these larger values of ϵ ρ .
The minimum value of χ 2 for each data set, compared to the χ 2 value for the standard ΛCDM parameter values, can be characterized by ∆χ 2 = χ 2 min −χ 2 ΛCDM .It is given by −2810.4(Pantheon+), −1.1 (Pantheon), and −94.5 (Union 2.1).The discrepancy for Union and Pantheon+ arises from the dependence on H 0 in the χ 2 min , which is calibrated using Cepheids.The superficial appearance that the ghost model gives a better fit than ΛCDM is a reflection of the Hubble tension between the supernovae and the larger cosmological data sets.This tension is hidden for the Pantheon analysis since H 0 is marginalized out.We will confirm a mild preference for the ghost model over ΛCDM in Section IV.2.
An interesting question is whether Ω g can explain all of the dark energy, replacing the cosmological constant.Fig. 6 suggests that this is indeed possible, using the Union and Pantheon data sets, by taking Ω g ∼ 0.6.However it is in tension with the Pantheon+ upper bound Ω g < 0.45, which requires Ω Λ ≳ 0.2.In general, the capability of Ω g to function as dark energy can be understood through the effective EOS of the ghost fluid whose redshift dependence is plotted in Fig. 7 (right).Notably, it violates the weak energy condition since w eff < −1.3 both at early and late times.The early-time limit w eff → −1.5 can be derived analytically from the exact solution (31) that will be introduced below, in the regime where back-reaction of the ghost fluid is negligible.The late time limit w eff ∼ = −1.3depends weakly on the value of m νs /Λ.It is striking that the ghost fluid (in conjunction with the produced positive-energy particles) can function as dark energy even in the absence of a classical condensate, or indeed a potential V (ϕ), which was the original motivation for phantom fields in cosmology.

IV. CMB AND LARGE SCALE STRUCTURE
Even if their combined contributions to the energy density of the Universe are negligible, the creation of massless ghosts and positive-energy particles can leave an imprint on the CMB and matter power spectra.In particular, if both species are massless, their energy densities cancel exactly, but they can nevertheless alter the cosmological evolution and the formation of structure through their perturbations, depending on the production rate.
As noted above, the expansion history of the Universe is modified if the ghosts and ν s are not degenerate in mass.We estimate that the results for m νs > 0 will revert to those of the massless case when m νs ≲ 2×10 −5 Λ.This was found by comparing the effects on the CMB and matter power spectra, including the linear order perturbations, and either including or neglecting the zerothorder contributions of the new species.In the following sections we will study the two qualitatively distinct regimes, depending on the scale of m νs .

IV.1. Massless νs case
If ν s is sufficiently light, its contribution to the total energy density of the Universe cancels that of the presumed massless ghosts and there is no appreciable effect at the background level.The first nontrivial contribution arises at linear order in the cosmological perturbations.We derive the relevant equations in appendix B and find that the perturbations for ϕ and ν s do not cancel each other, mainly because of the different phase-space distribution between the two species, shown in Fig. 3; see also Eqs. (B29) and (B30).
The perturbations were implemented in a modified version of the Boltzmann code CAMB [23].To the six parameters comprising the ΛCDM model,10 we add two new ones: Γ ρ , which controls the amplitude of the background energy densities of the dark species, and Γ ρ /Λ 4 , which determines the strength of the first-order perturbed collision term.We assume vanishing initial conditions for the ϕ and ν s perturbations in the early radiation era (a ≪ 10 −5 ).
For the background energy densities, the Boltzmann equations ( 22) can be rewritten in terms of the scale factor a and solved analytically, since the net density of the new species vanishes in the Friedmann equation (18).The solutions are The effects of the two new species on the CMB temperature and linear matter power spectra are shown in Fig. 8.The main modification on the former is an overall decrease of the late Integrated Sachs-Wolfe (ISW) effect, visible at large angular scales (small ℓ), which is due to the increase of the gravitational potential along the path of photons between the time of last scattering and today.This makes the photons lose energy as they travel through gravitational potential wells that are growing in time.
More precisely, the late ISW contribution to the CMB anisotropy in Fourier space is given by [29] where τ * and τ 0 are the conformal times at last scattering and today, respectively, Ψ and Φ are the potential and curvature in conformal Newtonian gauge [30], the prime indicates the conformal time derivative, and j ℓ are the spherical Bessel functions.The Newtonian curvature Φ and potential Ψ for the mode k are set by the total density contrast δ tot , velocity θ tot and shear σ tot perturbations [31], 11 Eqs.(31) are valid whenever back-reaction of the new species on the Hubble-parameter evolution can be neglected, and they apply for any value of mν s .One should keep in mind that Γρ is generally a function of mν s (see Eqs. ( 23) and ( 10)), which in turn affects the value of wν s as described by Eq. (24).
where H = aH(a) is the conformal Hubble parameter.The terms in the square brackets in Eq. ( 33) are gauge invariant and hence they can be evaluated in the frame comoving with the cold DM (CDM), corresponding to synchronous gauge.In a flat ΛCDM universe, the late ISW contribution is zero during matter domination because both Ψ and Φ remain constant, but as dark energy begins to dominate, the curvature Φ decays, reducing the magnitude of the potential Ψ.This effect is enhanced if (i) there is extra radiation pressure, (ii) the Universe expands faster than in ΛCDM, or (iii) a nonzero spatial curvature is present [32][33][34].
In the present case, one can understand the decrease of the late ISW effect directly from Eqs. (33).Fig. 9 (left) shows the large-scale evolution of the synchronousgauge physical perturbations for sterile neutrinos ν s and ghosts ϕ, which can be treated as a single effective fluid since both species are relativistic.The ghost nature as a negative-energy particle manifests itself through the opposite sign in the physical perturbations (ρ ϕ F ϕ,ℓ ) with respect to those for ν s , having negative density contrast and velocity, but positive shear stress perturbations, which is due to ρ ϕ ≲ 0. Positive physical shear stress is a generic feature of phantom cosmological species and the main effect is to drive clustering [35].This is exactly the opposite of what happens for normal species like SM neutrinos, where the presence of the shear stress leads to erasure of fluctuations that might initially be present in the fluid, suppressing further growth.
Due to the larger phase space available for ν s , 12 which increases its first-order perturbed collision term, the relative perturbations for ν s are expected to dominate over those for ϕ at high hierarchy multipoles (see discussion in section (a) of Appendix B).Since the relative shear σ and higher hierarchy multipoles for a single relativistic species are generally negative at large scales, as occurs for SM massless neutrinos and photons, we expect σ ϕ to be more negative than σ νs .As a consequence, the physical shear for the effective ν s + ϕ fluid, defined by ρ νs σ νs + ρ ϕ σ ϕ , becomes positive, as borne out in Fig. 9 (left) with the dot-dashed green curve.The effect of the new species is therefore to decrease the value of the potential Ψ in Eq. ( 33), favoring matter clustering.This is shown with the red dashed curve in Fig. 9 (right).
On the other hand, the density contrast δ and the velocity perturbation θ for the effective single fluid turn out to be positive and negative, respectively, which makes it harder to predict a priori the net contribution of the new species to Φ.The evolution of the latter is shown with the green dot-dashed curve in Fig. 9 (right).The net effect of the ν s + ϕ fluid is to reduce Ψ − Φ, which generally leads to smaller C ISW ℓ in Eq. ( 32) than in ΛCDM.A similar trend occurs in phantom dark energy models with non-negligible anisotropic stress [36].
For larger values of Γ ρ and Γ ρ /Λ 4 , at the largest angular scales, the shear perturbation σ increases faster than the lower multipole perturbations δ and θ, driving (Ψ − Φ) to smaller values on shorter time scales.This results in a larger contribution to the late ISW effect in the CMB temperature power spectrum because of the square in Eq. (32), which is visible in Fig. 8 (left) around ℓ ∼ 2 − 3. The divergent nature of the perturbations for large values of Γ ρ /Λ 4 or Γ ρ is not surprising, since theories with ghosts are generically subject to instabilities.The growth of the instability is controlled by the ghost momentum cutoff Λ and the mediator mass M i , and the analysis carried out below will exclude parameter values Red and blue regions correspond to the 95% C.L. cosmological bounds for the vector and scalar mediator models, respectively.Violet indicates the gamma-ray excluded region, Eq. ( 3). Green region is inconsistent with the technical assumption f 0 i ≲ 0.1 in Eq. ( 13).
leading to excessive growth. 13We note that the effect of weak gravitational lensing on the CMB power spectra (e.g., Ref. [37]), in the presence of ghosts and sterile neutrinos, does not deviate appreciably from the ΛCDM scenario because the Weyl potential (Ψ − Φ)/2 in the two cases is essentially identical at small angular scales (large ℓ).
The modifications to the matter power spectrum in our scenario are also negligible for allowed values of Γ ρ and Γ ρ /Λ 4 , and the main effect is a small rise of the spectrum at large scales, as shown in Fig. 8 (right).The small positive contribution to P (k) is due to its dependence on the ratio between the matter density perturbation δρ m and the background matter density ρ m .In fact, neither ν s or ϕ contribute directly to δρ m nor ρ m , except through the metric perturbations appearing in the equation governing the evolution of δρ m .The increase in the matter power spectrum signals that the ghost and and sterile neutrino effective fluid favors a mildly more efficient matter clustering than ΛCDM at large angular scales, which was already noticed from the CMB power spectrum. 13Conceivably, higher order perturbations could cut off the instabilities.In this work we constrain the effect of the first order perturbations, such that the higher-order contributions, proportional to higher powers of Γρ/Λ 4 ≪ 1, are negligible.

IV.1.1. Method and Results
To constrain the new parameters Γ ρ and Γ ρ /Λ 4 , in conjunction with the ΛCDM ones, we ran the publicly available MCMC code CosmoMC [38].Convergence of the MCMC chains was monitored through the Gelman-Rubin parameter R − 1 [39], by requiring R − 1 ≲ 0.01.The cosmological data sets we considered are: • Planck: the full range of the CMB measurements from Planck 2018, which include the temperature and polarization anisotropies, as well as their crosscorrelations [22]; • Lensing: the 2018 Planck measurements of the CMB lensing potential power spectrum, reconstructed from the CMB temperature four-point function [40]; • BAO: the distance measurements of the baryon acoustic oscillations given by the 6dFGS [41], SDSS-MGS [42], and BOSS DR12 [43] surveys; • DES: the first-year of the Dark Energy Survey galaxy clustering and cosmic shear measurements [28,44,45]; • Pantheon: distance measurements of (uncalibrated) type Ia supernovae from the Pantheon sample, comprising 1048 data points in the redshift range We did not include the Halofit modelling of the nonlinear part of the matter power spectrum because additional modifications of the standard nonlinear clustering model might be needed.Since this is beyond the scope of the present work, we consider only linear observables in our analysis.Fig. 10 shows the constraints on Γ ρ and Γ ρ /Λ 4 , for the vector and scalar mediator models, as derived from the CosmoMC analysis, in addition to the COMPTEL bound of Eq. (3).Regarding the cosmological limits, one observes two features: the models in agreement with the data lie in the region of Γ ρ /Λ 4 ≲ 10 −38 MeV, and for sufficiently small Γ ρ /Λ 4 , the dependence on Γ ρ disappears.We qualitatively explain such features in appendix B, since they are closely related to the form of the perturbation equations.The constraints on the scalar mediator model are slightly stronger than those of the vector model, due to differences in the perturbation equations for the two interactions, which arise from their energy dependence; see, e.g., Eq. ( 9).
The constraints on Γ ρ /Λ 4 versus Γ ρ can be mapped onto the plane of the ghost momentum cutoff Λ versus the mediator mass scale M i , using Eq.(23).The resulting limits are shown in Fig. 11 (left) and they are essentially independent of the mediator type.Also shown is the consistency condition for neglecting the Pauli-blocking and Bose-Einstein stimulated-emission factors in the collision term of Eq. ( 12).This is inferred by demanding  25 GeV for both the vector (i = v) and scalar (i = s) mediator cases.The 95% C.L. bounds on M i versus Λ can be approximately expressed by where Λ ≡ Λ/MeV ∈ (10 −7 , 10 2 ).We found that the two new parameters in the massless ν s case leave the standard six ΛCDM parameters almost unchanged compared to their Planck 2018 best-fit values, without ameliorating any current cosmological tensions, such as the H 0 and S 8 tensions [46,47].This will be considered in more detail in Section IV.3 below.The correlations between these parameters are shown in appendix C. Based on the foregoing discussion, the constraints on Γ ρ and Γ ρ /Λ 4 are expected to be dominated by the CMB power spectra rather than large-scale-structure observables.We verified this by excluding BAO, DES and Pantheon data sets in a separate CosmoMC run, finding that the results in Fig. 11 (left) were only slightly modified.
From the bounds in Fig. 10, one can estimate how large the abundance of massless sterile neutrinos can be today, while remaining consistent with cosmological data.This is shown in Fig. 11 (right), displaying Ω νs versus the ghost momentum cutoff Λ for both mediator models.The surprisingly large values of Ω νs ≫ 1 are viable because of their exact cancellation by Ω ϕ , so that the constraints arise only at the level of the perturbations.The maximum values of Ω νs are reached by saturating the gamma-ray constraint, Eq. ( 3), leading to for Λ ∼ 19 MeV.Despite the large numbers, these ν s gases are still dilute, by many orders of magnitude, relative to a degenerate Fermi gas with Fermi energy Λ.
Thus the occupation number f νs (p) ≪ 1, even for this extreme case.

IV.2. Massive νs case
In section III.5, we made a preliminary study of the effects of producing massive neutrinos plus massless ghosts, comparing to supernova data but not to the other cosmological data sets.Here we undertake the more complete analysis, implementing the energy-density evolution of ϕ and ν s at the background level, given by Eq. ( 26), in a modified version of the CAMB code.It is embedded in CosmoMC to scan the full parameter space.
In addition to the six parameters of ΛCDM, we considered Ω g ≡ Ω νs (a 0 ) + Ω ϕ (a 0 ), the net contribution of the ghost plus ν s densities today, and the equation of state for sterile neutrinos w νs as the independent new-physics parameters, even though they are related to each other via Eqs.(10,(23)(24).This allows for simultaneous study of the scalar and vector mediator models, which differ only by how Ω g and w νs are related to the microphysical parameters m νs , Λ, M i .
The linear perturbations for the two new dark species were included in an approximate form, given by Eqs.(B38) in appendix B. A more exact form for the perturbations for massive sterile neutrinos is technically challenging to derive, and is left for a future study.We expect that the bounds on Ω g and w νs derived here could only be moderately strengthened by such a study, since the main effect is through the background evolution rather than the perturbations, as discussed in ap- 12. Like Fig. 8, but for mν s > 0. We assumed Ωg ≃ 0.6, wν s ≃ 0.02, ΩΛ ≃ 0.1, and spatial flatness.

pendix B.
Figure 12 shows the impact of massless ghosts and massive sterile neutrinos on the CMB temperature anisotropies and on the matter power spectrum, for Ω g = 0.6 and w νs = 0.02.To highlight the effects of the new parameters, those of ΛCDM are fixed at their Planck 2018 best-fit values.(Note that Ω Λ is derived from the flatness assumption and is not a fundamental parameter in ΛCDM.)The features in Figure 12 differ qualitatively from those where m νs = 0, Fig. 8.In the CMB spectrum, there are two main effects: an overall shift of the CMB acoustic peaks toward smaller angular scales (larger ℓ) and an increase of the late ISW contribution at large scales.

IV.2.1. Acoustic peak shift
The peak shift comes from the decrease of the darkenergy abundance Ω Λ [34] (a derived parameter in ΛCDM), which is required to preserve spatial flatness.If one simultaneously keeps the other ΛCDM parameters fixed to their Planck best-fit values, 14 the Hubble rate H(a) in Eq. ( 18) decreases relative to the standard scenario because the density of the new species becomes significant only at late times, as shown in Fig. 13 (top).Since the age of the Universe is computed as t = d ln a/H(a), reducing H(a) results in an older Universe, and the thermal history, including photon decoupling, has to occur earlier in the past relative to ΛCDM.This shifts the acoustic CMB peaks to smaller scales, as seen in Fig. 12 (left). 14In our model, the preferred values of H 0 and Ωm change by relatively smaller amounts than Ω Λ and have a subdominant effect on the peak shift.In detail, the location of the first acoustic peak is given by ℓ s ≃ π/θ * s , where [34] θ is the angular size of the sound horizon at the time of last scattering, is the sound horizon with c s (z) the sound speed of the photon-baryon fluid, and is the comoving angular diameter distance between the last-scattering surface at redshift z * and today.A decrease of Ω Λ does not affect the sound horizon r * s , since it depends only on the physics before the time of photon decoupling, when the dark energy and the new species were negligible, but it does affect D * A by reducing H(z).Therefore, the comoving angular diameter distance increases with respect to ΛCDM, reducing the value of θ * s and shifting the location of the first CMB peak to larger ℓ.

IV.2.2. Late ISW enhancement
The second feature in the CMB temperature anisotropy spectrum, when nonrelativistic ν s are included, is an increase in the late ISW contribution, which can be explained by the interplay of two competing phenomena: the behaviour of the cosmological perturbations of ν s and ϕ, and the decrease of Ω Λ .
Fig. 14 (left) shows the large-scale evolution of the synchronous-gauge physical perturbations for the ghosts and sterile neutrinos.Similarly to the massless ν s case, the physical shear perturbation ρ ϕ σ ϕ for ghosts (blue dot-dashed line) dominates over that for sterile neutrinos (red dot-dashed curve), leading to an overall positive contribution to the total shear (green dot-dashed line) and a consequent decrease of the Newtonian potential Ψ, according to Eq. (33).The impact on the Newtonian curvature Φ from the perturbations of the new species is small, because the positive contribution of the density contrast δ, which is dominated by ν s , is largely cancelled by the negative contribution of the velocity perturbation θ, dominated by ϕ, as shown in Fig. 14 (left).The net effect on (Ψ − Φ), whose conformal time derivative enters the late ISW contribution in Eq. (32), is to drive it to more negative values on a shorter timescale, leading to an increase in the late ISW effect, as shown in Fig. 14 (right).
The situation is complicated by the decrease of Ω Λ in the background density evolution, which tends to counteract the late ISW effect.The dark energy impacts the decay of the gravitational wells which photons traverse.In ΛCDM, the Newtonian curvature Φ and potential Ψ, which parameterize the depth of the gravitational wells, remain constant during matter domination, but they start to decay when dark energy takes over, illustrated by the thin curves in Fig. 14 (right).This is due to the gravitational wells becoming shallower when the Universe starts accelerating, enhancing the late ISW effect.Conversely, the effect is diminished if Ω Λ is reduced [34], attenuating the increase of the late ISW contribution arising from the ϕ and ν s perturbations.
Which of these competing effects wins is determined by the relative values of Ω Λ and Ω g , due to the ν s contribution to the latter.If Ω Λ ≫ Ω g , the ϕ + ν s contribution becomes negligible and the late ISW effect is reduced.In general, we see that for relativistic sterile neutrinos, with w νs ≳ 0.1, the late ISW effect in the CMB temperature anisotropy spectrum decreases, similarly to the massless ν s case.

IV.2.3. Spatial curvature
The previous discussion assumes that the Universe is spatially flat, which forces Ω Λ to decrease with respect to the Planck prediction, whenever massive sterile neutrinos are produced.The presence of a positive curvature (Ω k < 0) might reduce the impact of the new species on the CMB spectrum since a nonzero curvature provides opposite effects to those seen here [34].A closed Universe could also reduce some discrepancies between the ΛCDM model and the Planck 2018 data [48], although it is still unclear whether they might be due to internal systematics [49].We leave the study of phantom fluids in a non-flat Universe for a future work.

IV.2.4. Matter power spectrum
As concerns the matter power spectrum, shown in Fig. 12 (right), the vacuum decay into massless ϕ and nonrelativistic massive ν s leads to a suppression at all scales.The dominant effect is the free-streaming of the new species, which occurs immediately after they are produced, causing a damping of their density fluctuations compared to standard nonrelativistic species, such as CDM.This is particularly important for sterile neutrinos, since they contribute to the total nonrelativistic matter and hence to the matter power spectrum.The density contrast for ν s is smaller than that for CDM, δ νs < δ cdm , due to the collision term (see appendix B), whereas its energy density ρ νs can be of the same order as ρ m , but in any case it leads to δρ νs ≲ δρ cdm .Since the matter power spectrum is defined in terms of δρ m /ρ m to which both ν s and CDM contribute, the net effect is a decrease of P (k) relative to ΛCDM.This differs from massive SM neutrinos, where the suppression occurs only at scales smaller than the free-streaming length set by the time of relativistic-to-nonrelativistic transition [50].Unlike SM neutrinos, the equation of state w νs for ν s from vacuum decay is approximately time-independent, depending only upon the value of m νs .
However, the free-streaming of ν s is not the only effect on P (k).Sterile neutrinos and ghosts contribute to the Friedmann equation and their abundance increases with time.If the ΛCDM parameters and H 0 are fixed at their Planck values, increasing Ω g requires reducing Ω Λ to keep the total dark energy unchanged.The net effect of the time-evolution of Ω g (a) + Ω Λ on the Friedmann equation is to lower the value of H(a) (Fig. 13, upper plot) and to delay the transition from matter-domination to the epoch when dark energy, including the ϕ + ν s fluid, started to dominate (Fig. 13, bottom plot).This allowed more matter to fall into gravitational wells and thereby enhance the matter power spectrum.
Similarly, the perturbations of the new species, and in particular the nonzero shear for ν s and ϕ, tend to reduce the decay of the Newtonian gravitational potential Ψ and curvature Φ, leading to an increase of the matter clustering.This is a direct gravitational back-reaction effect at the level of perturbations, which adds to the background effect discussed above, leading to an increase of P (k) primarily affecting large angular scales, as in the case of massless ν s .
For nonrelativistic massive ν s , the free-streaming effect dominates over the background and back-reaction ones, as shown in Fig. 12 (right), producing less matter clustering.On the other hand, the latter two effects dominate in the case of relativistic or mildly relativistic massive ν s , because such sterile neutrinos do not contribute to the total nonrelativistic matter and hence do not directly affect the matter power spectrum.

IV.2.5. Results
Similarly to the massless ν s case, the MCMC code CosmoMC was used to constrain the two new parameters Ω g and w νs , and to study their correlations with the ΛCDM ones.We adopted the same convergence requirement and data sets as used for the m νs = 0 case, as described in subsection IV.1.1.
Fig. 15 shows the correlations between the new parameters and those derived from ΛCDM that are significantly correlated.The latter include the total matter abun-dance Ω m , which is separate from the ν s contribution, the dark energy abundance Ω Λ , the Hubble constant H 0 , and the clustering amplitude parameter S 8 = σ 8 Ω m /0.3.The extended version of the correlation parameter plot shown in Fig. 15 is discussed in appendix C, including a table with the 68% C.L. limits.
Increasing Ω g , defined as the present-day energy fractions of ϕ plus ν s , leads to an increase in H 0 , which can be intuitively understood by the phantom nature of the effective ϕ + ν s fluid.Its equation of state w eff was found to be < −1.3 (Fig. 7), more negative than that of the cosmological constant.Therefore, increasing the relative abundance of the phantom fluid causes the Universe to expand faster today than in the ΛCDM model, leading to a larger value of H 0 .
Such accelerated expansion entails weaker growth of large scale structure, which is parameterized by S 8 = σ 8 Ω m /0.3, where σ 8 is the linear-theory standard deviation of matter density fluctuations in a sphere of radius 8 h −1 Mpc.Namely, Fig. 15 shows a slow decrease in the parameter S 8 with increasing H 0 , which is most pronounced at small values of w νs .This confirms the observation made in subsection IV.2.4,where the decrease of the amplitude of the matter power spectrum (and hence σ 8 and S 8 ) was observed for nonrelativistic massive ν s .We will quantitatively discuss the implications of those findings for H 0 and S 8 in the next subsection, devoted to possible resolutions of known cosmological tensions.
Another striking correlation is the decrease of both Ω Λ and Ω m as the present-day abundance of the new species increases.The former effect derives from the flatness condition Ω Λ + Ω m + Ω g = 1.Although Ω m is also a derived parameter in ΛCDM, it is determined by Ω b h 2 and Ω c h 2 , which are constrained by the relative heights of the CMB peaks [51].Since our model affects only the late-time evolution and hence does not alter photon decoupling, whose physics impacts the peak height, the value of (Ω b +Ω c )h 2 remains nearly unchanged relative to ΛCDM.This is verified in the extended correlation plot in appendix C. Consequently, the decrease in Ω m with increasing Ω g comes mostly from the growth of H 0 , as explained above.
For constraining Ω g and w νs , which are taken as independent input parameters, the second-from-top left plot of Fig. 15 displays their allowed values.It shows that Ω g is roughly independent of w νs , which can be understood from the definition of Ω g as the net amount of ρ νs + ρ ϕ that has been produced, whose value already incorporates the effect of w νs .In the limit of m νs → 0, w νs → 1/3, with Γ ρ fixed, it is true that Ω g would be driven to zero, but in the Ω g -w νs parametrization, Γ ρ is not fixed and the former are independent of each other.In fact, w νs does not correlate strongly with the other ΛCDM parameters either.Its main effect is on the linear perturbations.At the more important level of the homogeneous energy densities, w νs has a relatively small impact.
The 68% and 95% C.L. exclusion plots for Ω g and w νs corresponding to Fig. 15 are shown in Fig. 16 (left).We also map them onto the plane of w νs and ϵ ρ = Γ ρ /(3H 0 ρ crit ), or directly Γ ρ , using Eq. ( 26) and Fig. 15.The resulting limits, shown in Fig. 16 (right), resemble those obtained from the supernova analysis of section III.5, Fig. 7 (left).As expected, large values of Γ ρ are allowed as w νs → 1/3, since Ω g can remain fixed in (the vacuum decay rate) and mν s /Λ, for vector (red) and scalar (blue) mediator models.Thick (thin) curves show the 68% (95%) C.L. exclusion limits.Right: 95% C.L. bounds on the parameters Mi (the mediator mass scale, with i = s, v), versus the ghost momentum cutoff Λ, corresponding to the left plot.Solid lines correspond to mν s /Λ = 10 −2 , and dashed to the values that give the strongest constraints, namely mν s /Λ ≃ 0.53 (0.36) for the vector (scalar) mediator case.that limit.
Using Eqs.(A6), (23), and (24), one can constrain Γ 0 ρ ≡ Γ ρ | mν s =0 versus m νs /Λ.This is shown in Fig. 17 (left) for the vector and scalar mediator models.It does not correspond to independent bounds on the microphysical parameters Λ and M i , since cosmological observables are sensitive to the derived combination Γ 0 ρ .However we can use Eqs.( 23) and (A6), together with the bounds on Γ 0 ρ , to limit M i as a function of Λ, for fixed values of m νs /Λ, as shown in Fig. 17  Analytically, the new physics mass scales are bounded from below as for neutrino masses m νs /Λ ≳ 0.01.The strongest limits occur for particular values of m νs /Λ ≃ 0.53 (0.36) for the vector (scalar) mediator respectively.As expected, these bounds are much stronger than those in the massless ν s scenario, (section IV.1.1),resulting mainly from modifications at the background level, as opposed to perturbations.They depend weakly on m νs /Λ, since Γ 0 ρ in Fig. 17 (left) changes by at most few orders of magnitude from varying m νs .

IV.3. Cosmological tensions
In recent years, discrepancies have emerged among various observational probes regarding the values of certain cosmological parameters.The most notable one is the Hubble tension, which refers to a mismatch between early-time and late-time estimations of the Hubble constant H 0 [47,[52][53][54].This tension has crossed the critical 5σ level between the two most precise measurements: Planck finds H 0 = 67.66 ± 0.42 km/s/Mpc [22] from fitting the CMB data within the ΛCDM framework, whereas SH0ES finds H 0 = 73.04 ± 1.04 km/s/Mpc [55] based on the distance-ladder approach. 15A milder tension also exists concerning the large-scale structure of the Universe [44,[57][58][59][60][61][62][63].Estimations of the amplitude of late-time matter fluctuations S 8 from galaxy redshift surveys are consistently lower than those inferred from CMB data, which points to S 8 = 0.825 ± 0.011 [22], assuming the ΛCDM model.The most recent late-time measurement S 8 = 0.790 +0.018 −0.014 , from the DES Y3+KIDS-1000 surveys [64], agrees with the Planck value within the ∼ 2σ level.
This motivates us to reconsider our main results for the parameter distributions, Fig. 15, with respect to the cosmological tensions.In particular, we want to determine whether there exist models within the ΛCDMϕν s scenario that can ameliorate the tensions, while remaining compatible with all the cosmological data.Fig. 18 shows the same correlation plot between the parameters as in Fig. 15, restricted to those models that resolve the H 0 and S 8 tensions within the 1σ level.Interestingly, both tensions can be resolved in a narrow region of the parameter space with w νs ∼ 0.3 and Ω g ∼ 0.3 − 0.5, forcing both Ω Λ and Ω m to decrease from their Planck 2018 best-fit values.In our MCMC chains, only about 0.015% of the models can address both tensions, while keeping a total χ 2 comparable to that of the ΛCDM model.
However, this does not constitute a successful resolution of the H 0 + S 8 tensions, because it comes at the expense of creating new ones.In Table II

V. CONCLUSIONS
In this work, we have reconsidered the cosmological consequences of a phantom field, a scalar with the wrong-sign kinetic term.Such a theory must be a lowmomentum effective description, valid below some cutoff Λ.We refined a previous estimated upper limit to show that Λ ≲ 19 MeV, from observational constraints on diffuse gamma rays.At higher scales, phantom models must revert to a theory with normal-sign kinetic terms.
Unlike the great majority of phantom dark energy studies, we assumed the vacuum expectation value of the field vanishes, so that its potential makes no contribution to the Hubble expansion.Instead, the spontaneous breakdown of the vacuum into phantom particles plus normal ones produces a fluid that behaves as dark energy, as long as there is a mismatch between the masses of the two kinds of particles.We showed that the equation of state of this exotic "ghost fluid" always violates the null energy condition (NEC), since it lies in the interval −1.3 ≲ w eff < −1.5 at all times (see Fig. 7, right).Therefore, even if at the present time the net contribution of the ghost fluid to the energy budget of the Universe is small, Ω g ≪ 1, eventually it will come to dominate over the cosmological constant and lead to a "big rip" singularity in the future [67].We found that the NEC-violating equation of state is at odds with combined cosmological data if one tries to make the ghost fluid account for all of the dark energy.Instead, there is an upper bound Ω g ≲ 0.2, requiring it to be subdominant to the contribution from cosmological constant, 0.5 ≲ Ω Λ ≤ 0.7.
Such large values of Ω g cannot be attained if the phantoms couple with only gravitational strength to other particles.Instead, we allowed for stronger interactions, mediating the production of phantoms plus light dark matter (with mass < Λ).It leads to a larger Hubble rate at late times, in the case of massless phantoms, since the negative ghost energy redshifts faster than the posi-tive dark matter density.This constant rate of creation of energy density is reminiscent of the old steady-state cosmology [68,69], with the difference that the present mechanism is rooted in quantum field theory, rather than an ad hoc assumption.Naively, one might think it could resolve the Hubble tension, but we find that this exacerbates the S 8 tension, as expected on general grounds for late-time mechanisms [70,71] (see however Ref. [72]).The best-fit phantom model does soften the H 0 tension, reducing it to ∼ 3σ, and it provides a mild improvement relative to ΛCDM in fitting the data sets.
A special case is where massless phantoms couple to massless dark radiation with strength exceeding that of gravity.Then their energy densities cancel exactly, and the magnitude of each can exceed usual bounds on a new light species by many orders of magnitude, since there is no impact on the Hubble expansion.Nevertheless, linearorder cosmological perturbations grow in this scenario, and provide constraints via the ISW effect.In the extreme case where the cutoff Λ saturates the gamma-ray bound, and the new-physics scale coupling ghosts to dark radiation is M ∼ 10 10 GeV, the dark radiation could contribute as much as Ω νs ∼ 10 26 (see Fig. 10), being cancelled by the negative contribution of the ghosts.
In a companion paper, we considered a possible observable effect of such a large density of dark radiation, whose energy is of order the cutoff Λ: in the presence of an additional interaction coupling it to electrons, it could be seen in direct detection experiments [73].It would be interesting to identify other potentially observable consequences of such a large population of energetic particles in the late Universe.For example, there could be a mismatch between the preferred reference frame of the ghost fluid and that of the CMB, which for simplicity we neglected in the present study. where The two Heaviside functions restrict the ghost momenta to lie below the cutoff.The sterile neutrinos (with momenta p 3 , p 4 ) are assumed to be Dirac fermions and the ghosts (momenta p 1 , p 2 ) are complex scalars.For real scalar ghosts, an extra factor of 1/2 should be included.Using rotational invariance to make ⃗ p 1 parallel to ẑ and ⃗ p 2 to lie in the x-ẑ plane, the measure becomes where c jk ≡ cos θ jk for j, k = 1, 2, 3, and we used Using the remaining delta function to solve for φ 3 gives δ(. . . where c 23 = c 12 c 13 + s 12 s 13 cos φ 3 , with s jk ≡ sin θ jk .Here δ(g(x)) = xs δ(x − x s )/|g ′ (x s )|, where x s are the roots of g(x).In this case the single root is with g ′ (cos φ 3 ) = (p 2 p 3 /E 4 )s 12 s 13 .Using dc jk = −s jk dθ jk , the decay rate integrals reduce to where x j ≡ p i /Λ, y j ≡ E j /Λ and there is a factor of 2 for the two values of φ 3 satisfying Eq. (A5).The upper limit on E 3 comes from energy conservation with p 4 = 0, the maximum energy allowed by kinematics.
For the scalar mediator, |M n | 2 in Eq. ( 9) simplifies to and the result for the corresponding decay rate is Similarly, for the vector mediator, the matrix element (A6) is leading to a vacuum decay rate density contrast.The phase-space distribution evolves according to the Boltzmann equation, which is generically given in terms of the comoving momentum by where ε i = q 2 i + a 2 m 2 i = aE i is the comoving energy and we neglected the term proportional to the derivative of f i with respect to the direction qi , it being a second-order quantity.Putting Eq. (B1) into Eq.(B2) gives at the background level and at the linear-perturbation level C is the perturbed collision term.On the right-hand side of Eq. (B4) we used Eq.(B3) to replace the time-derivative of f 0 i with the zeroth-order collision term.Eq. (B3) reduces to Eq. ( 11) in terms of the physical time t and momentum p i , and identifying It is common practice to decompose Ψ i into Legendre polynomials P ℓ ( k • qi ), motivated by left-hand side of Eq. (B4) depending only on q i , k = | ⃗ k| and k • qi in Fourier space.Hence [31] Ψ This allows the Fourier version of Eq. (B4) to be expressed as a hierarchy of equations for the multipole moments Ψ i,ℓ (k, q i , τ ) [31,77] where the collision-term multipoles are [77] C and C,ℓ with Ω ⃗ k the solid angle subtended by the wave vector ⃗ k.

(B9)
Expanding f i as in Eq. (B1) and neglecting the Pauli-blocking and Bose-Einstein stimulated-emission factors, the background contribution and the linear-order perturbation of the collision term are respectively and C,ℓ where the dependence on ( ⃗ k, τ ) is suppressed and the symmetry q 1 ↔ q 2 was used to combine two terms in the last line.Similar expressions can be derived for the ghosts ϕ.Subleading contributions were kept in Eq. (B10), arising from the terms f 0 i ≪ 1, that will be relevant in the collision multipoles (B7).Integrating Eq. (B10), weighted by the sterile neutrino energy ε 3 , and neglecting subleading contributions gives the numerical fit 16   Γ ρ ≡ 1 a which appears in Eq. (22).The same result arises from the integral of C weighted by the ghost energy q 1 .The computation of the perturbed collision term in Eq. (B11) in the case of massive sterile neutrinos is nontrivial, requiring not only the integral of Ψ i,ℓ (defined in Eq. (B5)) and solution of the multipole hierarchy (B6), but also depending on two free parameters: the Λ and am νs .It simplifies in the relativistic limit, q ≫ am νs , approximated by m νs ∼ 0 in the preceding expressions.This is the most interesting case for the perturbation analysis, since then the background contributions to the energy density from ghosts and sterile neutrinos cancel and modifications of the CMB and matter power spectrum can arise only through the perturbations.We consider this scenario in detail in the following subsection, B.1, and derive an adequate approximate description for the more general case in subsection B.2.

Relativistic sterile neutrinos
In the relativistic limit, it is a common simplification to integrate over ⃗ q i in (B6), weighted by the energy and the background phase-space distribution [31].The equations can then be expressed in terms of the relative perturbation variables: the density contrast δ i , the divergence of the velocity field θ i and the anisotropic stress σ i , where i = ϕ, ν s .They are defined in terms of F i,ℓ ≡ dq q 3 f 0 i (q)Ψ i,ℓ (q) dq q 3 f 0 i (q) , (B13) 16 See footnote 3 by δ i = F i,0 , θ i = 3kF i,1 /4 and σ i = F i,2 /2.Then the perturbation equations (B6) reduce to [31,77] with collision multipoles C,ℓ The numerators can be evaluated by integrating Eqs.(B10) and (B11).In terms of x j ≡ q j /Λ, the collision terms for ν s are C,ℓ with cos φ 2 and cos φ 1 given by cos and c jk = cos θ jk , s jk = sin θ jk .Analogously, the collision terms for ghosts (with overall − signs from their negative energy) are  III.Numerical coefficients for the polynomials in Eq. (B24) used to fit the kernels in Eq. ( B18) and (B22).Kernels with higher ℓ than those displayed are negligible.Kernels with subscripts i = 1, 2 refer to ghost quantities, whereas those with subscripts i = 3, 4 refer to sterile neutrinos.K3 and W1 do not have the ℓ subscript because they refer to the background contribution. where and The kernels in Eq. (B18) and (B22) are fitted by polynomials with coefficients a j given in table III for the K and W kernels.The kernels associated with ghost particles (subscripts i = 1, 2) have only the two leading coefficients nonvanishing; therefore, the integrals over ghost momenta in Eqs.(B17) and (B21) take the form the first term of which is proportional to F ϕ,ℓ , while the second one is a higher-order correction.To estimate the latter, we assume Ψ ϕ,ℓ (x) depends weakly on x so that [a 0 + 0.66 a 1 ], scalar [a 0 + 0.69 a 1 ], vector , (B26) using f 0 ϕ (x) defined in Eq. ( 14); see Fig. 3.The factor 1/Λ 4 comes from the normalization of the background phasespace distribution according to the number density, so that f 0 i (x) ∝ Γ ρ /(Λ 4 H 0 ), while ρ i (a) ∝ Γ ρ /H 0 as in Eq. ( 31).For ν s , the analogous procedure leads to The series [a 0 + . . .] in Eqs.(B26,B27) differ from each other because of the larger phase space available for ν s relative to ϕ. Eqs.(B16) and (B20) can be similarly estimated, but more simply since all the functions are now known.Therefore, Eqs.(B16), (B17), (B20) and (B21) reduce to where a i,ℓ and b i,ℓ are given in table IV.There one notices that the ℓ = 0, 1 terms of Eq. (B28) for ν s are equal and opposite in sign to those for ϕ, as well as the background terms.This follows from energy and momentum conservation both at the zeroth and linear order in the perturbations [77].Moreover, the second term in the first two equations in (B28) is subleading compared to Γ ρ , scaling as Γ 2 ρ /(H 0 Λ 4 ).Although it can be neglected in the background Boltzmann equations, it is needed in the perturbation equations since it is of the same order as the perturbed collision term (see Eqs. (B15) and (B28)).
a. Qualitative features of the perturbed massless equations Eqs. (B29) and (B30) can be further simplified since ρ νs (a) = −ρ ϕ (a) for massless sterile neutrinos.Then the last term on the right-hand sides, which arise from the time-dependence of the background phase-space distributions f 0 i , are negative assuming F i,ℓ > 0 and therefore provide damping of the perturbations for all values of ℓ.Such damping is typical of particle decay processes (e.g., Ref. [78]), hence we will refer to it as the "decay term" below.Its effect is larger at early times because ρ i (a) evolves faster than a.The damping is independent of Γ ρ since ρ i ∝ Γ ρ /H 0 , according to Eq. (31).
On the other hand, the second to last term of the right-hand sides of Eqs.(B29) and (B30), which we will refer to as the "collision term," scales as a, so that it grows at late times.Being proportional to Γ ρ /Λ 4 ∼ Λ (Λ/M i ) 4 , it is generally smaller than the last term.We expect that deviations in the CMB or matter power spectra with respect to ΛCDM will be apparent only when the collision term starts to dominate over the decay term, otherwise Eqs.(B29) and (B30) become identical, giving F νs,ℓ ∼ F ϕ,ℓ , implying ρ νs F νs,ℓ + ρ ϕ F ϕ,ℓ ∼ 0. However, the collision term in the high-ℓ equations leads to exponential instabilities because of the positive sign in front of the relative perturbation F i,ℓ .17Hence, it cannot be too much larger than the decay term while remaining consistent with the CMB and matter power spectrum.This criterion gives the estimate Γ ρ /Λ 4 ≲ H 0 ∼ 10 −39 MeV, in rough agreement with the quantitative result of Section IV.
This reasoning also explains why Γ ρ becomes unconstrained at sufficiently small values of Γ ρ /Λ 4 (see Fig. 10).In that regime the collision term is small compared to the decay term, and the perturbations become irrelevant.Moreover, deviations between ΛCDMϕν s and ΛCDM are more pronounced at small k (large scales), because the Liouville term (the first term) of the perturbation equations above becomes negligible, being proportional to k.This is borne out in Fig. 8.
It is noteworthy that F νs,ℓ is less suppressed than F ϕ,ℓ , at ℓ > 2, where the ϕ and ν s equations decouple.This is because (a ϕ,ℓ − a ϕ,0 ) + b ϕ,0 < (b νs,ℓ − b νs,0 ) + a νs,0 in the comparison of their collision terms; hence F νs,ℓ ≳ F ϕ,ℓ for ℓ ≳ 2. 18 The physical origin for this can be traced to the larger phase space available for sterile neutrinos than for ghosts, which increases its collision term.For the ℓ = 0, 1 perturbations, the situation is more complicated since the ν s and ϕ equations are coupled.
The perturbations for the scalar mediator have a modestly larger effect on cosmological observables relative to the vector mediator.This is ultimately due to the kinematical dependences of the matrix elements (9), which result in the coefficients for the scalar model in table IV being larger than the vector ones.This trend is confirmed by the CosmoMC results in Fig. 10.

Massive sterile neutrinos
A detailed analysis of the cosmological perturbations for the most general case where ν s can be nonrelativistic, is beyond the scope of the present paper, since in that case the perturbations play a subdominant role compared to the background contributions to the energy density.Nevertheless, we develop an approximate treatment in this appendix.
The new collision term given by Eq. (B35) is complicated in the general m νs > 0 case.However, based on the results we obtained for m νs = 0, some inferences can be made.In particular, the first term of Eq. (B35) is proportional to Γ ρ /Λ 4 , since it depends linearly on the background phase-space distribution f 0 i (q) with i = ϕ, ν s , whereas the second term is exactly given by Γ ρ δ νs /ρ νs and is thus proportional to H 0 .In order to be consistent with type Ia supernovae data, Fig. 7 implies that Γ ρ ≲ O(10 −70 ) MeV 5 .For such small values, the constraints on Γ ρ and Γ ρ /Λ 4 for m νs = 0, Fig. 10, become almost independent of Γ ρ /Λ 4 , suggesting that the first term of Eq. (B35) is negligible compared to the second.Physically, this means that redshifting of the new species, described by the Liouville terms and the second part of the collision term, is more important than the growth from particle production, when m νs > 0. In this case the perturbations of the new species decouple from each other.
By neglecting the first piece of the new collision term in Eq. (B35), we expect to find conservative bounds for the parameter Γ ρ , or quantities derived from it such as Ω g , that might be moderately strengthened by a more exact treatment.The bounds derived for m νs > 0 come primarily from the background contributions of the new species, rather their perturbations.
Carrying out the same procedure as for δ νs to higher ℓ, the perturbation equations for ϕ and ν s are δ ′ νs ∼ −( assuming a barotropic fluid, where we neglected the equations for the multipoles with ℓ > 2 since they will become negligible for nonrelativistic ν s .Similar equations for the anisotropic stress σ i were derived in Refs.[29,36,79] to approximate the closing of the ℓ hierarchy at the quadrupole, valid for both massless and massive particles.We take the viscosity parameter in Ref. [79] to be c 2 i,vis ≃ w i .These equations encapsulate the free-streaming effect of the new species, which cannot be neglected; Ref. [79] showed that doing so in the SM with massless neutrinos induces an error of ∼ 10% in the CMB power spectrum.Different prescriptions for the shear stress equation can be found in Ref. [80], which provide similar numerical accuracy. The perturbation equations derived above become the same for ϕ and ν s in the limit w νs → 1/3, since the part of the collision term sensitive to their phase space differences has been neglected.Hence the physical perturbations of the new species cancel each other in that limit, (ρ νs δ νs + ρ ϕ δ ϕ ) → 0. Similarly to the m νs = 0 treatment, the equations are solved within CAMB up to the ℓ max = 45.

FIG. 2 .
FIG. 2. Diagrams of vacuum decay mediated by a scalar Φ or a vector boson Z ′ .

FIG. 5 .
FIG.5.Left: Present-day equation of state of sterile neutrinos wν s as a function of the sterile-neutrino mass to cutoff ratio.The numerical results are shown with solid lines, while green dotted curves correspond to the numerical-fit functions(24).Right: Scale-factor evolution of wν s for different choices of mν s .The solid curves are for the vector mediator model, whereas the dashed curves are for the scalar model.The band thicknesses and Ωm and ΩΛ are as in Fig.4.

− 1 .FIG. 6 .FIG. 7 .
FIG. 6. Left: 95% C.L. allowed regions for the ghost plus sterile neutrino contribution to the present composition of the Universe, for the three supernovae data sets considered.The triangles indicate the best-fit values for Ωg.Right: Most likely values of ΩΛ (solid), Ωm (dashed) and H0 (dot-dashed) as a function of Ωg.The 95% C.L. regions of the latter are indicated as dotted vertical lines.

FIG. 9 .
FIG.8.For the special case of mν s = 0. Left: Prediction of the ΛCDMϕνs model (black curve) for the CMB temperature autocorrelation spectrum, compared to the ΛCDM best-fit model (blue curve).The Planck 2018 data points are taken from Ref.[22].Right: Predictions for the linear matter power spectrum.The data points, used only for display purposes, are given by Ref.[24], which contains data from Planck 2018[25], eBOSS DR14[26], SDSS DR7[27] and DES YR1[28].For both plots, the vector mediator is assumed, we used Γρ ≃ 10 −69 MeV 5 , Γρ/Λ 4 ≃ 10 −41 MeV and fixed the six ΛCDM parameters and H0 to their Planck 2018 best-fit values.Residuals between the model predictions and the data points are shown at the bottom.

FIG. 11 .
FIG.11.For the special case of mν s = 0. Left: Mapping of excluded regions in Fig.10to the microphysics parameters, the mediator mass Mi (i = v, s) versus the ghost momentum cutoff Λ. Right: Conversion of the same constraints to the present-day sterile neutrino abundance Ων s as a function of the ghost momentum cutoff Λ.The color coding is the same as in Fig.10

FIG. 13 .
FIG.13.Late-time evolution of (Top) the Hubble parameter H(a) normalized to its present-day value H0, and (Bottom) the energy densities of the main species normalized to the critical density of the Universe today.Parameter values are as in Fig.12.Red curve is the total dark energy density, ΩΛ + Ωg.

8
FIG. 15.Correlations between the ΛCDMϕνs model parameters, as inferred from CosmoMC for massive νs, valid independently of the mediator model used.The data sets Planck + Lensing + BAO + DES + Pantheon are the same as described in subsection IV.1.1.The darker and lighter green shaded regions correspond to the 68% and 95% C.L. intervals, respectively.

8 solution of S 8
FIG. 18. Correlations between the ΛCDMϕνs parameters for the models addressing the H0 and S8 tensions, selected from the full MCMC scan shown in Fig.15.Gray regions solve only the S8 tension, cyan addresses only the H0 tension, and blue indicates models ameliorating both.Darker and lighter regions correspond to the 68% and 95% C.L. intervals, respectively.
[65,66]are the goodness of the fit for these models (3rd column) to ΛCDM (2nd column).Although the fit to Dark Energy Survey (DES) data improves, that to Planck, BAO and Pantheon degrades, resulting in a significantly worse global fit, with ∆χ 2 ≃ 26.Nevertheless, we find that the ΛCDMϕν s model can ameliorate the H 0 tension to the ∼ 3σ level, with H 0 ≃ Mpc, while simultaneously addressing the S 8 tension.This is shown in the third column of TableII, for the best-fit ΛCDMϕν s model.Coincidentally, this model predicts Ω g ∼ 0.15, which is the same value that was preferred by the Pantheon type Ia supernovae in Fig.6.The net improvement of ∆χ 2 = −8.8relative to ΛCDM provides a mild preference of the ghost model, based on the Akaike and Bayesian Information Criteria[65,66].
TABLE II.Cosmological parameters for ΛCDM (2nd column), for the preferred ΛCDMϕνs model solving both the H0 and S8 tensions (3rd column), and for the best-fit ΛCDMϕνs model found in the full MCMC scan, Fig. 15 (last column).

TABLE V .
68% C.L. parameter intervals for the ΛCDMϕνs model, as inferred from CosmoMC in the massless νs scenario, for the vector and scalar mediators.The first 6 parameters in bold are for the ΛCDM model, whereas the following two are associated to the new physics.The remaining ones are derived from the cosmological model.

TABLE VI .
Similar to table V, but for the massive νs scenario.The results are independent of the mediator model chosen.