Probing Dark Energy and Modifications of Gravity with Ground-Based Millimeter-Wavelength Line Intensity Mapping

Line intensity mapping (LIM) can provide a powerful means to constrain the theory of gravity and the nature of dark energy at low and high redshifts by mapping the large-scale structure (LSS) over many redshift epochs. In this paper, we investigate the potential of the next generation ground-based millimeter-wavelength LIM surveys in constraining several models beyond $\Lambda$CDM, involving either a dynamic dark energy component or modifications of the theory of gravity. Limiting ourselves to two-point clustering statistics, we consider the measurements of auto-spectra of several CO rotational lines (from J=2-1 to J=6-5) and the [CII] fine structure line in the redshift range of $0.25<z<12$. We consider different models beyond $\Lambda$CDM, each one with different signatures and peculiarities. Among them, we focus on Jordan-Brans-Dicke and axion-driven early dark energy models as examples of well-studied scalar-tensor theories acting at late and early times respectively. Additionally, we consider three phenomenological models based on an effective description of gravity at cosmological scales. We show that LIM surveys deployable within a decade (with $\sim 10^8$ spectrometer hours) have the potential to improve upon the current bounds on all considered models significantly. The level of improvements range from a factor of a few to an order of magnitude.


INTRODUCTION
Deciphering the physics responsible for the observed current accelerated expansion of the Universe -and whether it is driven by a cosmological constant or a dynamical dark energy (DE) component, or results from modifications to Einstein's theory of gravity (MG) -is one of the key open questions in modern cosmology.Shedding light on this question is the driving science case for several upcoming galaxy surveys, such as Euclid (Amendola et al. 2018) and DESI (Levi et al. 2019), which will map the large-scale structure (LSS) at redshifts z 2 by measuring shape and clustering of galaxies.In the standard cosmological model (among others), DE is a sub-dominant component of the Universe's energy density at redshifts z 1.On the other hand, models with an increased density of DE density at higher redshifts have been also considered in the literature.Examples include models with early-DE (EDE) component * emilio.bellini@ts.infn.itprior to recombination (e.g., Karwal & Kamionkowski (2016)) that has attracted significant attention in recent years in the context of the Hubble tension (e.g., Poulin et al. (2019); Ivanov et al. (2020); Hill et al. (2020); Simon et al. (2022)).In addition, there are also models exhibiting a "tracking behavior" at latetimes: the equation of state of DE w(z) tracks the dominant component of the universe in that particular epoch (w ≃ 0 during matter domination and w ≃ −1 for z 1), so that their energy density at high-redshifts is non-negligible (e.g., Urena-Lopez & Matos (2000); Bassett et al. (2002); Rakhi & Indulekha (2009)).Beyond constraining possible modifications to the theory of gravity that can explain the current accelerated expansion, testing general relativity (GR) on cosmological scales is of significant interest and is highly complementary to constraints from Solar System, and binary systems (Berti et al. 2015).Modifications to GR and a dynamic DE component can both modify the background expansion of the Universe and the growth of structure.
Therefore, statistical properties of the LSS enable us to constrain the nature of gravity and DE.
LIM is emerging as a powerful probe of the LSS.In contrast to galaxy surveys that map the LSS by resolving individual galaxies, LIM relies on measuring cumulative spectral-line emission from ensemble of galaxies or intergalactic medium (Pritchard & Loeb 2012;Kovetz et al. 2017).The measurements of spatial fluctuations in the intensity of the spectral lines together with their precisely measured frequencies enable us to measure the expansion history and growth of structure over extended redshift epochs and a wide range of scales, largely inaccessible to traditional galaxy surveys.In addition to 21cm line of hyperfine transition of neutral Hydrogen, LIM at mm/submm wavelengths, targeting spectral emission lines with rest-frame wavelengths in the far infrared (e.g., rotational lines of carbon monoxide CO(J→J-1) (Righi et al. 2008;Lidz et al. 2011;Breysse et al. 2014;Mashian et al. 2015;Li et al. 2016;Padmanabhan 2018) and fine structure of ionized carbon [CII]) (Gong et al. 2011a;Silva et al. 2015;Pullen et al. 2018;Padmanabhan 2019), has recently attracted a growing attention both in the context of galaxy/star formation (Gong et al. 2011b;Kovetz et al. 2020;Bernal & Kovetz 2022) and cosmology (Karkare & Bernal et al. 2021;Bernal et al. 2021;Schaan & White 2021a,b;Scott et al. 2022).The wide array of pathfinder experiments such as CCAT-Prime (Chapman et al. 2022), COMAP (Cleary et al. 2021), CONCERTO (Ade et al. 2020a), COPSS (Keating et al. 2016), EX-CLAIM (Ade et al. 2020b), mmIME (Keating et al. 2020), SPT-SLIM (Karkare et al. 2022b), TIME (Crites et al. 2014a), are paving the way for future wide-field surveys capable of providing high-precision constraints on cosmology by providing robust detection of the line intensity power spectrum and constraining the astrophysical model dependencies of the signal.Establishing the science cases for future surveys and determining optimized survey strategies to achieve relevant theoretical thresholds is of uttermost importance.
Building upon our previous works on investigating the potential of mm-wavelength LIM in constraining primordial non-Gaussianity (Moradinezhad Dizgah et al. 2019;Moradinezhad Dizgah & Keating 2019), and properties of neutrinos and light relics (Moradinezhad Dizgah et al. 2022a), in this paper, we study the prospects for constraining DE and MG.Contrary to upcoming spectroscopic galaxy surveys mapping the LSS at redshifts z 2, the broad redshift coverage of LIM observations enable us to map the background expansion and growth of structure both in the matter and DE domination epochs.As such, LIM surveys offer the possibility of testing the nature of DE or MG in the regime that has not been probed before.Furthermore, mapping the LSS on ultra large scales allows searching for horizon-scale imprints of possible modifications of GR (see e.g.Sailer et al. (2021); Berti et al. (2022); Casas et al. (2023); Scott et al. (2022) for recent studies on constraining beyond ΛCDM models with 21cm and mm-wavelength LIM).
For this work, we consider a ground-based survey capable of probing six rotational lines of CO and [CII] over the redshift range of 0 < z < 12, assuming the same survey characterizations as those used in Moradinezhad Dizgah et al. (2022a); Karkare et al. (2022a).We analyze several beyond ΛCDM models, falling into two broad categories of covariant and effective description models.More specifically, we consider the Jordan-Brans-Dicke (JBD) theory of gravity and EDE driven by an ultra-light axion Poulin et al. (2019) as examples of covariant models.On the other hand, within the effective description, we consider the Chevallier-Polarski-Linder (CPL) parameterization for the DE equation of state and two models for the effective description of luminal Horndeski theories.Lastly, in a hybrid approach, we consider an effective description of ShiftSymmetric Horndeski models imposing theoretical priors that ensure the consistency with the the underlying covariant theory.This selection of models enables us to explore phenomenologically rather broad regions of the parameter space of viable models to determine the potential of future LIM surveys in constraining beyond ΛCDM scenarios.
In our forecasts, we use limFisher package, which will be publicly released upon publication of this paper.This code is written in C language and is designed to perform forecasts for LIM power spectrum.It utilizes CLASS Boltzmann solver (Blas et al. 2011) 1 , as well as its extensions for beyond ΛCDM models; In particular, we use the publicly available CLASS EDE code (Hill et al. 2020) 2 for EDE model, and hi class code (Zumalacárregui et al. 2017;Bellini et al. 2018Bellini et al. , 2020) ) 3 for other models, including JBD, shift-symmetric Horndeski, and effective description of DE models.We only account for the imprints of the considered models at the linear level, as implemented in the Boltzmann codes, consistently neglecting non-linear scales.
The rest of the paper is organized as follows; in §2, we review the theoretical models that we will consider, summarizing the model parameters and current bounds, as well as the imprints on the LSS.In §3, we describe the model of the power spectrum of line intensity fluctuations that we use in our forecasts.In §4, we summarize the Fisher forecast methodology, and in §5, we describe the survey specifications assumed in the forecasts.In §6, we present the results of the forecasts, and finally in §7, we draw our conclusions.

REVIEW OF THEORETICAL MODELS BEYOND ΛCDM
In this section, we first lay out the theoretical framework for selecting the models to study, and then review the main features of each of the models that we consider and their connection to the general framework, as well as the imprints of the model on the background evolution and the growth of structures.
We perform our analysis following two complementary approaches.To study alternative theories of gravity, it is common either to specify a covariant Lagrangian, in the following referred to as covariant approach or to build an effective description that encapsulates the effects of the additional degree of freedom on the evolution of perturbations on a FLRW background.The two approaches offer different advantages, and we use both to explore a meaningful portion of the parameter space of viable cosmological models.
Our starting point for the covariant approach is the following Lagrangian where L m is the matter Lagrangian, G 2, 3 are functions of a scalar field φ, and its canonical kinetic term X ≡ −φ ;µ φ ;µ /2, and G 4 is a function of the scalar field alone.This describes the sub-class of Horndeski models (Horndeski 1974;Deffayet et al. 2009;Kobayashi et al. 2011) with the additional requirement that the speed of gravitational waves (GW) has to be equal to the speed light (Baker et al. 2017;Creminelli & Vernizzi 2017;Ezquiaga & Zumalacárregui 2017;Sakstein & Jain 2017).This is why the G 4 function depends only on the scalar field and not on its derivatives; a derivative coupling, allowed in the full Horndeski framework, would modify the speed of gravitational waves.Our choice of this subclass of Horndeski models is driven by the nearsimultaneous observation of the GW signal (GW17081) and the gamma-ray burst (GRB 170817A) emitted from the merger of a binary neutron stars system, which established that the GW speed has to be equal to the speed of light within one part in 10 15 (Abbott et al. 2017).Upon fixing the functional form of G 2, 3, 4 in Eq. ( 1) and setting the initial conditions for the scalar field, it is possible to solve the background expansion history and the perturbations consistently.However, for non-standard choices of G 2, 3, 4 , it is difficult to predict the evolution of the background and the perturbations.
An alternative approach that is not specific to a particular (covariant) Lagrangian is to utilize the "effective description" of dark energy, often denoted as Effective Field Theory for DE (EFTofDE).Here the information content of Eq. ( 1) can be compressed into a few meaningful functions of time (Gubitosi et al. 2013;Gleyzes et al. 2013;Bellini & Sawicki 2014).Up to linear order in perturbation theory (PT), the dynamics of Eq. ( 1) can be fully described in terms of four 4 functions of time 5 where w ≡ p DE /ρ DE is the equation of state of DE that contributes to the expansion history of the universe, while the remaining three functions affect only the evolution of the perturbations.In particular, α K is the standard kinetic term of perfect fluids, α B is the kinetic term arising from the coupling of the derivatives of the scalar field with the ones of the metric, and M 2 * is an effective Planck mass squared characterizing the change in the strength of gravity6 .One of the key advantages of this approach is that each one of these functions has a clear physical meaning, and they describe particular classes of models.In this scenario, one has the freedom of parameterizing these functions of time, keeping in mind that in ΛCDM they reduce to The two approaches to modeling alternative theories of gravity are to some extent equivalent, but each one has its own peculiarities and is preferable in specific situations.The covariant approach is fully consistent, i.e. for a single theory it is possible to describe gravitational interactions in any regime.Still, it is non-trivial to jump from one model to the other, and one needs some clever intuition to choose G 2, 3, 4 .On the other hand, the effective description is ideal for investigating general classes of models at once.It is closer to observations but can lead to non-physical theories.Indeed, while any choice 4 Note that here we are considering a sub-class of Horndeski.The full Horndeski theory requires five functions.
The derivative of the effective Planck mass M 2 * is commonly written as α M ≡ dlnM 2 * /dlna.
of the G 2, 3, 4 functions has a corresponding time evolution of {w, α K , α B , M 2 * }, it is not guaranteed that any time evolution of {w, α K , α B , M 2 * } is the representation of a particular choice of G 2, 3, 4 .The choice of approach depends on the needs of the analysis.

Jordan-Brans-Dicke
First, we study the constraints on the Jordan-Brans-Dicke (JBD) theory of gravity (Brans & Dicke 1961), the workhorse of modified gravity models, extensively studied in the literature, e.g.Lima & Ferreira (2016); Umiltà et al. (2015); Ballardini et al. (2016); Joudaki et al. (2022).The JBD theory adds a new scalar d.o.f. with respect to GR, that is non-minimally coupled to the metric.Following the covariant approach, JBD is recovered by specifying where V (φ) is a generic potential and ω BD is the coupling constant of the theory.The GR limit is recovered when ω BD → ∞.Therefore, the accelerated expansion of the Universe can be obtained for a sufficiently flat potential and a slow-rolling scalar field.We restrict our analysis to a constant potential, i.e.V (φ) = Λ (see, e.g., Ballardini et al. (2016) for constraints on power-law potential).Since Λ can be obtained by imposing that in a flat universe, the sum of the fractional densities has to be one, the free parameters of the theory are {ω BD , IC}, where IC stands for initial conditions of the scalar field.The initial value for the scalar field's first derivative is irrelevant since the scalar field evolves quickly to an attractor solution irrespective of the value of the first derivative.It is possible to recast the scalar field as an effective Newton's constant using Then, the only two parameters that will be varied are The strongest bounds on JBD gravity come from astrophysical probes.In particular, the stellar triple system PSR J0337+1715, where the inner pulsar-white dwarf binary is in orbit with another dwarf, was shown to provide the strongest bound on the market, i.e. ω BD > 1.4 × 10 5 (95% CL) (Archibald et al. 2018;Voisin et al. 2020).The current cosmological constraints are much weaker, i.e. ω BD > 2230 and Geff = 0.996 ± 0.029 (95% CL), from KiDSx2dFLenS, BOSS full shape, Planck18 and Pantheon data (Joudaki et al. 2022).This is because the JBD model has no screening mechanism, which implies that the fifth force introduced by the scalar field at cosmological scales is not suppressed at small scales.Therefore, in the quest for the ultimate theory of gravity, JBD becomes less interesting than other competitors.However, it is still possible to use JBD to infer when cosmological probes will have the same constraining power as astrophysical probes.Furthermore, the JBD gravity is one of the few models that has survived the recent constraints from the neutron star merger that constrained the speed of gravitational waves.It is also worth noting that the constraints on the JBD gravity can, to some extent, give insights into the long-wavelength features of the generalized scalar-tensor theories (Avilez & Skordis 2014).
At the background level in JBD gravity, the scalar field begins to evolve during the radiation-dominated era as a power-law of the scale factor.The expansion is slowed down during matter domination, and at latetimes, the constant gravitational potential of the scalar field sources the accelerated expansion of the universe.The scalar field evolution depends on the value of ω BD ; the larger it is, the more frozen is it.Therefore, we expect strong degeneracies between the Geff , the presentday Hubble parameter, H 0 , the total matter overdensity, Ω m , the amplitude of fluctuations S 8 = σ 8 √ Ω m , as well as a relatively mild correlation with the total sum of neutrino masses M ν (see Joudaki et al. (2022) for recent analysis).

Early Dark Energy driven by ultra-light axions
The models in which a new component acting as DE dominates the energy density of the universe briefly before recombination and rapidly decays after the last scattering have recently attracted attention in the literature as a solution to alleviate the Hubble tension between the early and late universe, as proposed originally in Karwal & Kamionkowski (2016).An EDE component increases the expansion rate, thus the value of H 0 , and acts to shrink the sound horizon at the last scattering.We consider the phenomenological model of EDE proposed in Poulin et al. (2019), which has been used in several recent analyses of the cosmic microwave background (CMB) and LSS data (see, e.g., Ivanov et al. (2020); Hill et al. (2020); Sailer et al. (2021)) and got the silver medal in the H 0 Olympics (Schöneberg et al. 2022).This model envisages an ultra-light axion, much lighter than fuzzy dark matter axions and displaced from the minimum of its potential at early times, as an EDE component, and can be described in the covariant approach with 9) Here, m is the axion mass, f is the axion decay constant, n is an exponent that generalizes the well-motivated axion-like potential (recovered by setting n = 1) and controls the rate of dilution after the field becomes dynamical, and Λ is a cosmological constant driving the late-times accelerated expansion of the universe.Given that the dynamics are shown to be relatively insensitive to changes of 2 n 4.5 (Smith et al. 2020;Agrawal et al. 2019), as done in previous analyses of Hill et al. (2020) and Schöneberg et al. (2022), we fix n = 3.Early on, due to Hubble friction, the axion stays in its displaced position and acts as an additional DE component.Once the Hubble parameter drops below the axion's mass (H m), the field starts rolling down the potential and oscillates at its minimum.To have this EDE component dissipate fast and become a subdominant energy component, near the minimum, the potential should be steeper than quadratic.The physics of this model is, therefore, governed by three parameters, where f EDE ≡ max(ρ EDE /3M 2 pl H 2 ) t is the maximal fractional contribution of the EDE density to the total energy density of the universe, z c is the critical redshift at which this maximum is reached, and θ i = φ i /f is the normalized initial value of the scalar field such that −π < θ i < +π7 .Once f EDE and z c are fixed, the role of θ i is to set the dynamics of perturbations right around z c via EDE sound speed, c s .The ΛCDM limit is recovered by setting the axion mass to zero in Eq. ( 9).This class of EDE models has been recently constrained using various combinations of datasets (e.g., Poulin et al. (2019); Smith et al. (2020); Hill et al. (2020); Schöneberg et al. (2022); Simon et al. (2022)).We consider the results from Simon et al. (2022), and use the best-fit values as our fiducial model and their quoted 1σ errors as a baseline for comparison with our results, We note that the above constraints are obtained from joint analysis of Planck (temperature, polarization, and lensing) data and BOSS full-shape galaxy power spectrum, imposing prior on h 0 from SH0ES local distanceladder data.Since our forecasts only include Planck priors (from temperature and polarizations data) on ΛCDM and not on the EDE model, the comparison of LIM results with current constraints should be considered only a first exploration of the potential of LIM in constraining EDE.
To illustrate the potential of LIM data in constraining the EDE model, it is helpful to review the imprints of an EDE component on the clustering of matter and the biased tracers.The faster expansion rate before recombination in the EDE model implies a larger acoustic sound horizon at drag time, which shifts the peaks of Baryon Acoustic Oscillations (BAO) to smaller wavenumbers compared to the ΛCDM.The presence of an EDE field also affects the dynamics of the perturbations within the comoving horizon during the epochs in which EDE is relevant, slightly suppressing the growth of perturbations around the recombination.To have an EDE model consistent with the primary CMB power spectra, some of the ΛCDM parameters must shift; Ω c is increased to compensate for the less efficient growth of perturbations, and n s is increased since the EDE is only relevant in a short time period and its effect on perturbations is scaledependent.These parameter shifts alter the matter power spectrum, the effect (especially on smaller scales) being more prominent at higher redshifts.This redshift dependence is also encoded in the growth rate of the structure, which is measured via the redshift-space distortions (RSD) of galaxy clustering (Hill et al. 2020).In addition to modifications of the matter power spectrum, the halo/galaxy abundance and clustering can also be altered in the EDE scenario.Using N-body simulations, Klypin et al. (2021) showed that the halo mass function in the EDE model is enhanced compared to ΛCDM model, i.e., the EDE model predicts more halos, with the enhancement being substantially more significant at higher redshifts and higher mass halos.Furthermore, halos of a given mass were shown to form earlier and have higher concentrations8 .
Apart from mapping a larger number of modes and hence having a lower cosmic variance, access to the wide range of spatial scales and redshifts makes LIM particularly powerful for constraining the EDE scenario.In terms of the scales, measuring the line power spectrum from horizon to nonlinear scales allows for constraints on the scale-dependent effect of the EDE (Hill et al. 2020).Access to large scales is crucial for models of EDE with lower critical redshifts since the growth suppression is pushed to later times, and thus affects correspondingly larger scales.For example, the EDE that peaks after recombination (not a viable resolution to Hubble tension) suppresses the linear matter power spectrum in a scale-independent way on small scales and can only be probed via its effect on large scales (Sailer et al. 2021).Since the nonlinear gravitational evolution tends to reduce the difference between the EDE and ΛCDM models, probing higher redshifts is advantageous.Moreover, being sensitive to the halo mass function, the line signal picks up the imprint of the EDE on the halo abundance, which was shown to be substantial at high redshifts (Klypin et al. 2021).

Standard CPL parametrization of Dark Energy
A simple and widely used description for DE models is the CPL parametrization (Chevallier & Polarski 2001;Linder 2003).It falls in the effective description category, and it assumes the presence of a DE fluid with the equation of state where w 0 is the value of the DE equation of state today, and w a modulates its time dependence, with a being the scale factor.This fluid is a minimal extension of the standard cosmological model, which is recovered by setting w 0 = −1 and w a = 0.The role of w is only to modify the expansion rate of the universe, and it (mildly) affects the evolution of the perturbations only through a modified Hubble parameter (in the next Sections, we will present present models that have intrinsic additional freedom for the perturbations).
Recent constraints from CMB, BAO, and Supernovae (Aghanim et al. 2020) show that these parameters are compatible with the standard ΛCDM model, i.e.
In our forecasts we set the fiducial values of w 0 = −1 and w a = 0 corresponding to ΛCDM.LIM can provide tight constraints on the CPL model by measuring the expansion history over a wide redshift range, thus, is expected to tightly constrain redshift-dependence of DE equation of state, i.e., the w a parameter.

Shift-symmetric Horndeski Models
To build a very general, yet tractable, MG theory from Eq. ( 1), it is possible to assume that the scalar field is slowly evolving on a cosmological background and Taylor expand the G 2, 3, 4 functions for small X.
To further compress the parameter space of the theory, we place constraints only on shift-symmetric models, i.e. models invariant under φ → φ + c (where c is some constant).Shift-symmetric theories are interesting because radiative corrections are parametrically suppressed around (quasi) de Sitter backgrounds (Luty et al. 2003;Nicolis et al. 2009).Following the covariant approach, a shift-symmetric version of Eq. ( 1) expanded up to quadratic order in X looks like where {c 01 , c 02 , d 01 , d 02 } are parameters of the theory, and commonly it is assumed that Λ 4 2 = M 2 Pl H 2 0 and Λ 3 3 = M Pl H 2 0 to ensure that all interactions have O(1) contributions to present-day cosmological background evolution.This theory is a rather minimal, yet rich four-parameter extension of ΛCDM.However, instead of constraining the parameters of this theory, we follow an alternative and more powerful approach proposed in García-García et al. (2020); Traykova et al. (2021), in which the parameter space of Eqs.(19-21) is projected into the effective description 9 .In other words, they investigate how the effective description functions evolve with time if they are obtained consistently from Eqs. (19-21), instead of parameterizing them directly (as in the spirit of the effective description formalism).This provides theoretical priors on the effective description functions that ensure we are dealing with physical theories.The theoretical priors have to be interpreted as the region in the parameter space where a "covariant theory" will likely fall.A possible tension between data driven constraints and theoretical priors would not indicate the data have preference on non-physical theories, but rather that only few of the physical theories are in agreement with data.
As anticipated, following the results of Traykova et al. (2021), we can safely ignore Eqs.(19-21), and study their García-García et al. ( 2020) is an in-depth treatment of quintessence models, while Traykova et al. (2021) shows results for the model considered here.
proposed parameters, i.e., {w 0 , w a , αB , m}, which come from simple parametrizations of the effective description functions It is clear that αB represents the magnitude of α B today, the parameter m modifies its dependence on the Hubble parameter H.The ΛCDM limit is obtained fixing w 0 = −1, w a = 0 and αB = 0, in which case the value of m becomes irrelevant.
Combining observational constraints from CMB, BAO, RSD and Supernovae Type IA data with theoretical priors, and assuming that the scalar field is entirely responsible for the late time acceleration (case Λ = 0 in Traykova et al. ( 2021))10 , in Traykova et al. (2021) they obtained the following constraints on the four model parameters which we will use to set the fiducial values for our forecasts and as a comparison point for results.
As already mentioned, one of the advantages of the effective description is that different physical effects can be encapsulated by independent functions of time.These functions can thus be considered as orthogonal and treated separately.As described in §2.3, the DE equation of state, w, affects primarily the expansion rate of the universe, which in turn modifies the growth of structures.Here, at the level of linear perturbations we have one more ingredient, which does not modify the expansion history, but introduces new scale dependence in the matter power spectrum.In particular, its amplitude is responsible for lowering the amplitude of the power spectrum at large scales (k < 10 −3 h/Mpc) and enhancing it at small scales11 .This effect is redshift dependent, being less evident at higher redshifts, and in general depends on the time evolution of α B (regulated by the parameter m).

Effective Description of Luminal Horndeski Models
The last class of models we will investigate belongs to the effective description class of luminal Horndeski theories.Starting from Eq. ( 2), we parameterize the four EFT functions as The kineticity α K is largely unconstrained by cosmological data (Bellini et al. 2016), due to the building blocks of α K being those of simple DE models, with at most first derivatives in the scalar field Lagrangian.As a consequence, in the sub-horizon limit, where terms with higher number of spatial derivatives become dominant, the DE fluid becomes almost negligible.We additionally fix α B in our forecast.Taking into account that its effect has been explored in §2.4, and to avoid excessively enlarging the parameter space, we focus on the remaining d.o.f. of luminal Horndeski theories, i.e., M 2 * .It is also important to remember that by construction the effective description allows to separate the different physical effects, which makes it possible to study the effects of one single function at time and then draw conclusions on the general model12 .
We keep α K and α B fixed to the evolution of the scale factor while varying the parameters related to the expansion history and Planck mass, as is similarly done in Raveri et al. (2017), where This parametrization is a power law expansion around the present epoch up to first-order in the expansion history and second-order in the effective Planck mass.It allows a smooth GR limit (i.e.where Eq. ( 33) can be expressed as {−1, 0, 0, 0, 0}), and can have rich dynamics.While the vast majority of the models proposed in the literature as late-times DE/MG kick in at z ∼ 1, we explore the potential of LIM experiments spanning a wide range of redshifts.Thus, we choose two ad-hoc models (defined in Table 1), designed to have signatures at high-redshift, i.e., 2 z 10.In particular, the effective Planck mass of model1 increases smoothly with time, and its deviations w.r.t. 1 are ∼ 10% at z = 10.model2 is even more extreme, with deviations ∼ 25% at z = 10, but also a feature at z ≃ 1 after which M 2 * starts decreasing.
Two type of information can be derived from forecasts on this model: (i) maximize the constraining power of LIM experiments on DE/MG parameters, and (ii) understand to what extent a time varying effective Planck mass can be constrained in this regime.The effective Planck mass M 2 * , together with its derivative α M , enhances the growth of structures at large scales (k ≃ 10 −3 Mpc −1 h), while keeping the amplitude of the matter power spectrum almost constant at smaller scales.The effect is redshift dependent, and only models with high-redshift signatures will have the matter power spectrum modified at high-redshift.

THE POWER SPECTRUM OF LINE INTENSITY FLUCTUATIONS
We use a simple model for the line intensity signal,in which line intensity is calculated via scaling relationships between host halo mass and the galaxies CO or [CII] luminosities.The details of the model and relevant references to previous works on modeling line power spectra can be found e.g., in Moradinezhad Dizgah et al. (2022a).Here we just give the final form of the line power spectrum used in our forecasts.
The total observed power spectrum of fluctuations in a given line has three contributions: clustering, shot and instrumental noise, Here µ is the cosine of the angle between wavevector k and the line-of-sight direction.The clustering contribution (typically in units of µK 2 Mpc −3 h 3 ) is anisotropic due to the RSD and the Alcock-Paczynski (AP) effect.Assuming linear bias, on large scales the line power spectrum is given by where T is the mean brightness temperature of the line, b line is the linear bias of the line, P 0 is the linear matter power spectrum, H is the Hubble parameter, D A is the comoving angular diameter distance, β = f /b line , and f = dlnD/dlna is the logarithmic growth rate of structure with D being the growth factor and a the scale factor.The wavenumbers and angles in the true and reference cosmologies (characterizing the AP effect) are related as The exponent σ v in Eq. ( 35) includes both finger-of-god (FoG) effect and the intrinsic line width of individual emitters, both of which result in smearing of power on small scales.
The mean brightness temperature of the line (typically in units of µK) at redshift z is given by where M min and M max are the minimum and maximum masses of the halos that host galaxies emitting in a given line, c is the speed of light, k B is the Boltzmann factor, ν obs is the observed frequency of the line, dn/dM is the halo mass function, L line is the specific luminosity of the line-emitting galaxy located in a halo of mass M at redshift z, and D L is the luminosity distance.The terms dl/dθ and dl/dν reflect the conversion from units of comoving lengths, l, to those of the observed specific intensity: frequency, ν, and angular size, θ.The parameter p 1,σ accounts for scatter in the relations between SFR and specific luminosity with halo mass (Li et al. 2016;Keating et al. 2016).
In the Poisson limit, the shot-noise contribution arising from the discrete nature of line emitters is given by We note that to obtain more realistic forecasts, the line power spectrum model used in this work should be improved to include nonlinearities in matter distribution, RSD and line biasing.Furthermore, one should account for deviations from Poisson shot noise (Moradinezhad Dizgah et al. 2022b) and marginalize over these corrections in addition to the full set of bias parameters to obtain cosmological parameters.Including more free parameters allow us to consider smaller scales in the forecasts, but at a fixed small-scale cutoff, would result weakening of the constraints compared to the linear model (see e.g., Sailer et al. (2021)).We leave a forecasts using a more complete model to future work.

ANALYSIS METHODOLOGY
We use Fisher forecast formalism to determine the constraining power of a ground-based mm-wavelength intensity mapping survey, probing CO rotational lines and [CII] line.
Our analysis framework is the same as what we used in our previous paper, Moradinezhad Dizgah et al. (2022a).Therefore, we only summarize the main features here and refer the reader to that work for further details.
Binning the survey in redshift bins of 0.1 dex, and neglecting correlations between redshift bins, for each emission line x, the total Fisher matrix is the sum of Fisher matrices of individual redshift bins, with the Fisher matrix in the i th redshift bin given by where λ are the model parameters that are varied, V i is the volume of i th redshift bin extended between z min and z max and is proportional to sky coverage of the survey, f sky , while Var −1 [P x line ] is the variance of the line power spectrum.We include the noise due to interloper lines from lower and higher redshifts in the variance of the line of interest, these interloper contributionstypically dominated by low-redshift CO -are assumed to be removable in the power spectrum domain (e.g.Cheng et al. 2016).
We set the fiducial ΛCDM parameter values to those from Planck 2018 data 13 (Aghanim et al. 2020): ln(10 10 A s ) = 3.0447, n s = 0.96589, h = 0.6732, Ω cdm = 0.2654, Ω b = 0.04945.The fiducial values of free parameters of each of the theoretical models considered are described in section 2. We assume three degenerate massive neutrino species for all forecasts and fix their total mass to M ν = 0.06 MeV.For the finger-of-god contribution, we set the fiducial value of σ FOG,0 = 250 km/s.We fix the line bias and mean brightness temperature to their theoretical values, accounting for their cosmology dependence.For the combined results of LIM and Planck, we use the available Planck 2018 Monte Carlo Markov chain for ΛCDM 14 to compute the parameter covariances and the Fisher matrix, marginalizing over optical depth and assume the CMB and LIM Fisher M , c (2)  1.
In addition to marginalizing over nonlinear biases and corrections to the Poisson shot noise, to obtain a more realistic forecast, the line mean brightness temperature, T , should also be marginalized over.While in general considering a more extended parameter space degrades the constraints at a fixed small-scale cutoff, it also can offer advantages; At linear level, the linear bias and T are highly degenerate, with RSD only partially breaking the perfect degeneracy between these two parameters in isotropic power spectrum (Barkana & Loeb 2005;McQuinn et al. 2006;Sailer et al. 2022).Considering nonlinearities can further alleviate this degeneracy.As mentioned earlier, we leave consideration of more extended model of LIM power spectrum to a future work.

SURVEY SPECIFICATIONS
In this section, we briefly summarize the survey specifications that we assume in our forecasts and refer to Moradinezhad Dizgah et al. (2022a) for a detailed de-scription of the hypothetical next-generation mm-wave LIM survey and how those specifications translate into sensitivity estimates used in our forecasts.
The instrument noise contributing to the measured power spectrum, P N , is given in terms of per-voxel noise, σ N , of the original image cube as V vox is the volume of individual voxels within the image cube, which can be further expressed as Here Ω vox and δν are the solid angle and bandwidth covered by a single voxel, respectively.The per-voxel noise is parameterized by spectrometer hours, τ sh , as a proxy for the "level-of-effort" of an experiment.If the survey area is Ω s , such that the number of independent pointings is given by Ω s /Ω vox , then we can express Eq. ( 41) as where σ NET is the noise-equivalent temperature (NET) of the detector (units of K• √ s), and η opt is the optical efficiency of the instrument.In computing σ NET , we assume that the spectrometer performance is close to background-limited, with the primary contributors to the instrument noise being atmospheric emission, with secondary contributions from the telescope emissivity, and various astrophysical backgrounds (including Galactic dust and the CMB).Further details on the simulated spectrometer performance can be found in Moradinezhad Dizgah et al. (2022a).
We define an effective instrumental noise, to account for attenuation of the signal due to the finite resolution of the instrument (α max ) and the finite redshift and angular coverage of the survey (α min ).
The two attenuation factors are defined in terms of the largest and smallest recoverable modes in parallel and perpendicular to line-of-sight directions.The former is determined by the redshift and angular coverage of the survey, while the latter is set by the spectral resolution of the instrument and the diameter of the aperture (see Moradinezhad Dizgah & Keating (2019); Moradinezhad Dizgah et al. (2022a) for the full description and explicit expressions).
For our analysis, we consider ground-based observations from an accessible observing site with excellent proven mm-wave observing conditions, such as the South Pole or the Atacama desert.Each spectrometer is assumed to be sensitive to the entire frequency range, i.e., each spectrometer can measure all individual spectral channels simultaneously.
We vary spectrometer hours over a wide range, starting with first-detection experiments and extending to larger-scale surveys that could be fielded in the next ten years.We set the lower bound comparable to the current-generation instruments, which feature ∼ 50 spectrometers (Crites et al. 2014b) and are capable of completing surveys of order 10 5 spectrometer hours.As described in Moradinezhad Dizgah et al. ( 2022a), for each model, we calculate the parameter constraints while varying f sky and spectrometer hours to determine the optimal survey strategy.We compute the constraints both with and without the addition of the Planck priors on ΛCDM parameters (only the former ones are shown in §6).
For our analysis, we tabulate and report the results with and without interlopers -the latter representing the optimistic view that the interloper emission can be (near-)perfectly removed by leveraging the information contained with multiple overlapping lines (e.g.Cheng et al. 2020), and the former representing the pessimistic view that no suppression of interloper "noise" would be possible.

RESULTS
In this section, we present the forecasted constraints for the models discussed in §2.For each model, we first describe the results of survey optimization by showing the marginalized constraints as a function of spectrometer hours and show the corresponding sky coverage required to achieve those constraints.In these plots, we show also a gray horizontal band, representing current constraints on each parameter.It is important to stress here that, when comparing the current constraints on the same parameter but for different models (e.g., w 0 from §2.3 and §2.4) these can be different for two reason: (i) changes in the number of parameters changes, but most importantly (ii) the quoted current constraints come from different analysis which use different datasets.To illustrate the parameter degeneracies and the redshift-dependence of the constraints for each model, we show the constraints on each parameter as a function of the maximum redshift, and the marginalized constraints on parameter pairs, assuming fixed spectrometer hours and sky coverage.Before discussing the results of individual models, let us summarize the main trends that are evident for all of them.Considering the dependence on survey parameters, the parameter constraint improves with increasing spectrometer hours and sky coverage.Accounting for interlopers degrades the constraints.For some parameters (which we discuss below), the degradation is most significant at higher values of spectrometer hours, causing saturation of the constraining power.This implies that in those cases, the noise from the interloper dominates the thermal noise determined by the number of spectrometer hours.Furthermore, the interloper noise pushes the required sky coverage to higher values.For all models, the redshifts corresponding to the peak of star formation rate at which the amplitude of the line power spectrum is the largest (z ∼ 2) provide considerable constraining power.However, the information content of the signal at different redshifts differs in the models considered, which we will discuss in what follows.
In table 2 and 3 of Appendix A, we show the 1-σ marginalized constraints for all the considered models for five stages of LIM surveys defined in terms of an approximate number of spectrometer hours.For each model and each stage, we quote the constraints with (top row) and without (bottom row) interlopes.

Jordan-Brans-Dicke
In Fig. 1, we show the 1-σ constraints on parameters of the JBD model as a function of spectrometer hours, assuming the fiducial value of ω fid BD = 800 .The solid (dashed) blue lines show the constraints, including (neglecting) the interloper noise, and the red lines correspond to the required sky coverage to obtain those constraints.The plateau in the red lines corresponds to the maximum sky coverage of f sky ≃ 0.5 that we have imposed.The grey bands are the current constraint from Joudaki et al. (2022).
For the comparison with the current bounds, it should be kept in mind that the JBD constraints (in particular ω −1 BD ) depends on the choice of fiducial values.Our choice of ω fid BD = 800 is comparable to the peak of the 1d posterior of w −1 BD in Joudaki et al. (2022).While only upper bounds on ω BD are reported in Joudaki et al. (2022), the gray band we show here loosely corresponds to σ(ω −1 BD ) ∼ 0.002.For 10 8 spectrometer hours, the extra noise due to interlopers degrades the constraints on ω −1 BD and Geff by less than a factor of 2. At higher values of spectrometer hours, not removing the interlopers causes saturation of the constraining power as the interloper noise dominates the thermal noise.Furthermore, the unsubtracted interloper noise increases the required sky coverage by nearly a factor of 2; without interlopers the maximum sky coverage of 50% is saturated at ∼ 2 × 10 8 spectrometer hours, while with the interloper noise, it is reached at ∼ 5 × 10 7 .Fig. 2 shows the 1-and 2-σ marginalized constraints on parameter pairs to illustrate parameter degeneracies.Blue contours are from LIM only, while the red ones include Planck priors.As expected from the discussion in §2, for LSS probes, there are strong degeneracies between Geff , h, and Ω c (and to a lesser extent with Ω b ) due to their similar effects on expansion rate and growth of structure.Imposing Planck priors improves the LIM-only results, both by tightening the constraints on ΛCDM parameters and by alleviating parameter degeneracies.We caution that while here we imposed Planck ΛCDM priors, a more careful joint analysis of CMB and LIM, assuming the same model, is necessary to quantify the information content of the individual and combined probes more accurately.
Fig. 3 shows the constraints as a function of maximum redshift z max (top panel) and as a function of the center of each redshift bin (bottom panel), assuming a fiducial value of ω BD = 800 and imposing Planck ΛCDM priors.For both parameters, the constraints improve with increasing z max up to z max ∼ 2.5, after which a near plateau is reached.At z ∼ 6 when the [CII] signal also contributes, we see a decrease of ∼ 15% on σ(w −1 BD ), while σ( Geff ) is only lowered by 3%.At high redshifts, the survey we have considered here only access to high-J rotational lines of CO, which are expected to be much fainter (in brightness temperature units) then their low-J counterparts.As a result, the constraints at z ∼ 6 are primarily driven by [CII].
Comparing the constraints for fiducial values of ω fid BD = 800 and ω fid BD = 10 5 , we find that the latter are tighter.To understand this behaviour, it is important to stress that the larger the value of ω BD is, the lesser the scalar field is evolving.In other words, in the limit ω BD → ∞ the scalar field is completely frozen and reduces to a cosmological constant.Therefore, requiring a value of Geff different than unity today would imply almost the same value at earlier times.On the other hand, it is easier to accommodate deviations from Geff = 1 for lower values of ω BD , since an evolving scalar field would go back (in the past) to its tracker value (φ = 1) rapidly.This is reflected also in the constraints of the other parameters.As an intuitive example, one could roughly think that the constraints on the parameters of the ω fid BD = 10 5 model are the ones on the ω fid BD = 800 model, but fixing Geff = 1, rather than marginalizing over it.Therefore, in comparing our results with previous forecasts, e.g., those of Alonso et al. (2017), the difference in the assumed fiducial values, should be kept in mind.

Early Dark Energy driven by ultra-light axions
In Fig. 4, we show the 1-σ constraints on the three parameters of the EDE model as a function of spectrome- ter hours (in blue), and the corresponding sky coverage (in red).Again solid (dashed) lines show the results with (without) interloper noise, and the grey bands are the current constraints from Simon et al. (2022).As described in §2.2, in comparison of our results with existing bounds, caution should be taken since not having Planck constraints on EDE parameters and using only ΛCDM priors, makes our LIM+Planck constraints sig-nificantly more pessimistic.Keeping this in mind, we see that the constraints on f EDE can be improved significantly over the current bounds.On the other hand, lowering the uncertainties for log 10 (z c ) and in particular θ i is more challenging, requiring large number of spectrometer hours and larger sky coverage.In addition to the point above about the comparison with the current bounds, one can argue that the poor constraint on θ i is also driven by the fact that the signature left by a change in its value, i.e. the normalized initial value of the scalar field, is erased by the proceeding evolution, mainly due to the EDE feature at recombination.For f EDE and log 10 (z c ), removing the interloper noise would allow a factor of ∼ (4 − 5) improvement over the current bounds for ∼ 10 8 spectrometer hour we considered, while the constraints on θ i (without full Planck constraints) would be a factor of ∼ 2 worse.In Fig. 5, we show the 1-and 2-σ marginalized constraints on parameter pairs.There are several degeneracies between EDE and ΛCDM parameters in LIM-only results in addition to those among the EDE parameters.The most notable ones are between f EDE , and Ω b , A s , and n s , and between log 10 (z c ), A s , and n s .Adding Planck priors on ΛCDM improves the constraints, in particular by tightening the constraints on Ω b .In Fig. 6, we show the 1-σ constraints as a function of maximum redshift z max .Again, the trend is similar to that of the JBD model; constraints improve with increasing redshift up to z max 2.5.After that, the constraints almost completely plateaus.Without imposing Planck priors (not shown in this plot), the constraint on log 10 (z c ) improves by ∼ 14% at z max ∼ 6, at which [CII] signal is also included.Let us emphasize that in interpreting the information content of LIM measurements at different redshifts some care should be taken.First, we have only included the effects of EDE on liner matter power spectrum, and modeled the line power spectrum at linear level, and also did not include any modification to halo mass function due to and EDE (Klypin et al. 2021).Second, for the survey specification that we have assumed, the thermal noise rapidly increases as we go to higher redshifts beyond z ≃ 2.5.Therefore, the fact that our results do not show a significant improvement beyond z ≃ 2.5 should not be taken as a strong statement about the constraining power of LIM measurements at high redshifts.

Standard CPL parametrization of Dark Energy
In Fig. 7, we show the results of survey optimization for the CPL model, imposing ΛCDM Planck priors  and including (solid) or neglecting (dashed) interlopers.Similar trends to other models is observed, i.e. improvement of constraints as a function of spectrometer hours, and a near power-law increase of the required sky coverage.The considered survey can provide significantly tighter constraints in the w 0 − w a plane than the current bounds from the combination of CMB, BAO, and Supernovae (Aghanim et al. 2020).For ∼ 10 8 spectrometer hours, and including the interloper noise, we obtain σ(w 0 ) ∼ 0.016 and σ(w a ) = 0.034.As shown in Fig. 7, even a relatively modest survey covering ∼ 10% of the sky with ∼ 10 7 spectrometer hours can improve the current bounds on w 0 and w a by factors of about 2 and 3, respectively.Fig. 8 shows the 2D marginalized constraints on CPL and ΛCDM parameters with (red) and without (blue) Planck priors.There are strong degenracies between w 0 and w a , as well as between each of them and h, Ω b , and Ω c .Again here, imposing Planck priors on ΛCDM improves the constraints on CPL parameters, primarily by tightening ΛCDM constraints.In actual joint analysis of the two data sets varying the same model, it is expected breaking parameter degeneracies will also play essential role (see e.g., Moradinezhad Dizgah et al. 2022a).
In Fig. 9, we show the 1-σ constraints on w 0 and w a as a function of maximum redshift.We also show the constraint for the model with w a = 0 in dashed line to highlight the information gain with increasing z max for the CPL model.While for a non-dynamic DE, there is no information gain at z 1.8, for the CPL model, the constraints keep on improving (only at 5 − 9% level in the redshift range of 2.5 < z max < 6).The addition of [CII] signal at z ∼ 6 improves constraints on both w 0 and w a by about 20% and 25%, respectively.
Our results are comparable to the forecasts obtained for Euclid-like experiments.Blanchard et al. (2020) found that the 1-σ errors on w 0 vary from 0.02 (optimistic scenario) to 0.04 (pessimistic scenario).This range of values improves current constraints and it is at the level of ours for a LIM experiment with ∼ 10 8 spectrometer hours.The situation is a bit more favorable to LIM-experiments for w a , where their quoted constraints are 0.1 − 0.2.From the bottom panel of Fig. 7 we notice that the next generation of LIM surveys improves easily these numbers, having a redshift excursion much higher than LSS surveys.

Shift-symmetric Horndeski Models
In Fig. 10, we show the results of survey optimization for the shift-symmetric Horndeski model, imposing the Planck priors and including (solid) or neglecting (dashed) interlopers.Comparing with the constraints on w 0 in this model w.r.t. the corresponding one for the pure CPL parametrization (shown in the top plot of Fig. 7), we can see that the expected constraints with interlopers (solid blue lines) are essentially the same.This indicates that opening up the parameter space with αB and m does not degrades the constraining power of LIM experiments for the effective constant equation of state, as one would naively expect.This clearly suggests that the new parameter have orthogonal effects when considering LIM observables.It is also important to stress that what we quote here as "current constraint" for w 0 is lower than the corresponding constrain for the CPL model in is due to the fact that in obtaining these constraints different datasets are used.Similar observation can be made for w a , with the only remarkable difference at low spectrometer hours its here are slightly worse than in the previous analysis.This suggests that the constraining power of current LIM surveys leaves some that broaden the allowed space.Considering αB , LIM alone will be to current constraints for 4 × 10 8 spectrometer hours.With the parameter the situation is a bit different.First, the minimum spectrometer hours needed to improve the current constraints are 7 × 10 8 ; second, there is an almost-plateau at high spectrometer hours, which signals a diminishing return in pushing the survey specifications to their limits in constraining m.Fig. 11 shows the 1-and 2-σ marginalized constraints on the shift-symmetric Horndeski and |LambdaCDM parameter pairs, fixing the number of spectrometer hours to ∼ 10 8 .The LIM-only (blue) contours display strong degeneracies between w 0 and h, Ω b and Ω c as in Fig. 8. Also, w a has some mild to large degeneracy with the same parameters.The degeneracy between w 0 and (α B , m) is removed when adding Planck priors, which tightly constrain h and Ω b .The constraints on w 0 do not degrade with respect to CPL because the additional parameters are only mildly degenerate with w 0 when imposing Planck priors.Finally, it is worth noticing a strong degeneracy between αB and m, which is due to the fact that the time modulation parameterized by m can be reabsorbed into the overall amplitude αB for small redshift variations.
Fig. 12 shows the marginalized 1σ constraints as a function of redshift.The improvements of the con-straints for dynamical-DE parameters, w 0 and w a is significant up to z ≃ 3.5.Higher redshifts have little to none constraining power, similarly to the other models.It is worth noticing that here the constraints saturate at higher redshifts w.r.t.those of Fig. 9.This is because we chose a different fiducial model.Instead of a pure Λ behaviour, here we have an evolving DE component, which leaves clear imprints at higher redshifts.The saturation of the constraining power for αB and m is reached at lower, but still large, redshift (z ≃ 2.5).It is easy to show that the time evolution of α B (t) goes to zero more slowly than the most commonly considered models of modified gravity/DE which are designed to have low redshift signatures.

Effective Description of Luminal Horndeski Models
The effective description of DE is, to some extent (as discussed in §2) the most general class of models we are investigating.On top of that, our fiducial models were designed ad-hoc to have signatures beyond z ≃ 1 as dynamical-DE models.It is important to stress that the analysis of these models is meant to be illustrative of the potential of LIM surveys to probe general EFTofDE models, rather than providing constraints on a specific model.Given our chosen parametrization, shown in Eqs. ( 23), constraints on this class of models provide insights into evolution of dark energy equation of state, as well as time evolution of Planck mass.As described in §2, both considered models are designed to have signatures at high-redshift, i.e., 2 z 10.In particular, the effective Planck mass of model1 increases smoothly with time, deviating from unity at the level of ∼ 10% at z ≃ 10. model2 is even more extreme, with deviations ∼ 25% at z ≃ 10, but also a feature at z ≃ 1 after which Planck mass starts decreasing.Fig. 13 shows the results of survey optimization; the marginalized 1-σ constraints as a function of spectrometer hours, together with the required sky coverage15 .For all parameters, similar to all other models considered, the constraints improve as we increase the spectrometer hours.However, contrary to previous models, here removing the interlopers seems to be less important and we do not see saturation of the constraints in the range of spectrometer hours considered.This suggests that it should be possible (even if unrealistic with the next generation of surveys) to extract more information on these parameter increasing the observation time.One notable difference between model1 and model2 (not shown here) is the role of the interloper lines; perfect removal of the interlopers does not considerably improve the constraints in model2, while it does for model1.
In Fig. 14, we show the 1-and 2-σ marginalized constraints for model1 (top panel) and model2 (bottom (1) M .Fig. 15 shows the marginalized 1-σ constraints as a function of redshift for model1 (left) and model2 (right).While the perturbation parameters are better constrained in model2, the constraints on background parameters are tighter in model1.The redshift dependencies of the constraints on w 0 and w a are similar in both model, mostly saturating at z ≃ 1.5.For model2 (right), the constraints on c (2) M improve sharply (and by ∼2 orders of magnitude) as a function of redshift.To interpret this result, one should take into account that the fiducial value for c (2) M in model1 was zero, while in model1 we chose a value different from zero.This makes M and c (1) M (in the right plot) lie nearly completely on top of each other.Here, interloper noise is included, Planck ΛCDM priors are imposed, and the spectrometer hours is fixed to ∼ 10 8 .small variations around the fiducial model much more important in the second model.

CONCLUSION
Elucidating the nature of the dark energy component driving the current accelerated expansion and searching for possible modifications of gravity are two of major science goals for upcoming optical galaxy surveys.By offering the possibility of probing the LSS over a wide range of redshifts and scales, line intensity mapping can provide significant improvements in our understanding of the theory of gravity and the component of the Universe's energy density acting as DE, in the regime largely beyond the reach of stage IV galaxy surveys.
In this paper, we performed Fisher forecasts to determine the potential of a hypothetical future multi-stage ground-based LIM survey, targeting several rotational lines of CO (from J=2-1 to J=6-5) and the fine structure line of [CII], in constraining several beyond ΛCDM models.The considered models include two classes; the first class are consistently described by a specific covariant Lagrangian for a scalar field such as Jordan-Brans-Dicke gravity and (pre-recombination) early DE model.The second class are described in terms of effective description of the evolution of the background and perturbations, without defining a Lagrangian.An hybrid approach has been used for shift-symmetric Horndeski models.i.e., instead of constraining the parameters appearing in the Lagrangian, we used the corresponding effective description.For each model, we varied two survey specifications, the number of spectrometer hours and the fraction of sky covered by the survey, and optimized the survey strategy to obtain the best parameter constraints minimizing the survey cost.We considered optimistic and pessimistic scenarios for removal of line interlopers vs. accounting for them as a source of (anisotropic) noise, and obtained the parameter constraints with and without imposing Planck priors on ΛCDM parameters.
Our results show that for all models, the considered LIM survey can improve upon the current bound by a factor of a few to an order of magnitude.The parameter uncertainties decrease as one increase the number of spectrometer hours and required sky coverage grows nearly linearly with decreasing 1-σ errors.The interloper noise increases the required sky coverage to compensate for additional noise contribution.Including the interloper lines degrade the constraints on all models, but affects some parameters more significantly than others, causing a saturation of the parameter constraints even when increasing the number of spectrometer hours.Imposing Planck priors, thus tightening constraints on ΛCDM parameters, improves on beyond ΛCDM parameters, especially those most degenerate with them.
Among the models considered, for the JBD, the axiondriven EDE, and the CPL models, the constraints improve with increasing z max up to z max = 2.5, after which a near plateau is hit.Addition of [CII] signal at z ≃ 6 can improve the constraints on some of the parameters by at most (15-25)%.For shift-symmetric Horndeski model, the constraints on DE equation of state keeps on improving up to higher redshifts (z max ≃ 3.5) and plateaus afterwards.For our selected effective description of luminal Horndeski models, the impact of interlopers doesn't lead to saturation of the constraining power at high spectrometer hours, and instead just degrade the constraint by nearly constant factors.The interlopers affect the sky required sky coverage for best constraints on parameters describing the time evolution of Planck mass more significantly than the dark energy equation of state parameters.
It is worth noting that for several science goals considered in §6, a millimeter-wave LIM experiment begins to produce constraints that are competitive with current measurements at 10 7 spectrometer-hours, concurrent with previous work focused on LIM-enabled constraints on neutrino properties (Moradinezhad Dizgah et al. 2022a).With a dedicated multi-year survey (with total observing time of 10 4 hours), such a survey could be completed with a broadband camera of several dozen pixels -only a few times larger than LIM-focused instruments currently being deployed (e.g.Karkare et al. 2022b).As demonstrated, such surveys would optimally cover a survey area of a few thousand square degrees, well within the reach of current survey telescopes (e.g., SPT).
There are several directions in which this work can be extended.First, here we have assumed that the line bias and mean brightness temperature are perfectly predicted by the theory.In reality, given the uncertainty on the theoretical predictions of these parameters (e.g.due to uncertainty in the model of line luminosity), one should marginalize over them.Second, since the variation of bias and mean brightness temperature is expected to degrade the constraints due to their degeneracies with cosmological parameters, studying how much the cross-correlations between different emission line signals, as well as between LIM and galaxy clustering and the CMB lensing can ameliorate these degeneracies is of great interest.On a more general ground, quantifying the potential improvements of the inferred constraints resulting from including both auto and cross-correlations is highly relevant.Third, in order to to go beyond the Fisher forecasts, and in the absence of current wide-field LIM survey, performing likelihood analysis of mock LIM data, imposing Planck priors on the specific beyond ΛCDM models considered is necessary.While here we have used Planck priors on ΛCDM parameters, joint analysis of LIM and CMB primary anisotropies assuming the same model for the two is essential.Fourth, we have assumed a linear model of the line power spectrum, and limited our analyses to relatively large scales.Including nonlinearities and promoting the model to include one-loop corrections, and marginalizing over a complete set of line biases is a natural next step.The improved model will allow us to go to smaller scales, and potentially help in breaking the degeneracies between mean brightness temperature and linear line bias.Lastly, accurate characterization of foregrounds and observational systematics are crucial in obtaining more realistic forecasts.

Figure 1 .BDFigure 2 .Figure 3 .
Figure1.Blue lines show the 1-σ constraints on JBD model parameters as a function of spectrometer hours, including (solid) and neglecting (dashed) interlopers.The red lines are the corresponding required sky coverage, and the grey bands are the existing constraints from a combination of KiDSx2dFLenS, BOSS full shape, Planck18 and Pantheon data(Joudaki et al. 2022).

Figure 4 .Figure 5 .Figure 6 .
Figure4.Blue solid (dashed) lines show the marginalized 1-σ marginalized constraints on EDE driven by ultra-light axions, including (neglecting) line interloper noise.The red lines show the corresponding required sky coverage, and the grey band shows the current bounds from Planck (temperature, polarization, and lensing) data and BOSS full-shape galaxy power spectrum, imposing prior on h0 from SH0ES local distance-ladder data(Simon et al. 2022).The ΛCDM Planck priors are applied.

Figure 7 .
Figure 7. Blue solid (dashed) lines show the marginalized 1-σ constrained and optimized survey configuration for CPL model, including (neglecting) interlopers, and imposing Planck priors on ΛCDM parameters.Red lines are the corresponding required sky coverage, and the grey bands are the current bounds from CMB, BAO, and Supernovae data.

Figure 8 .Figure 9 .
Figure8.Marginalized 1-and 2-σ constraints on CPL and ΛCDM parameter pairs from LIM alone (blue) and when imposing Planck priors on ΛCDM parameters.Interloper noise is included here.The sky fraction is assumed to be f sky = 0.5, and the spectrometer hours is fixed to ∼ 10 8 .

Figure 10 .
Figure 10.Blue solid (dashed) lines show marginalized 1-σ constraints on background (top) and perturbation (bottom) parameters of the shift-symmetric Horndeski model, including (neglecting) interloper noise.Planck priors on ΛCDM parameters are imposed.Red lines are the corresponding required sky coverage and grey bands are the current bounds from CMB, BAO, RSD and Supernovae Type IA data jointly with theoretical priors of Traykova et al. (2021).

Figure 11 Figure 12 .
Figure11.1-and 2-σ marginalized constraints on ShiftSym-Horndeski vs. ΛCDM parameters from LIM alone (blue) and when imposing Planck priors on ΛCDM parameters.Interloper noise is included here.The sky fraction is assumed to be f sky = 0.5, and the spectrometer hours is fixed to ∼ 10 8 .

Figure 13 .
Figure 13.Blue solid (dashed) lines show marginalized 1-σ constraints on background (top row) and perturbation (bottom row) parameters of the Effective-DE model1, including (neglecting) interloper noise.Red lines are the corresponding required sky coverage.Planck ΛCDM priors are imposed.

Figure 14 .Figure 15 .
Figure14.Marginalized 1-and 2-σ constraints on Effective-DE models and ΛCDM parameter pairs with (red) and without (blue) Planck ΛCDM priors.The top plot shows the constraints for model1 while the bottom one correspond to model2.Interloper noise is included here, and the spectrometer hours is fixed to ∼ 10 8 .

Table 1 .
Summary of the models considered in the Fisher forecasts, and their corresponding free parameters.
matrices are independent so can be summed to obtain the joint constraints.A more consistent forecast including Planck and LIM data requires computing the Fisher matrices of the two observables for the same beyond ΛCDM model.The choices of minimum and maximum scales, [k min , k max ], and the modeling of interlopers are the same as MoradinezhadDizgah et al. (2022a).The list of the models that we consider in our analysis is given in Table