BAO scale inference from biased tracers using the EFT likelihood

The physical scale corresponding to baryon acoustic oscillations (BAO), the size of the sound horizon at recombination, is precisely determined by CMB experiments. Measuring the apparent size of the BAO scale imprinted in the clustering of galaxies gives us a direct estimate of the angular-diameter distance and the Hubble parameter as a function of redshift. The BAO feature is damped by non-linear structure formation, which reduces the precision with which we can infer the BAO scale from standard galaxy clustering analysis methods. Many methods to undo this damping via the so-called BAO reconstruction have so far been proposed; however, they all rely on backward modeling. In this paper, we present the first results of BAO inference from rest-frame halo catalogs using forward modeling combined with the EFT likelihood, in the case where the initial phases of the density field are fixed. We show that the remaining systematic bias is less than 2% when we consider cutoff values of $\Lambda \leq 0.25 \,h\,{\rm Mpc}^{-1}$ for all halo samples considered, and below 1% and consistent with zero for all but the most highly biased samples. We also demonstrate that, when compared to the standard power spectrum likelihood approach under the same assumption of fixed phases, the 1$\sigma$ errors associated to the field level inference of the BAO scale are 1.1 to 3.3 times smaller, depending on the value of the cutoff and the halo sample. Our analysis therefore unveils another promising feature of using field-level inference for high-precision cosmology.


Introduction
Baryon acoustic oscillations (BAO) are an oscillatory feature in the matter power spectrum. The same feature is visible in the correlation function as a bump located at the characteristic BAO scale r s . The origin of the BAO can be found in the early Universe era when photons and baryons were tightly coupled by Compton scattering, forming the baryon-photon fluid. During this era, the gravitational force acting on the baryon perturbations was balanced by the radiation pressure resulting in acoustic oscillations of the baryon-photon fluid [1]. As the Universe expands and cools down, photons decouple. Traces of these sound waves remain visible as the acoustic oscillations in the CMB temperature anisotropies with the characteristic scale of the sound horizon at recombination, r s . Essentially the same scale is imprinted in the acoustic density perturbations in the baryon distribution. Since baryons are coupled to dark matter gravitationally, and both jointly evolve under gravitational evolution after decoupling, the imprint of these early-time oscillations is visible at fixed comoving scale in the late-time clustering of matter. Given that the size of the sound horizon at recombination has been well measured through CMB experiments, determining its apparent size in the late-time matter distribution allows us to estimate the angular-diameter distance and the Hubble parameter as a function of redshift. For a more detailed review of the BAO method see [2,3].
Before we can apply this method however, we have to face the problem that 1) matter evolved nonlinearly; 2) we do not directly observe the evolved matter density field, but rather biased tracers of this field such as galaxies, galaxy clusters, quasars and others. The distribution of such objects at low redshifts is affected by the highly non-linear structure formation, both of the matter distribution itself and of the formation of the tracers themselves (see [4] and [5] for reviews on these topics, respectively). For the BAO feature specificially, nonlinear structure formation shifts and broadens the peak in the correlation function, and equivalently dampens the oscillations in the power spectrum on small scales [6][7][8]. These effects reduce the precision with which the BAO can be measured from galaxy clustering [9] by relying only on information available in the power spectrum. As it was shown in [6], the dominant source of the broadening comes from bulk flows which are induced by large-scale modes, meaning that these effects can potentially be reversed. Therefore, there has been much interest in the BAO reconstruction methods [10][11][12][13][14][15][16]. These methods start by smoothing the galaxy density field to filter out high-k non-linearity, then this density is used to estimate the displacement field. Finally, the estimated displacement field is used to take tracers back to their estimated initial positions. Thus, all these methods rely on a reverse or backward modeling approach. In addition, they have to make assumptions about galaxy bias and the cosmological model to infer the displacement field, rather than inferring all parameters jointly with the BAO scale.
In this paper we instead use the forward model approach to test how well we can constrain the BAO scale by starting from the initial conditions. Forward modeling has gained a lot of momentum in the past few years [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32]. One of the main advantages of this approach is in the fact that it does not rely on the correlation functions, instead it exploits the amplitudes and the phases of the tracer field directly. This is done by writing down a joint posterior for the initial density field, cosmological parameters and nuisance parameters (bias parameters and stochastic amplitudes). A crucial ingredient in this posterior is the likelihood function of observing a tracer field δ h given the evolved matter density field. A likelihood function in the context of the effective field theory (EFT) of large scale structures [33,34] has been derived in [35][36][37]. A natural part of every EFT theory is a cutoff scale Λ which corresponds to the maximum wavenumber of modes included in the calculations. The precise value of Λ is arbitrary and its role is similar to the one of k max in standard power spectrum analyses. Crucially, the results of all measurements should be independent of Λ. A natural upper limit for Λ in the case of the EFT of LSS is the nonlinearity scale (Λ NL ≈ 0.25 h Mpc −1 at z = 0), where perturbation theory of LSS breaks down. If this condition is satisfied, a controlled inference of the cosmological parameters and initial conditions can be performed. The result of applying this EFT likelihood to σ 8 inference has been presented in [38,39]. A significant feature of the EFT likelihood is that it allows us to constrain the parameter of interest at the field level. Furthermore, Ref. [35] has shown that the forward model combined with the EFT likelihood naturally includes the BAO reconstruction. In this paper, we follow up and perform an unbiased inference of the isotropic BAO scale, which we will refer to simply as BAO scale, from rest-frame halo catalogs using the EFT likelihood. Our tests are based on simulations, which enables us to fix the initial phases of the linear density field to their correct values (to avoid any possible misunderstanding, fixing the phases of the initial density field means both fixing its amplitude and phase in each grid voxel). This removes cosmic variance as much as possible and reduces the size of the error bars.
This paper is organized in the following way. In Sec. 2 we give a short summary of the most important properties of the EFT likelihood. In Sec. 3 we give an overview of the method used for the inference of the BAO scale in the EFT based approach. In Sec. 4, we use the field-level EFT likelihood to find the BAO scale value. To gauge what are the improvements from the field level likelihood over standard, power-spectrum-based approaches, we also determine the BAO scale value using a likelihood constructed from the power spectrum (Sec. 5). For the predicted power spectrum we do not perform any additional BAO reconstruction; instead, we use the deterministic halo field found from the forward model. This means that the EFT likelihood is still used to constrain the bias parameters in the construction of the predicted power spectrum, but is not used to constrain the BAO scale. On the other hand, our power spectrum covariance takes into account that the phases are fixed to the ground truth, i.e. there is no cosmic variance. Thus, we perform a fair comparison between both methods. We conclude in Sec. 6.

The EFT Likelihood
We begin with a brief review of the EFT likelihood. Throughout this paper we refer to the tracers considered as halos, simply because we are working with halo catalogs from simulations. However, since the EFT approach only assumes that the formation of the tracer is spatially local, all the results are equally applicable to galaxies or any other cosmological tracer. With δ h we will denote the observed fractional number density perturbation of a given halo sample. In the rest frame of a halo, this density field is given by where τ is the conformal time, n h (x, τ ) denotes the comoving rest-frame halo density andn h (τ ) is its position-independent mean. In [35] a joint posterior for the initial density field δ in , cosmological parameters θ and nuisance parameters (bias parameters b O and stochastic amplitudes σ a ), P (δ in , θ, b O , σ a |δ h ), was introduced. All the important physics of halo formation is contained within the likelihood P (δ h |δ in , θ, b O , σ a ) giving the probability of observing the halo density δ h given the initial conditions, cosmological and nuisance parameters. Once the initial density field δ in has been specified, there are three important parts of P (δ h |δ in , θ, b O , σ a ) we need to focus on. Those are the deterministic forward model for matter, the bias relation and the conditional likelihood for finding a measured halo density field given the matter density field and bias parameters. In the following we summarize the most important information about each of them. The deterministic forward model for matter δ = δ fwd [δ in ] used in this paper is third-order Lagrangian perturbation theory (3LPT), as explained below. All initial perturbations with wavenumber k > Λ, where Λ is the initial cut-off, are set to zero. Further, we use the Lagrangian bias expansion where b L O and O L are the Lagrangian bias coefficients and operators, respectively, and q is the Lagrangian coordinate marking the initial position of the particle as τ − → 0. In the Lagrangian bias expansion, we first construct the bias operators and then displace them to Eulerian frame. This can be done conveniently in the same step as the LPT calculations. The relationship between the Lagrangian position q and and the final Eulerian position of the matter particles x is given through the displacement vector s, In LPT, we treat the components of the displacement tensor as small parameters, which allows us to write This expansion of the displacement tensor is related to the expansion in powers of δ performed in Eulerian Perturbation Theory (EPT). This can be seen from the following consideration. From Eq. (2.3), we find the Jacobian J ij to be where we have introduced the Lagrangian deformation tensor Then using the continuity relation together with Eq. (2.5), we find the relationship between the deformation tensor and the density field From Eq. (2.7), it is clear that at first order in perturbations The basis of the Lagrangian set of operators can be conveniently expressed in terms of the symmetric part of the Lagrangian deformation tensor [40] (2.9) The antisymmetric part, ∂ q,[i s j] , appears from third order in perturbations, but does not need to be included in the bias expansion as it is redundant [40]. We work by treating the components of the deformation tensor as small parameters and find the Lagrangian operators at each perturbative order by taking all the scalar contractions of M (n) ij . We do not need to include tr[M (n) ] for n > 1 since those can always be expressed in terms of scalars constructed using the lower order operators. Eq.
(2.8) is the starting point in a recursion relation that can be used to construct the tensors M (n) at all orders [41]. In our calculations, we will be using the operators up to third order in perturbations which are listed here according to their perturbative order [42] (2.10) The corresponding bias coefficients are (2.11) With these ingredients, we are finally able to construct the Lagrangian bias expansion presented in Eq. (2.2). The set of Eulerian operators in turn is then obtained by displacing each of the Lagrangian operators via Eq. (2.3). We define a grid of 512 3 cells in which O L (q i ) is set as the weight (or mass) of a "particle" at position q i . With the aim of preventing noise generation on large scales, we then deposit the particle mass at its Eulerian position x i using a could-in-cell scheme, such that the total mass is guaranteed to be conserved. This displacement technique is performed for all Lagrangian bias operators (up to the desired order) with the exception of tr[M (1) ]. For the latter field, we instead displace a unity-weight field to obtain the Eulerian density, which is associated to the well-known Eulerian bias parameter b δ commonly called b 1 . Thus, we will hereafter use b 1 ≡ b δ following standard convention. More details about the implementation of this procedure can be found in [39].
Note that the operator ∇ 2 x δ is not derived from the recursion relations arising from Eq. (2.8), but contains two more spatial derivatives. This operator is the leading higher-derivative operator which accounts for the non-locality of halo formation. The coefficients of higher-derivative operators are thus related to the spatial scale R * which quantifies the size of the spatial region involved in the process of halo formation, and their contribution to δ h,det is suppressed by powers of k 2 R 2 * on large scales.
Finally, we turn to the conditional likelihood which provides the probability for finding a measured halo density field given the matter density field and the bias parameters. The EFT likelihood which we use in this paper has the following form: The parametrization of σ 2 is chosen in such a way to ensure that σ 2 is positive definite, We can interpret σ ε as the amplitude of halo stochasticity in the large-scale limit (k − → 0). σ ε,2 is the leading scale-dependent correction to halo stochasticity and it captures the backreaction of small physical scales in real space. Since we found that σ ε,2 has a negligible contribution to our results, we set it to zero throughout the paper. A distinct feature of this likelihood is the existence of a hard cutoff Λ which marks a boundary above which all the modes k are integrated out. In other words, this cutoff ensures that we are focusing only on large scales. On such large scales, the central limit theorem guarantees that the noise fields for k < Λ can be approximated to leading order as independent Gaussian degrees of freedom, making the resulting likelihood normal with diagonal covariance (see Fig. 8 of [38] for an explicit demonstration). In our analysis we will be calculating the EFT likelihood for different cutoffs Λ to see how the results of the BAO scale inference are influenced by different modes included in the calculation. Since the number of modes scales as k 3 , and constraints are thus dominated by modes close to the cutoff Λ, the results for different Λ values will be essentially independent.
For the results in Sec. 4, we use the EFT likelihood to infer the value of the BAO scale. In this application, the likelihood is marginalized over all bias parameters as described in [36]. The EFT likelihood is also used in Sec. 5 to find the values of the bias parameters b O to be used in the deterministic halo power spectrum. To find the value of a parameter b O , we marginalize the EFT likelihood over all remaining bias parameters and then run the MINUIT minimizer [43] to find the best fit for the b O parameter. We repeat this procedure for all the fields appearing in Eq. (2.10) (see also [44,45] for a recent study of halo and galaxy bias using this method).

Method
The BAO is normally used to infer the angular diameter distance to a given observed redshift by comparing the predicted scale of the BAO feature to the data for a given assumed distance. This kind of approach is not suitable in our case since we are working with simulations on a cubic box with periodic boundary conditions. To change the distances inside such a box, we would have to introduce a window function, and would not be able to keep the initial conditions fixed to the ground truth (since changing the fiducial distance amounts to changing the comoving volume of the data as well). In order to avoid these significant complications, we adopt a different approach, essentially rescaling the predicted comoving sound horizon.

Approximating the power spectrum
We want to constrain the BAO scale r s just from the information available in the oscillatory part of the power spectrum, without referring to its broad-band part. This is because the broad-band power spectrum depends on other cosmological parameters as well. One possible way to constrain r s from the power spectrum is by varying the baryon density ω b and checking which value agrees the best with the data. However, varying the value of ω b changes not only the oscillatory part of the power spectrum, but also its broad band. Therefore, this approach is not suitable.
Instead, we approximate the linear matter power spectrum as where A and k D are constants and r fid is the fiducial BAO scale. Through this equation we separate the broad band part of the power spectrum, described with the function P L,sm (k), from its oscillatory feature. In the oscillatory feature we recognize the contribution sin(kβr fid ) describing the baryon acoustic oscillations and the exponential envelope corresponding to the primordial photon diffusion, or Silk damping. The Silk damping term absorbs all the physics that is not captured within the fluid approximation to the baryon-photon system before recombination. Finally, we introduced the factor β as By changing β, we are changing the size of the BAO scale r s to match the data while keeping the distances fixed. Most importantly, changing β will result in changes in the oscillatory part of the power spectrum while keeping its overall shape intact. Notice that, since the BAO scale was imprinted in the power spectrum during the early Universe, varying it in the initial (linear) density field is the physically correct approach. The function P L,sm (k) can be written in the form where N is a normalisation constant that is proportional to the primordial normalization A s times the growth factor squared, k p is the pivot scale and T (k) is the "no-wiggle" transfer function which  we take from Ref. [1]. We found the value of N , k D and A by fitting Eq. (3.3) to the linear power spectrum produced by the CLASS code [46]. Fig. 1a shows the ratio of the CLASS power spectrum to P L,sm (k). We can clearly see the damped oscillation in the BAO range which indicates that P L,sm really does describe the smooth power spectrum with no BAO wiggles. Fig. 1b shows the ratio of the CLASS power spectrum to the power spectrum approximated by Eq. (3.1). While we do see some residual wiggles in the plot, we also notice that they are suppressed at high k where most of the constraints come from. Therefore, we can conclude that Eq. (3.1) is a good approximation of the linear power spectrum and we can use it for the BAO scale inference. Given the known fiducial power spectrum, i.e. the power spectrum from which the initial conditions of the N-body simulations were drawn, it is easy to find the power spectrum with a different BAO scale, using Eq. (3.1). We introduce the factor f (k, β) as Notice that f (k, 1) = 1. From f (k, β) it is straightforward to find the relationship between the fiducial and rescaled linear density fields To recap, δ β is the linear matter density field with all fiducial phases but for which the BAO scale is of the size r s = βr fid . Throughout the paper we will be using different δ β as the initial fields for our forward model.

Profile likelihood
All numerical results presented here were obtained for a spatially flat ΛCDM cosmology with parameters Ω m = 0.3, Ω Λ = 0.7, h = 0.7, n s = 0.967 and a box with size L = 2000h −1 Mpc. We use four halo mass bins in the mass range 10 12.5 h −1 M ⊙ -10 14.5 h −1 M ⊙ . We present results on two simulation realizations, "run 1" and "run 2", which differ in their initial phases. In Tab. 3.5, we present the number density of halos in run 1 at different redshifts.
z As mentioned earlier, we do not sample the initial density field; it is instead fixed to the exact initial conditions used in the N-body simulations within which the halos are identified. To get the initial density with different BAO scales, we apply Eq. (3.5) for a set of values {β i }. The default set spans the range [0.8, 1.02]; in all cases, we make sure that the maximum-a-posteriori (MAP) value of β is safely within the range. Fixing the initial phases not only saves the computational time, but it also minimizes the cosmic variance as much as possible resulting in smaller error bars for the inferred value of β.

Mass rangen
To find the MAP estimate for β, which we denote asβ, we use the profile likelihood [47]. For a probability distribution P (β, σ ε |δ h ) (recall that the bias coefficients are analytically marginalized over) and parameter β, the profile likelihood is defined as where the parameter σ ε has been profiled out. For a fixed Λ, halo sample, redshift and β i we maximize the profile likelihood using the MINUIT minimizer [43]. In this way we obtain a set {β i , −2 ln P prof (β i )} which is nicely fit by a parabola for all halo samples and all cutoffs. An example of this parabola for two different cutoffs is shown in Fig.  2, where the elements of the set {β i , −2 ln P prof (β i )} are represented with orange dots and the blue line corresponds to the parabolic fit. The MAP valueβ is located at the minimum of the best fit parabola. The estimated 68% confidence-level error onβ is given by the inverse square root curvature of the parabolic fit.

Field-Level Results
In this section, we show the results of applying the EFT likelihood to the halo catalogs. We start by comparing the results for two different bias orders-second and third order-at fixed redshift z = 0. Fig. 3 shows the deviation of the MAP valuesβ from 1 as a function of Λ for different halo mass ranges. For all of halo mass bins except the highest one,β is consistent with being unbiased within the error bar obtained from the profile likelihood. Moreover,β is moving closer to 1 as Λ is increased, consistent with the shrinking error bar as more k modes are being included in the likelihood and forward model. We notice that the MAP valuesβ are closer to 1 for the third order bias expansion than in the case of second order, for every halo sample. This indicates that the systematic error in β in the 3rd order bias case is reduced, as expected if one is in the converging regime of the EFT. Therefore, in the rest of the paper, we focus only on the 3rd order bias expansion. Fig. 4 depicts the value of the 1σ error bar, σ F (β), for the field-level inference (as emphasized by the subscript F ) as a function of Λ for the 3rd bias order for both runs 1 and 2. We see that σ F (β) is smoothly decreasing with increasing Λ. Since our initial conditions are exactly the ones used in the halo simulations, the statistical uncertainty σ F (β) is only sourced by the halo stochasticity which appears in the EFT likelihood. Note that we do expect the σ F (β) results to change once we start sampling the initial phases instead of keeping them fixed.
We also notice that the σ F (β) values do not change much between the different halo samples. This trend can be understood by inspecting how the numerator and denominator of Eq. (2.12) change with halo mass. On the one hand, more massive halos are rarer, and hence have larger noise (stochasticity), i.e. larger σ 2 ε in the denominator. On the other hand, higher-mass halos are more biased, and hence show a stronger clustering signal. Hence the numerator also increases with halo mass. As a consequence, the ratio of both quantities is actually roughly constant, so that we get a similar σ F (β) for all halo bins considered. Notice that this result only applies to the fixed-phase study done here.  Next, we look at the results found at different redshifts. Fig. 5 and Fig. 6 show the value ofβ as  a function of Λ at redshifts z = 0.0, 0.5, 1.0 for run 1 and run 2, respectively. For run 1, these results are also summarized in Tab. 2 for a fixed cutoff, Λ = 0.16 h Mpc −1 . We find that the remaining systematic bias is very low across all the redshifts and mass ranges. Even for this cutoff, the bias is less than 1% at z = 0 for all mass ranges and the results generally keep improving with growing Λ. If we look across all redshifts, it is clear that the remaining bias is still below 2%, and in fact consistent with zero, for most of the cases. It goes over 2% only for the most massive halos which are very rare at higher redshifts. Those halo samples would most likely benefit from from going to higher bias orders in the bias expansion. From Fig. 5 and Fig. 6 we notice that the remaining systematic bias inβ is increasing with growing redshift. This occurrence is counter-intuitive, since from perturbation theory we would expect a better performance at higher redshifts where perturbation theory extends to higher wavenumbers. A similar trend was noticed with inference of σ 8 from the halo catalogues described in [38]. It was found there that this trend is caused by the higher-order bias terms. Although the higher-order bias terms are suppressed by powers of the normalized growth factor D norm (z) = D(z)/D(0) at higher redshifts, it is possible that the increase in their coefficients with redshift more than compensates for this suppression. To check if this was the case for us as well, we use the test suggested in [38] which was based on the assumption from [42] that the higher order bias terms can be approximated as being a function of (b 1 − 1)D norm (z). Results are shown in Fig. 7, where we plotted the |β − 1| values for all halo mass bins and redshifts against (b 1 − 1)D norm (z). There is a hint of a correlation between |β − 1| and (b 1 − 1)D norm (z), although all but one points are consistent withβ = 1 within one sigma.
Let us also comment on the limits of the cutoff we are using. For matter, the EFT is under perturbative control for Λ 0.25 h Mpc −1 at z = 0. For highly biased tracers, the cutoff is reduced due to the growing size of bias parameters at higher orders. Thus, we are going beyond that limit, and not all of our values of Λ are strictly under perturbative control. However, because the BAO is an oscillatory feature, while higher-order corrections are expected to be smooth functions of k, the BAO inference seems to be still robust at these high k. We leave a more systematic investigation of this important issue to future work.      Table 2. Summary of the results found using the field-level EFT likelihood at the cutoff Λ = 0.16 h Mpc −1 for different redshifts and halo mass bins.

Comparing the field-level results to the power spectrum approach
Having presented the results of constraining the BAO scale using the EFT likelihood, we now turn to comparing these results to a more traditional BAO inference approach based on the power spectrum.

Power spectrum likelihood
Care is needed in order to ensure that the comparison we are making is valid, since in the EFT approach we use fixed phases in the matter density field. Therefore, we adopt the following Gaussian likelihood for the halo power spectrum: is the power spectrum of the deterministic halo field found using the same forward model as in EFT case for a fixed β value; P h (k) is the measured halo power spectrum, m k is the number of modes in a wavenumber bin, and P ε is the noise spectrum. Notice that the covariance appearing in the numerator of the likelihood, Var fix [P h (k)], is modified to reflect the fact that we are using fixed phases. The derivation of the power spectrum covariance for fixed phases can be found in Appendix A and its final form is given in Eq. (A.18). It is also important to note that we are not performing any additional BAO reconstruction on the halo data, but comparing the halo power spectrum directly with the theory predictions from the full forward model. Therefore the comparison we are making is at the level of likelihoods: the EFT likelihood is performing at the level of the field, while the likelihood in Eq. (5.1) compresses the data to the power spectrum in bins of k. Both likelihoods however consistently assume fixed initial conditions.
To find the best fit for β, we use the following procedure for different values β i . We start by finding the initial matter fields with the BAO scales r s = β i r fid using Eq. (3.5) as in the field-level likelihood calculations. Once we have the linear matter density field δ in (k, β i ) , we use the 3LPT forward model to generate the evolved matter field δ(k, β) = δ fwd [δ in (k, β i )], where we set all modes with k > Λ to zero. For the bias operators, we use the same bias model as described in Sec. 2. The MAP for the bias parameters is found by maximizing the EFT likelihood. We keep one bias parameter free at a time and marginalize over all other bias coefficients. Once we found the MAP value for that parameter, we move on and repeat the procedure for the remaining ones. This gives us the deterministic halo field whose power spectrum P det (k|δ in , β, {b O }) is straightforward to measure in the same k bins as the halo sample.
We now turn to the determination of P ε . Ideally, one would fit for this together with β and the bias parameters. In our simplified analysis, we only fit P ε , and use the same noise spectrum value P ε across all Λ and β values. This value is found for Λ = 0.2 hMpc −1 and β = 1.00 by fitting the difference P h (k) − P det (k|δ in , β = 1, {b O }) to a constant, using w = 1/σ w as the weight where σ w = |P h − P det |/ 2/m k . Fitting the noise separately from bias terms and β leaves us with some uncertainties in its estimate. We roughly estimate this uncertainty by repeating the same analysis for run 2 halo samples, resulting in values of P ε that differ by around 20%, which results in a corresponding 20% shift in the 1σ error forβ. We conclude that our results for the latter carry an uncertainty of ∼ 20%. This is sufficient for the approximate comparison we are aiming for in this paper. We aim to improve this in future work.
Finally, by inserting P det (k|δ in , β i , {b O }) and P ε in Eq. (5.1), we find the likelihood value for each β i . Repeating this procedure at fixed halo sample, redshift and Λ, leads to a set {β i , −2 ln P prof (β i )}, which is nicely fit by a parabola. An example of this parabola fit is shown in Fig. 8.β and σ P S (β), the value of the 1σ error bar for the power spectrum inference, are found as the location of the minimum and the inverse square root of the parabolic fit, respectively.

Results
We now turn to the results for MAP.β is found using the likelihood given in Eq. (5.1). The residual values ofβ as a function of Λ at the three different redshifts are shown in Fig. 9. For the most massive  Table 3. MAP values of β for cutoff Λ = 0.16 h Mpc −1 inferred from the power spectrum likelihood, at different redshifts for different halo mass ranges.
halo range log 10 (M/h −1 M ⊙ ) = 14.0 − 14.5, we show results only at redshift zero. For this halo range at higher redshifts, the set {β i , −2 ln P prof (β i )} does not yield a well-defined maximum. We also exclude all the samples for which the MINUIT algorithm does not converge for the bias coefficients due to a poor signal to noise ratio. The quantitative results are summarized in Tab. 3. We see that, for most of the samples, the residual bias inβ is between 1.20% and 2.3%. The MAP values of the linear bias parameter b 1 are also listed in the table. We notice that, for a fixed mass bin, b 1 is increasing with halo mass and redshift as is expected. In Fig. 10, we show σ P S (β) as a function of Λ at redshift z = 0. While for the field-level likelihood σ F (β) reduces about 2.4 times from Λ = 0.1 h Mpc −1 to Λ = 0.2 h Mpc −1 , here we do not see such a trend. Instead, σ P S (β) stays fairly constant across all Λ for the power spectrum likelihood. This is presumably because the field-level likelihood can still make use of the phase information at wavenumbers for which the power spectrum likelihood is already dominated by the noise P ε .
The most interesting result is Fig. 11, which compares the error onβ from the power spectrum approach, σ P S (β), to the one from the field level approach, σ F (β). This ratio is shown for three different halos mass ranges at three different redshifts. For smaller cutoffs, both likelihoods give similar results, which is the expected result if the data (δ h ) are well approximated as a Gaussian random field. However, as Λ grows, the EFT likelihood starts outperforming the power spectrum  based likelihood. At the highest Λ considered, the σ F (β) value is around 2.5 times smaller than σ P S (β). The field-level EFT likelihood performs better because it operates at the level of the field. This means that it includes not only all the information coming from the power spectrum, but also information from the from N-point functions of arbitrarily high orders. Concretely in the case of the BAO, the field-level likelihood knows about the bulk flow field, and can thus compare the expected BAO scale at a given location with the data. The power spectrum on the other hand is averaged over all locations, and thus suffers from the damping of the BAO peak [7,8]. Thus, the fact that the field-level likelihood outperforms the power spectrum based one comes as no surprise. Finally, note that we have fixed the bias coefficients in the theory prediction for the power spectrum to the values obtained from the field-level likelihood. In practice, those would have to be marginalized over in a power spectrum analysis.   Figure 11. Ratio of the uncertainty on the BAO scale inferred from the power spectrum likelihood, σP S (β), to that from the field-level likelihood, σF (β), as a function of cutoff for different redshifts. Each panel corresponds to a different halo mass range.

Summary and Conclusions
In this paper we compared the inference of the BAO scale from the halo catalogs using an LPT-based forward model combined with the EFT likelihood with the standard approach which compresses the data to the power spectrum. The forward model uses a combination of 3LPT expansion for the matter field and a third-order bias expansion. Our results were expressed in the terms of the parameter β, defined as the ratio of the measured value of the BAO scale to its fiducial value.
The field-level inference results are summarised in Fig. 5 and Tab. 2. From these it is clear that the remaining systematic error inβ is at most ∼ 2% for all samples. If we ignore the most biased sample (log 10 (M/h −1 M ⊙ ) > 14.0 at z = 0.5 and z = 1.0), at Λ = 0.16 h Mpc −1 , the remaining bias in β is less than 1% for all remaining samples, which is remarkably low. Moreover,β is statistically consistent with being unbiased for all halo samples except the highest mass one. It is also interesting to notice that the bias inβ is under control for all halo mass ranges, even for Λ = 0.25 h Mpc −1 , which is close to the nonlinear scale. For the lighter halos, log 10 (M/h −1 M ⊙ ) < 13.5, this even applies to While we consider halos here, the EFT approach is equally applicable to galaxies. This is confirmed by the results of [44], who demonstrated an unbiased inference of the linear power spectrum normalization σ 8 on fully hydrodynamical simulated galaxies.
In order to assess the performance of the field-level inference of the BAO scale, we compare it to the more traditional approach of the BAO inference from the power spectrum. For this we utilized the likelihood defined in Eq. (5.1), where the theory model for the power spectrum is based on the same forward model as in the EFT likelihood approach (in particular, the field-level likelihood was still used to find the best-fit values for the bias parameters). We modified the covariance in this likelihood to reflect the fact that we are using fixed phases, so that we can compare the two approaches on the same footing. Results found using this likelihood are shown in Fig. 9 and Tab. 3. For a fixed Λ = 0.16 h Mpc −1 , the remaining systematic errorβ is between 1.2% and 2.3% for those halo samples that yielded converged profile likelihoods. Fig. 11 shows the relative performance of the field-level and power spectrum based likelihoods. Across all halos samples, we notice a similar trend. For smaller cutoffs, both likelihoods show similar performance. However, for Λ > 0.12 h Mpc −1 , the field-level likelihood gives better results across all the halo masses and redshifts. For the highest cutoff value considered in both likelihoods, Λ = 0.2 h Mpc −1 , the error on the BAO scale inferred from the power spectrum is between 2.47 − 3.3 times larger than that obtained from the field-level likelihood, depending on the halo sample. Since the field-level likelihood contains all the information that would come from the higher order correlation functions, including the precise bulk-flow field, while the information available in the likelihood from Eq. (5.1) are only those from the power spectrum, a better performance of the EFT likelihood was to be expected.
In future work we will investigate how well we can constrain the BAO from the EFT likelihood in the cases when the initial conditions are not fixed, but sampled. This will allow for a realistic comparison of the constraining power on the BAO scale that can be obtained from the field-level inference as compared to that based on the galaxy power spectrum.

A Power spectrum covariance for fixed phases
In this section we derive the power spectrum covariance in the case where the initial phases are fixed. Inside a thin shell bin of magnitude k, which we keep fixed throughout, the prediction for the halo power spectrum can be written as with the sum running over all the modes q inside the bin of magnitude k.