An Optimally Weighted Estimator of the Linear Power Spectrum Disentangling the Growth of Density Perturbations Across Galaxy Surveys

Measuring the clustering of galaxies from surveys allows us to estimate the power spectrum of matter density fluctuations, thus constraining cosmological models. This requires careful modelling of observational effects to avoid misinterpretation of data. In particular, signals coming from different distances encode information from different epochs. This is known as"light-cone effect"and is going to have a higher impact as upcoming galaxy surveys probe larger redshift ranges. Generalising the method by Feldman et al. (1994), I define a minimum-variance estimator of the linear power spectrum at a fixed time, properly taking into account the light-cone effect. An analytic expression for the estimator is provided, and that is consistent with the findings of previous works in the literature. I test the method within the context of the halo model, assuming Planck 2014 cosmological parameters. I show that the estimator presented recovers the fiducial linear power spectrum at present time within 5% accuracy up to $k \sim 0.80\;h\,\mathrm{Mpc}^{-1}$ and within 10% up to $k \sim 0.94\;h\,\mathrm{Mpc}^{-1}$, well into the non-linear regime of the growth of density perturbations. As such, the method could be useful in the analysis of the data from future large-scale surveys, like Euclid.


Introduction
Over the last few decades, more and more precise measurements of the Cosmic Microwave Background (CMB) have shown that the Universe was almost isotropic at the epoch of last scattering, presenting fluctuations in the temperature of the CMB of order 10 −4 K [2][3][4]. This remarkable discovery has led to the description of the Universe as initially homogeneous, except for tiny perturbations in the density of matter, which were set by inflation. Such seed fluctuations have then grown due to gravitational instability, eventually giving rise to the large-scale structure observable today.
The evolution of photons, dark matter and baryons is governed by the Einstein-Boltzmann equations. As long as matter density perturbations are small, they can be linearised around the solution relative to a perfectly homogeneous and isotropic Universe. In this regime the growth of fluctuations is relatively straightforward, since their amplitude is simply enhanced by a time-dependent growth function, set by the cosmological parameters. If the fluctuations are not much smaller than the mean density, one needs to perturb Einstein-Boltzmann equations to higher orders, or to abandon perturbation theory in favour of semianalytic or numerical methods.
Within or outside the regime of applicability of linear theory, the choice of the cosmological model impacts the observed large-scale structure. So, by observing structures today, one can infer how likely a certain model is, and constrain its parameters. Galaxy surveys play a major role in this respect. The distribution of galaxies is thought to be a biased tracer of matter and its number density is described as a Poisson sampling of the underlying matter distribution. Consequently, mapping the galaxy distribution tells us also how matter is distributed. More importantly, the statistical properties of galaxy clustering represent a formidable source of information.
Two widely used statistics are the two-point correlation function and the power spectrum. They are essentially the same statistics, expressed in real and Fourier space, respectively. In fact, the correlation function is defined as ξ(r − r ) = δ m (r)δ m (r ) , (1.1) where δ m (x) is the density fluctuations field. The power spectrum is nothing but its Fourier transform 1 P (k) = d 3 r e ik·(r−r ) ξ(r − r ) . (1.2) Both correlation function and power spectrum are affected by cosmological parameters. This is the reason why these functions are so important in cosmological studies. They can be predicted from observations and then used to constrain models. A systematic study of the statistical properties of clustering was initiated by [5,6] and [7]. These pioneering works informed an entire research field, focusing on the understanding of the two-point correlation function (e.g. [8][9][10][11][12][13][14][15][16][17][18][19], and references therein). The two-point correlation function contains full information if the distribution of the amplitudes of the primordial perturbations is Gaussian, otherwise it contains only partial information, but still represents a good starting point. Formally, studying the power spectrum instead of the twopoint correlation function is absolutely equivalent, the former being the Fourier transform of the latter. However, this is not necessarily true when considering the estimates of these two quantities, which are obtained from finite and noisy observational data. The statistical analysis of galaxy surveys through the power spectrum represented indeed a flourishing area of research (e.g. [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34], and references therein). An advantage of studying the power spectrum with respect to the two-point correlation function is that it is always positive, so it is easier to detect any non-physical result. Besides, its shape does not depend on the knowledge of the mean number density of galaxies, which affects the k = 0 mode only. Also, the power spectrum is the quantity which is generally provided by theories studying the early Universe (e.g. inflation). [1] It is of utmost importance to develop solid methods to infer the power spectrum from galaxy surveys. A fundamental work in this respect was published by [1], who defined an estimator of the power spectrum of galaxy clustering, given a certain galaxy distribution. The estimator is defined such that its variance divided by the power spectrum is minimised. I shall refer to [1] as "FKP" throughout this work.
The FKP method is subject to some approximations and limitations. In particular, with galaxy surveys becoming wider and deeper, providing the correct estimate of such quantity is not trivial at all. One must indeed take into account various observational effects, like the redshift evolution of the luminosity function clustering amplitude and bias (e.g. [12,33,[35][36][37]), the contribution of modes larger than the volume of the survey (e.g. [38]), and generalrelativistic effects (e.g. [39,40]). Various extensions of FKP have been realised, incorporating for example the light-cone effect [41] and redshift space distortions [42]. More sophisticated analyses followed, based for example on non-Gaussian contributions to likelihood analysis (e.g. [43]), Bayesian or hybrid approaches (e.g. [44,45]), N-body simulations (e.g. [46]), and multipole-analysis of the power spectrum [47,48].
In this paper, I shall focus on the light-cone effect. With this term I mean the fact that one is able to observe only one's own past light-cone. So, if the survey is particularly deep, one may receive signals coming from epochs spanning a significant lapse of time. Consequently, one is not really able to map the galaxy distribution, and thus compute its power spectrum, at a certain fixed time. A possible way to deal with this issue is looking for an estimator of the power spectrum, somehow averaged over a certain time lapse. This approach was adopted by [41], who was able to generalise the FKP method including, among others, the light-cone effect.
In the current paper, a different approach is attempted. I define an estimator of the power spectrum at a certain fixed time, including the light-cone effect in the FKP method, under certain assumptions on the behaviour of the correlation function. Section 2 contains the main result of this work. In that section, I present my generalisation of FKP, and in section 3 how I implemented it in a test-case to determine its accuracy. In sections 4 and 5 I present and discuss the results, respectively. In section 6 I state my conclusions and give an overview on the possible perspectives of this work. Unless otherwise indicated, all units mentioned in the text are comoving.

Description of the Method
In this section I generalise the FKP method, considering the growth of density perturbations across a galaxy survey. To do that, I will follow the logic and notation of the original FKP paper. In section 2.1 I introduce the formalism necessary to address the problem, and show that, within linear perturbation theory, generalsing FKP is relatively straightforward. There is no original result in subsection 2.1. It serves only as an introduction to the subsequent subsections, in which I consider non-linear growth of pertrubations. That is the most interesting case, since the applicability of linear theory is limited to large scales only. The method described in subsections 2.2-2.6 represents original work.

Growth of Perturbations
When performing observations, one actually receives light signals coming from a range of distances and thus from different epochs. The measured overdensities in two points r 1 and r 2 are then δ(r 1 , t 1 ) and δ(r 2 , t 2 ). In contrast, one is interested in determining δ(r, t 0 ), i.e. the overdensity field at different positions, at the same reference time t 0 (for example, today). It is then necessary to develop a formalism to connect the overdensity at the generic time t to the overdensity at t 0 .
Making the time-dependence of the overdensity field explicit, the number density of galaxies in a survey is given by wheren(r, t) is the expected number density of galaxies, given the luminosity and angular selection criteria of the survey. Since time is a function of the distance r, the explicit time dependence in equation (2.1) may be suppressed. Within linear theory, equation (2.1) can be written as where D(r) is the growth function [49], and δ 0 (r) is the overdensity at the reference time t 0 = t(r 0 ). I define the corresponding weighted galaxy fluctuation field as where w(r) the weight function and n s (r) is the synthetic catalogue. The constant α is a rescaling factor to normalise the total number of galaxies in the synthetic catalogue to the real catalogue. Since the synthetic catalogue has in general many more galaxies than the real one, α is usually small. In the limit of infinitely many synthetic galaxies, α tends to zero. The optimal estimator of the linear power spectrum at fixed time can be immediately obtained noting the analogy with the work published by Percival, Verde and Peacock in 2004 (hereafter, PVP) [37]. Their work focused on an extension of the FKP method considering a set of galaxies of a certain luminosity, from a Poisson sampling of a linearly biased density field. In PVP, the density of galaxies with luminosity L is given by where b(r, L) is the bias. Comparing the definition of the galaxy density in equations (2.4) and (2.2), it is clear that the linear growth of density perturbations can be formally treated as a linear bias. Therefore, it is sufficient to follow PVP replacing b(r, L) with D(r)/D(r 0 ). Following FKP, the estimator P 0 (k) of the linear power spectrum at fixed time t 0 is defined as where P shot is the shot noise. The resulting optimal weight function is .
(2.6) Equation (2.6) is consistent with the expression found in PVP. They are not exactly the same because I adopted a different definition of the galaxy density fluctuations in (2.3). In the context of this manuscript, that definition is preferable, since it allows for a more straightforward comparison with the main results of this work, which are presented in subsection 2.6. Nevertheless, it can be verified that, adopting the same definition of F (k) as in PVP, the expressions for the optimal weight function are fully consistent. Assuming linear growth of perturbations makes the estimate of the fixed-time linear power spectrum from galaxy surveys particularly straightforward. However, this assumption restricts the validity of the method at large scales only. It is then interesting to develop a formalism to estimate the power spectrum at fixed time, considering non-linear growth of perturbations across the galaxy survey. In the following subsections, I shall generalise the FKP method to the mildly non-linear regime. I shall call P (k) the fully, non-linear power spectrum, and ∆ 2 (k) = k 3 P (k)/2π 2 the dimensionless non-linear power spectrum. In the remainder of the paper, I shall denote with ∆ 2 L (k) the dimensionless power spectrum predicted by linear theory. When referring to ∆ 2 (k) and ∆ 2 L (k), I shall often omit "dimensionless". That should not be confusing, as it should be clear from the context whether I am referring to P (k) or ∆ 2 (k).
There have been various attempts to get an approximate analytic expression for nonlinear power spectrum. The most straightforward way to accomplish that is perhaps perturbation theory, where the equations for the evolution of the matter density field are perturbed beyond first order (see [50] for a review). A different idea to address the issue is expressing the non-linear power spectrum (or correlation function) as a function of its linear counterpart, evaluated at a different scale [51][52][53][54][55]. Another approach is assuming a power law for the shape of the primordial power spectrum and the time evolution of the scale factor, and then determine the time evolution of the scale at which perturbations become non linear. It can be shown that, under certain conditions, the power spectrum evolves in time according to a self-similar solution [56][57][58][59]. A different widely used analytic method to describe nonlinear clustering is the "halo model" ( [60][61][62][63][64][65]; see [66] for a review). The basic idea consists in splitting the power spectrum into two terms. One of them, the "two-halo" term, models large-scale correlations in the density of matter residing in different halos. The other one, the "one-halo" term, represents the contribution to the power spectrum due to correlations of the density field within the same halo. The parameters of the analytic functional form of the power spectrum provided by the halo model can be determined fitting results of N-body cosmological simulations [49]. This semi-analytic method is called "halofit" model [67]. In this work, I will consider the revised version of the halo model provided by [67].

Revised Halofit Model
I shall now briefly review the revised halofit model by [67]. Let us define where ∆ 2 Q (k) and ∆ 2 H (k) are the two-halo and one-halo terms, respectively. Unlike in linear theory, they cannot be expressed as a product of a function of k only and a function of time (or, equivalently, distance or redshift) only. Specifically, In the equation above, y = k/k σ , with k σ defined as σ 2 (k −1 σ ) = 1, where The effective spectral index n ef f and the curvature C are defined as All other parameters are determined as polynomial functions of n ef f and C. The coefficients of the polynomials are obtained by fitting results of N-body simulations. The details can be found in [67]. The halofit model can be visualised in figure 1. The solid red line is the linear power spectrum generated through CAMB ( [68,69]) with the cosmological parameters given by [2] (Ω m = 0.315, Ω Λ = 1 − Ω m = 0.683, Ω b = 0.049, h = 0.673, A s = 2.215 × 10 −9 , n s = 0.96). The dashed cyan and green lines are the corresponding two-halo and one-halo terms, respectively, with the parametrisation given by [67]. The solid black line is the sum of the two terms, i.e. the non-linear power spectrum. One can clearly see that the two-halo term dominates at scales k 0.02 h Mpc −1 , where the power spectrum is very well described by linear theory. Moreover, the two-halo term approximates the linear power spectrum up to k ∼ 0.2 h Mpc −1 . On the contrary, for k 0.6 h Mpc −1 , the power spectrum strongly deviates from linearity and is very well approximated by the one-halo term. In the range of modes between these two limits it is necessary to consider both terms for an accurate description of the non-linear power spectrum. The time evolution of the power spectrum is shown in figure 2. In the left panel, I show the linear (orange lines) and non linear (blue lines) power spectra at redshifts z = 5, z = 3 and z = 0, represented with dot-dashed, dashed and solid lines respectively. The linear power spectrum has always the same shape, only the amplitude being affected by time evolution. This is a direct consequence of linear theory, according to which density perturbations are simply rescaled by the growth factor, δ(z) = (D(z)/D(z 0 ))δ(z 0 ). At z = 5, the non linear power spectrum is barely distinguishable from the linear one. This means that the scales in the dynamic range shown (k 20 h Mpc −1 ) are still growing linearly. As the redshift decreases, larger and larger scales enter the horizon, so one can see a departure of the non- linear power spectrum from the one predicted by linear theory. The point discriminating between the two regimes moves towards larger scales. The right panel of figure 2 shows the time evolution of the two-halo (cyan lines) and one-halo (red lines) terms. The line styles have the same meaning as in the left panel. The two-halo term evolves linearly at small k and becomes smaller and smaller at large k as time evolves. On the contrary, the time growth of the one-halo term is predominant at large k. Although the one-halo term does present an evolution at small k as well, it is always negligible with respect to the two-halo term.
In developing the generalisation of FKP, I shall assume that the power spectrum can be written as in equation (2.7). I shall then disentangle time-evolution effects, providing an estimator for the liner power spectrum at fixed time.

Computing Weighted Density Fluctuations
Following the FKP treatment, I define where A(k) is the normalisation function (2.13) I have explicitly included the possible k-dependence of the weight function. Taking the expectation value of |F (k)| 2 and exploiting the same relations used in the FKP treatment, I obtain The second term is the shot noise, which I shall write as P shot (k) from now on. The first term contains the two-point correlation function. One can assume that in a realistic situation the correlation function approaches zero very rapidly as the distance between r and r increases. Within this assumption and the distant observer approximation, one can write [41] ξ where r * = (r + r )/2. In other words, the past light-cone effect is encoded by the average distance between the two points considered. One can now plug this relation in the expression for |F (k)| 2 . It yields where T (k) and H(k) contain only the two-halo and one-halo terms, respectively. In the next two sections I shall treat each term separately. The goal is to find a strategy to obtain an estimator for the linear power spectrum at fixed time ∆ 2 L0 (k).

Two-halo Term
The two-term proposed by [67] is given by equation (2.8). Rather than developing an ad hoc method for the revised halofit model proposed by [67], I shall consider a wider class of models. Specifically, I assume the relation between ∆ 2 Q (k) and ∆ 2 L (k) to be I have not assumed that ∆ 2 Q is a function of the norm of k only, to keep the treatment as general as possible. I had to include a dependence of ∆ 2 Q on r to make sense of the time evolution of the perturbations. In other words, equation (2.19) says that the mapping can be written as a product of a function of ∆ 2 L (k) and k only, and a positive function of k and r only. The mapping expressed in equation (2.8) falls into this category.
If the two-halo term has the shape described in equation (2.19), it is possible to re-write T in (2.18) simply as a polynomial function of the linear power spectrum at fixed time. To make notation lighter, let us set x = ∆ 2 L 0 (k), so that the function f in (2.19) can now be written as f (D(r)/D(r 0 )) 2 x, k . Keeping r and k fixed, a different shape of the power spectrum will yield a different value of x. For any fixed k, one can then Taylor expand the function f in equation (2.19) around x = 0, till a certain order m. For x 1, the expansion recovers the linear regime. Forcing slightly the aforementioned approximation according to which the correlation function goes to zero rapidly as r − r increases, one can assume r ≈ r ≈ r * . Equation (2.16) then simplifies as If the functions G j (k, k ) are narrow in k-space, that is, peaked around k and such that their width is ∼ 1/L, where L is the depth of the survey, then one can define and write (2.23)

One-Halo Term
The one-halo term, as defined by [67], cannot be factorised as in equation (2.19). Moreover, it is not a function of ∆ 2 L 0 only, which makes it hard to disentangle the dependence on r and expand it in a Taylor series as it has been done for the two-halo term. The one-halo term depends on ∆ 2 L 0 only through k σ . Since the linear power spectrum depends on the time t(r) considered, the r-dependence of the one-halo term is also encapsulated in k σ . It is then necessary to introduce an approximation to drop the dependence on r. For this purpose, I propose the following ansatz. I replace ∆ 2 H (k, r) in equation (2.17) with its average over r, weighted with the number density of galaxies, that is .
In other words, ∆ 2 H (k) can be seen an "effective one-halo term". The one-halo term ∆ 2 H (k, r) contributes to ∆ 2 H (k) at any value of r. The contribution is bigger where the number density of galaxies is larger.
Although applied to a different context, this approach is similar to considering a spaceaveraged estimator of the power spectrum, as in [41] (discussed at the end of subsection 2.6). However, in this work only the one-halo term is averaged, and not the entire power spectrum.
Adopting the same approximations that in the previous paragraph led to equation (2.23), one obtains Following what has been done for the G j (k, k ) functions when computing T (k), one can assume that also L(k, k ) is narrow in k-space. Under such approximation:

Estimator and Optimal Weighting
Using equations (2.23) and (2.27), equation (2.16) can be re-written as (2.28) I define x to be the estimator of the dimensionless linear power spectrum at present time. It can be determined by solving the above equation, for every value of k of interest. The estimator can then be averaged over a shell of volume V k in Fourier space. I shall indicate the averaged estimator of the linear power spectrum at present time as P L 0 (k), and ∆ 2 L 0 (k) its dimensionless counterpart. Equation (2.28) is the main result of this work.
The variance of the estimator is given by σ 2 P L 0 (k) = ( P (k) − P (k)) 2 . To compute the variance, I shall follow the same reasoning and approximations adopted in FKP, using the identity ( P L 0 (k) − P L 0 (k))( P L 0 (k ) − P L 0 (k )) = | F (k)F * (k ) | 2 . (2.29) Considering k = k + δk, with |δk| k, one obtains The variance of the estimator, averaged in a spherical shell of volume V k in Fourier space, is then given by (2.34) -10 -In the limit α → 0, it can be shown that To obtain the optimal weight function w 0 (r), I require the above expression to be stationary upon arbitrary variations δw(r) of the weight function. The result is = P L 0 (k) 1 +n(r)P (k, r) .
(2.36) Expanding the denominator the weight function in the above expression only to first order, equation (2.6) is recovered. As a caveat, I recall that the identity (2.29) holds only if F (k) is Gaussian distributed. In the strongly non-linear regime the distribution of density fluctuations is not Gaussian, so in principle the above weight function may not be valid for any k. Nevertheless, the assumption of Gaussianity is valid at large scales [1], so the "real" optimal weight function should coincide with equation (2.36) at these scales. Therefore, even if (2.36) is not the optimal weight function, one expects it to be close to the actual one, by continuity. Furthermore, I will show in section 3 that the method is accurate up to the mildly non-linear regime (∆ 2 (k) 1), so in practice there is no reason to worry about the strongly non-linear regime (∆ 2 (k) 1). The expression of the weight function at the right hand side of the last equality in equation (2.36) is consistent with the findings by [41]. However, the results in [41] are only formally identical to the ones presented here. Indeed, the estimator defined by [41] is This quantity does not estimate the power spectrum at fixed time, but the power spectrum integrated over a certain range in r. Thus, it is different from P 0 (k). In [41] no particular functional shape describing the growth of perturbations is assumed. In addition to the light-cone, also other effects, like redshift space distortions and bias, are included in [41]. Consequently, the result in [41] is more general, but provides an estimator of the power spectrum integrated over a certain redshift range. On the other hand, the method presented in this manuscript is more model-dependent, but provides an estimator of the power spectrum at fixed time, disentangling the growth of perturbations. Both methods are correct, as long as one clearly bears in mind which estimator was defined, and is coherent with that definition.

Implementation of the Method
I shall now outline the steps of the code that has been implemented to test the method, underscoring the approximations adopted. First of all, I choose the cosmological parameters given by [2] (Ω m = 0.315, Ω Λ = 1 − Ω m = 0.683, Ω b = 0.049, h = 0.673, A s = 2.215 × 10 −9 , n s = 0.96) as the fiducial cosmological model. I then generate the corresponding linear power spectrum at redshift z = 0 using software CAMB [68,69]. I assume it to be the true power spectrum. Density fluctuations are assumed to be an isotropic random field, so the power spectrum depends only on the norm of k.
Then, the code computes the redshift-distance relation for the chosen cosmology and thus is able to obtain the growth function D(r) using where γ is a constant depending on the theory of gravitation considered. For General Relativity, its value is 0.55 [70]. The growth function is plotted in the left panel of figure 3.
Functions σ 2 (R) and k σ (r) are then computed, as well as all the quantities needed to define the one-and two-halo terms according to [67]. In this particular model, the function f (∆ 2 L (k), k) defined in section 2.4 actually depends only on the norm of k and only through ∆ 2 L (k). When one considers time evolution, there is of course an additional r dependence. Regarding the function g(k, r) defined in section 2.4, the r-dependence is encoded by the function k σ (r). Now one has to define the number density and the weight function. Regarding the number density, I assume a realistic shape given by [71] n(z) = Cz 2 exp − z z I chose z = 1 and C = 400. The corresponding shape of n(z) can be seen in the right panel of figure 3. The dependence of the number density on r is then easily found using the redshiftdistance relation previously computed. I then adopted the weight function minimising the variance.
One can compute |F (k)| 2 − P shot (k) using (2.28). Since all quantities depend now only on the norm of r and k, the expressions for G j (k, k ) and I j (k) are simplified accordingly.
Using all the expressions computed so far, equation (2.28) is solved for different values of k and at different orders of the expansion, restricting the solutions to real numbers. In this way, one can see if the equation admits real solutions and, if so, how many. In case more than one real solution is possible, one can try to identify the correct one through physical arguments. For example, I found out that one of the solutions of equation (2.28) was negative for all k. The power spectrum is by definition positive, so one could clearly identify the physical solution. One can use such criterion at any order. Even if there are two positive solutions for all k at a certain order m, one can still understand which is the correct one by considering that, at small k, it should not deviate appreciably from the solution at order m − 1. I can then investigate the validity of the method with further analysis, comparing the solutions at different orders with the power spectrum provided by CAMB.
Before showing the results, let us point out one last remark. In all the quantities considered, integrals over r and k are extended over the whole three-dimensional space. Of course, it is not possible to integrate till infinity, so cut-off values must be set. The output of CAMB provides P (k) till a certain wavenumber. The power spectrum tends to zero at large k. I then extrapolate P (k) to k max , defined so that P (k max ) = 0. In other words, instead of considering a power spectrum which tends asymptotically to zero as k tends to infinity, I assume that it is identically zero from a sufficiently large value k max till infinity. Similarly, also σ 2 (R) tends to zero at large R. I extrapolate again and define r max such that σ 2 (r max ) = 0. With the cut-off values considered, I am able to span a range in which the growth of perturbations is significant. For example, at z = 3, where the chosen number density of galaxies is not yet small, the growth factor is about 3.

Results
In figure 4 I show the input linear power spectrum at present (solid red line), together with the solutions of equation (2.28) at orders 1, 2, 3 and 5 (dashed black, green, purple and orange lines, respectively). 2 The solutions at all orders considered reproduce very well the input linear power spectrum in the linear regime, as expected. At first order, one can recover ∆ 2 L 0 up to k ∼ 0.1 h Mpc −1 . Increasing the order, the solutions approach the input linear power spectrum at higher k, being able to recover it even in the non linear regime. At order 5 the method seems to approximate ∆ 2 L 0 even up to k ∼ 1 h Mpc −1 . The goodness of the solution obtained at each order is quantified by computing the ratio between ∆ 2 L 0 and that solution. I then looked for the value of k up to which the ratio deviates from unity within ±0.01, ±0.05 and ±0.1. In other words, I sought the range of k in which the prescribed solution reproduces ∆ 2 L 0 within 1%, 5% and 10% accuracy, respectively. In the following, I shall denote the upper bound of such range as k lim . The best result is obtained at fifth order. The ratio between the input power and this solution is plotted in figure 5. At small k, this is very close to 1. As one moves towards smaller scales, the ratio fluctuates, eventually becoming larger and larger. The values of k lim corresponding to the aforementioned levels of accuracy are 0.025 h Mpc in this work is able to recover the input linear power spectrum with 5% in a range of modes where non linearities are already significant. The analysis explained above has been done for all orders considered. The results can be seen in figure 6. The plot shows k lim as a function of the order of the solution, for 1% (black crosses), 5% (black circles) and 10% (black squares) accuracy. Once again, the shaded grey area indicates the dynamic range where the linear regime is valid. The points at order 4 are missing because there is no real solution beyond k ∼ 0.44 h Mpc −1 . I verified that, beyond fifth order, the range of validity of the method for each level of accuracy does not increase. Therefore, the solution at order 5 gives the best estimator of ∆ 2 L 0 .

Discussion
The method presented in this work is able to reproduce the fixed-time linear power spectrum of matter density fluctuations in a wide dynamic range. Its range of validity has been determined by considering three levels of accuracy, namely 1%, 5% and 10%.
The most restrictive threshold has been chosen because the future Euclid survey is expected to estimate cosmological parameters (which influence the shape of the power spectrum) with a precision of a few percent [70]. However, considering possible errors arising from data 1% Accuracy 5% Accuracy 10% Accuracy analysis pipelines based on N-body simulations, 1% accuracy may be too optimistic. Indeed, simulations showed that baryonic physics may have a non-negligible impact on the matter power spectrum, so ignoring the role of baryons when analysing Euclid-like data could yield systematic errors in the inference of cosmological parameters. Different feedback models predict discrepancies up to tens of percent (see e.g. [72]), so even a modestly accurate method to infer the matter power spectrum from data could be able to discriminate among different feedback models. Even though at present observational data show no evidence for significant feedback effects on the matter power spectrum at k 2 h Mpc −1 [73], it may still be sensible to expect an accuracy worse than 1% for Euclid. Furthermore, the accuracy of the halofit model is only around 5% for k ∼ 1 h Mpc −1 at 0 < z < 3 [67], so it makes sense to consider a looser requirement of accuracy for the method presented in this paper.
For all the aforementioned reasons, a more sensible accuracy threshold to determine the range of validity of the method would be 5-10%. I recall that the method presented 1% Accuracy 5% Accuracy 10% Accuracy Figure 6. Upper limit of the range in which the estimated linear power spectrum recovers the assumed true linear power spectrum within 1%, 5% and 10% accuracy (black crosses, circles and squares, respectively), as a funciton of the order of the Taylor expansion used to compute the estimator. The results plotted refer to redshift z = 0. The grey shaded area represents the region where the true linear power spectrum at z = 0 is smaller than 1, i.e. where linear growth is a legitimate approximation. The values of k lim are positive by definition. The y-axis starts from negative values only to make the plot more readable.
here reproduces the input power spectrum up to k ∼ 0.80 hMpc −1 , already well into the non linear regime, with an accuracy comparable to, or better than, the underlying halofit model. Relaxing the demands on the accuracy down to 10%, one can recover the input power spectrum up to almost 1 h Mpc −1 . These are non trivial results, especially if one considers that these performances could make the method applicable to future large-scale surveys like Euclid. The method allows estimating the linear power spectrum at a fixed time. The estimator is not integrated over a redshift range, so the time evolution of density perturbations has been disentangled. The attempt by [41] of dealing with the estimation of the power spectrum in the presence of the light-cone effect is more general, as it does not assume one particular model for the growth of density perturbations. Also, further effects, like redshift space distortions and bias, are considered. On the other hand, it does not completely disentangle time evolution. Thus, in that respect the method presented here is more effective. The two methods are then somewhat complementary to each other. The optimal weight functions obtained in the two cases are consistent.
When implemented, the method presented here is computationally very cheap. The estimation of ∆ 2 L 0 is reduced to the calculation of a finite number of integrals and the solution of a polynomial equation. Given sufficient computing time, it is always possible to find the solution. For the levels of accuracy tested, the method achieves its widest range of validity already at order 5. I also tested the procedure with a less realistic trend for the number density of galaxies (i.e. a constant) and still the results are consistent with the ones presented in this manuscript.
The method has been validated for the halo model, so it is model-dependent. This issue could be circumvented using forward modeling and statistical inference. One could assume cosmological parameters and a value for γ in equation (3.1). Consequently, one would have the shape of the growth function and the linear power spectrum at present time, obtaining a prediction for |F (k)| 2 −P shot (k). In a real survey, this quantity is obtained from data, so one can get the the likelihood and then, exploiting Bayes theorem, the probability distribution of the parameters, given the data.
There is still some merit in finding an estimator of P (k) with the strategy outlined in this work, though. An analogy with the CMB may be helpful to go through this point. In the case of the CMB, an enormous amount of data is collected (e.g. [2]). The information given by the temperature map is compressed into the angular power spectrum, in a model independent way. 3 This statistic actually encodes the information which is necessary to infer the probability of a certain model, given the data. It is not obvious that a similar compression of the data could be done in the case of the power spectrum of matter density fluctuations as well. The method described in this work represents an attempt of doing that. Not surprisingly, it turned out that I had to assume both the growth function and a model describing non linearities. Hopefully, in future work one could be able to relax the assumptions about the latter, making the technique less model-dependent. Regarding the growth function, the situation is different, since it does depend on the cosmological model. One could first of all check how sensitive the method is to a change of the functional form of the growth function. If it turns out that it is dependent on just one parameter, e.g. the derivative of the growth function at some intermediate redshift, one could consider a family of power spectra, each one corresponding to a different value of such parameter, thus providing an estimator for each one of them. So, further insight is needed also in this direction. Another perspective is to generalise the implementation of the method dropping the various approximations described in section 3.

Conclusions and Perspectives
I extended the method by [1], incorporating the light-cone effect. I defined a minimumvariance estimator for the linear power spectrum of matter density fluctuations at fixed time and provided and analytic expression for the weight function. This has been done considering non-linear growth of density perturbations, within the context of the halofit model.
The weight function minimising the variance of the estimator has the same shape as in the original FKP paper, except for the additional dependence of the power spectrum on distance, as shown in equation (2.36). The result presented in this work is formally identical to [41], but has a different physical meaning. Nevertheless, the two methods are consistent. With respect to Yamamoto's method, the technique presented in this work completely disentangles the time evolution of density perturbations, providing an estimator of the linear power spectrum at fixed time. On the other hand, [41] takes into account more effects, such as redshift space distortions and bias.
The estimator of the fixed-time linear power spectrum defined in this work is one of the solutions of (2.28), which is the main result of the paper. Equation (2.28) relies on a Taylor expansion, meaning that one obtains a different solution depending on the order considered. It is noteworthy that the method reduces to the computation of a finite number of integrals and the solution of a polynomial equation. Given sufficient computing time, a solution can always be found.
The method has been tested for the ΛCDM cosmology with the parameters given by [2], assuming the revised halofit model by [67] and a shape for number density of galaxies as in [71]. At order 5 of the Taylor expansion in equation (2.28), the estimator of the linear power spectrum at present time reproduces the input power spectrum with 1%, 5% and 10% accuracy up to k ∼ 0.025 h Mpc −1 , k ∼ 0.80 h Mpc −1 and k ∼ 0.94 h Mpc −1 , respectively. The range of validity does not increase at higher orders.
It is remarkable that the method is capable of deriving the fixed-time linear power spectrum within 5% − 10% well above k ∼ 0.20 h Mpc −1 , where linear theory breaks down [74]. This means that the method it potentially useful for analysing the data of future surveys, like Euclid. Indeed, although the expected accuracy of Euclid is 1% [70], ignoring baryonic physics in the data analysis may introduce systematic errors, degrading the accuracy. For example, different feedback models could introduce discrepancies in the predicted matter power spectrum up to the order of tens of percent (e.g. [72], but see also [73]).
Despite its simplicity and accuracy over a significant dynamic range, the method presented in this work is model dependent. It assumes that the growth of density perturbations is set by General Relativity and that non linearities in the power spectrum are described by the halo model. Of course, it would be desirable to relax the dependence of the method on any of these assumptions. The former could be addressed by obtaining a family of estimators for different values of the parameter γ in (3.1) and compare each one of them with data, to discern which value gives the best fit. Dealing with the latter hypothesis would probably require modeling non linearities with less parameters.
As an alternative method, one could consider a forward modelling approach. One could parametrise the growth function and the linear power spectrum at present time to predict the left hand side of (2.28). This would be estimated from data of surveys, so that one could apply Bayes theorem to obtain the probability distribution of the parameters, given the data. One could then investigate to which extent the data could be compressed in a model-independent way, although it is not obvious that it can be done.
There are also other ways to improve the method presented in this work. First of all, one could include other effects, like redshift space distortions and bias. Moreover, the technique has been implemented exploiting the fact that the correlation function drops rapidly to zero as the distance between points increases. Nevertheless, [18] suggest that this approximation does not hold well at large scales. All the aforementioned assumptions merit further study.