Searching for Light Relics with Large-Scale Structure

Light thermal relics of the hot big bang, often quantified by the parameter $N_\mathrm{eff}$, are one of the primary targets of cosmological measurements. At present, the energy density in such relics is constrained to be less than ten percent of the total energy density in radiation. Upcoming cosmic microwave background (CMB) experiments, however, have the potential to measure the radiation density at the one-percent level, which is close to well-motivated theoretical targets. In this paper, we explore to what degree the CMB observations can be enhanced by future large-scale structure surveys. We carefully isolate the information encoded in the shape of the galaxy power spectrum and in the spectrum of baryon acoustic oscillations (BAO). We find that measurements of the shape of the power spectrum can significantly improve on current and near-term CMB experiments. We also show that the phase shift of the BAO spectrum induced by relic neutrinos can be detected at high significance in future experiments.


Introduction
Future cosmological observations have the potential to measure the radiation density of the early universe at the subpercent level. This order of magnitude improvement over current constraints would provide a new window into the very early universe and allow us to search for extra light particles with very weak couplings to the Standard Model. Small changes to the radiation density of the early universe lead to well-understood changes in the anisotropy spectrum of the cosmic microwave background (CMB) [1][2][3][4]. The same effects also create imprints in the initial conditions for the clustering of matter and, hence, may be observable in the late universe. It is therefore natural to ask how much the constraints on extra relativistic species can be improved by including future observations of the large-scale structure (LSS) of the universe.
Within the Standard Model (SM) of particle physics, neutrinos make a significant contribution to the radiation density of the early universe. The cosmic neutrino background (CνB) was created about one second after the Big Bang, when the expansion rate of the universe dropped below the weak interaction scale. Shortly after neutrino decoupling, electrons and positrons annihilated, transferring their entropy to photons, but not to the neutrinos. This slightly reduced the energy density of the neutrinos relative to that of the photons. Nevertheless, 41 % of the total radiation density of the universe is still expected to be in the form of cosmic neutrinos. The gravitational effect of the CνB has recently been observed in the damping [2] and the phase shift [3,5] of the CMB anisotropy spectrum.
An interesting consequence of many proposals for physics beyond the Standard Model (BSM) are extra light particles [6], such as axions [7][8][9], axion-like particles (ALPs) [10], dark photons [11,12] and light sterile neutrinos [13]. These particles are often so weakly coupled to the SM that they escape detection in terrestrial experiments. However, in astrophysics and cosmology, we have access to high-density environments which can overcome the small cross sections and allow a significant production of the extra species. For example, new light particles can be produced in the interior of stars [14]. The absence of an anomalous extra cooling over the lifetime of stars puts some of the best current constraints on weakly coupled species. A similar argument can be applied to cosmology [15][16][17]. The high densities of the early universe allow these particles to have been in thermal equilibrium with the SM and can therefore make a significant contribution to the total radiation density of the universe. New particles that are more weakly coupled than neutrinos would have decoupled before the QCD phase transition. Their contribution to the final radiation density is then suppressed, explaining why these particles have not been detected yet. In this paper, we will explore the sensitivity of future cosmological observations to this type of BSM physics.
The search for light thermal relics has been adopted as one of the main science targets of the next generation of CMB experiments, such as the CMB-S4 mission [18]. Through improved measurements of small-scale anisotropies and polarization, future CMB observations will be extremely sensitive to the damping and the phase shift of the anisotropy spectrum. In this work, we explore the additional constraining power provided by current and future LSS experiments, such as (e)BOSS [19,20], DES [21], DESI [22], LSST [23] and Euclid [24]. It was established in [25][26][27] that these surveys carry information about relativistic species. We will examine how this infor-mation is encoded in both the shape of the matter power spectrum and the spectrum of baryon acoustic oscillations (BAO). We find that measurements of the shape of the power spectrum can significantly improve on the current CMB constraints, although the largest improvements are subject to the usual challenge of modeling the power spectrum. The peak locations of the BAO spectrum carry additional information about light relics that is robust to corrections to the overall shape of the power spectrum [28], such as those arising from nonlinear gravitational evolution [29][30][31]. We will explore in detail how this information can be isolated in the BAO spectrum. This protected information may play a useful role in elucidating apparent discrepancies between CMB and low-redshift measurements, and be a valuable tool in the search for exotic physics in the dark sector.
The outline of the paper is as follows. In Section 2, we present the theoretical motivation for a precise measurement of the radiation density in the early universe, focusing on the effects of extra light species on the spectrum of acoustic oscillations. We highlight that these effects are imprinted in both the CMB and BAO spectra. In Section 3, we forecast CMB and LSS constraints on the number of relativistic species, N eff , for a number of future observations. In Section 4, we isolate the information encoded in the phase shift of the BAO spectrum and study the prospects for extracting this information in upcoming surveys. Our conclusions are presented in Section 5.
A series of appendices contain technical details of our analysis: In Appendix A, we describe our CMB forecasts and present results for a range of experimental configurations. In Appendix B, we provide details of our LSS forecasts. We define the specifications for the galaxy surveys used in this work and present results for a range of data combinations and cosmologies. In Appendix C, we outline our method for extracting the broadband spectrum and the phase shift. Finally, in Appendix D, we show a few of the convergence tests that we performed to establish the stability of our numerical analysis.

Cosmological Signatures of Light Relics
It is rather remarkable that all current cosmological data (e.g. [32][33][34][35]) is fit by a simple sixparameter model-the ΛCDM model. In this section, we introduce the standard cosmological model and its extension to include extra relativistic species. We review the imprints that light particles leave on the cosmic microwave background and the large-scale structure of the universe. We will pay particular attention to the unique signature that these particles leave on the spectrum of acoustic oscillations. In the next section, we will quantify the level of constraints on extra light species to be expected from future cosmological observations.

The Standard Model
The ΛCDM model includes two parameters characterizing the initial conditions, namely the amplitude A s and the tilt n s of the spectrum of primordial curvature perturbations. The remaining four parameters are associated with the geometry and composition of the universe: The matter content of the universe is described by the physical baryon and dark matter densities, ω b ≡ Ω b h 2 and ω c ≡ Ω c h 2 , where h is the reduced Hubble constant h ≡ H 0 / 100 km s −1 Mpc −1 . Instead of the Hubble constant H 0 , we use the angular size of the sound horizon at decoupling, θ s ≡ r s (z * )/D A (z * ), where r s is the physical sound horizon and D A is the angular diameter distance, both evaluated at the redshift of decoupling, z * . The parameter θ s receives a contribution from the dark energy density Ω Λ . The standard six-parameter model is completed by the optical depth τ . In Table 1, we list the fiducial values of the ΛCDM parameters, based on the Planck best-fit cosmology [33].
Parameter Fiducial Value Description Physical dark matter density ω c ≡ Ω c h 2 100 θ s 1.04112 100 × angular size of the sound horizon at decoupling τ 0.066 Optical depth due to reionization ln(10 10 A s ) 3.064 Log of scalar amplitude (at pivot scale k 0 = 0.05 Mpc −1 ) n s 0.9667 Scalar spectral index (at pivot scale k 0 = 0.05 Mpc −1 ) Effective number of (free-streaming) relativistic species Y p 0.2478 Primordial helium fraction Table 1: Parameters of the reference cosmological model and their fiducial values based on [33].
In this work, we are interested in future measurements of the radiation density of the universe. The contribution from photons, ρ γ , is fixed by the measured value of the CMB temperature. In addition, the Standard Model of particle physics predicts a contribution from neutrinos. The expected radiation density from each neutrino species is The three neutrino species of the Standard Model (and their antiparticles) therefore contribute a significant amount to the total radiation density in the early universe: ρ ν /ρ r = i ρ ν i /ρ r ≈ 41 %.
Although neutrinos decoupled at early times, their gravitational effects are still relevant and have recently been observed in the CMB [3,5].

Extra Light Relics
Physics beyond the Standard Model may add an extra radiation density ρ X to the early universe. 1 It is conventional to measure this radiation density relative to the density ρ ν i of a single SM neutrino species: 2) and define N eff = 3.046 + ∆N eff as the effective number of neutrinos, although ρ X may have nothing to do with neutrinos. Current measurements of the CMB anisotropies and the light  Figure 1: Contributions of a single thermally-decoupled Goldstone boson, Weyl fermion or massless gauge boson to the effective number of neutrinos, ∆N eff , as a function of its decoupling temperature T dec . The drop in ∆N eff around 150 MeV is due to the QCD phase transition, where we employed the lattice QCD calculation of [37]. element abundances find [33,36] which is consistent with the SM prediction of N eff = 3.046. We expect that future cosmological observations will improve these constraints by up to an order of magnitude. Any non-zero value for ∆N eff would indicate physics beyond the standard models of particle physics and/or cosmology.
A natural source for ∆N eff = 0 are extra relativistic particles. Figure 1 shows the contribution to ∆N eff from a single thermally-decoupled species as a function of the decoupling temperature T dec and the spin of the particle. The plot assumes that the extra species was in thermal equilibrium at some point in the history of the universe and that the number of relativistic degrees of freedom at decoupling was not significantly larger than the SM value. We also assumed no significant entropy production after decoupling. We see that decoupling after the QCD phase transition produces a contribution to N eff that is comparable to that of a single neutrino species, which is in tension with current observations. Decoupling before the QCD phase transition, however, creates an abundance that is smaller by an order of magnitude and hence still consistent with current limits. Future observations will therefore give us access to particles that are more weakly coupled than neutrinos. The exclusion of the minimal thermal abundance ∆N eff = 0.027 would have important consequences for BSM physics [15][16][17]. We find it intriguing that this threshold seems to be within reach of future CMB and LSS observations. In this paper, we will quantify this expectation. N eff Figure 2: Variation of the CMB power spectrum as a function of N eff . The spectra have been rescaled, so that the fiducial spectrum for N eff = 3.046 is undamped, i.e. the exponential Silk damping was removed. Following [5], the physical baryon density ω b , the scale factor at matterradiation equality a eq ≡ ω m /ω r and the angular size of the sound horizon θ s are held fixed in all panels. The dominant effect in the first panel is the variation of the damping scale θ D . In the second panel, we fixed θ D by adjusting the primordial helium fraction Y p . The dominant variation is now the amplitude perturbation δA. In the third panel, the spectra are normalized at the fourth peak. The remaining variation is the phase shift φ (see the zoom-in in the fourth panel).

Phases of New Physics
Keeping the acoustic scale θ s fixed (e.g. by adjusting the Hubble constant H 0 ), an increase in the radiation density of the early universe reduces the mean free path of fluctuations in the photon-baryon fluid and increases the damping of small-scale fluctuations [2] (see Fig. 2). The constraint in (2.3) is mostly derived from measurements of the CMB damping tail [33,38]. However, the damping tail is also affected by changes to the primordial helium fraction, Y p , which induces a variation in the free electron fraction and hence the mean free path of photons. In our forecasts, we will both fix Y p to the value demanded by BBN consistency (ΛCDM+N eff ) and vary it (ΛCDM+N eff +Y p ) to explore the degeneracy with N eff . The part of N eff that is associated with free-streaming relativistic particles leads to a characteristic phase shift in the CMB spectrum [1,3] (see Fig. 2), which helps to break the degeneracy between N eff and Y p . The phase shift associated with SM neutrinos has recently been measured in the Planck spectrum [3,5].
In this work, we pay particular attention to the information about N eff contained in the BAO spectrum. To isolate the BAO signal, we split the power spectrum into a smooth ('nowiggle') part and an oscillatory ('wiggle') part, P (k) ≡ P nw (k) + P w (k) .
(2.5)   Figure 3: Variation of the matter power spectrum P (k) (top) and the BAO spectrum P w (k)/P nw (k) (bottom) as a function of N eff . The physical baryon density ω b and the physical sound horizon at the drag epoch, r s , are held fixed in all panels of the BAO spectrum. In the second BAO panel, we fixed the scale factor at matter-radiation equality, a eq ≡ ω m /ω r . The variation in the BAO amplitude δA is then the dominant contribution. In the third BAO panel, the spectra are normalized at the fourth peak and the bottom panel shows a zoom-in illustrating the remaining phase shift.
Our method for performing this separation is described in Appendix C. We will demonstrate that the most robust information about N eff lives in P w (k). In particular, it was shown in [28] that the phase of the BAO spectrum is immune to the effects of nonlinear gravitational evolution. In Figure 3, we show the dependence of the phase of the BAO spectrum on the number of relativistic species N eff . We claim that this information is preserved after nonlinear corrections are taken into account.

Future Constraints on Light Species
We have argued that measuring the radiation density at the percent level provides an interesting window into early universe cosmology and beyond the Standard Model particle physics. In this section, we will further quantify the constraining power of future cosmological observations. We will consider two types of forecasts based on P (k) and P w (k). We will refer to these as 'P (k)forecasts' and 'BAO-forecasts', respectively.

Fisher Methodology
We will use standard Fisher information theory to forecast the constraints of future observations. While Fisher forecasts have to be used with care, they provide useful guidance for the sensitivities and design of future experiments. In this section, we recall the basic elements of the Fisher methodology and its application to galaxy surveys [39,40]. The relatively standard Fisher forecasting of CMB observations is summarized in Appendix A. Further details on the LSS forecasting can be found in Appendix B.
Given a likelihood function L( θ ) for the model parameters θ ≡ {ω b , ω c , θ s , τ, A s , n s , N eff , Y p }, we define the Fisher matrix as the average curvature of the log-likelihood around the fiducial point in parameter space, where the expectation value denotes an average over all possible realizations of the data. If the likelihood is Gaussian, then the inverse Fisher matrix gives the covariance matrix. This means that F −1/2 ii is the error on the parameter θ i , when all other parameters θ j =i are known, while σ(θ i ) = (F −1 )

1/2
ii is the error on θ i after marginalizing over the other parameters. More generally, the Cramér-Rao bound, gives a lower limit on the marginalized constraints.
The Fisher matrix for a galaxy survey is [41] where P g (k, µ) is the anisotropic galaxy power spectrum, µ is the cosine between the wavevector k and the line-of-sight, and V eff is the effective survey volume, In the second equality, we have assumed that the comoving number density of galaxies is independent of position, n g ( r ) ≈n g = const, and introduced the actual survey volume V . To derive the constraints from independent redshift bins, we take V to be the volume within each bin and add the corresponding Fisher matrices. The minimum wavenumber accessible in a survey is given by the volume of the survey 2 as k min = 2π [3V /(4π)] −1/3 .

Modeling the Power Spectrum
In §2.3, we introduced the linear matter power spectrum P lin (k), and separated it into its smooth and oscillatory parts. In order to obtain semi-realistic constraints on most parameters of the cosmological model, it is often sufficient to model the observed galaxy power spectrum as P g (k) ≈ b 2 P lin (k), where b is the linear biasing parameter. However, the constraints on extra relativistic species are particularly sensitive to the way degeneracies are broken and to the nonlinear damping of the oscillatory feature, so we need to be more careful in the modeling of the signal [42][43][44]. Moreover, since observations only determine the angular positions and redshifts of objects, we need to take into account the corresponding redshift space distortions (RSD) and geometric projection effects.
Our model for the observed galaxy power spectrum is the following remapping of the linear matter power spectrum: (3.5) All functions in this expression have an implicit redshift dependence. We now define the different elements of (3.5): • O(k, µ): This function encodes the BAO signal and can be written as where O lin (k ) ≡ P w lin (k )/P nw lin (k ) is the normalized wiggle spectrum evaluated at the rescaled wavenumbers [45] This rescaling reflects the fact that the wavenumbers k cannot be measured directly, but instead have to be derived from the measured angles and redshifts using the angular diameter distance D fid A (z) and Hubble rate H fid (z) of a fiducial cosmology. This is often referred to as anisotropic geometric effects. In the limit of spherically-averaged clustering measurements, these become isotropic and k = k/q, where q = q , with the radial BAO dilation given by D V ∝ (D 2 A /H) 1/3 . To model uncertainties in the BAO extraction, we have introduced two free functions B(k) and A(k) in (3.6), which we take to be smooth polynomials in k (see §3.1.2). Ultimately, we will marginalized over these polynomials to remove any information that is not robust to the BAO signal itself.
• b(z): The bias of the target galaxies (e.g. luminous red galaxies, emission line galaxies or quasars) sets the overall amplitude of the signal in each redshift bin. We will make the common assumption that b(z) ∝ 1/D 1 (z), where D 1 (z) is the linear growth function. This means that the bias is larger at high redshifts, which implies that the galaxy power spectrum may get significant corrections from nonlinear biasing even at high redshifts.
• F (k, µ): This function characterizes the effect of redshift space distortions. Following [46], we write F (k, µ) = 1 where β ≡ f /b, with the linear growth rate f ≡ d ln D 1 / d ln a. The factors of q i account for differences in the cosmic volume in different cosmologies. Projection effects on the angle to the line-of-sight are included as [45] µ (k, µ) = µ/ µ 2 + (1 − µ 2 )Q 2 , (3.9) where Q ≡ q /q ⊥ , which becomes unity in the isotropic case. BAO reconstruction removes redshift space distortions on large scales, which we have modeled by adding the factor where the value of Σ s depends on the experimental specifications, in particular the noise levels. In our baseline forecasts, we take Σ s → ∞, i.e. R ≡ 1, but we comment on finite values of Σ s in §3.2.1.
• D(k, µ): This function models the nonlinear damping of the BAO signal [29,47] where the damping scales perpendicular and parallel to the line-of-sight are given by with σ 8 being the amplitude of (linear) matter fluctuations at a scale of 8 h −1 Mpc. We account for BAO reconstruction by decreasing these damping scales by an appropriate factor, e.g. 0.5 for 50 % reconstruction. Following [25,48], we include the degradation in the reconstruction due to shot noise using a reconstruction multiplier r(x), i.e. Σ i → r(x)Σ i . We obtain r(x) by interpolating over the table with r(x < 0.2) = 1.0 and r(x > 10.0) = 0.5, which depends on the number densityn g via x ≡n g P g (k 0 , µ 0 )/0.1734 evaluated at k 0 = 0.14 h Mpc −1 and µ 0 = 0.6. This means that we assume 50 % reconstruction at high number densities and no reconstruction for low densities.
• P nw (k, µ): The linear no-wiggle spectrum P nw lin (k, µ) is determined from the linear power spectrum using the method described in Appendix C. Nonlinear corrections to this spectrum can be parameterized as whereB(k) andÃ(k) are smooth functions (see §3. 1.2). For the purpose of our BAOforecasts,Ã(k) andB(k) are degenerate with A(k) and B(k) in (3.6) and it is therefore consistent to use the linear spectrum.
• Z(k, µ): For photometric surveys, we take the uncertainty in the redshift determination of the targets into account through the following function: where Σ z = c (1 + z) σ z0 /H(z) is given in terms of the root-mean-square redshift error σ z0 [49,50]. The redshift error, which depends on the experimental specifications, reduces the effective resolution for modes along the line-of-sight. We neglect this effect for spectroscopic surveys.
When evaluating the derivatives in the Fisher matrix (3.3), the parameters b(z), β, R(k), D(k, µ) and Z(k, µ) are always computed using the fiducial cosmology. We are assuming that, after accounting for modeling uncertainties, no relevant cosmological information can be recovered from these functions.

Accounting for Broadband Effects
Nonlinear evolution and biasing can change the shape of the power spectrum at high wavenumbers in a way that cannot be modeled from first principles. We account for this uncertainty by marginalizing over polynomials in k in both the P (k)-and BAO-forecasts. In particular, the functions introduced in (3.14) are defined as As indicated, we allow independent polynomials in each redshift bin centered around z i . The coefficientsã n,i andb m,i are included in the list of parameters θ i . Derivatives with respect to these parameters are determined analytically, using the fiducial valuesb 0,i = 1 andã n,i =b m =0,i = 0. A more careful treatment would replace this polynomial model with a perturbative model for the dark matter and biasing, and would marginalize over the bias parameters. In practice, this has been shown to give qualitatively similar forecasts [51]. Our marginalization procedure is therefore sufficient to illustrate the sensitivity of our forecasts to broadband information.
Our BAO-forecasts will marginalize over the 'broadband corrections' in (3.6), with A(k) and B(k) defined as in (3.16). 3 At the level of the Fisher matrix, marginalizing over a polynomial and an exponential are equivalent. As a result, the function B(k) captures the uncertainty in the damping scales Σ and Σ ⊥ in (3.10). This implies that our marginalization procedure will eliminate any cosmological information associated with the nonlinear damping of the power spectrum, leaving the distinct information contained in the oscillating part of the spectrum O lin (k (k, µ)). This type of procedure is used in the analysis of BAO data to correct for errors made in the modeling of P nw (k), see e.g. [52].
We will choose various levels of marginalization in our forecasts. This will help to distinguish the information encoded in the smooth shape of the spectrum, P nw (k), from that contained in the frequency and phase of the BAO spectrum, P w (k). In addition, these marginalizations also give a sense for the level of robustness of each type of information when accounting for the various uncertainties in modeling the data of a realistic galaxy survey.

Extracting the BAO Signal
In describing the power spectrum, we introduced the idea of marginalizing over polynomials to remove the information in P g (k) that is thought to be degenerate with nonlinear evolution and galaxy biasing. The BAO spectrum is known to be robust to these effects and should therefore survive any such treatment. In principle, the BAO signal could be isolated with sufficient marginalization. However, in practice, it is more useful to extract the information associated with the BAO signal before any marginalization. The robustness of the BAO spectrum to nonlinearities means we can be more aggressive with our choice of k max and less cautious with our marginalization. Consequently, it is convenient to treat the BAO signal and the broadband information independently.
The observed BAO spectrum is defined by where D(k, µ) and O(k, µ) were introduced in (3.5). To derive the new Fisher matrix for the BAO spectrum directly, we first write the derivatives of P g (k, µ) as We then drop the term proportional to ∂ θ i P nw g since it is degenerate with the marginalization over the broadband corrections. For the same reason, we write we do not act with the derivatives on the functions D(k, µ) and bF (k, µ). The derivative in (3.18) therefore becomes While the derivatives that we have dropped are non-zero, the marginalization procedure described above is designed to remove them and the forecasts for cosmic parameters should consequently be the same. Removing this information by hand (and marginalizing) ensures that our BAOforecasts do not include these broadband effects, as we will show in Fig. 5. The resulting Fisher matrix is then given by We note that this Fisher matrix depends on P nw g (k, µ) only through V eff (k, µ), which determines the signal-to-noise. For photometric surveys, we replace V eff (k, µ) → Z(k, µ) 2 V eff (k, µ) to account for the redshift error and the associated reduction of power along the line-of-sight. In principle, we should model P nw g (k, µ) using the nonlinear (galaxy) power spectrum, given that we will work close to the nonlinear regime. However, nonlinear evolution also correlates the modes and produces a non-Gaussian covariance matrix. Since most of the surveys under consideration in this paper are limited by shot noise, using the nonlinear power spectrum without taking into account the associated mode coupling in the covariance would artificially increase the number of signal-dominated modes. To be consistent with the use of a Gaussian covariance, our forecasts will therefore use the linear broadband spectrum.

Summary of Results
We are now ready to forecast the constraints of current and future CMB and LSS observations on the effective number of relativistic species N eff . Unless stated otherwise, our baseline analysis assumes a ΛCDM+N eff cosmology in which the primordial helium fraction Y p is fixed by consistency with BBN. At the end of the section, we will also present results with Y p as a free parameter. We will further dissect the information content of the BAO spectrum in the next section.
In Appendix A, we present detailed forecasts for current and future CMB experiments. The expected 1σ constraints for representative versions of the Planck satellite, a near-term CMB-S3 experiment and a future CMB-S4 mission are σ(N eff ) = 0.18, 0.054, 0.030, respectively. In §A.3, we show how these constraints depend on variations of the experimental configurations. We would like to know how much these CMB constraints would improve with the addition of LSS data.
We will give the results of two types of forecasts based on P (k) and P w (k). Our P (k)forecasts apply the Fisher matrix (3.3) with k max = 0.2 h Mpc −1 and marginalize over b m≤1 . To be conservative about nonlinear biasing, we do not increase k max at large redshifts, despite the (near-)linearity of the matter power spectrum. Our BAO-forecasts use the Fisher matrix (3.20) with k max = 0.5 h Mpc −1 and marginalize over a n≤4 , b m≤3 . We will also show how these forecasts depend on the choice of k max and the level of marginalization.

Constraints from Planned Surveys
A number of galaxy surveys are expected to take place over the next decade. The power of these surveys to constrain N eff is most sensitive to the survey volume, the number densities of galaxies and the redshift errors (spectroscopic versus photometric). The precise specifications of the surveys used in our analysis are given in Appendix B, where we also present more detailed forecasts for the full set of parameters.
Baseline results In Table 2 Table 3: Forecasted 1σ constraints on N eff for various combinations of current and future CMB and LSS experiments using BAO-forecasts with k max = 0.5 h Mpc −1 .
of current and future CMB and LSS experiments using the full P (k)-forecast. In Table 3, we compare these results to the same experiments using our BAO-forecasts. At BOSS levels of sensitivity and number densities, the BAO feature makes the most significant impact on constraints, particularly when combined with a CMB experiment like Planck. In contrast, with the larger volume and redshift range of DESI, the broadband shape carries most of the information and can lead to a significant improvement in the constraint on N eff both for Planck and a typical CMB-S3 experiment. Finally, photometric redshift surveys like DES and LSST generally perform worse than spectroscopic surveys because they are effectively two-dimensional for the scales of interest. However, the employed redshift error is conservative and we do not take the full potential of these surveys into account as we are only considering observations of galaxy clustering and have not included weak gravitational lensing measurements, for instance. We expect the constraints to improve with these additional LSS observables, but quantifying this is beyond the scope of this work.
Sensitivity to k max The broadband signal is sensitive to nonlinear effects and we should therefore understand how sensitive these results are to the choice of k max . In particular, we have chosen k max = 0.2 h Mpc −1 in Table 2, but the usable range of scales is uncertain. Figure 4 shows how the constraints vary as a function of the maximal wavenumbers included in the analysis, k max , for both the P (k)-and BAO-forecasts. For the BAO-forecasts, we see a clear plateau for k max > 0.2 h Mpc −1 . This behavior is due to the damping of the oscillations at higher k relative to the smooth power spectrum. Cosmic variance is ultimately determined by the amplitude of the smooth power spectrum and one cannot recover the high-k oscillations even by lowering the shot noise. In contrast, the P (k)-forecasts show improvements out to k max > 0.3 h Mpc −1 .
Given that the BAO spectrum is robust to nonlinear evolution, it is natural to consider an optimal combination of the P (k) and BAO spectra that uses all the available information. This means using P (k) up to a certain k max and adding BAO-only information for larger k. The k max of the P (k) analysis then becomes the k min of the BAO analysis to avoid double counting the information. Results for this optimal combination are shown as the dotted line in Fig. 4.
Sensitivity to marginalization High-redshift galaxy surveys benefit significantly from measuring highly biased objects. These large biases can offset the growth function, b(z)D 1 (z) = const, and keep the amplitude of the galaxy power spectrum effectively fixed at high redshift. This boost is important for maintaining a signal above the shot noise, which we have assumed is redshift-  independent. As a consequence, high-redshift and low-redshift galaxy power spectra are equally sensitive to uncertainties in the biasing coefficients. This is particularly significant when determining the largest wavenumbers that carry useful cosmological information. While taking k max > 0.2 h Mpc −1 is appealing to maximize the constraints on N eff , we must also marginalize over successively more bias parameters. Figure 5 shows how the results depend on the marginalization scheme. While both the P (k)-and BAO-constraints degrade significantly when going from no marginalization to a few bias parameters, the BAO-forecasts quickly become robust to the marginalization. In contrast, the P (k)-forecasts weaken notably with additional biasing, but always lie below the BAO-only results, as one would expect. This confirms the intuition that the information that is primarily driving the constraints derived from P (k) is present in the no-wiggle power spectrum, P nw (k), instead of the BAO spectrum.
It is instructive to compare the results of our BAO-forecasts with those of a standard BAO analysis. Specifically, it is conventional to use the BAO signal to constrain only q i , i =⊥, , defined in (3.7) and derive parameter constraints from them. 4 These derived limits on N eff are shown as Planck+DESI-P (k) Figure 5: Forecasts for BOSS and DESI combined with Planck as a function of the largest Fourier modes used in the forecast, k max , using various levels of both additive and multiplicative marginalization, cf. the a i and b i -terms in (3.16). We have varied the number of parameters in the marginalization from none (/ a) to five (a n≤4 ) and none ( / b) to six (b m≤5 ), respectively. The dashed line shows the constraints from a standard isotropic BAO analysis for comparison.
the dashed lines in Fig. 5. The fact that the standard BAO constraints are slightly weaker than those of our full BAO-forecasts, even after marginalization, suggests there is information in the BAO spectrum beyond the BAO scale. We will explore this further in Section 4.
Degeneracy with Y p To explore possible degeneracies between the effective number of relativistic species N eff and the primordial helium fraction Y p , we now consider a ΛCDM+N eff +Y p cosmology. In Tables 4 and 5, we present the 1σ constraints on N eff and Y p for various combinations of current and future CMB and LSS experiments using P (k)-forecasts and BAO-forecasts, respectively. As expected, the CMB-only constraint on N eff become worse due to the well-known degeneracy between N eff and Y p in the CMB damping tail. When broadband information is included, we find significant improvements in the constraints on both N eff and Y p . However, this improvement cannot be attributed to the phase shift as we see only modest improvements in our BAO-forecasts. The broadband shape of the matter distribution is sensitive to the expansion history and to free-streaming neutrinos, but is not significantly affected by Y p . As a result, the broadband information in P (k) can break CMB degeneracies even without the phase shift information.   Comments on reconstruction In our baseline forecasts, we took R ≡ 1 in (3.8), which is equivalent to taking Σ s → ∞. A few comments are in order regarding the effect of a finite Σ s . As discussed in [53], the optimal smoothing scale Σ s used in the BAO reconstruction depends on the noise levels of the experiment. Having said that, we have found only small changes in our results when going from Σ s = ∞ to finite Σ s . The constraints quoted in Tables 2 to 5 are basically unaffected, except for DESI and Euclid in the P (k)-forecasts, where the impact is also mild. Changing Σ s from 30 h −1 Mpc to 15 h −1 Mpc and 10 h −1 Mpc, the constraint on N eff slightly weakens from 0.090 to 0.093 and 0.096 for Planck+DESI (0.082, 0.086 and 0.090 for Planck+Euclid) in ΛCDM+N eff compared to the quoted 0.087 (0.079) in Table 2. In practice, this roughly 10 % effect has to be compared to the impact on the reconstruction efficiency.

Designer's Guide for Future Surveys
One of the main benefits of a Fisher forecast is that it can inform the design of future experiments. For spectroscopic surveys, the basic parameters are the total number of objects, N g , the maximal redshift, z max , and the sky area in square degrees, Ω. From these, we derive the survey volume, V , and the comoving number density,n g . 5 In this section, we will explore how the constraints on N eff depend on these parameters.
Most of the survey characteristics are encoded in the effective survey volume, 6 V eff , cf. (3.4) and (3.20). The dependence of V eff on the survey parameters is somewhat non-trivial. Increasing V (by increasing z max and/or Ω), at fixed N g , will also reducen g . For signal-dominated modes,n g P g 1, this effect is not important and the effective volume scales approximately as V eff ∝ V . However, forn g P g 1, the shot noise is important and the reduction in the comoving density is more important than the increase in the volume, so that the effective volume scales as V eff ∝ V −1 . This means that we will only benefit from an increase in the volume as long as the As mentioned before, the increased linearity of the matter distribution at high redshifts is undermined by the larger biasing. As a result, the main benefit of large z max is the increased survey volume and hence the total number of modes. Unfortunately, the survey volume only grows slowly with redshift for z > 2 and the resulting improvements in parameters is relatively modest for large increases in z max . The situation is slightly different for the BAO spectrum as the nonlinear damping factor D(k, µ) depends on the clustering of the matter directly and is therefore less important at high redshifts. However, the BAO signal alone has a relatively modest effect on N eff forecasts in general and the change to the damping factor consequently does not make a visible difference in our forecasts.
In the top panel of Fig. 6, we present P (k)-forecasts for N eff for a variety of survey configurations, assuming Y p is fixed by BBN consistency. We see that the largest improvement comes from increasing N g from 10 7 to 10 8 . As we increase the number of objects further, we reach the cosmic variance limit for all modes of interest. We see that an optimistic future survey combined with a near-term CMB experiment can provide constraints that are comparable to (or slightly stronger than) those projected for CMB-S4 alone. Having said that, it does not appear that one can push the measurement of N eff well beyond the CMB-S4 target. Moreover, as in the case of the planned experiments, the improvements from the BAO signal alone are rather small.
The value of LSS becomes more significant as we expand the space of parameters. The bottom panel of Fig. 6 shows P (k)-forecasts for ΛCDM+N eff +Y p . We again see that the most significant jump in sensitivity arises when N g increases from 10 7 to 10 8 . We note that a factor of two improvement in σ(N eff ) over CMB-S4 seems possible. We also see that the P (k)-forecasts for N eff marginalized over Y p are competitive with CMB-only forecasts with Y p held fixed. In this sense, P (k) adds robustness to the measurement of N eff under broader extensions of ΛCDM. The improvement in Y p is slightly weaker, but shows the same general trend.
The range of accessible modes in near-term galaxy surveys is limited by their reliance on highly biased objects, but more futuristic surveys may not have the same limitations. Future surveys can also have high signal-to-noise beyond k = 0.2 h Mpc −1 , making it worth to consider the impact of increasing k max . In Figure 7, we show the potential reach of two representative surveys. The first, denoted "Future", is characterized by N g = 10 8 , f sky = 0.5 and z max = 3, which is roughly the same as a spectroscopic follow-up to LSST. The second, denoted "CVL", is cosmic variance  limited for all k ≤ k max over f sky = 0.5 and z max = 6. In principle, a 21 cm intensity mapping survey could achieve similar performance [27]. We see that σ(N eff ) ∼ 0.015 is achievable through the measurement of P (k) in either survey for k max = 0.5 h Mpc −1 , although the improvement with CVL is more robust to marginalization.

Measurements of the Phase Shift
In the previous section, we showed how much the combination of future CMB and LSS measurements can improve the sensitivity to extra relativistic species. The dominant source of improvement came from the broadband shape of the power spectra, P nw (k), rather than the BAO spectrum, P w (k). Nevertheless, the shift of the acoustic peaks is a particularly robust signature of free-streaming, relativistic species [28] and is therefore an interesting observable in its own right. In this section, we will isolate the signal coming from the phase shift and forecast our ability to measure it in future surveys. Measuring the BAO phase shift provides an independent test of pre-recombination physics in a low-redshift observable. This could be used to shed light on possible discrepancies between low-and high-redshift measurements or as a discovery channel for exotic new physics.

Isolating the Phase Shift
The BAO feature in Fourier space can be written as where the parameter α represents changes in the BAO scale r s , and the amplitude modulation A(k) and the phase shift φ(k) encode a number of physical effects that alter the time evolution of the baryons. While α and A(k) are implicit functions of redshift, φ(k) is redshift independent. Relativistic species are the unique source of a constant shift in the locations of the BAO peaks in the limit of large wavenumbers, i.e. φ(k → ∞) = φ ∞ [1,3]. In practice, however, the measurement of the BAO spectrum occurs over a relatively small range of scales with a small number of (damped) acoustic oscillations. On these scales, the k-dependence of the shift can be relevant. Furthermore, additional k-dependent shifts from other cosmological parameters may also have to be taken into account [4].
To measure the phase shift φ(k), we will construct a template for the k-dependence as a function of the relevant parameters. For small variations around their fiducial values, it is a good approximation to treat the shifts arising from each cosmological parameter independently. By varying one parameter at a time and measuring the change in the peak locations, we can construct a template φ(k) = i β i ( θ )f i (k). For ΛCDM+N eff , the parameters A s , n s , and τ do not affect the evolution of the baryons prior to recombination and, therefore, do not change the phase of the oscillations. The parameters ω b and θ s do alter the BAO spectrum, but are effectively negligible for any realistic parameter range. The shifts induced by ω c and N eff , on the other hand, can be significant.
The parameter that is most independent of N eff is not the dark matter density ω c , but the scale factor at the time of matter-radiation equality, a eq . Since CMB data essentially fixes a eq , our template model can be reduced to namely the shift induced by changing N eff at fixed a eq . This is the same choice made by Follin et al. [5] in their CMB measurement of the phase shift. Fixing a eq also reproduces the expected constant phase shift at large wavenumbers. The template for the phase shift at fixed ω c , in contrast, does not approach a constant at large wavenumbers, which implies that the change of a eq to maintain constant ω c is introducing a phase shift of comparable size to the constant shift induced by varying N eff . For our applications, this additional shift plays no role, but it could be useful in future investigations.
We describe the measurement of the phase shift and the construction of the template in Appendix C. In short, we determine the shift in the locations of the peaks/troughs and zeros of the BAO spectrum compared to the fiducial cosmology with N eff = 3.046 and sample 100 different cosmologies with varying N eff at fixed a eq . It is convenient to normalize the template f (k) such that β = 0 and 1 for N eff = 0 and 3.046, respectively. In Figure 8, we illustrate how the peaks/troughs and zeros of the BAO spectrum change in response to this variation in N eff . We see that the phase shift created by N eff approaches a constant at large wavenumbers in line with physical expectations. The measurement of the phase shift is challenging because it requires a very accurate model of the no-wiggle spectrum P nw (k) across a wide range of cosmological parameters. Errors in P nw (k) effectively change the functions A(k) and B(k) in (3.6) and lead to errors in the measurement of the BAO peaks and zeros, respectively. The small size of the phase shift in Fig. 8 only exacerbates this problem. Fortunately, while the template is difficult to generate, our forecasts using the template are very stable. Furthermore, the template is well approximated by a simple fitting function, where φ ∞ = 0.227, k = 0.0324 h Mpc −1 and ξ = 0.872 were obtained by a weighted fitting procedure. From the analytic treatment at high wavenumbers k, we expect φ ∞ = 0.191π fid + O( 2 fid ) ≈ 0.245 to linear order [1,3], where (N eff ) = N eff /(a ν + N eff ) is a measure of the excess radiation density, (ρ r − ρ γ )/ρ r , with a ν ≈ 4.40 as introduced in (2.1). This approximation overestimates the value obtained using the fitting formula by about 8 %, which is consistent with the expected corrections from higher orders in fid ≈ 0.41. Around k ∼ 0.1 h Mpc −1 , where BOSS and DESI have the largest signal-to-noise ratio, the relative difference is almost 50 %, which makes it evident that the offset from the analytic approximation has to be taken into account in an analysis such as the one proposed below, whereas the precise shape of the template plays a sub-dominant role. We also note that this template is basically independent of changes to the BAO scale r s , for example due to changes in the dark matter density.
We use the measured phase template to write the BAO spectrum in terms of the spectrum in the fiducial cosmology: where α ≡ α(z i ) takes an independent value in each redshift bin centered around z i and β is a single parameter for the entire survey. A measurement of α(z i ) and β can then be translated into constraints on cosmological parameters using where the parameters q and D V were introduced in §3.1.1. With this normalization, the largest possible phase shift due to N eff is given by β(N eff → ∞) = 2.45.
In §4.3, we will show that the forecasts produced using only these templates are in agreement with the forecasts using the full BAO spectrum. From a measurement of β > 0, one gets a constraint on N eff that is only associated to the size of the phase shift. This approach is analogous to the template-based measurement of the phase shift in the CMB by Follin et al. [5]. The measurement of N eff from the phase alone ignores the effects of N eff on α, but has the advantage that any detection is unambiguously 7 a measurement of free-streaming relativistic particles.
We will also be interested in the measurement of β when a prior on α is included, e.g. from the CMB. 8 In a given cosmological model, the parameter α is fully determined by the set of cosmological parameters, α = α( θ ). Since the α(z i ) inferred from the CMB are correlated between the n redshift bins of a galaxy survey and n is in general larger than the number of cosmological parameters, we compute the n-dimensional inverse covariance matrix according to C −1 α = A T F A, where F is the Fisher matrix and A is the pseudo-inverse of ∇ θ α. We use the CMB Fisher matrices for the ΛCDM+N eff cosmology as in Section 3. We can then impose the α(z i )-prior on the redshift-binned likelihood function L(α, β; z i ) according to L(β) ∝ z i dα i z i L(α i , β; z i ) π(α 1 , . . . , α n ), where α i ≡ α(z i ) and π is the n-dimensional Gaussian prior with covariance matrix C α . The observed posterior distribution of α(z i ) could also be constructed by evaluating α(z i ) for each point in a given CMB Markov chain.

Constraints from Planned and Future Surveys
We will now show how well the phase shift can be measured in planned galaxy surveys. It is useful to first understand the parameter space α-β without imposing a prior on α. Both parameters affect the locations of the acoustic peaks and are therefore quite degenerate. We will use likelihood-based forecasts to ensure accuracy. We will confirm that the posterior distributions 9 of α and β are Gaussian, while the constraints on N eff derived from this parameterization are significantly non-Gaussian. This suggests that a Fisher matrix forecast in terms of α and β would be more reliable than one that starts directly from N eff . 7 We have explicitly checked that our template gives an unbiased measurement of β. In particular, we have verified that we reproduce β ≈ 0 for a cosmology with N eff = 0. 8 We also indirectly use the CMB data to constrain other cosmological parameters, in particular the scale factor at matter-radiation equality aeq, so that we can ignore any additional phase shifts not associated with N eff . 9 Since we assume flat priors for the parameters, we can identify the posteriors with the likelihoods.
We define the phase shift relative to the fiducial model with N eff = 3.046. The broadband spectrum for the fiducial model can be isolated by using the method in Appendix C or through the use of a fitting function along the lines of [52]. These methods generate the BAO spectrum O fid (k) and hence O(k) via (4.4). We compute the log-likelihood using the same noise and modeling as in the Fisher matrix (3.20).

Planned surveys
Forecasts for the one-and two-dimensional posteriors are shown in Fig. 9 for both BOSS and DESI. We see that for both surveys the posterior distributions are Gaussian. The best-fit Gaussian for BOSS and DESI has σ(β) = 1.3 and 0.47, respectively, which corresponds to a rejection of β = 0 at 77 % and 98 % confidence. Clearly, BOSS cannot exclude β = 0 (and hence N eff = 0) without any prior information from the CMB. Since the weakness of the constraint on β is driven by the degeneracy with α (see the left panel in Fig. 9), we expect to get significant improvements in the constraints on β after imposing a CMB prior on α. Inspection of the two-dimensional contours already shows that we will sizeably limit the range of β. The posterior distribution with the prior from Planck (CMB-S4) is shown in the right panel of Fig. 9. For BOSS, we find σ(β) = 0.76 (0.50) which implies that β > 0 at 81 % (95 %) confidence. Evidence for this signature of free-streaming neutrinos has been seen in existing data [54]. For DESI, we should find strong evidence for a phase shift with σ(β) = 0.30 (0.22) which excludes β = 0 at 3.5σ (4.6σ).
To translate these results into constraints on N eff , we use the relationship between β and N eff given in (4.6). This map is nonlinear over the measured range of β and we therefore anticipate the posterior to be non-Gaussian. The derived N eff -posteriors in Fig. 10 indeed show a highly  non-Gaussian distribution. As anticipated from the β-posterior for BOSS, the constraints on N eff are relatively weak without imposing a Planck prior on α.
We also see that the constraining power is significantly weaker at bounding large values of N eff than small ones. This asymmetry is simply a reflection of the fact that increasing N eff does not produce proportionally larger phases shifts. This asymmetry was also seen in the CMB constraints of Follin et al. [5], likely for the same reason. Recall that we have an upper limit on the phase shift of β < 2.45, which is saturated for N eff → ∞. In practice, this means that for N eff a ν ≈ 4.40, we will have an equal likelihood 10 for every value of N eff because they produce identical spectra. As a result, a flat prior on N eff (rather than β) will lead to ill-defined results because the integral ∞ dN eff L(N eff ) will diverge. On the other hand, for highly-significant detections of β > 0, a flat prior over any reasonable range of N eff will produce stable results. We are not quite in this regime with BOSS, which is why we will only quote constraints on β. Table 6 shows the projected constraints on β for a variety of planned surveys with and without priors from the CMB. We see that roughly a factor of three improvement can be achieved in spectroscopic surveys going from BOSS to Euclid. Both DESI and Euclid should have sufficient sensitivity to reach a more than 5σ exclusion of β = 0 when imposing a Planck prior. As before, galaxy clustering measurements in photometric surveys do not lead to competitive constraints as they are effectively two-dimensional on the relevant scales.  Table 6: Forecasted 1σ constraints on the amplitude of the phase shift β for current and future LSS experiments. We also show the constraints on β after imposing a redshift-dependent prior on the BAO parameter α from Planck and CMB-S4.

Future surveys
Given the robustness of the phase shift as a probe of light relics, a high-significance detection of the phase shift in LSS would be a valuable piece of cosmological information. We have seen that current and planned surveys can detect the phase shift, but are not expected to produce constraints on N eff that are competitive with those from the CMB. It is natural to ask if future surveys can reach this level of sensitivity.
Like the measurement of the BAO scale, the measurement of the phase requires large signalto-noise for 0.1 h Mpc −1 k 0.3 h Mpc −1 . As long as the number density is sufficiently large to keep the shot noise below cosmic variance, we gain primarily by increasing z max to achieve larger survey volumes. At larger levels of the shot noise, we only measure a few peak locations well which increases the degeneracy between α and β. Figure 11 shows results for a variety of possible  Figure 11: Future constraints on the amplitude of the phase shift β as a function of z max and N g , assuming f sky = 0.5. The dashed and dotted lines indicate that a CMB prior on α has been imposed using Planck and CMB-S4, respectively. The corresponding 1σ lower limit on N eff , which is N eff = 3.046 − σ − (N eff ), is indicated by the right axis.
survey configurations. As before, the constraints on β can be mapped into constraints on N eff using (4.6). We see that with 10 8 objects and z max > 3, we consistently obtain σ(N eff ) < 0.5 (1.0) with (without) a prior on α included.
To put these results into context, the measurement of Follin et al. [5] of N φ eff = 2.3 +1.1 −0.4 from the Planck TT spectrum is comparable to a survey with N g = 10 9 objects out to redshift z max = 3. Follin et al. also forecasted σ(N φ eff ) = 0.41 for Planck TT+TE+EE which is near the sensitivity of future LSS surveys when increasing the redshift range to z max = 6. Reaching this level of sensitivity will be extremely challenging with an optical survey, but could potentially be achieved with 21cm intensity mapping [27].

Comparison to Parameter-Based Approach
It is instructive to compare the results of our template-based forecasts to a more direct parameterbased approach to isolating the phase shift. In the parameter-based approach, we define two new parametersθ s andÑ eff that play the role of θ s and N eff in the BAO signal, but are taken to be independent of the same parameters in the CMB. We will then fix all remaining cosmological parameters in the BAO spectrum using the CMB, except ω c which we traded for a eq . As with our template extraction, holding a eq fixed ensures that the phase shift approaches a constant at large wavenumbers, whose value is determined byÑ eff . Beside measuring the phase of the BAO signal, the parameterÑ eff also contributes to the scale parameter α and could therefore be constrained by the standard BAO-scale measurement if all the other cosmological parameters are fixed to their Planck best-fit values. Introducing the additional parameterθ s gives enough freedom to remove this effect and any constraint onÑ eff must be coming from the phase shift alone. This is analogous to isolating the phase shift in the CMB by marginalizing over Y p or any other parameters that are degenerate with the N eff -induced change to the damping tail [3]. We will confirm this expectation in our forecasts.
Typically, the advantage of the parameter-based approach is that it is easy to implement. However, in this case, we found it more difficult to set up reliably. The phase shift ultimately controls the breaking of the degeneracy betweenθ s andÑ eff and, as we discussed in §4.1, P nw (k) must therefore be determined sufficiently accurately to not produce errors in this shift. To compute the likelihood directly, we must re-compute P nw (k) for every value of the cosmological parameters. Producing stable results for the BAO spectrum across a wide range of parameters can be very computationally expensive and technically challenging. Simpler and faster methods can work well near the fiducial cosmology (such as the use of a fitting function), but often produce noisy results as the parameters vary significantly and typically underestimate the likelihood as we depart from the fiducial cosmology (and, hence, overestimate the constraining power).
Despite the challenge presented by a parameter-based approach, it has the advantage that it should capture all of the cosmological information available. It is therefore useful to compare the results of the parameter-based and template-based approaches to see if the template is missing information. Fortunately, we will see that the posterior distributions ofÑ eff andθ s can be largely reproduced as a derived consequence of the template-based forecasts. From the previous subsections, we should anticipate that the posteriors forÑ eff andθ s will be non-Gaussian, and will therefore require the calculation of the likelihood forÑ eff andθ s directly and not only the Fisher  Figure 12: BOSS (left) and DESI (right) two-dimensional 1σ and 2σ contours forÑ eff andθ s , determined ('directly') from the likelihood for the BAO spectrum for each value of the parameters and derived ('from (α, β)') from the redshift-binned likelihood for α and β. We find good agreement between both methods, suggesting that the two-dimensional parameterization is capturing most of the relevant information. The dashed lines indicate the fiducial values.
matrix. We will follow the same approach as described in §4.2. Computing the full likelihood is quite involved, which is the reason why we will assume that the CMB data fixes the other cosmological parameters to their fiducial values, except forÑ eff andθ s .
Results of the likelihood analysis in terms of these parameters for both BOSS and DESI are shown in Fig. 12. We see that the results are similar, which establishes that our templates are capturing most of the information available in the BAO spectrum. This is an important observation because it allows us to simplify the analysis to a two-parameter template without much loss of information. In fact, the conclusion that these likelihoods are the same is not easily reproduced with any method of BAO extraction, but requires the robustness and stability of a method such as the one we use. Given instead our phase shift template, one can reliably compute Fisher matrices or likelihoods for α and β, and derive the implications for cosmological parameters from them. Future surveys such as DESI show somewhat larger differences between the two methods, which suggests that more information could potentially be extracted by using additional and/or alternative templates.
The doubling of cosmological parameters to treat the CMB and LSS independently, like in the case ofÑ eff andθ s , has useful conceptual advantages even if we derive constraints on these parameters from the posterior of α and β. Growing tensions between the CMB and certain low-z measurements have garnered much attention, but lack a compelling explanation. Measuringθ s andÑ eff in the BAO spectrum may provide a new perspective on this issue without the need for a CMB anchor.

Conclusions
Large-scale structure surveys are an untapped resource in the search for light relics of the hot big bang. The growing statistical power of these surveys will make them competitive with the CMB in terms of the constraints they will provide on a broad range of cosmological parameters. Moreover, the combination of CMB and LSS observations will allow powerful and robust tests of the physical laws that determined the structure and evolution of the early universe.
In this paper, we have explored the potential impact of LSS surveys on measurements of the parameter N eff . We have found that the dominant statistical impact of future surveys lies in the shape of the galaxy power spectrum. The distribution of dark matter in the universe is altered through the gravitational influence of the free-streaming radiation, leading to changes in the shape of the power spectrum that can be detected at high significance. A summary of the reach of selected planned and future surveys is given in Fig. 13. We see that BOSS and DESI can extend results significantly beyond the current CMB constraints. Futuristic surveys combined with a future CMB-S4 mission could achieve σ(N eff ) ∼ 0.015, which is close to reaching the target of ∆N eff = 0.027 at a significance of 2σ.
Future LSS surveys will also be able to detect the coherent shift in the peak locations of the BAO spectrum. This would be an intriguing measurement as this phase shift is a highly robust and unambiguous probe of light relics and the cosmic neutrino background [28]. Moreover, it is sensitive to extensions of ΛCDM without requiring the CMB as an anchor. Improved measurements of the phase shift may therefore play a useful role in elucidating apparent discrepancies between CMB and low-redshift measurements of the Hubble parameter H 0 [55].
Future CMB and LSS observations could have a significant impact on our understanding of fundamental physics. In this paper, we have argued that, in the case of N eff , these observations can play complimentary roles by both increasing the raw sensitivity and adding to the robustness of the measurement. We have also shown that the BAO spectrum carries more accessible cosmological information than only the acoustic scale. A broader exploration will likely reveal more targets that benefit from this complementarity.

A Forecasting CMB Constraints
Forecasting the sensitivities of future CMB observations is by now a standard exercise; see e.g. [18,[60][61][62]. For completeness, this appendix collects the basic elements of our CMB Fisher analysis, as well as the specifications of the CMB experiments that were used in our analysis.

A.1 Fisher Matrix
The Fisher matrix for CMB experiments can be written as The covariance matrix C XY l for each multipole l and X = ab, Y = cd, with a, b, c, d = T, E, B, is defined by where C X l are the theoretical CMB power spectra, N X l are the (Gaussian) noise spectra of a given experiment and f sky is the effective sky fraction that is used in the cosmological analysis. We employ perfectly delensed power spectra and omit the lensing convergence for simplicity as it is sufficient for our purposes. We however comment on the effects of these assumptions below. The noise power spectra are where ∆X = ∆T, ∆P are the map sensitivities for temperature and polarization spectra, respectively, and θ b is the beam width (taken to be the full width at half maximum). Note that we set N T E l ≡ 0 as we assume the noise in temperature and polarization to be uncorrelated. For a multi-frequency experiment, the noise spectrum (A.3) applies for each frequency channel separately. The effective noise after combining all channels is where N X,ν l are the noise power spectra for the separate frequency channels ν.

A.2 Experimental Specifications
Our specifications for the Planck satellite are collected in Table 7. The adopted configuration is the same as that used in the CMB-S4 Science Book [18]. For the low-l data, we use the unlensed TT spectrum with l min = 2, l max = 29 and f sky = 0.8. We do not include low-l polarization data, but instead impose a Gaussian prior on the optical depth, with σ(τ ) = 0.01. For the high-l data, we use the unlensed TT, TE, EE spectra with l min = 30, l max = 2500 and f sky = 0.44. Since the low-l and high-l modes are independent, we simply add the corresponding Fisher matrices.  Table 7: Specifications for the Planck-like experiment used in [62] and in the CMB-S4 Science Book [18]. The dashes in the first two columns for ∆P indicate that those frequency channels are not sensitive to polarization.
We parameterize future CMB experiments in terms of a single effective frequency with noise level ∆T , beam width θ b and sky fraction f sky . We will present constraints as a function of these three parameters. We take θ b = 3 , ∆T = 5 µK arcmin and f sky = 0.3 as the fiducial configuration of a CMB-S3-like experiment. For a representative CMB-S4 mission, we adopt the same configuration as in the CMB-S4 Science Book [18]: θ b = 2 , ∆T = 1 µK arcmin and f sky = 0.4. For both experiments, we use unlensed temperature and polarization spectra with l min = 30, l T max = 3000, l P max = 5000. We add the low-l Planck data as described above, include high-l Planck data with f sky = 0.3 and f sky = 0.2 for CMB-S3 and CMB-S4, respectively, and impose the same Gaussian prior on the optical depth τ as for Planck.
Unlike the CMB-S4 Science Book, we do not include delensing of the T-and E-modes. For N eff forecasts, this was shown to have a negligible impact [63], while using unlensed spectra overestimates the constraining power of the CMB by roughly 30 % for N eff +Y p . We are primarily interested in the improvement in parameters from adding LSS data, which should be robust to these relatively small differences. We also ignore the lensing convergence as it does not impact the constraints on these parameters.

A.3 Future Constraints
As a point of reference, we present constraints derived from CMB observations alone. In Table 8 Table 9: Forecasted sensitivities of Planck, CMB-S3 and CMB-S4 for the parameters of ΛCDM+Y p and ΛCDM+N eff +Y p .
and CMB-S4. In Table 9, we display how these constraints vary when we allow the primordial helium fraction Y p to be an additional free parameter. The differences in the forecasted sensitivities for Planck compared to the constraints published in [33] can be attributed entirely to the improvement in σ(τ ) which arises from the imposed prior on the optical depth τ . The forecast of N eff for CMB-S3 is a rough estimate and will be subject to the precise specifications of the respective experiment. While the precise design of CMB-S4 is also undetermined at this point, σ(N eff ) = 0.03 is a primary science target and is therefore more likely to be a reliable estimate of the expected performance.
In Figure 14, we demonstrate how the constraints on N eff depend on the sky fraction f sky , for three different values of θ b and fixed noise level ∆T = 1 µK arcmin. When varying the total sky fraction, we also appropriately change the contribution of the included high-l Planck data. In Figure 15, we illustrate the constraint on N eff as a function of the beam size θ b and the noise   level ∆T , for fixed sky fraction f sky = 0.4. Comparing Figure 15 to the equivalent figure in the CMB-S4 Science Book [18] (Fig. 22), we see that the difference between the two forecasts is ∆σ(N eff ) ≈ 0.002. This can be attributed to the effects of imperfect delensing and is completely negligible for our purposes.

B Forecasting LSS Constraints
In this appendix, we collect the specific information regarding the planned LSS surveys which we used in our Fisher and likelihood forecasts. We also provide the full set of constraints on all of the cosmological parameters and cosmologies that are studied in this paper.

B.1 Survey Specifications
Below, we provide the experimental specifications for the galaxy surveys used in our forecasts. We have slightly simplified the details compared to [25], for example. In particular, we group different types of tracers (e.g. luminous red galaxies, emission line galaxies or quasars) into a single effective number density and bias. We find our results to be fairly insensitive to many of these details and well approximated by a fixed number of objects distributed with a constant comoving number density over the same redshift range.
The employed parametrization of the spectroscopic redshift surveys BOSS, eBOSS, DESI and Euclid are provided in Tables 10 to 13. For eBOSS, we combine BOSS and the two eBOSS configurations of Table 11 into one survey neglecting the small overlap. We effectively treat each redshift bin with mean redshiftz as an independent three-dimensional survey. Our Fisher matrix is the sum of the Fisher matrices associated with each bin, F = z Fz. We translated the survey specifications used in [25] into three numbers per redshift bin: the linear galaxy bias b, the comoving number density of galaxiesn g and the bin volume V . This is sufficient to fully specify the Fisher matrix in each bin. The spherical bin volume is given by where f sky is the sky fraction, d c (z) is the comoving distance to redshift z, and z min =z − ∆z/2 and z max =z + ∆z/2 are the minimum and maximum redshift of the respective bin. Throughout this paper, we use redshift bins of width ∆z = 0.1.
For DES, we employ a survey area of Ω = 5000 deg 2 and a redshift coverage of 0.1 ≤ z ≤ 2.0, while we take 20 000 deg 2 and 0.1 ≤ z ≤ 3.5 for LSST. This results in approximately 1.4 × 10 8 and 5.9 × 10 8 objects in a total survey volume of about 24 h −3 Gpc 3 and 215 h −3 Gpc 3 for the two surveys, respectively. We neglect the spectroscopic redshift error as it is expected to be comparable to (or smaller than) the longitudinal damping scale Σ , but use a conservative root-mean-square photometric redshift error of σ z0 = 0.05 for both DES and LSST. Finally, we reiterate that, by   Table 11: Basic specifications for eBOSS derived from [25]. The redshift range is covered twice, first showing the survey covering Ω = 1500 deg 2 that will include emission line galaxies (resulting in roughly 3.8 × 10 5 objects in a volume of about 8.0 h −3 Gpc 3 ), and then the survey with Ω = 6000 deg 2 that will not (resulting in roughly 7.2 × 10 5 objects in a volume of about 32 h −3 Gpc 3 ).
considering galaxy clustering alone, we only take a subset of the cosmological observables into account, in particular for photometric surveys, and we therefore expect to underestimate the full power of these experiments.

B.2 Future Constraints
Using these specifications, we generated forecasts for all of the cosmological parameters discussed in the main text in combination with the Fisher matrices for Planck, CMB-S3 and CMB-S4. We include both P (k)-and BAO-forecasts for ΛCDM (Table 14), ΛCDM+N eff (Table 15), ΛCDM+Y p (Table 16), and ΛCDM+N eff +Y p (   bin. Since we marginalize over galaxy bias, our forecasts show no improvements beyond the CMB for ln(10 10 A s ) and τ . We therefore do not include these two parameters in the following tables.
Apart from the improvements in the constraints on N eff and Y p , which we already discussed in §3.2.1, we see that mainly ω b and ω c benefit from combining the discussed LSS surveys with CMB experiments. The sensitivities may be enhanced by factors of three (two) and more compared to Planck (CMB-S3). We note that the DESI specifications of Table 12 are slightly more optimistic overall than what was considered in [22] resulting in roughly the same BAO-forecasts and up to about 5 % better P (k)-forecasts.
Comparing our forecasts with the ones obtained from the BAO scale alone (combined with Planck), we see that the BOSS analysis for ΛCDM is nearly optimal, but improvements on the constraints of more than 10 % can be achieved in extended cosmologies. For instance, the constraints on ω b , n s and N eff improve by 3 % or more, and ω c in ΛCDM+Y p even by 12 %. For DESI, the obtained sensitivities can generally be increased by a larger amount, e.g. up to 15 % for ω b and n s in ΛCDM+N eff , and for ω c in ΛCDM+Y p .     Table 14, but for an extended ΛCDM+Y p cosmology. The constraints on ln(10 10 A s ) and τ are the same as in Table 9 for all combinations.  Table 17: As in Table 16, but for an extended ΛCDM+N eff +Y p cosmology.

C Broadband and Phase Shift Extraction
In this appendix, we describe our implementation of a robust broadband extraction method and the computation of the phase shift template.

C.1 Broadband Extraction
The split of the matter power spectrum into a broadband ('no-wiggle') part and an oscillatory ('wiggle') part, P (k) = P nw (k) + P w (k), is not uniquely defined, but depends on the method that is being used. In the following, we describe our method for extracting the broadband spectrum which is robust and stable over a very large parameter space.
Computationally it is easier to identify a bump over a smooth background than to find the zeros of oscillations on top of a smooth background. This suggests that it is convenient to sine transform the matter power spectrum to discrete real space where the oscillations map to a localized bump. We then remove this bump and inverse transform back to Fourier space.
An algorithm for the discrete spectral method was outlined in §A.1 of [64]. Concretely, the relevant steps of our implementation are: 8. Provide P nw (k) and P w (k): In order to cut off numerical noise at low and high wavenumbers, perform two cuts at k 1 and k 2 , where k 1 = 3 · 2 −n and the value of k 2 is found as the trough of |P (k) − P nw (k)|/P nw (k) following the smallest maximum (before the oscillation amplitude increases again due to the numerical artefacts intrinsic to the procedure). The reliably extracted no-wiggle spectrum P nw (k) is then valid for k ∈ [k 1 , k 2 ]. In practice, choose n and k max large enough initially so that k 1,2 are outside the range of wavenumbers of interest, e.g. those covered by a survey. The wiggle spectrum in this range is then given by P w (k) = P (k) − P nw (k).
Examples of the broadband extraction using this procedure are shown in Fig. 16. We see that the extraction method is unbiased, i.e. the resulting wiggle spectrum both oscillates around zero and asymptotes to zero for large wavenumbers. In addition, it is robust and stable over a large parameter space at small computation time (depending on n). Since the position of the first BAO peak is close to the peak of the matter power spectrum, it is sensitive to how exactly the baryonic bump is removed. However, we have checked that the computed constraints on cosmological parameters are insensitive to this uncertainty. The same holds for varying the parameters n and k max with fixed shifts in step 5 as long as k 1,2 are outside the range of wavenumbers of interest. Extracted broadband spectrum P nw (k) compared to the full power spectrum P (k) for N eff = 3.046 (left) and 10 (right). The spectra are rescaled by k 3/2 to exaggerate any oscillations. Bottom: Extracted BAO spectrum P w (k)/P nw (k) for N eff = 3.046 and 10 with linear (left) and logarithmic (right) k-axis.

C.2 Phase Shift Measurement
In the following, we describe our method for computing the phase shift template used in the likelihood analysis of Section 4.
First, we compute the BAO spectrum using CLASS and the broadband extraction method detailed above for a given value of N eff . In practice, we set the primordial helium fraction Y p to the fiducial value, but the final template is independent of this choice. As discussed in §4.1, we keep the time of matter-radiation equality fixed at its fiducial value by changing the dark matter density ω c according to where a ν is defined in (2.1). We then fit the following envelope function to the maxima of the absolute value of the BAO spectrum: The parameters A i , a i , κ i , with i = d, e, are fitting parameters, while k e is the location of the peak of P nw (k) and k d is the wavenumber associated with the mean squared diffusion distance. These fitting functions are motivated by the modeling in [5,65]. We define the 'undamped spectrum' as O(k) ≡ a(k) −1 P w (k)/P nw (k) .
Before we can measure the phase shift, we have to match the sound horizon at the drag epoch, r s , to that in the fiducial cosmology to remove the change to the BAO frequency induced by N eff . By rescaling the wavenumbers as k → r fid s /r s k, we fix r s k to the fiducial model for all wavenumbers k. For convenience, we also normalize the spectrum such that the amplitude of the fourth peak is the same as in the fiducial cosmology.
Finally, we can extract the phase shift as the shift of the peaks/troughs and zeros of O(k) relative to the fiducial cosmology, δk * = k * − k fid * . To obtain the template f (k), we sample 100 cosmologies with varying N eff ∈ [0, 3.33], 11 and define where β(N eff ) is the normalization introduced in (4.6). The bars in Fig. 8 indicate the locations of the peaks/troughs/zeros of the fiducial spectrum O(k) and their length shows the standard deviation in these measurements which is generally small.

D Convergence and Stability Tests
One of the motivations for including our full list of forecasts in Appendix B is to make the results reproducible. It is therefore also important that we explain how the numerical derivatives were computed in the Fisher matrix, including the employed step sizes. In this appendix, we provide this information and demonstrate that the step sizes are appropriate for the convergence and stability of our calculations.
The numerical derivatives in (3.3) and (3.20) are computed using a symmetric difference quotient or two-point stencil, f (θ) = [f (θ + h) − f (θ − h)]/(2h), with fiducial parameter value θ and absolute step size h. For each parameter, we choose the step sizes given in Table 18   In Figures 17 and 18, we show that our results are converged for both the P (k)-and BAOforecasts. The results in these figures (as in the rest of this paper) use CLASS with a high-accuracy setting. We have also checked that the forecasted constraints are converged when employing the standard accuracy setting, but note that the results are slightly less stable to changes away from these values. For the P (k)-forecasts, we see that a sufficiently small step size is needed, but a further decrease in the step size still leads to converged results. The BAO-forecasts, by contrast, show islands of convergence where performance decreases both when the step size is increased and when it is decreased. This behavior is more noticeable using the standard accuracy setting of CLASS, but likely reflects the fact that the BAO feature is itself a small effect and small step sizes are therefore more likely to produce effects comparable to numeric or modeling errors.  σ(N eff ) Figure 18: As in Figure 17, but for the BAO-forecasts of DESI.