Understanding the Reconstruction of the Biased Tracer
Abstract
Recent development in the reconstruction of the largescale structure (LSS) has seen significant improvement in restoring the linear baryonic acoustic oscillation (BAO) from at least the nonlinear matter field. This outstanding performance is achieved by iteratively solving the MongeAmpere equation of the mass conservation. However, this technique also relies on several assumptions that are not valid in reality, namely the longitudinal displacement, the absence of shellcrossing and homogeneous initial condition. In particular, the conservation equation of the tracers comprises the biasing information that breaks down the last assumption. Consequently, direct reconstruction would entangle the nonlinear displacement with complicated bias parameters and further affect the BAO. In this paper, we formulate a theoretical model describing the reconstructed biased map by matching the tracer overdensity with an auxiliary fluid with vanishing initial perturbation. Regarding the performance of the reconstruction algorithm, we show that even though the shot noise is still the most significant limiting factors in a realistic survey, inappropriate treatment of the bias could also shift the reconstructed frame and therefore broaden the BAO peak. We suggest that, in principle, this biasrelated BAO smearing effect could be used to independently selfcalibrate the bias parameters.
Subject headings:
cosmology: largescale structure of universe1. Introduction
The extraction of cosmological information from the abundant galaxies survey data has long been one of the primary objectives of the LSS studies. At the large scale where the density perturbations are still linear, this process is quite straightforward from measuring the twopoint function of the galaxies distribution. However, as the nonlinearity of the structure formation gradually takes effect at smaller scales, the information starts to leak into higherorder statistics (Rimes and Hamilton, 2005, 2006). One wellknown example is the baryonic acoustic oscillation (BAO) peak which has been used as a standard ruler to measure the expansion history of our Universe (Eisenstein, 2003; Blake and Glazebrook, 2003; Hu and Haiman, 2003; Seo and Eisenstein, 2003; Eisenstein et al., 2005). This sharp feature was smeared out when particles/galaxies started to move out of their original location (Crocce and Scoccimarro, 2008) and eventually causing us to lose the constraining power on dark energy (Seo and Eisenstein, 2007; Ngan et al., 2012).
The question is whether we can undo the structure formation and recover the primordial information from the data. Two different strategies emerged to tackle this problem, one is the forward modeling, and the other one is the backward reconstruction. Both approaches have made significant progress recently. The first approach samples the vast parameter space of initial fluctuations and compares its forward evolution, using simplified Lagrangian model, FastPM simulation, or the full Nbody nonlinear dynamics, against observations. It usually involves finding the maximum a posteriori solution using Gibbs or/and Hamiltonian Markov chain Monte Carlo sampling (Jasche and Wandelt, 2013; Wang et al., 2009, 2013, 2014; Feng et al., 2018), which are usually computationally expensive.
The backward reconstruction, on the other hand, directly operates on the observed data, expecting to recover the initial fluctuation. Earlier examples include the logarithmic transformation and Gaussianization(Weinberg, 1992; Neyrinck et al., 2009, 2011; Wang et al., 2011), both of which attempt to Gaussianize the onepoint probability function. However, as local transformations, they are incompatible with the process of the structure formation, therefore do not genuinely reproduce the initial condition. Eisenstein et al. (2007) focused on the nonlinear BAO degradation instead and demonstrated a simple yet powerful method to reduce the smearing of the signal. Recently, various improved algorithms, including the isobaric reconstruction (Zhu et al., 2017; Yu et al., 2017; Wang et al., 2017; Pan et al., 2017; Zhu et al., 2018) and other iterative solutions (Schmittfull et al., 2017; Shi et al., 2018; Hada and Eisenstein, 2018), have been shown to be able to almost entirely remove the degradation, and perfectly recover the linear BAO signal. Despite their divergent technical details, all these algorithms solve the displacement potential from the following mass conservation equation
(1) 
Here and are Lagrangian and Eulerian coordinates of particles, and the displacement vector and potential are defined as , where is the Kronecker delta. By solving equation (1), the reconstruction algorithm eliminates this nonlinear coordinates transformation and produces the field on a grid that is close to the Lagrangian system, namely the isobaric frame.
Despite their excellent performance, these algorithms have at least made the following three assumptions

the displacement field is longitudinal,

there is no shellcrossing,

the initial perturbation is negligible, i.e. .
The first assumption guarantees that the problem is solvable, i.e., one unknown solved from one equation; assumption (2) necessary for the uniqueness since otherwise, any particle permutation would produce a new solution; finally the uniform initial distribution makes sure no extra information would be needed.
Unsurprisingly, these assumptions introduce complications and errors to the real application of the algorithm. For example, the neglected transverse component, which appears until the third order in the Lagrangian perturbation theory, add extra contributions to the reconstructed potential . Moreover, the reconstruction result is most likely meaningless below the shellcrossing scale. However, these higherorder smallscale errors have limited effects on the cosmological constraint at the BAO scale. The third assumption breaks down for any biased samples.
Of course, one could proceed and perform the reconstruction regardlessly. As demonstrated by Yu et al. (2017), the performance, characterized by the crosscorrelation with the initial condition, is notably lower than that of the matter field. While the shot noise undoubtedly plays a significant role, it is unclear how the clustering bias would affect the reconstruction. In practice, we would like to understand this question because it might lead to different planning strategies for the future surveys.
Theoretically, we can study a noiseless biasing model, i.e., a (deterministic) mapping from matter overdensity to the protohalo field . For given bias function or its Taylor series, we would like to understand in details the properties of the reconstructed field. As first observed by Yu et al. (2017), the reconstructed halo field has the same largescale bias as the halo overdensity itself. It is unsurprising since the reconstruction algorithm does not differentiate an input dark matter map from the biased tracer. However, the question remains for smallscale, nonlinear and scaledependent biases. How can we extend this simple intuition to describe more complicated biasing model? Moreover, is there any practically feasible way of correcting their effects?
The reason why the reconstruction could restore the linear BAO lies in its ability to recover the Lagrangian frame in which the BAO features were defined. Because of those algorithmic assumptions and numerical errors, the isobaric frame deviates from even for matter reconstruction, even though the difference is quite small. With the clustering bias, however, it is not difficult to imagine a much larger frame shift and possibly some extra smearing of the BAO.
In this paper, we would like to address these questions. In section 2, we investigate the consequences of the direct reconstruction without any preprocessing of the bias. We formulate a theoretical model describing the reconstructed biased map by matching the tracer overdensity with an auxiliary fluid with vanishing initial perturbation. We focus on the frameshift and BAO smearing in section 3 and then discuss the possible bias selfcalibration in section 4. Finally, we conclude in section 5.
We sometimes denote the argument of a function as subscripts, e.g. , where the comma are used as separators. The derivative, on the other hand, uses semicolon, .
2. Theoretical Modeling of the Isobaric Reconstruction
As demonstrated by Wang et al. (2017), at large scale the isobaric reconstruction produces the longitudinal component of the nonlinear displacement field, and the equation (1) is only solvable assuming the homogeneity of the initial density distribution. In this section, we discuss the additional effect introduced by the clustering bias and the consequences of the direct reconstruction without any preprocessing of the map.
By definition, the reconstruction of a general map could be described by equation (1), with the substitution that . Conceptually, this is equivalent to introducing an extra auxiliary fluid with a uniform initial distribution whose displacement is described by no matter how was formed at the first place.
2.1. Revised Newtonian Dynamics
Assuming our auxiliary fluid also follows the Newtonian dynamics, the comoving Eulerian position of a fluid element follows the trajectory accelerated by the gravitational force
(2) 
where denotes the gradient in Eulerian space, is the Lagrangian derivative with respect to the conformal time , , and is the Newton gravitational potential which obeys the Poisson equation
(3) 
For the biased reconstruction, the gravitational potential here is instead determined by the galaxies overdensity, i.e. , where is the number density. By definition, these auxiliary fluid element starts from the isobaric frame, i.e. , and could be described by along its movement, where is our reconstructed displacement potential. Since we only need the longitudinal part of equation (2), combining with equation (3), one simply has
(4) 
Here we have defined the nonlinear operator . We could further rewrite the above equation with the Jacobian matrix between and the isobaric frame ,
(5) 
where is the Kronecker delta. Following Matsubara (2015), one could expand above equation and have
(6) 
where is the LeviCivita symbol, is the determinant of the Jacobian matrix . We have also assumed the galaxy overdensity is some nonlinear/nonlocal function of the matter density , and here is the the Jacobian from Lagrangian to Eulerian coordinate . Let us now only keep the linear bias, denoted as , we have
Remember that both sides of the equation are evaluated at even though and are explicitly defined in and respectively, with the assumption that neither nor is singular. Substituting the definition of , we have
(8) 
To proceed, we need the dynamical equation of as well. At the linear order, however, the homogeneous part of equation (2.1) is simply
(9) 
Here is just the linear density perturbation, which is the solution of the equation . Hence, it is straightforward to see that the has the linear solution . Notice that one property of this approach is that the specific growth function (as a function of ) is not relevant as long as it gives correct values at the boundary, i.e. at the initial and final time.
2.2. Density Matching in Configuration Space
One advantage of the above dynamical approach is that it provides a complete framework to derive the recurrence relation of the revised displacement in all order, e.g., following Matsubara (2015). However, since the exact growth history is irrelevant, such derivation seems unnecessarily complicated. A much more straightforward and more intuitive approach is also available. Since all the reconstruction algorithm does was to find a longitudinal displacement field that could reproduce density map as accurate as possible, one could then define by matching the density field
(10)  
For a given model of the biased tracer , we could then obtain by solving this equation. Here, is the nth order bias parameter and is the reconstruction error, i.e., the density residual. Notice that the density on the righthand side of equation (10) could depend on many other ingredients as well, including the transverse component of the displacement field , which only contribute after the second loop order, therefore, we will neglect in this paper.
We do not have to provide any mathematical proof of the existence and uniqueness of above equation since the density residual is unknown and will never be zero in practice. Even with the assumption that , as we will do in the following of the paper, one could show that this is always solvable at least perturbatively. For example, in the linear Eulerian bias model, equation (10) then reduces to
(11) 
Therefore, the reconstructed field is also biased by , which agrees with our previous intuition and the conclusion of Yu et al. (2017).
Since the Eulerian bias model is not very convenient in describing the movement of galaxies^{1}^{1}1One then has to add scaledependent biases, we will work within the Lagrangian approach in the rest of this section. To proceed, we now consider the biased clustering model that tracers also obeys the same following conservation equation as matters
(12) 
where is the velocity divergence. We have assumed that the tracers comove with the dark matter, i.e., , which is right at large scale because of the equivalence principle. Combining these two equations, we have the following equation
(13) 
with a straightforward solution (Desjacques et al., 2016)
(14) 
Here is the formation time of these tracers, and is their moving Eulerian coordinates. Given the as a function of the real matter displacement potential and combining equations (14), (10) and (1), we further have a general expression of the reconstructed potential
(15) 
Notice that the potential field and are defined in two different coordinates in general, i.e. and respectively. Again, assuming the nonsingular mapping and , this equation is actually evaluated in the Eulerian space , and . Now expanding this equation perturbatively and neglecting all nonlinear terms, we then obtain a Poisson equation for the reconstructed field
(16)  
Here, and are first order perturbations of the tracer and matter at . For simplicity, let us assume that the number density of the tracer was linearly biased with respect to by at the formation time (Desjacques et al., 2016)
(17) 
This then leads to the linear equation of the reconstructed displacement potential
(18)  
Here we have once again used the fact that at the linear order . Following the usual convention, we have also taken the limit while keep the Lagrangian bias fixed. Since the potential field is already a first order quantity, we are safe to drop the coordinates difference of the Laplacian, i.e., , so that the linear solution of the reconstructed field is simply
(19) 
Since the linear Eulerian bias , this agrees with equation (10).
A more important aspect of deriving equation (19) in the Lagrangian framework is that the reconstructed isobaric frame will deviate from the primordial coordinate by , which at the linear order equals
(20)  
It then gives that . As will be shown in the next section, this coordinate shift could have a major consequence on the smearing of the baryonic acoustic oscillation signals.
2.3. Higher Orders
One could extend above calculation to higher orders. Instead of configuration space, we present in appendix A a systematic derivation in the Fourier space. Similar to the usual Lagrangian perturbation expansion, one could define the reconstructed displacement vector
(21)  
with the modified LPT kernels
(22)  
where are the original LPT kernels in equation (A). Once again, the first order agrees with equation (19) although here the linear Lagrangian bias could also be scaledependent. At the second order, however, the kernel not only depends on and , but also lower order and as well. This statement is true in general that depends on all lower order () LPT kernels and bias parameters.
3. The Frame Shift and BAO Damping
As mentioned previously, the reconstruction of the biased tracers inevitably induces a coordinate shift from the Lagrangian frame to isobaric frame . In this section, we will particularly focus on its effect on the baryonic acoustic oscillation.
To proceed, we recall that the BAO features are defined in the primordial Lagrangian coordinates, and one consequence of the nonlinear gravitational evolution of the matter field is the BAO smearing caused by the pairwise displacement of particles (Crocce and Scoccimarro, 2006a, b, 2008) , where the nonlinear power spectrum of the density fluctuation could, in general, be expressed as a combination of two contributions
(23)  
namely the propagator part which is proportional to and the socalled modecoupling terms . Here , the nonlinear propagator , where is the variance of the displacement (Crocce and Scoccimarro, 2006a). The BAO signal, defined as ratio between and the nowiggle power spectrum ,
(24) 
will be smeared out by the variance of this displacement .
Similarly, the coordinate shift induces an extra smearing on the baryonic acoustic oscillation. To see this, we could express our reconstructed displacement potential in Fourier space as
(25)  
assuming the existence of the onetoone mapping . At the leading order, the power spectrum of the reconstructed field could be approximated as
(26)  
where is the linear power spectrum and we have neglected all higher order contributions. Consequently, without correcting the bias, one receives an extra BAO smearing proportional to , i.e.
(27) 
We demonstrate this effect in figure (1). To eliminate the complications from the shot noise, we perform the reconstruction on the matter field. In this idealized example, we simply consider the matter field as a linearly biased tracer with the . Therefore, any incorrect estimation of the bias, i.e., reconstructing where , would introduce this BAO damping effect after the reconstruction. Assuming the bias estimator differs from the true value by , i.e. , the coordinate shift
(28)  
so we have when .
In Figure (1), we demonstrate the reconstructed BAO after rescaling the density map with corresponding equivalent shift parameter and respectively. The BAO wiggles are derived from two otherwise identical simulations, particularly the random seed, with different transfer functions, i.e., with vs. without BAO. The simulations shown in this figure have box with particles. As shown in the upper panel of the figure, the larger the is (solid lines), the more suppressed the BAO wiggles appear to be after the reconstruction. In the same plot, we also illustrate the simple analytic solution, i.e., equation (27), as dashed lines. For smaller shifting parameters or , this formula describes these curves reasonably well, but it starts to deviate for larger . The error of the model is shown in the lower panel of the figure.
We then perform the reconstruction on the real halo samples. Due to the shot noise contamination, the BAO wiggles, of both halo field and its reconstruction, are averaged over realizations of nbody simulations with box and particles. The catalog comprises the largest halos with corresponding comoving number density and respectively.
The result is presented in Figure (2). As shown in each panel, the red dashed line is the direct reconstruction from halo number density field; the solid black line corresponds to the linearly ’debiased’, i.e., , reconstruction; and the green dashed line is the halo reconstruction multiplied by the inverse of the exponential damping model equation (27). Meanwhile, the BAO of the original halo field and the matter reconstructions are shown in blue and red dotted lines respectively.
From the figure, we can see that the reconstruction produces more and clearer BAO wiggles compared to the halo map itself in every tested sample, regardless of any biasrelated preprocessing details. The shot noise is still the most significant limiting factor regarding the performance of the reconstruction. A halo sample with number density will provide enough information for the algorithm to recover almost all BAO signals up to , which is enough cosmologically since the Fisher information on the sound horizon scale saturates after (Wang et al., 2017). Moreover, the clustering bias of such sample is usually very close to one, which further help to simplify the analysis. For example, as shown in the top panel, the BAO of the reconstructed and are almost identical.
The frameshift effect is noticeable for highbiased samples. In the middle panel of Figure (2), where and linear bias , one could see the improvement of the linearly ‘debiased’ sample (blacksolid) compared to the direct reconstruction (reddashed). At least for the first four peaks (), this improvement is consistent with our analytic solution (greendashed). For the sample with lower number density, however, such simple estimation starts to fail, likely caused by the combination of shot noise and larger bias deviation .
4. Selfcalibration of the Clustering Bias
4.1. Linear Bias Calibration
Because of the biasrelated coordinate shift and wiggle smearing, in principle, one could use the sharpness of BAO peak to calibrate or even constrain the bias parameters. In Figure (3), we replot the BAO wiggles of those two high bias samples, demonstrating this idea. Since for both of the samples, one could consider the direct reconstruction (red dashed) as an underestimated bias reconstruction. We also include the wiggles from the reconstruction of an overestimated bias (yellow solid). Compared to the correct (black solid), both of them suffer extra damping, as expected. From equation (28), for , whereas the direct reconstruction has a shift parameter , which equals and for two samples respectively. So, both situations produce a similar amount of damping, which is indeed the case from the figure.
In practice, one could constrain the linear bias by repeatly adjusting this parameter before the reconstruction. The best fit bias would be the one producing the least amount of the BAO damping. In the perspective of the configuration space correlation function, we mainly use the width information of the BAO peak which was not utilized before. This method is independent of other types of the bias measurement.
It is then interesting to see how stringent the constraint could be. For this purpose, we calculated a twoparameter Fisher matrix estimation. The Fisher matrix is expressed as (Seo and Eisenstein, 2003)
(29) 
where the effective volume
(30) 
where is the comoving density of the sample, is the survey volume, which we assume . Here we are only interested in two parameters: the sound horizon scale and the bias deviation . Particularly, we choose to use for its nontrivial derivative ^{2}^{2}2The same derivative would be zero for at the fiducial value .
(31) 
The twodimensional constraint is shown in Figure (4). We only display the upper half of the counter as . The comoving number density assumed here are , , and respectively. As mentioned previously, the constraining power on the bias deviation originates from the sharpness of the BAO peak whereas the sound horizon scale from the peak location. Hence there is essentially no degeneracy among these two parameters. For a reasonably sampled survey, e.g. , the 1 constraint is for and for . They then reduces to about half for a perfect shotnoiseless survey to and respectively.
4.2. Scaledependent and Nonlinear Biases
So far we have only discussed the selfcalibration and constraints on the scaleindependent linear bias. For the recovery of the linear BAO, one is further encouraged to construct a matter density estimator as accurate as possible. For example, one could apply a dependent Wiener filter so that
(32) 
with . In this paper, however, we did not proceed along this direction due to technical limitations of our reconstruction solver at this moment. In practice, this filter is clearly sample dependent and need to be carefully parametrized. The caveat is then a subtle balance between the model accuracy and the number of parameters needed.
It is also interesting to see what the higherorder corrections are if we only partially remove, say, the linear bias. To proceed, let us divide the map by some scaledependent bias parameter before the reconstruction. From the derivation in the Appendix, this corresponds to defining the modified LPT kernel such that
(33) 
where is the Eulerian kernel constructed from corresponding LPT kernel . It is straightforward to calculate the reconstructed kernels as
(34)  
When and assuming all biases are local, we could then remove the bias effect at linear order
(35) 
but the second order kernel
(36)  
is still quite complicated, and the biasrelated corrections do not vanish in general. Therefore, unless we know in advance the nonlinear and nonlocal bias parameters accurate enough, any limited corrections would inevitably cause some complex residuals at higher orders. We defer a more detailed study in the future.
5. Discussion and Conclusion
In this paper, we attempt to achieve a better understanding of the isobaric reconstruction of the biased tracer. Compared to the matter field, the performance is largely limited by the shot noise and clustering bias. Still, for all samples with various number density we have tested, the reconstruction always improves the BAO signal regardlessly. Particularly, for a reasonable spectroscopic survey with , the reconstructed BAO will be washed out by the shot noise at . To some extent, that will only partially affect the dark energy constraint since the information gain saturates around (Wang et al., 2017). On the other hand, the postreconstruction BAO degradation at low k is likely to be caused mainly by the clustering bias.
We have demonstrated that even with a simple linear debiasing, one is already able to sharpen the BAO signal. The improvement is consistent with simple Gaussian damping model. It is possible that a more sophisticated nonlinear debiasing scheme might improve the reconstruction furthermore. Alternatively, assuming all matter exists in dark matter halos, and assuming we have a reasonably accurate estimation of the cluster mass, one could in principle construct the massweighted density field as a proxy to the underlying matter distribution.
Sine any incorrect bias estimator would further smear the BAO wiggles, we proposed a selfcalibration scheme to constraint the linear bias. A simple twodimensional Fisher prediction showed that the constraints on the bias and sound horizon scale are orthogonal with each other, which is because they independently use the sharpness and the scale of the BAO peak respectively.
Acknowledgements: XW would like to thank productive discussion with Marcel Schmittfull and Matias Zaldarriaga.
Appendix A Nonlinear Reconstruction Kernels
In this appendix, we will apply the density matching technique to derive the higher order LPT kernels of the reconstructed field. We start by considering the following conservation equation (Matsubara, 2011)
(A1) 
Here is the displacement vector and we have introduced the Lagrangian nonlinear bias function with the property that , so it has the perturbative expression , the Lagrangian bias here could also be nonlocal, so generally we have as a function of in Fourier space. Transforming equation (A1) to Fourier space and expand the bias function perturbatively, one obtains (Matsubara, 2011)
where the nonlinear displacement field has its own perturbative expansion
(A3)  
where is the th LPT kernel. The first two are
(A4) 
where .
On the other hand, one could expand directly in Eulerian space
(A5)  
Here we have neglected some intermediate steps that expresses the as a function of nonlinear matter density and condense all processes into the compact kernel . Further combining equation (A) and (A5) produces a perturbative relation between these two types of kernels (Matsubara, 2011)
(A6)  
When there is no bias, these kernels then reduces to the Eulerian matter perturbation kernels, which usually denoted as , with a similar relation (Matsubara, 2011)
(A7) 
Since the algorithm assumes the input map to be matter field whose initial distribution is homogeneous, the reconstruction of the biased map is equivalent of finding the effective LPT kernels such that the corresponding satisfy the following relation
(A8) 
where is a function of as shown in equation (A). To the second order, it is straightforward to derive the effective LPT kernels
(A9)  
The relation becomes tedious at higher and we stop at the second order. Finally, the reconstructed displacement field of biased tracer could be described by
(A10)  
References
 Rimes and Hamilton (2005) C. D. Rimes and A. J. S. Hamilton, MNRAS 360, L82 (2005), astroph/0502081 .
 Rimes and Hamilton (2006) C. D. Rimes and A. J. S. Hamilton, MNRAS 371, 1205 (2006), astroph/0511418 .
 Eisenstein (2003) D. Eisenstein, ArXiv Astrophysics eprints (2003), astroph/0301623 .
 Blake and Glazebrook (2003) C. Blake and K. Glazebrook, ApJ 594, 665 (2003), astroph/0301632 .
 Hu and Haiman (2003) W. Hu and Z. Haiman, Phys. Rev. D 68, 063004 (2003), astroph/0306053 .
 Seo and Eisenstein (2003) H.J. Seo and D. J. Eisenstein, ApJ 598, 720 (2003), astroph/0307460 .
 Eisenstein et al. (2005) D. J. Eisenstein, I. Zehavi, and e. a. Hogg, ApJ 633, 560 (2005), astroph/0501171 .
 Crocce and Scoccimarro (2008) M. Crocce and R. Scoccimarro, Phys. Rev. D 77, 023533 (2008), arXiv:0704.2783 .
 Seo and Eisenstein (2007) H.J. Seo and D. J. Eisenstein, ApJ 665, 14 (2007), astroph/0701079 .
 Ngan et al. (2012) W. Ngan, J. HarnoisDéraps, U.L. Pen, P. McDonald, and I. MacDonald, MNRAS 419, 2949 (2012).
 Jasche and Wandelt (2013) J. Jasche and B. D. Wandelt, MNRAS 432, 894 (2013), arXiv:1203.3639 .
 Wang et al. (2009) H. Wang, H. J. Mo, Y. P. Jing, Y. Guo, F. C. van den Bosch, and X. Yang, MNRAS 394, 398 (2009), arXiv:0803.1213 .
 Wang et al. (2013) H. Wang, H. J. Mo, X. Yang, and F. C. van den Bosch, ApJ 772, 63 (2013), arXiv:1301.1348 .
 Wang et al. (2014) H. Wang, H. J. Mo, X. Yang, Y. P. Jing, and W. P. Lin, ApJ 794, 94 (2014), arXiv:1407.3451 .
 Feng et al. (2018) Y. Feng, U. Seljak, and M. Zaldarriaga, ArXiv eprints (2018), arXiv:1804.09687 .
 Weinberg (1992) D. H. Weinberg, MNRAS 254, 315 (1992).
 Neyrinck et al. (2009) M. C. Neyrinck, I. Szapudi, and A. S. Szalay, ApJ 698, L90 (2009), arXiv:0903.4693 [astroph.CO] .
 Neyrinck et al. (2011) M. C. Neyrinck, I. Szapudi, and A. S. Szalay, ApJ 731, 116 (2011), arXiv:1009.5680 [astroph.CO] .
 Wang et al. (2011) X. Wang, M. Neyrinck, I. Szapudi, A. Szalay, X. Chen, J. Lesgourgues, A. Riotto, and M. Sloth, ApJ 735, 32 (2011), arXiv:1103.2166 .
 Eisenstein et al. (2007) D. J. Eisenstein, H.J. Seo, E. Sirko, and D. N. Spergel, ApJ 664, 675 (2007), astroph/0604362 .
 Zhu et al. (2017) H.M. Zhu, Y. Yu, U.L. Pen, X. Chen, and H.R. Yu, Phys. Rev. D 96, 123502 (2017), arXiv:1611.09638 .
 Yu et al. (2017) Y. Yu, H.M. Zhu, and U.L. Pen, ApJ 847, 110 (2017), arXiv:1703.08301 .
 Wang et al. (2017) X. Wang, H.R. Yu, H.M. Zhu, Y. Yu, Q. Pan, and U.L. Pen, ApJ 841, L29 (2017), arXiv:1703.09742 .
 Pan et al. (2017) Q. Pan, U.L. Pen, D. Inman, and H.R. Yu, MNRAS 469, 1968 (2017), arXiv:1611.10013 .
 Zhu et al. (2018) H.M. Zhu, Y. Yu, and U.L. Pen, Phys. Rev. D 97, 043502 (2018), arXiv:1711.03218 .
 Schmittfull et al. (2017) M. Schmittfull, T. Baldauf, and M. Zaldarriaga, Phys. Rev. D 96, 023505 (2017), arXiv:1704.06634 .
 Shi et al. (2018) Y. Shi, M. Cautun, and B. Li, Phys. Rev. D 97, 023505 (2018), arXiv:1709.06350 .
 Hada and Eisenstein (2018) R. Hada and D. J. Eisenstein, MNRAS 478, 1866 (2018).
 Matsubara (2015) T. Matsubara, Phys. Rev. D 92, 023534 (2015), arXiv:1505.01481 .
 Desjacques et al. (2016) V. Desjacques, D. Jeong, and F. Schmidt, ArXiv eprints (2016), arXiv:1611.09787 .
 Crocce and Scoccimarro (2006a) M. Crocce and R. Scoccimarro, Phys. Rev. D 73, 063519 (2006a), astroph/0509418 .
 Crocce and Scoccimarro (2006b) M. Crocce and R. Scoccimarro, Phys. Rev. D 73, 063520 (2006b), astroph/0509419 .
 Matsubara (2011) T. Matsubara, Phys. Rev. D 83, 083518 (2011), arXiv:1102.4619 [astroph.CO] .