Cosmology Inference from Biased Tracers using the EFT-based Likelihood

TL;DR

This paper applies an EFT-based likelihood approach to galaxy clustering data to infer cosmological parameters, demonstrating accurate recovery of the power spectrum normalization $\sigma_8$ with potential for improved robustness against small-scale complexities.

Contribution

It introduces a novel EFT-based likelihood method that utilizes high-order information for cosmological inference from biased tracers, improving robustness and accuracy.

Findings

Successfully recovered $\sigma_8$ within 95% CL across various samples and redshifts.

Incorporated bias terms up to second order and leading higher-derivative terms.

Identified potential sources of systematic bias for future refinement.

Abstract

The effective-field-theory (EFT) approach to the clustering of galaxies and other biased tracers allows for an isolation of the cosmological information that is protected by symmetries, in particular the equivalence principle, and thus is robust to the complicated dynamics of dark matter, gas, and stars on small scales. All existing implementations proceed by making predictions for the lowest-order $n$-point functions of biased tracers, as well as their covariance, and comparing with measurements. Recently, we presented an EFT-based expression for the conditional probability of the density field of a biased tracer given the matter density field, which in principle combines information from arbitrarily high order $n$-point functions. Here, we report results based on this likelihood by applying it to halo catalogs in real space, specifically on the inference of the power spectrum normalization $\sigma_8$. We include bias terms up to second order as well as the leading higher-derivative term. For a cutoff value of $\Lambda = 0.1 h\,{\rm Mpc}^{-1}$, we recover the ground-truth value of $\sigma_8$ to within 95% CL for different halo samples and redshifts. We discuss possible sources for the remaining systematic bias in $\sigma_8$ as well as future developments.

Figures

Tables

Table 1. Table 1 : The halo samples used in our numerical tests. Throughout, masses M ≡ M 200 ​ m 𝑀 subscript 𝑀 200 𝑚 M\equiv M_{200m} are spherical-overdensity masses with respect to 200 times the background matter density.

Redshift	Mass range ${\log }_{10}M\left[h^{-1}{M}_{\odot }\right]$	$N_{\text{halo}}$ (run 1)	$N_{\text{halo}}$ (run 2)	$\overline{n}\left[{\left(h^{-1}\text{Mpc}\right)}^{-3}\right]$
0	$[13,13.5]$	2807757	2803575	$3.5\times {10}^{-4}$
0	$[13.5,14]$	919856	918460	$1.1\times {10}^{-4}$
1	$[13,13.5]$	1507600	1506411	$1.9\times {10}^{-4}$
1	$[13.5,14]$	301409	302182	$3.8\times {10}^{-5}$

Table 2. Table 2 : Summary of results for Λ = 0.1 ​ h ​ Mpc − 1 Λ 0.1 ℎ superscript Mpc 1 \Lambda=0.1\,h\,\text{Mpc}^{-1} and k max = Λ / 2 subscript 𝑘 max Λ 2 k_{\rm max}=\Lambda/2 , with the likelihood settings described in the text. For the best-fit scaled σ 8 subscript 𝜎 8 \sigma_{\mathrm{8}} estimate α ^ ^ 𝛼 \hat{\alpha} , results from run 1 and run 2 are shown individually with estimated 68% CL error bars. b 1 subscript 𝑏 1 b_{1} and stochastic amplitudes are reported for the fiducial σ 8 = σ 8 fid subscript 𝜎 8 superscript subscript 𝜎 8 fid \sigma_{\mathrm{8}}=\sigma_{\mathrm{8}}^{\text{fid}} and averaged over both runs. The stochastic variance σ ε 2 superscript subscript 𝜎 𝜀 2 \sigma_{\varepsilon}^{2} is scaled to the Poisson expectation for the given halo sample, as described in Appendix  B . The last column shows the ratio of the higher-derivative stochastic amplitude to the lowest-order one, indicating the scale associated with the expansion of σ 2 ​ ( k ) superscript 𝜎 2 𝑘 \sigma^{2}(k) in k 𝑘 k .

Redshift	Mass range ${\log }_{10}M\left[h^{-1}{M}_{\odot }\right]$	$\hat{\alpha }$ (run 1)	$\hat{\alpha }$ (run 2)	$b_1$	${\sigma }_{\epsilon }^2$ $[$Poisson$]$	${\sigma }_{\epsilon ,2}/{\sigma }_{\epsilon }$ $\left[{\left(h^{-1}\text{Mpc}\right)}^2\right]$
0	[13.0-13.5]	$0.96\pm 0.02$	$0.90\pm 0.02$	$1.20$	$1.11$	$-25.5$
0	[13.5-14.0]	$1.02\pm 0.03$	$1.01\pm 0.03$	$1.61$	$0.96$	$-11.7$
1	[13.0-13.5]	$1.05\pm 0.04$	$1.16\pm 0.04$	$2.36$	$0.93$	$1.3$
1	[13.5-14.0]	$1.05\pm 0.05$	$0.93\pm 0.05$	$3.49$	$0.89$	$10.5$

Equations

$${\vec{\delta }}_{h,\det }[\vec{\delta },\{b_O\}]=\sum _O{b}_O\vec{O}[\vec{\delta }],$$

$$O\in \{\delta ,{\nabla }^2\delta ,{\delta }^2,K^2\},$$

$$K^2\equiv (K_{ij}{)}^2={\left(\left[\frac{{\partial }_i{\partial }_j}{{\nabla }^2}-\frac{1}{3}{\delta }_{ij}\right]\delta \right)}^2$$

$${\mathcal{B}}_O=\{b_1,c_{{\nabla }^2\delta },b_2/2,b_{K^2}\}.$$

$$\ln P\left(\vec{{\delta}}_{h}\Big{|}\vec{{\delta}},\{b_{O}\}\right)=-\sum_{\bm{k}\neq 0}^{k_{\rm max}}\left[\frac{1}{2}\ln\sigma^{2}(k)+\frac{1}{2\sigma^{2}(k)}\left|{\delta}_{h}(\bm{k})-{\delta}_{h,\rm det}[\vec{{\delta}},\{b_{O}\}](\bm{k})\right|^{2}\right]\,.$$

$${\sigma }^2(k)={\left({\sigma }_{\epsilon }+k^2[{\sigma }_{\epsilon ,2}+b_1{\sigma }_{\epsilon {\epsilon }_m,2}]\right)}^2.$$

$${\delta }_{h,\det }(\mathbf{k})=\mu (\mathbf{k})+\sum _{O\in {\mathcal{O}}_{\text{marg}}}{b}_{O}O(\mathbf{k}),\mu (\mathbf{k})=\sum _{O\in {\mathcal{O}}_{\text{all}}\setminus {\mathcal{O}}_{\text{marg}}}{b}_{O}O(\mathbf{k}),$$

$$\displaystyle P\left(\vec{{\delta}}_{h}\Big{|}\vec{{\delta}},\{b_{O}\}\right)=\>$$

$$\displaystyle\times\exp\bigg{\{}-\sum_{\bm{k}\neq 0}^{k_{\rm max}}\bigg{[}\frac{\left|{\delta}_{h}(\bm{k})-\mu(\bm{k})\right|^{2}}{2\sigma^{2}(k)}-2\sum_{O\in\mathcal{O}_{\text{marg}}}b_{O}\frac{\Re\>[{\delta}_{h}(\bm{k})-\mu(\bm{k})]O^{*}(\bm{k})}{2\sigma^{2}(k)}$$

$$\displaystyle\hskip 71.13188pt+\sum_{O,O^{\prime}\in\mathcal{O}_{\text{marg}}}b_{O}b_{O^{\prime}}\frac{O(\bm{k})O^{\prime*}(\bm{k})}{2\sigma^{2}(k)}\bigg{]}\bigg{\}}\,,$$

$$P\left(\vec{{\delta}}_{h}\Big{|}\vec{{\delta}},\{b_{O}\}\right)=\frac{1}{{\cal N}}\exp\left[-\frac{1}{2}\sum_{\bm{k}\neq 0}^{k_{\rm max}}\ln\sigma^{2}(k)\right]\\ \times\exp\left\{-\frac{1}{2}C+\sum_{O\in\mathcal{O}_{\text{marg}}}b_{O}B_{O}-\frac{1}{2}\sum_{O,O^{\prime}\in\mathcal{O}_{\text{marg}}}b_{O}b_{O^{\prime}}A_{OO^{\prime}}\right\}\,,$$

$$C$$

$$B_O$$

$$A_{O{O}^{\prime }}$$

$$\int d^n\vec{x}\exp \left[-\frac{1}{2}{\vec{x}}^{\top }A\vec{x}+{\vec{B}}^{\top }\cdot \vec{x}\right]=\frac{(2\pi {)}^{n/2}}{|A{|}^{1/2}}\exp \left[\frac{1}{2}{\vec{B}}^{\top }{A}^{-1}\vec{B}\right].$$

$$\displaystyle P\left(\vec{{\delta}}_{h}\Big{|}\vec{{\delta}},\{b_{O}\}_{\text{unmarg}}\right)$$

$$\displaystyle=\frac{(2\pi)^{n_{\text{marg}}/2}}{{\cal N}}\Big{|}A_{OO^{\prime}}\Big{|}^{-1/2}\exp\left[-\frac{1}{2}\sum_{\bm{k}\neq 0}^{k_{\rm max}}\ln\sigma^{2}(k)\right]$$

$$\times \exp \left\{-\frac{1}{2}C(\{b_O\})+\frac{1}{2}\sum _{O,O^{\prime }\in {\mathcal{O}}_{\text{marg}}}{B}_O(\{b_O\})(A^{-1}{)}_{O{O}^{\prime }}{B}_{O^{\prime }}(\{b_O\})\right\}.$$

$$-2\ln P({\delta }_h|\delta )$$

$$=-2\ln P_{\text{fid}}({\delta }_h|\delta )-2\sum _O{B}_O^{\text{err}}{b}_O+\sum _{O,O^{\prime }}{A}_{O{O}^{\prime }}^{\text{err}}{b}_O{b}_{O^{\prime }},$$

$$B_O^{\text{err}}$$

$$A_{O{O}^{\prime }}^{\text{err}}$$

$$-2\ln P_{\text{fid}}({\delta }_h|\delta )$$

$$A_{O{O}^{\prime }}$$

$$B_O$$

$$\Delta \mathbf{b}\equiv \mathbf{b}-{\mathbf{b}}^{\text{fid}}=\left({\mathbf{A}}^{-1}{\mathbf{A}}^{\text{err}}\right)\left({\mathbf{b}}^{\text{fid}}\right)-{\mathbf{A}}^{-1}\left({\mathbf{B}}^{\text{err}}\right).$$

$$\Delta b\Big{|}_{\text{loops}}\sim\frac{\sum_{\bm{k}\neq 0}^{k_{\rm max}}\sigma^{-2}(k)P_{\text{2-loop}}(k)}{\sum_{\bm{k}\neq 0}^{k_{\rm max}}\sigma^{-2}(k)P_{\text{L}}(k)}\,.$$

$$P_{\text{2-loop}}(k)\sim {⟨[{\delta }^3](\mathbf{k})[{\delta }^3]({\mathbf{k}}^{\prime })⟩}^{\prime }=6\int d^{3}r[{\xi }_{\mathrm{L}}(r){]}^3{e}^{i\mathbf{k}\cdot \mathbf{r}},$$

$$\Delta b\Big{|}_{\text{stoch.}}\sim\frac{|{\varepsilon}(\bm{k})|}{|{\delta}(\bm{k})|}\Big{|}_{k_{\rm max}}\sim\sqrt{\frac{P_{\varepsilon}}{P_{\text{L}}(k_{\rm max})}}\stackrel{{\scriptstyle\mbox{\small Poisson}}}{{=}}\Big{[}\bar{n}P_{\text{L}}(k_{\rm max})\Big{]}^{-1/2}\,,$$

$$P^{\mathrm{prof}}(a)=\underset{\{b_i\}}{\mathrm{argmax}}[P(a,\{b_i\})],$$

$${\sigma }_8$$

$$\hat{\alpha }\equiv \frac{{\hat{\sigma }}_8}{{\sigma }_8^{\text{fid}}}.$$

$$[\delta ](\mathbf{k})$$

$$\displaystyle\Bigl{[}{\delta}^{2}\Bigr{]}(\bm{k})$$

$$\displaystyle\Bigl{[}K^{2}\Bigr{]}(\bm{k})$$

$${\Sigma }_{1-3}^2(k)=4{\int }_{\mathbf{p}}{W}_{\Lambda }(\mathbf{p})W_{\Lambda }(\mathbf{k}-\mathbf{p})F_2(-\mathbf{k},\mathbf{p})P_{\text{L}}(p).$$

Full text

1 Introduction

State-of-the-art approaches for the analysis of large-scale structure (LSS) data typically make use of summary statistics like the two-point correlation function to compare theoretical models to observational data [e.g., 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. However, as has been realized early on [15], the galaxy distribution is non-Gaussian, and even an infinite hierarchy of higher order correlation functions may be insufficient to fully capture all cosmological information encoded in LSS data sets [16, 17].

Alternative approaches have since been developed that take a different, more ambitious avenue to cosmological signal inference. Instead of focusing on summary statistics, they aim directly at reconstructing the three-dimensional underlying matter density field from observations of astrophysical tracers like galaxies [18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]. In an important step, Ref. [28] improved on previous implementations by introducing a physical forward model for the matter field as a means to more precisely predict the intricate statistical properties of the evolved density field. Starting from a set of simple Gaussian initial conditions at high redshift as probed by cosmic microwave background radiation experiments [e.g., 31, 32, 33, 34, 35, 36, 37], nonlinear effects of gravitational collapse are taken into account via approximate analytical or purely numerical methods to compute the corresponding evolved density field at low redshift that can then be compared to observations. Owing to the Gaussian nature of the initial density field, the two-point correlation function becomes a lossless form of data compression in this regime and fully captures all information relevant to cosmology. As a by-product of extending the reconstruction to the initial density field, a fully probabilistic description of all possible histories of the structure formation process compatible with the analyzed data set becomes available, offering interesting opportunities to implement other cosmological tests [e.g., 38].

Unsurprisingly, these sophisticated approaches are faced with numerous challenges that are both theoretical and numerical in nature. For example, exploring an extremely high-dimensional parameter space typical of the reconstruction of a three-dimensional field necessitates the use of sophisticated sampling algorithms [39], [e.g., 40, 24]. Further, since a comprehensive theoretical understanding of the formation of actual observed tracers such as galaxies is still lacking, an effective model for galaxy biasing is necessary to connect any underlying dark matter density field to actual observables like galaxy positions and redshifts [41, 42, 43, 44, 45, 46, 47]. Finally, a statistical comparison of observational data to theoretical predictions requires the use of a likelihood function that quantifies the probability to observe the given data realization under the model hypothesis.

Numerous functional forms for this likelihood have been proposed to compare model predictions to observations [e.g., 48, 49, 19, 50, 51]. Based on the effective field theory (EFT) approach [52, 53], however, only recently a rigorous description for the LSS likelihood has been developed [54, 55]. Most notably, and in contrast to previous attempts at modeling the statistics of LSS data, this likelihood is formulated in Fourier space. In this paper, we use numerical simulations to provide an in-depth analysis of the EFT LSS likelihood, using dark matter halos in different mass bins as tracers. In particular, we test whether it allows for unbiased cosmological parameter inference from reconstructed matter density fields.

Since they are based on simulations, our tests use the correct phases of the initial density field. Hence, we eliminate an important source of uncertainty that would normally result in noticeably larger error bars on cosmological parameters in real surveys. Our tests are thus more stringent to pass compared to a realistic application, where the initial density field is unknown and has to be simultaneously inferred from the data.

The paper is organized as follows. In Sec. 2 we briefly review the statistical framework to analyze LSS data developed in [54], which forms the basis of our analysis. We also present an analytical marginalization scheme which greatly reduces the number of free parameters that need to be varied. We then introduce a set of simulations used to assess the performance of the likelihood in Sec. 3. After presenting our implementation in Sec. 4, we turn to the results in Sec. 5. We summarize our findings in Sec. 6. The appendices present further details and a derivation of the Poisson expectation for stochasticity.

2 Method

We start by giving an overview of the main aspects of the methodology used in our analysis, and then introduce a useful extension, analytical marginalization, that aims at obtaining numerical results more efficiently.

2.1 Recap: Fourier-space likelihood

In this section, we provide a short summary of the results presented in [54]. Our main objective is to obtain a probabilistic reconstruction of the initial density field $\vec{{\delta}}_{\text{in}}$, cosmological parameters, and additional nuisance parameters necessary to capture the uncertainties of the process of structure formation. To do so, we need to model the effects of various aspects of the underlying physical processes. More precisely, we have to specify a prior characterizing the statistical properties of the initial density perturbations, and provide a forward model, describing the gravitational collapse of matter overdensities over time. Additionally, we need to include a bias model, connecting the observed tracers to the matter density field, and, finally, a likelihood, describing the stochastic aspects involved in this process and therefore allowing for a quantitative comparison of model predictions to observations.

We begin by reviewing the basic notation introduced in [54] that we will make use of in the remainder of this paper. For simplicity and in view of the tests presented in Sec. 2.4, we will refer to the tracers of the matter density field as halos. Then, given the evolved matter density field $\vec{{\delta}}$, at any order in perturbation theory, the deterministic halo density, $\vec{{\delta}}_{h,\mathrm{det}}$, can be expressed as linear sum over a finite set of operators, $\vec{O}$,

[TABLE]

where the $b_{O}$ are the associated bias parameters (see [47] for a review). Here and throughout, the vector notation denotes fields, in our case discretized on a uniform cubic grid. In the following, we will use the discrete Fourier transform, so that fields in Fourier space are dimensionless as well. Following Ref. [54], we use the bias expansion up to second order in perturbations, i.e., we restrict ourselves to the following set of operators:

[TABLE]

where ${\delta}$ is the fractional matter density perturbation, and

[TABLE]

is the tidal field squared. The corresponding bias parameters are denoted as

[TABLE]

The coefficient of $\nabla^{2}\delta$ is denoted as $c_{\nabla^{2}{\delta}}$ to emphasize that it also absorbs other higher-order contributions, as discussed in [54]. All operators in Eq. (2.2) are constructed from the density field filtered with a sharp low-pass filter in Fourier space which only retains Fourier modes at $|\bm{k}|<\Lambda$, where $\Lambda$ functions as a cutoff for the theoretical description. The sharp cut in Fourier space naturally arises in effective field theory approaches which integrate out initial density perturbations above a cutoff $\Lambda$ [55], although this does not mean that other filter shapes are excluded in principle. However, we derived in [54] that a sharp-$k$ filter is a condition for obtaining an unbiased maximum-likelihood point of the Fourier-space likelihood in the form used here.

Moreover, the operators in Eq. (2.1) are renormalized, which in this case essentially means that the quadratic operators ${\delta}^{2}$ and $K^{2}$ are made to be orthogonal to ${\delta}$. We refer to Appendix A for details on the exact definition and renormalization procedure. It is worth noting that the renormalization is not essential to obtain an unbiased $\sigma_{\mathrm{8}}$ estimate, but allows for a comparison of the resulting bias parameters with measurements from $n$-point functions.

Notice that Eq. (2.1) involves an expansion in two small parameters, essentially orders of perturbations and spatial derivatives (see Sec. 4.1 of [47] for a detailed discussion). Here, we assume that both small parameters are comparable, which leads us to include terms up to second order in perturbations as well as the leading higher-derivative operator $\nabla^{2}{\delta}$ in Eq. (2.2). The reasoning behind this is discussed in greater detail in [54].

As derived in [54], for an observed halo field, $\vec{{\delta}}_{h}$, we can then compute the conditional probability for the halo density field given the matter density in Fourier space,

[TABLE]

In the following, we will refer to this conditional probability simply as “likelihood,” since it is the part of the overall likelihood of biased tracers that is relevant for the study presented in this paper. Here, $\sigma^{2}(k)$ is a dimensionless variance that is analytic in $k$, specifically a power series in $k^{2}$. For our numerical implementation, we include terms up to order $k^{4}$, and parametrize $\sigma^{2}(k)$ as follows:

[TABLE]

The parametrization is chosen so that $\sigma^{2}(k)$ is positive definite. $\sigma_{\varepsilon}^{2}$ can be interpreted as the amplitude of halo stochasticity in the large-scale limit ($k\to 0$). In Appendix B we derive the expectation for $\sigma_{\varepsilon}^{2}$ for a Poisson process, which will be useful for the interpretation of numerical results on this parameter. $\sigma_{{\varepsilon},2}^{2}$ is the leading scale-dependent correction to the halo stochasticity (see Sec. 2.7 of [47] for a discussion). Finally, $\sigma_{{\varepsilon}{\varepsilon}_{m},2}$ captures the cross-correlation of stochasticity in the halo and matter fields; the stochasticity in the matter field ${\varepsilon}_{m}(\bm{k})$ is constrained by mass- and momentum conservation to be of order $k^{2}$ on large scales.

The closed-form expression for the likelihood in Eq. (2.5) affords a straightforward way to derive maximum likelihood estimates for bias parameters, and, as shown in [54], for the cosmological parameter $\sigma_{8}$. This framework was used in [54] to demonstrate unbiased parameter estimation from LSS data without the explicit use of conventional summary statistics. Here, we will obtain results using the field-level likelihood, rather than the analytical maximum-likelihood point discussed in [54].

The degeneracy between $b_{1}$ and $\sigma_{\mathrm{8}}$, which is perfect in linear theory, is broken when including nonlinear information. In particular, the fact that the displacement term contained in the second-order matter density is also multiplied by $b_{1}$, coupled with the fact that the second-order matter density scales differently with $\sigma_{\mathrm{8}}$ than the linear-order one. Thus, fundamentally, the possibility of estimating $\sigma_{\mathrm{8}}$ in this way is due to the equivalence principle, which ensures that galaxies move on the same trajectories as matter on large scales, and thus requires that the second-order displacement term is multiplied by the same bias coefficient as the linear-order density field (see also Sec. 2 of [47]). At higher orders in perturbations, more such terms that are protected by the equivalence principle appear, and the EFT likelihood will consistently capture those as well once extended to higher order.

2.2 Marginalizing over bias parameters

In the approach summarized in the previous section, a number of nuisance parameters have been introduced to capture the uncertainties associated with some of the poorly understood physical processes describing halo and galaxy formation. In practical applications, these parameters would then be estimated alongside cosmological parameters in a statistical analysis of observational data. Interestingly, owing to the simple functional form of the likelihood in Eq. (2.5), it is easily possible to analytically marginalize over some of the bias parameters $b_{O}$. Here, we consider only those that do not appear in the variance $\sigma^{2}(k)$; given Eq. (2.6), this includes all bias parameters except $b_{1}$. While it might be possible to extend the analytical marginalization to parameters which appear in the variance, the marginalization performed here is sufficient for our purposes.

Let us thus write

[TABLE]

where $\mathcal{O}_{\text{marg}}$ denotes the subset of operators, whose bias parameters we wish to marginalize over (we denote the cardinality of this set as $n_{\text{marg}}$). We can then write the likelihood Eq. (2.5) as

[TABLE]

where ${\cal N}$ is a normalization constant which is independent of all parameters. This expression can be more compactly written as

[TABLE]

where

[TABLE]

Note that $A_{OO^{\prime}}$ is a Hermitian and positive-definite matrix. The former is immediately obvious from its definition. The latter follows from the fact that $A_{OO^{\prime}}$ is the zero-lag covariance matrix of a set of sharp-$k$-filtered real fields $O(\bm{x})$. Eq. (2.9) then allows us to perform the Gaussian integral over the $b_{O}$. Here, we will assume uninformative priors, although Gaussian priors can trivially be introduced by adding a prior covariance to $A_{OO^{\prime}}$. The result is111Using the well-known Gaussian integral identity

$\int d^{n}\vec{x}\exp\left[-\frac{1}{2}\vec{x}^{\top}A\vec{x}+\vec{B}^{\top}\cdot\vec{x}\right]=\frac{(2\pi)^{n/2}}{|A|^{1/2}}\exp\left[\frac{1}{2}\vec{B}^{\top}A^{-1}\vec{B}\right]\,.$

[TABLE]

We have thus reduced the parameter space from $\{b_{O}\}$ to $\{b_{O}\}_{\text{unmarg}}$. This marginalization applies to an arbitrary number of bias coefficients to be marginalized over. The price to pay is that we now need to compute the vector $B_{O}$ and the matrix $A_{OO^{\prime}}$. Further, we need the determinant of $A_{OO^{\prime}}$ and its inverse. Note, however, that $n_{\text{marg}}$ – and therefore the size of $B_{O}$ and $A_{OO^{\prime}}$ – will rarely become a large number in practical applications. More specifically, in the tests presented below, we marginalize over $c_{\nabla^{2}{\delta}},b_{2},b_{K^{2}}$ and so $n_{\text{marg}}=3$. Further, $A_{OO^{\prime}}$ is independent of the remaining, unmarginalized bias parameters. On the other hand, $A_{OO^{\prime}}$ does depend on the parameters entering the variance $\sigma^{2}(k)$, Eq. (2.6).

At this point, let us briefly comment on the relation of this analytic marginalization to other approaches presented in the recent literature. In particular, Refs. [56, 57] perform a Gram-Schmidt orthogonalization on the fields entering the bias expansion Eq. (2.1), arguing that lower-order operators thus become independent of higher-order fields and make their corresponding bias parameters more robust to higher-order corrections. This approach is directly related to the analytic marginalization pointed out here. To see this, consider the case where an orthogonalization has been performed on the operators. In particular, this implies that $\left\langle\mu(\bm{k})O(\bm{k})\right\rangle=0$ for all $O\in\mathcal{O}_{\text{marg}}$. This in turn renders $B_{O}$ independent of all unmarginalized bias parameters, so that it becomes a constant vector. Then, the marginalized likelihood reduces to the same form as the unmarginalized likelihood keeping only the terms involving unmarginalized bias parameters. In this sense, orthogonalization is equivalent to marginalizing over bias parameters. The computational cost of both approaches is expected to be essentially the same. However, the marginalization described here does offer the possibility of including prior information on the bias parameters that are marginalized over.

2.3 Estimating systematic errors

One of the crucial advantages of the rigorous perturbative approach pursued here is that it allows for an estimate of the systematic error due to imperfections in the likelihood. We can distinguish three principal sources of such systematic errors (see [55] for more details on Types 2 and 3 in particular):

Type 1:

Errors in the forward model for the matter density field (and correspondingly the operators constructed from it);

Type 2:

Higher-order bias terms neglected in the expansion in ${\delta}_{h,\rm det}$;

Type 3:

Higher-order contributions to the variance $\sigma^{2}(k)$ as well as in the form of the likelihood itself.

The most rigorous way to evaluate the size of these contributions is to include the set of leading higher-order terms that have been neglected in the forward model, bias expansion, and likelihood, and evaluate the shift in resulting parameter values. In case of the forward model (Type 1), this can be tested by using the density field from N-body simulations instead of 2LPT to construct the bias operators. This will be presented in Sec. 5. In case of the bias expansion (Type 2), this test is not too difficult either, since the coefficients can be marginalized over analytically, as shown in the previous section. We defer an implementation of the higher-order contributions to future work however.

Let us here approximately estimate the size and scaling with $k_{\rm max}$ of the systematic error of Type 2. Note that strictly speaking we have two cutoffs: the cutoff $\Lambda$ of the sharp-$k$ filter, and $k_{\rm max}<\Lambda$. In practice, one will choose $k_{\rm max}$ as a fixed fraction of $\Lambda$ (see Sec. 5); hence, it is sufficient to consider the dependence on $\Lambda$ here. For simplicity, we evaluate the systematic shift in the bias parameters $b_{O}$. As described in [54], one can similarly evaluate the shift in $\sigma_{8}$ by introducing scaled bias parameters $\beta_{O}$. We will also count higher-derivative terms as higher order in perturbations, which assumes that the scale controlling higher-derivative terms does not differ greatly from the nonlinear scale. Thus, “higher-order contributions” include higher-derivative contributions in what follows.

Incorporating both sources of error described above, the correct likelihood can be written as

[TABLE]

where

[TABLE]

Here, $O_{\text{err}}(\bm{k})$ denotes the error field in the operator $O$ due to deficiencies in the forward model, while ${\delta}_{h}^{\text{h.o.}}(\bm{k})$ denotes the higher-order bias contributions to the actual halo density field. Finally, $P_{\text{fid}}({\delta}_{h}|{\delta})$ stands for the fiducial likelihood, which differs from the correct one due to the systematic error terms. In the second line of Eq. (2.12), we have dropped an irrelevant constant term which does not depend on the parameters being varied. Let us write the fiducial likelihood in analogy to Eq. (2.9) as

[TABLE]

Under the assumption that the parameter shift $\Delta b_{O}$ due to the systematic errors is small, one can immediately solve for this shift based on the maximum-likelihood points of the correct and fiducial likelihoods. One can then estimate the expected amplitude of the shift by taking the expectation values of $\bm{A},\bm{B},\bm{A}^{\text{err}},\bm{B}^{\text{err}}$. This is closely analogous to the “Fisher bias.” Using bold-face to denote vectors in the $n_{O}$-dimensional vector space of operators considered, we obtain in matrix notation

[TABLE]

This expression involves the expectation value of the correlators in Eq. (2.13) and Eq. (2.14) which are straightforward to evaluate in perturbation theory. We begin by estimating at which order in perturbation theory the various correlators contribute.

First, the expectation values of $A_{OO^{\prime}}$ and $B_{O}$ are of order $P_{\text{L}}(k)+P_{\text{1-loop}}(k)$ (in case of $A_{{\delta}{\delta}}$ and $B_{\delta}$), or of order $P_{\text{1-loop}}(k)$ (all other elements). On the other hand, both $\bm{B}^{\text{err}}$ and $\bm{A}^{\text{err}}$ are of order $P_{\text{2-loop}}(k)$. To see this, notice that both $O^{\text{err}}$ and ${\delta}_{h}^{\text{h.o.}}$ are at least of cubic order in the linear density field. This means that all correlators which involve these error fields are two-loop contributions, apart from the cross-correlation with ${\delta}^{(1)}$, which is at 1-loop order. The latter however only appears in $B^{\text{err}}_{\delta}$, via $\left\langle{\delta}^{\text{err,(3)}}{\delta}_{h}^{(1)}\right\rangle$ and $\left\langle{\delta}_{h}^{(3)}{\delta}^{(1)}\right\rangle$. As we argue in App. C of [54], these particular 1-loop contributions are of very similar shape as that coming from the higher-derivative bias, and are thus largely absorbed by $c_{\nabla^{2}{\delta}}$.

Thus, without performing any detailed calculation, we can very roughly estimate that

[TABLE]

As an approximate estimate of the size of two-loop correlators, we will use the auto-correlation of $[{\delta}^{3}]$:

[TABLE]

where $\xi_{\rm L}(r)$ is the linear matter correlation function. We emphasize that this is a very rough estimate: in reality, $\Delta\bm{b}$ involves many different contributions with various order-unity coefficients, which could add up or partially cancel. The main prediction of Eq. (2.16) is the scaling with $k_{\rm max}$.

Finally, let us consider systematics of Type 3, i.e. higher-order terms in the likelihood itself. Similar to the bias expansion, these come in two forms: an expansion in powers of $k$, equivalent to spatial derivatives; and an expansion in powers of perturbations, in this case the error field ${\varepsilon}(\bm{k})$ whose variance is $\sigma_{\varepsilon}^{2}$. Beginning with the former, a naive counting following loop contributions to the power spectrum indicates that a term $Ck^{2}$, where $C$ is a constant and which corresponds to the term $\propto\sigma_{{\varepsilon}}\sigma_{{\varepsilon},2}$ in Eq. (2.6), is of 2-loop order (see, e.g., Sec. 4.1. of [47]). Hence, terms of order $C^{\prime}k^{4}$ should have a negligible impact. This is corroborated by our numerical results (Sec. 5). The second type of higher-order stochasticity corresponds to non-Gaussian corrections such as the stochastic three-point function $\left\langle{\varepsilon}{\varepsilon}{\varepsilon}\right\rangle$ as well as coupling between stochasticity and the long-wavelength perturbations. These are briefly discussed in Sec. 5 of [54] (see also [55]). A derivation of the precise form of these contributions to the likelihood requires more theoretical investigation, and is left for future work. However, we can guess the approximate magnitude of these contributions by relating them to the terms we have kept here, which are nonlinear in the long-wavelength modes that determine ${\delta}_{\Lambda}$. The higher-order stochastic contributions are expected to be suppressed at least by

[TABLE]

where in the last equality we have assumed the Poisson expectation, $P_{\varepsilon}=1/\bar{n}$ where $\bar{n}$ is the mean number density of halos, which is a reasonable first-order estimate for this purpose. While the proper result will involve a summation over $\bm{k}$ modes similar to Eq. (2.16), we conservatively evaluate the ratio at $k_{\rm max}$ here, as it is unclear what precise weighting should be employed for this type of higher-order contribution. Notice that Eq. (2.18) also approaches 0 as $k\to 0$, but depends sensitively on the abundance of halos. In particular, it becomes large for small number densities.

2.4 The $\bm{\sigma_{8}}$ profile likelihood

We now introduce the framework used in our numerical tests to obtain maximum likelihood values and confidence intervals for cosmological parameters. Below, we center our discussion around the normalization of the primordial power spectrum, described by $\sigma_{\mathrm{8}}$. The reasoning for choosing $\sigma_{\mathrm{8}}$ as parameter is that it can only be inferred by using information in the nonlinear density field, as mentioned in Sec. 2.1. An unbiased inference thus means that the specific part of the information content in the nonlinear density field that is robust has been properly isolated. In particular, nonlinear information is explicitly necessary in order to break the degeneracy with the bias parameter $b_{1}$, rendering it the most direct test of our nonlinear inference approach. Future work will consider other cosmological parameters as well.

Besides from being a function of $\sigma_{\mathrm{8}}$, the marginalized likelihood Eq. (2.2) also depends on bias and other nuisance parameters (including the entire set of Fourier modes of the three-dimensional matter density field). Since the probabilistic inference of the initial matter density field from tracers like a halo catalog is numerically very expensive [see, e.g., 28], we instead constrain it to the actual initial conditions used in the simulations, evolved to low redshifts using either second-order Lagrangian perturbation theory (2LPT), or the N-body code directly. This forward evolution is then performed for a set of discrete $\sigma_{\mathrm{8}}$ values around the fiducial $\sigma_{\mathrm{8}}=0.85$ (see Sec. 4). We can then maximize the likelihood to simultaneously obtain best-fit values for cosmological and the remaining nuisance parameters. As mentioned in Sec. 1, this is in fact the most stringent test possible for any systematic bias, since the absence of any flexibility in the phases means that there is less room for errors in the likelihood to be absorbed by changing the initial conditions.

On the other hand, by fixing the phases, the only way to obtain rigours error estimates would be to analyze a large ensemble of large-volume simulations. Since these are costly, we here resort to a different method, allowing us to obtain error estimates from the likelihood itself: the profile likelihood method [58] provides estimates of confidence intervals for individual parameters of multivariate distributions within a frequentist approach. For a probability distribution $P(a,\{b_{i}\})$, the profile likelihood for parameter $a$ is defined as

[TABLE]

where the additional set of parameters $\{b_{i}\}$ has been profiled out. Constructing a full profile likelihood for $\sigma_{\mathrm{8}}$ is still numerically expensive, since it formally requires a recomputation of the final matter density field each time the function argument is updated. To speed up the analysis, we instead interpolate the profile likelihood evaluated on a predefined grid in $\sigma_{\mathrm{8}}$ centered about the fiducial value of the simulation. The details of this procedure will be described in Sec. 4.

3 Simulations

All numerical tests presented below are based on the N-body simulations presented in [59]. They are generated using GADGET-2 [60] for a flat $\Lambda$CDM cosmology with parameters $\Omega_{\mathrm{m}}=0.3$, $n_{\mathrm{s}}=0.967$, $h=0.7$, and $\sigma_{\mathrm{8}}=0.85$, a box size of $L=2000\,h^{-1}\mathrm{Mpc}$, and $1536^{3}$ dark matter particles of mass $M_{\mathrm{part}}=1.8\times 10^{11}\,h^{-1}M_{\odot}$. Initial conditions for the N-body runs were obtained at redshift $z_{\mathrm{ini}}=99$ using second-order Lagrangian perturbation theory (2LPT) [61] with the 2LPTic algorithm [62, 63]. Ref. [59] also presented runs for two further values of $\sigma_{\mathrm{8}}$ bracketing the fiducial value in order to perform the derivative of the halo mass function with respect to $\sigma_{\mathrm{8}}$ (for studies related to the scale-dependent bias induced by primordial non-Gaussianity). We use these, in addition to simulations for 4 additional $\sigma_{\mathrm{8}}$ values which were performed specifically for this paper, as well.

Dark matter halos were subsequently identified at redshift $z=0$ as spherical overdensities [64, 65, 66] applying the Amiga Halo Finder algorithm [67, 68], where we chose an overdensity threshold of 200 times the background matter density. The halo samples considered here consist of two mass ranges each at two redshifts, and are summarized in Tab. 1. We only consider halos above $10^{13}\,h^{-1}M_{\odot}$, corresponding to a minimum of 55 member particles. Note that differences in the number densities of the halo samples imply differences in the expected parameters as well as errors on the inferred $\sigma_{\mathrm{8}}$, an important aspect in our validation of the inference framework.

4 Implementation

We now provide additional details of the setup and numerical implementation used in our tests. We take as given the halo catalogs described in the previous section, as well as a set of matter particles generated either via 2LPT or full N-body for a set of values. Since the matter density field is a function of $\sigma_{\mathrm{8}}$, our main parameter of interest, mapping the profile likelihood as a continuous function of $\sigma_{\mathrm{8}}$ would require to recompute the set of operators for each function evaluation. To expedite the analysis, we instead generate representations of the evolved matter density field for fixed initial phases and a discrete set of values for $\sigma_{\mathrm{8}}$ (as well as redshifts $z=0$ and $z=1$ of the halo samples) given by

[TABLE]

For a given halo sample at a given redshift, and a fixed value $\sigma_{\mathrm{8}}^{i}$, the steps for evaluating the $\sigma_{\mathrm{8}}$ profile likelihood are as follows:

The halos and matter particles are assigned to a $1024^{3}$ grid using a cloud-in-cell density assignment. The high resolution is chosen to avoid leakage of the assignment kernel to the low wavenumbers of interest. 2. 2.

A sharp-$k$ filter is applied to the matter and halo fields in Fourier space, such that modes with $|\bm{k}|>\Lambda$ are set to zero. The Fourier-space grids are subsequently restricted to $384^{3}$, chosen such that the Nyquist frequency of each grid remains above $3\Lambda$ for all values of $\Lambda$ considered here. 3. 3.

Quadratic and higher-derivative operators are constructed from the sharp-$k$ filtered matter density field and held in memory. The quadratic operators are renormalized following Appendix A. 4. 4.

The maximum of the likelihood over the parameter space spanned by the remaining bias and stochastic variance parameters is then found via function minimizer MINUIT [69] (in practice we minimize $-2\ln\mathcal{L}$, i.e. the pseudo-$\chi^{2}$). The operator fields do not need to be recomputed for each evaluation, as only their coefficients are varied.

More precisely, we employ the analytic marginalization described in Sec. 2.2 for the parameters $b_{2}$, $c_{\nabla^{2}{\delta}}$ and $b_{K^{2}}$, leaving only $b_{1}$ and the three stochastic amplitudes to be varied in the minimization. We have not found any significant impact of the term $\sigma_{{\varepsilon}{\varepsilon}_{m},2}$ in Eq. (2.6), but a significant degeneracy with $\sigma_{{\varepsilon},2}$. For this reason, we fix the former to zero in our default analysis. This leaves a three-dimensional parameter space to be searched in the minimization, which typically converges quickly. We have found the minimization robust to varying initialization points and number of successive MINUIT cycles. Our default choice for the maximum wavenumber in the likelihood $k_{\mathrm{max}}$ is $\Lambda/2$.

This procedure results in a set of values $\{\sigma_{\mathrm{8}}^{i},\,-2\ln P^{\text{prof}}(\sigma_{\mathrm{8}}^{i})\}_{i}$ which we find is fit well by a parabola in all cases (we disregard values of $\sigma_{\mathrm{8}}$ where the minimization failed to converge). The best-fit value $\hat{\sigma}_{\mathrm{8}}$ is given by the minimum point of the best-fit parabola, while the estimated $1\sigma$ error on $\hat{\sigma}_{\mathrm{8}}$ is given by the inverse square-root of the curvature of the parabolic fit. We emphasize that this error does not include any residual cosmic variance, and is essentially purely governed by the halo stochasticity which appears in the variance of the likelihood.

5 Results

In the following, we present results for the best-fit value $\hat{\sigma}_{\mathrm{8}}$. For convenience, we phrase these in terms of the ratio to the fiducial value, introducing

[TABLE]

First, Fig. 1 shows examples of the profile likelihood determined as described in the previous section, with the parabolic fit that is used to determine $\hat{\sigma}_{\mathrm{8}}$. All panels are for $\Lambda=0.1\,h\,\text{Mpc}^{-1}$. Clearly, the log-profile likelihood is well approximated by a parabola, so that we expect maximum point and curvature to yield unbiased estimates of the maximum and 68%-level confidence intervals with respect to a full scan of the profile likelihood. The results for all halo samples and $\Lambda=0.1\,h\,\text{Mpc}^{-1}$, our fiducial choice, are summarized in Tab. 2. We find that an unbiased value of $\hat{\sigma}_{\mathrm{8}}$ is recovered to within $\sim 2\sigma$ in most cases. Notice that the run-to-run variance is larger than the estimated error bars in several cases. This could be due to residual cosmic variance, which is not contained in the estimated error bars as discussed in the previous section, although possible issues with the minimizer in isolated cases also cannot be excluded. In order to investigate this, more realizations would be needed.

The remaining columns in Tab. 2 show the value of $b_{1}$ as well as the stochastic amplitudes, all corresponding to the maximum-likelihood point for the fiducial value $\sigma_{\mathrm{8}}^{\text{fid}}$. Recall that all other parameters are analytically marginalized over. The bias $b_{1}$ is essentially fixed by the cross-correlation of ${\delta}_{h}$ with ${\delta}$. Correspondingly, we find that the combination $b_{1}\sigma_{\mathrm{8}}$ is constant for all $\sigma_{\mathrm{8}}$ values to within several percent. The stochastic amplitude $\sigma_{\varepsilon}^{2}$ is scaled to the Poisson expectation following Appendix B. Values greater (less) than one thus correspond to super- (sub-)Poisson stochasticity. We do find evidence for a smaller stochasticity for the rarer halo samples, in agreement with previous findings [70]. The last column shows the ratio of the higher-order (in $k^{2}$) stochastic parameter to the leading-order one. This gives a rough indication for the spatial length scale squared associated with the scale-dependent stochasticity. We thus find this length scale to be of order $(1-5)\,h^{-1}\text{Mpc}$. Notice however that this parameter is expected to also absorb various higher-order contributions not explicitly included in the likelihood, as discussed in [54], so that one cannot robustly infer a physical length scale from this value.

Allowing $\Lambda$ to vary, we obtain the results shown in Fig. 2. In each case, we show results for both simulation runs. While the differences in $\hat{\alpha}$ from unity are broadly consistent with being residual stochasticity and cosmic variance for $\Lambda\leq 0.1\,h\,\text{Mpc}^{-1}$, this no longer holds for higher cutoff values which at $z=1$ should still be under good perturbative control. On the other hand, the results appear remarkably stable toward higher values of $\Lambda$ up to $\Lambda=0.2\,h\,\text{Mpc}^{-1}$ in many cases. Notice that the majority of the modes contributing to the profile likelihood with $\Lambda=0.14\,h\,\text{Mpc}^{-1}$, say, are not included in the likelihood with $\Lambda=0.1\,h\,\text{Mpc}^{-1}$, so that these are largely independent measurements.

Before turning to possible explanations for these trends, let us briefly comment on the choice of $k_{\text{max}}/\Lambda$. We do not find strong trends with this parameter. Increasing $k_{\text{max}}$ at fixed $\Lambda$ yields similar trends as increasing $\Lambda$ itself, which is shown in Fig. 2. For this reason, we fix $k_{\text{max}}=\Lambda/2$ throughout.

As a test of the systematics of Type 1, we use the N-body density field itself instead of 2LPT for the construction of the bias operators. The result is shown in Fig. 3. We only find minor shifts in $\hat{\sigma}_{\mathrm{8}}$. Indeed, the cross-correlation coefficient between the 2LPT and N-body density fields in Fourier space is better than 0.97 for all scales and redshifts considered here, so a large shift would be surprising. Although there is some improvement, we conclude that the 2LPT density field is not primarily responsible for the bias in $\hat{\sigma}_{\mathrm{8}}$ found.

We next turn to systematics of Type 3, specifically the possible impact of the chosen implementation of $\sigma^{2}(k)$ in terms of two free parameters (since we fix $\sigma_{{\varepsilon}{\varepsilon}_{m},2}$ to zero). Since the halo number density is smaller at higher redshift, and the stochasticity correspondingly larger, this could possibly explain the increased bias in $\hat{\sigma}_{\mathrm{8}}$ at $z=1$ compared to $z=0$. Performing the profile likelihood analysis with varying $\sigma_{{\varepsilon}{\varepsilon}_{m},2}$ on the one hand, and both the former parameter and $\sigma_{{\varepsilon},2}$ fixed to zero on the other, leads to sub-percent shifts in $\hat{\sigma}_{\mathrm{8}}$ (this can also be gleaned by the values shown in the right-most column of Tab. 2, which, when multiplied by $k_{\text{max}}^{2}$ indicate the upper bound on the fractional contribution of the term $\sigma_{{\varepsilon},2}k^{2}$ to the total variance). Thus, we conclude that the parametrization of $\sigma^{2}(k)$ is unlikely to be responsible for the bias in $\hat{\sigma}_{\mathrm{8}}$ as well.

This leaves two possible sources of systematic error: our Type-2 systematic, i.e. higher-order terms in the bias expansion (both in perturbations and derivatives); and other systematics of Type-3, namely higher-order corrections to the form of the likelihood itself. Both types of terms are expected to be largest for the rarest and most highly biased halo samples. It is worth nothing that higher-order bias contributions are not necessarily smaller at higher redshift for fixed halo mass, since higher-redshift samples are more biased. Further, higher-derivative terms, which are expected to be tied to the Lagrangian radius of halos, most likely do not decrease toward higher redshift (see for eample the results of [71]). Indeed, preliminary investigations show that incorporating higher-order terms in the derivative expansion, such as $(\nabla^{2})^{2}{\delta}$ and $\nabla^{2}({\delta}^{2})$, in the set of operators does have an impact on the profile likelihood. We leave a detailed investigation of the impact of higher-order bias terms to upcoming work.

The remaining Type-3 systematics are expected to be controlled by the ratio in Eq. (2.18); this turns out to be of order unity or larger for the halo samples considered here, implying that higher-order stochastic corrections could be significant. At this point, lacking an explicit expression for these terms, it is difficult to quantitatively evaluate their impact however.

Additional investigations have pointed to a likely cause for the bias in the inferred value of $\sigma_{8}$ which is related to a higher-order term that is enhanced on large scales. Specifically, the auto-correlation of the quadratic bias operators $A_{OO^{\prime}}$ contains a contribution from the connected matter four-point function (trispectrum), which, for the likelihood and second-order bias expansion used here, corresponds to an error term $A_{OO^{\prime}}^{\text{err}}$ (Sec. 2.3). For small $k$, the dominant term from the particular trispectrum configuration involved scales as $\left\langle{\delta}_{\Lambda}^{2}\right\rangle^{2}P_{\rm L}(k)$, where $P_{\rm L}(k)$ is the linear matter power spectrum. While this contribution is suppressed compared to the leading, disconnected term, the latter approaches a constant at small $k$, while the trispectrum contribution grows toward small $k$ (assuming $k\gtrsim 0.02\,h\,\text{Mpc}^{-1}$) due to the factor $P_{\rm L}(k)$. For this reason, it can bias the maximum-likelihood value for $\sigma_{8}$ even for very small values of $\Lambda$ which correspondingly push $k$ to small values (see the behavior of $\alpha$ for $\Lambda\lesssim 0.1\,h\,\text{Mpc}^{-1}$ in Fig. 2). We leave a more detailed investigation, and possible remedies, for future work.

6 Conclusions

We have presented the results of a first application of the effective-field-theory-based Fourier-space likelihood derived in [54] to halo catalogs. More precisely, the test case is to obtain an unbiased estimate of the amplitude of the linear power spectrum $\sigma_{\mathrm{8}}$ (or equivalently normalization of scalar perturbations $\mathcal{A}_{s}$) purely based on the nonlinear information in the halo density field that is protected by the equivalence principle. For this, we vary four bias parameters as well as two stochastic amplitudes. The reasoning for choosing $\sigma_{\mathrm{8}}$ as parameter is that it can only be inferred by using information in the nonlinear density field. An unbiased inference—which we have not completely achieved yet—would thus mean that the robust, protected nonlinear information content has been isolated. We expect that other parameters, such as the BAO scale or the matter power spectrum shape, will then also be unbiased, as they rely to a lesser degree on purely nonlinear information.

We further presented a method to analytically marginalize over bias parameters, which we apply to three of the four bias parameters in our implementation. We expect that this analytical marginalization will prove extremely powerful when going to higher order in the bias expansion (both in orders of perturbations and derivatives): when using this technique, the cost of finding the maximum-likelihood point, or more generally, sampling from the likelihood, only increases quadratically with the number of bias terms (since the cost is dominated by the evaluation of the matrix $A_{OO^{\prime}}$).222This scaling is based on the fact that the computational time of evaluating the likelihood is dominated by the determination of the quantities $A_{OO^{\prime}}$ and $B_{O}$ in the notation of Sec. 2.2. The cost of the matrix inversion is negligible for a realistic number of operators. On the other hand, the computational cost would grow much more rapidly if one where to explicitly vary all bias parameters.

Our results indicate that $\sigma_{\mathrm{8}}$ can be recovered with a systematic error under $\sim 10\%$ for a range of halo samples at different redshifts when using a cutoff value $\Lambda=0.1\,h\,\text{Mpc}^{-1}$. Our assessment is that the most likely explanation for the residual systematic bias in $\sigma_{\mathrm{8}}$ can be traced to a trispectrum contribution to the quadratic operator correlators. Precisely how this contribution can be taken into account via higher-order terms in the bias expansion or likelihood is left for future work. It would further be interesting to perform a joint inference from the different halo samples that we have considered separately here, which however requires a generalization of the likelihood to include the stochasticity cross-covariance between different tracers. We leave these developments to future work.

Our approach has some resemblance to what was recently presented in [57]. Instead of the sharp-$k$ filter employed here (necessary for an unbiased inference following [54]), the authors of [57] used a Gaussian filter. More importantly, they allowd for the bias coefficients to be free functions of $k$, $b_{O}\to b_{O}(k)$, which removes the information on $\sigma_{\mathrm{8}}$, and instead focused on the degree of cross-correlation of the field we call ${\delta}_{h,\rm det}$ here with the actual halo density field. It would be interesting to study the corresponding correlation coefficient for our, sharp-$k$-filtered field ${\delta}_{h,\rm det}$ given the best-fit bias parameters. We defer this to future work as well.

A natural next question is: what is the expected statistical error for the inferred $\sigma_{\mathrm{8}}$ in the realistic case when the phases of the linear density field are unknown? In order to determine this, one unfortunately has to marginalize over those phases, which requires an implementation of the EFT likelihood into a sampling framework along the lines of, e.g. [28]. We thus have to defer this question to future work as well. It is clear however that this uncertainty will be very sensitive to the value of the cutoff $\Lambda$. Regardless of the expected statistical precision, it is worth emphasizing that, by fixing the phases to their ground-truth values throughout, the test of unbiased cosmology inference presented here is the most stringent test possible.

Finally, in parallel to the empirical studies on halo samples, a more rigorous theoretical study of the EFT likelihood expansion should be performed [55]. This is essential in order to obtain proper estimates of the relative size of the different expansion parameters, and to consistently carry out the expansion to higher order. After all, one of the main advantages of the EFT likelihood approach is that it is very simple to systematically go to higher orders.

Acknowledgments

We thank Alexandre Barreira, Giovanni Cabass, Stefan Hilbert, Marcel Schmittfull, Marko Simonović, and Matias Zaldarriaga for helpful discussions. We further thank Titouan Lazeyras for supplying us with the N-body and 2LPT codes used to generate initial conditions for the simulations our results are based on. FE, MN, and FS acknowledge support from the Starting Grant (ERC-2015-STG 678652) “GrInflaGal” of the European Research Council. GL acknowledges financial support from the ILP LABEX (under reference ANR-10-LABX-63) which is financed by French state funds managed by the ANR within the Investissements d’Avenir programme under reference ANR-11-IDEX-0004-02. This work was supported by the ANR BIG4 project, grant ANR-16-CE23-0002 of the French Agence Nationale de la Recherche. This work is done within the Aquila Consortium.333https://aquila-consortium.org

Appendix A Operator correlators and renormalization

As argued in [54] (see also [57, 56]), nonlinear operators constructed out of the matter density field must be renormalized in order to suppress dependencies of their cross-correlations on small-scale modes [72]. Specifically, using square brackets to denote renormalized operators $[O]$, we found using a tree-level (leading-order) calculation [54]

[TABLE]

where

[TABLE]

Here, the $W_{\Lambda}(k)$ are sharp filter functions defined in Fourier space. Finally,

[TABLE]

Since this calculation is only valid at leading order, we use a numerical renormalization procedure in our likelihood implementation instead. Specifically, we measure

[TABLE]

on a linear grid in $k$ (we choose 100 bins between the fundamental and Nyquist frequencies). The same is done for the density field itself, yielding $P_{{\delta}{\delta}}(k)$. Then, for each mode $\bm{k}$, we renormalize through

[TABLE]

where a cubic-spline interpolation is used to obtain the power spectra at each value of $\bm{k}$. Fig. 4 shows the cross-correlation of ${\delta}$ and the two quadratic operators ${\delta}^{2}$ and $K^{2}$ before and after renormalization. For the $k$ values that matter most in the likelihood, $k\gtrsim 0.02\,h\,\text{Mpc}^{-1}$, the cross-correlation is removed to high accuracy by the renormalization procedure. Also shown is the tree-level perturbation-theory prediction for the correlator before renormalization, which matches the measurement reasonably well, although not perfectly even at low $k$, since modes near the cutoff $\Lambda$ contribute to this cross-correlation.

Fig. 5 shows the cross-correlation of the quadratic operators among each other. As argued in [54], the renormalization also removes the dominant higher-order (trispectrum) contribution to these correlators. Indeed, the leading perturbation-theory prediction matches the cross-correlation of the quadratic operators well. We have verified that the good agreement also holds for other values of redshift and $\Lambda$, and that the deviations from the perturbation-theory prediction show the expected scaling, with agreement improving toward higher redshift and for smaller values of the cutoff $\Lambda$. We conclude that the operator cross-correlations, which form the practical basis of the EFT likelihood as discussed in detail in [54], are well understood.

Appendix B Interpreting the variance $\bm{\sigma_{\varepsilon}^{2}}$

In this appendix we derive the expectation for the variance parameter $\sigma_{\varepsilon}^{2}$ for Poisson noise. Neglecting long-wavelength perturbations, let $n_{i}\equiv n(\bm{x}_{i})$ denote the number density of halos in the grid cell centered around $\bm{x}_{i}$. Assuming this is Poisson distributed, we obtain

[TABLE]

where $N_{h}$ is the total number of halos in the box. The noise in the fractional halo density perturbation ${\delta}_{h}$ is then given by ${\varepsilon}_{i}=n_{i}/\lambda$, where we neglect the subtraction of the mean here since it is irrelevant for modes of finite $\bm{k}$. The noise field obeys, under the Poisson assumption,

[TABLE]

Finally, its power spectrum is given by

[TABLE]

where $\bar{n}_{h}=N_{h}/L_{\text{box}}^{3}$ is the number density of halos. Notice that the value depends on the grid resolution adopted, which in our implementation is $N_{g}=384$. The values of $\sigma_{\varepsilon}^{2}$ given in Tab. 2 are divided by this Poisson expectation. Values greater (less) than one thus correspond to super- (sub-)Poisson stochasticity.