Post-collapse perturbation theory in 1D cosmology -- beyond shell-crossing
Atsushi Taruya, St\'ephane Colombi
TL;DR
This paper introduces a new perturbation theory approach for modeling the gravitational dynamics of large-scale structures in 1D cosmology after shell-crossing, improving accuracy over previous models.
Contribution
The authors develop a post-collapse perturbation theory with adaptive smoothing to accurately describe multi-stream flows after shell-crossing in 1D cosmology.
Findings
Accurately reproduces power spectrum at small scales
Effectively models phase-space structures post-shell-crossing
Outperforms Zel'dovich approximation in simulations
Abstract
We develop a new perturbation theory (PT) treatment that can describe gravitational dynamics of large-scale structure after shell-crossing in the one-dimensional cosmological case. Starting with cold initial conditions, the motion of matter distribution follows at early stages the single-stream regime, which can, in one dimension, be described exactly by the first-order Lagrangian perturbation, i.e. the Zel'dovich solution. However, the single-stream flow no longer holds after shell-crossing and a proper account of the multi-stream flow is essential for post-collapse dynamics. In this paper, extending previous work by Colombi (2015, MNRAS 446, 2902), we present a perturbative description for the multi-stream flow after shell-crossing in a cosmological setup. In addition, we introduce an adaptive smoothing scheme to deal with the bulk properties of phase-space structures. The filtering…
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Post-collapse perturbation theory in 1D cosmology – beyond shell-crossing
Atsushi Taruya1,2 and Stéphane Colombi1,3
1Center for Gravitational Physics, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan
2Kavli Institute for the Physics and Mathematics of the Universe, Todai Institutes for Advanced Study, The University of Tokyo,
Kashiwa, Chiba 277-8583, Japan (Kavli IPMU, WPI)
3Institut d’Astrophysique de Paris, CNRS UMR 7095 and UPMC, 98bis bd Arago, F-75014 Paris, France
(Accepted XXX. Received YYY; in original form ZZZ)
Abstract
We develop a new perturbation theory (PT) treatment that can describe gravitational dynamics of large-scale structure after shell-crossing in the one-dimensional cosmological case. Starting with cold initial conditions, the motion of matter distribution follows at early stages the single-stream regime, which can, in one dimension, be described exactly by the first-order Lagrangian perturbation, i.e. the Zel’dovich solution. However, the single-stream flow no longer holds after shell-crossing and a proper account of the multi-stream flow is essential for post-collapse dynamics. In this paper, extending previous work by Colombi (2015, MNRAS 446, 2902), we present a perturbative description for the multi-stream flow after shell-crossing in a cosmological setup. In addition, we introduce an adaptive smoothing scheme to deal with the bulk properties of phase-space structures. The filtering scales in this scheme are linked to the next-crossing time in the post-collapse region, estimated from our PT calculations. Our PT treatment combined with adaptive smoothing is illustrated in several cases. Predictions are compared to simulations and we find that post-collapse PT with adaptive smoothing reproduces the power spectrum and phase-space structures remarkably well even at small scales, where Zel’dovich solution substantially deviates from simulations.
keywords:
Large-scale structure – Cosmology – Vlasov-Poisson equation
††pubyear: 2015††pagerange: Post-collapse perturbation theory in 1D cosmology – beyond shell-crossing–F
1 Introduction
It is currently admitted that processes of structure formation in the Universe are mainly dominated at large scale by an invisible component called dark matter. Although the microscopic origin of dark matter is still unclear, it is macroscopically described as a self-gravitating collisionless fluid following the collisionless Boltzmann or Vlasov equation in a cosmological background,
[TABLE]
supplemented with the Poisson equation for the Newton potential ,
[TABLE]
where f(\mbox{\boldmathx},\mbox{\boldmathv},t) is the phase-space density at comoving position , peculiar velocity and time , is the expansion factor of the Universe, is the gravitational potential and is the average dark matter density.
In the standard picture of structure formation, dark matter was initially cold i.e. with a virtually null local velocity dispersion, so the six-dimensional phase-space distribution is effectively reduced to a three-dimensional hyper-surface and this remains true at all times thanks to the Hamiltonian nature of the system. At early times, dark matter thus follows the single-stream flow regime, with a velocity field {\overline{}\mbox{\boldmathv}} uniquely determined as a function of position, and its evolution is that of a pressure-less fluid with a phase-space distribution function given by
[TABLE]
where is the density contrast of the dark matter distribution, initially of very small magnitude, as well as {\overline{}\mbox{\boldmathv}}.
Substituting Eq. (3) into Eqs. (1) and (2) and taking the zeroth and first velocity moments yields the Eulerian formulation of large-scale structure dynamics
[TABLE]
Although the single-stream flow is, strictly speaking, valid only during the early phase of structure formation, the above equations have been shown in practice to describe accurately nonlinear mode-coupling in the weakly nonlinear regime and provide a solid basis for perturbation theory (PT) calculations to predict statistical quantities of large-scale structure, such as the power spectrum or the two-point correlation function of the matter distribution (see, e.g. Bernardeau et al., 2002).
With the advent of cosmological observations aiming at precisely mapping galaxy clustering at large scales, techniques have been developed to improve on the slow convergence of the perturbative expansion as well as to accelerate higher-order calculations (e.g., Crocce & Scoccimarro, 2006a, b, 2008; Valageas, 2007; Taruya & Hiramatsu, 2008; Matsubara, 2008; Bernardeau et al., 2008; Pietroni, 2008; Hiramatsu & Taruya, 2009; Taruya et al., 2009; Okamura et al., 2011; Matsubara, 2011; Bernardeau et al., 2012; Crocce et al., 2012; Bernardeau et al., 2014; Taruya et al., 2012; Valageas et al., 2013). Improved calculations involving the next-to-next-to-leading order called two-loop give a remarkable agreement with cosmological -body simulations at weakly nonlinear scales and have been applied to observations. However, this does not ensure the convergence of perturbative calculations including higher-order corrections, e.g., three-loop. In fact, a direct calculation at three-loop order suggests a very large UV contribution to the large-scale modes through nonlinear mode-coupling (Bernardeau et al., 2014; Blas et al., 2014) and indicates a break down of higher-order perturbative expansions even at large scales.
Deficiency of PT calculations has been also highlighted in recent numerical analyses. Nishimichi et al. (2016) directly measured the coupling between different scales in cosmological -body simulations and found that the actual contribution from small scales to the large-scale modes is suppressed, as opposed to the prediction of PT based on the single-stream approximation. These facts imply that the validity of the single-stream treatment is questionable even at large scales, and higher-order perturbative corrections need to be remedied with a proper account of small-scale dynamics, where the multi-stream flow contribution is important.
One way to account for small-scale dynamics consists of using some ansatz to summarize the main physical effects in the multi-stream regime, such as Burgers’ equation (e.g., Gurbatov et al., 1989; Bernardeau & Valageas, 2010) and various alternatives (e.g., Sahni & Coles, 1995), which roughly amount to adding a source term in the right hand side of equation (5). In particular, it was proposed recently to use effective fluid equations giving account of the non-vanishing stress tensor arising in Eq. (5) when calculating the first velocity moment of Vlasov equation in the multi-stream regime. This effective-field theory approach has attracted a lot of interest and has been studied in detail (e.g., Baumann et al., 2012; Carrasco et al., 2012; Hertzberg, 2014; Baldauf et al., 2015). The drawback of this approach is however that the parameters in the stress tensor characterizing the small-scale dynamics need to be calibrated with -body simulations to keep predictions of perturbative calculations under control. Furthermore, these parameters generally vary with cosmology and with time and no prediction with PT is really possible in this framework independently of -body simulations.
An alternative approach that we consider in this paper consists in going back to a more fundamental description, i.e. the Vlasov-Poisson system, Eqs. (1) and (2), and trying to follow accurately the phase-space distribution function in the multi-stream regime, which is essential for describing the formation of dark matter halos. For this purpose, it is useful to employ Lagrangian PT (e.g., Zel’dovich, 1970; Shandarin & Zeldovich, 1989; Bouchet et al., 1992; Buchert, 1992; Buchert & Ehlers, 1993; Bouchet et al., 1995; Bernardeau et al., 2002, and references therein), where the small parameter is the displacement field. We shall consider one-dimensional gravitational dynamics in standard (three-dimensional) cosmology. In this case, large-scale structure dynamics is described by the gravitational interaction of massive parallel infinite planes moving left and right along a fixed axis, while Hubble expansion is taking place as usual in all the directions following standard Friedman-Lemaître equations. Despite its simplicity, one dimensional dynamics defined as such still displays a rich physical content which somewhat shares the same features as 3D clustering. This is partly the reason why the 1D model has recently attracted much attention (e.g., Benhaiem et al., 2013; McQuinn & White, 2016; Vlah et al., 2016; Baldauf et al., 2016). In particular, the Zel’dovich solution provides an exact solution for the dynamics of massive sheets before shell-crossing (Zel’dovich, 1970; Shandarin & Zeldovich, 1989), and thus, starting with Zel’dovich solution, a tractable perturbative treatment of multi-stream flows is made possible based on a Lagrangian description. The analysis in the present paper is an extension of the method developed in Colombi (2015) to the cosmological setup. We will describe perturbatively post-collapse dynamics around the shell-crossing region and apply it to several cases including random initial conditions.
In addition, we shall present a novel regularization scheme reducing the impact of small-scale clustering and improving greatly our post-collapse PT predictions of large-scale structure statistics. The idea is to apply adaptive smoothing to initial density peaks and to better capture the bulk properties of phase-space structures in the post-collapse regions, where the interaction or merger of halos is supposed to be significant. The idea is similar to the peak-patch treatment proposed by Bond & Myers (1996) and subsequent works (see, e.g. Monaco et al., 2002), but we here implement it in the PT prescription in order to better describe the late-time post-collapse dynamics. Indeed, our perturbative approach will not allow us to follow post-collapse dynamics beyond next-crossing time, although an iterative prescription such as proposed by Colombi (2015) could make this possible but is out of the scope of the present work. Here, we will show instead that the predictions with adaptive smoothing reproduce simulations remarkably well even at nonlinear scales, where Zel’dovich solution significantly deviates from simulations.
Note finally that an analytical study in 1D of course represents only a first step toward a proper description of 6D phase-space dynamics. Apart from a few examples including self-similar solutions (Fillmore & Goldreich, 1984a, b; Bertschinger, 1985; Ryden, 1993; Lithwick & Dalal, 2011; Alard, 2013), little is indeed known analytically for the phase-space dynamics of Vlasov-Poisson systems, particularly in cosmology. Therefore, further development of analytical treatment in this framework seems indispensable, complementary to simulations, and even helpful to cross check simulation codes. This is all the more motivated by the fact that simulations in 6D phase space have now become available thanks to recent efforts on the development of pure Vlasov codes (Yoshikawa et al., 2013; Sousbie & Colombi, 2016; Hahn & Angulo, 2016).
This paper is organized as follows. In § 2, we begin by describing the basic setup of our calculations in one-dimensional cosmology. We then discuss in § 3 the analytic treatment beyond shell-crossing and develop post-collapse PT. § 4 introduces a novel filtering scheme to improve the perturbative description of post-collapse dynamics, suited for a system with interacting clusters. In § 5, analytic calculations in the post-collapse perturbative framework are tested against controlled -body experiments. Several cases are considered, including the single initial sine wave and random initial conditions. Finally, § 6 is devoted to discussion and conclusion.
2 1D cosmology
2.1 Basic setup
We consider an ensemble of massive parallel infinite planes moving along axis and interacting through gravitational force in the expanding Universe. The Lagrangian equations of motion of the planes are given by
[TABLE]
where and are respectively the comoving position and peculiar velocity of each plane, the gravitational potential, the average matter density, the density contrast and the expansion factor of the Universe.
To simplify the equations, it is useful to introduce the super-conformal time defined by (e.g., Doroshkevich et al., 1973; Martel & Shapiro, 1998),
[TABLE]
the new velocity and potential :
[TABLE]
Then, Eqs. (7)–(9) are simplified:
[TABLE]
where is the matter density parameter of the Universe and is the Hubble constant.
To deal with this system, especially for the dynamics after shell-crossing, it is useful to introduce in the cold case considered here the Lagrangian coordinate defining the initial position of the planes in the absence of perturbation
[TABLE]
and to express the subsequent position and velocity of each mass element as and . In particular, mass conservation implies
[TABLE]
This equation is valid in the single-stream regime, i.e. as long as remains monotonic as a function of .
With the new expressions above, the solution of the equations of motion can be formally written as:
[TABLE]
where and are the initial positions and velocities given at a starting time , which will be specified below.
In what follows, we consider the dynamics of the cosmological system in a finite-size box with periodic boundaries, . The solution of Poisson equation, Eq. (14), can be expressed in an integral form as:
[TABLE]
The derivation of this expression is presented in Appendix A. Then, the force exerted on a mass element at position is given by:
[TABLE]
where function represents the Heaviside step function. In the above, we used the fact that fluctuations averaged over space vanish, . In the limit , equation (20) of course converges to the well-known result in the infinite space.
2.2 Pre-collapse dynamics
In one-dimension and in the cold case, Zel’dovich approximation (Zel’dovich, 1970) is known to provide the exact solution for the dynamics of mass elements before shell-crossing. It can be explicitly written as
[TABLE]
Here, function corresponds to the linear growth factor satisfying the following equation:
[TABLE]
Note that in terms of cosmic time , Eq. (22) reduces to the standard form of linear evolution equation:
[TABLE]
The Zel’dovich solution in Eq. (21) contains an arbitrary function that we call displacement field. It is related to the linear density field given at a very early time ( or ) through
[TABLE]
Since Zel’dovich solution is exact until shell-crossing, we do not necessarily assume that the evolved density field is small. One may thus consider the situation where in a region around a Lagrangian coordinate , the density field becomes large and the region will undergo shell-crossing at time . Recalling the fact that in 1D, a shell-crossing point corresponds to an inflection point of the mapping from Lagrangian to Eulerian space, the actual conditions for shell-crossing are given by
[TABLE]
At the time of shell-crossing, , the solution (21) around the shell-crossing region can be expanded as follows:
[TABLE]
Using Eq. (24), the conditions (25) for shell-crossing can be rewritten as
[TABLE]
In other words, shell-crossing takes place at local density peaks. Hence, in the 1D case, the conditions for shell-crossing are equivalent to peak constraints and collapse take places exactly when the Eulerian linear density contrast at the peak positions becomes equal to unity.
3 Perturbative treatment of post-collapse dynamics
We are interested in the dynamics of mass elements after shell-crossing, when the Zel’dovich solution is no longer valid. In this section, extending the work of Colombi (2015), we develop perturbative calculations to deal with the multi-stream motion around the shell-crossing region.
3.1 Post-collapse PT
The basic formalism to treat post-collapse dynamics can be described as follows. Starting with the cold initial conditions in Sec. 2.2, we first follow pre-collapse dynamics with the exact Zel’dovich solution. Then, in regions undergoing shell-crossing, we switch to a perturbative treatment and compute the backreaction to the Zel’dovich flow based on an explicit functional form of the displacement field in the shell-crossing region. To be precise, we compute the force from Eq. (20) exerted at each position using the extrapolation of the Zel’dovich flow. Integrating the force over time, we obtain the correction for the velocity to the Zel’dovich motion from Eq. (18). Further integrating the corrected velocity over time, a correction for the position is obtained from Eq. (17). Throughout these calculations, we assume that the collapsing region is small, which allows us to express the phase-space configuration in terms of polynomial forms of low order of the Lagrangian coordinate.
Let us focus on a collapse point centered on Lagrangian position and perform a perturbative description of the multi-stream flow around . As we mentioned earlier, the conditions for shell-crossing are given by Eq. (25). Shortly after shell-crossing time , the multi-stream flow has just started growing in a small Lagrangian region around and the displacement field can be approximately described by a third-order polynomial in . Since the deviation from Zel’dovich flow is small, the motion of a mass element around may be expanded at third-order in for as
[TABLE]
with the time-dependent coefficients , and defined by
[TABLE]
where the prime denotes the derivative with respect to . Here, we used Eqs. (24) and (27). The above expressions imply that the shell-crossing point slightly moves from to and the shell-crossing structure develops as an “” shape, as shown in Fig. 1.
We now introduce the new Lagrangian coordinate, , for which the structure of the local density peak is symmetric. Since we expand the displacement field at third-order, the equation inside the shell-crossing region has three solutions. An example is shown in Fig. 1 for the equation with solutions . In general, the three ordered roots of the equation , when they exist, are related to each other through the following expressions:
[TABLE]
Here, the Lagrangian extent of the boundary of the multi-valued region, also related to the Lagrangian position of the left caustic in Fig. 1 through , is expressed as
[TABLE]
Because of the time dependence of the coefficients (see Eqs. 30 and 31), the position of the boundary gradually changes in time. The leading-order expression in powers of time reads
[TABLE]
with
[TABLE]
In other words, it also means that a given element of fluid of Lagrangian position will enter the multi-valued region at some time defined by and then subsequently coincide exactly with the Lagrangian position of the caustic at a time with . By inverting the relations and , we obtain
[TABLE]
and
[TABLE]
The quantities defined above thus play a crucial role to disentangle the single-stream regime from post-collapse dynamics, as well as the inner part of the “” shape of Fig. 1 from its outer part, the Lagrangian position of the caustic acting as a separator in the latter case. The expression of the force will indeed be different between the Zel’dovich single-flow regime and the multi-valued regime composed of an outer caustic region, , and an inner caustic region, .
One important remark for the subsequent calculations is that the boundaries coordinates and are assumed to be small and that the backreaction to the position and velocity inside the multi-valued region can be described perturbatively in a polynomial form of [partly allowing the fractional power of ]. Thus, the resultant expressions are, rigorously speaking, only valid for a short period after collapse time, but we shall see in practice, as also shown by Colombi (2015) in the non-cosmological case, that they remain impressively accurate even up to next-crossing time (see Sec. 5), which will allow us to set up the framework for a powerful self-adaptive scheme.
3.2 Computing the force in the multi-valued region
To derive the corrections to the motion, we first compute the force exerted on a mass element inside the multi-valued region, (see Fig. 1). Note again that the outer regions and are described by the Zel’dovich solution.
The force in the multi-valued region is computed using Eq. (20), dividing each integral of the right-hand-side into three contributions:
[TABLE]
Assuming that the collapse region, , is small enough, the contributions to the integrals from each domain can be computed analytically, based on the geometrical setup in Fig. 1. The detailed calculations are presented in Appendix B. Summing up all the contributions given in Eqs. (102), (103), (106) and (110), the force exerted on the mass element at inside the multi-valued region becomes
[TABLE]
with the functions and respectively defined by
[TABLE]
and
[TABLE]
where the quantities , , and are defined by Eqs. (29)–(31). Note that in deriving Eq. (40), we have assumed that the system follows Zel’dovich solution if . Since the resultant expressions are written in terms of the local quantities characterizing the density peak at position and the shell-crossing time , Eq. (40) is in fact still applicable to other shell-crossing regions possibly appearing elsewhere in the region , allowing one to generalize the result to a smooth random initial field with multiple peaks. One issue, discussed later, is then to treat mergers, i.e. the case when for instance two multi-stream regions start overlapping.
3.3 Corrections to the Zel’dovich flow: basic post-collapse PT results
Given the explicit expression for the force in the multi-stream region and using the formal solution given by Eqs. (17) and (18), we now compute corrections to the Zel’dovich flow that we write as follows:
[TABLE]
We noticed in previous section that, depending on the Lagrangian position of interest, the expression for the force is different and we have to divide the domain of the integrals in Eqs. (49) and (50) into several pieces:
: the position is located in the single-valued region (i.e., ) and the motion is still described by the Zel’dovich solution. We have
[TABLE]
: the position is in the multi-valued region and satisfies , i.e. lies in the outer part of the caustic. In addition to the Zel’dovich flow, the corrections arising from the multi-stream flow need to be added and we note them as follows:
[TABLE]
In this region, which correspond to the tails of the shape in Fig. 1, the system is globally expanding in phase space and gaining energy, at variance with the central part, which on the contrary, contracts. This can be understood from computing the variation of energy of a typical test particle during a fraction of orbit. This phenomenon was studied in the non-cosmological case by Colombi (2015) and the results should not be fundamentally different in the cosmological case studied here.
: this corresponds to , i.e. the position now lies in the inner part of the multi-stream region with respect to the caustic. Similarly to the above case, the backreaction to Zel’dovich flow needs to be computed, including at present both the multi-stream dynamics of the inner part and the incoming flow from the outer part. We may write
[TABLE]
Note that the corrections and partly come from outer part contributions, and . As mentioned in the previous point (ii), in this region, which corresponds to the central part of the on Fig. 1, the system is globally contracting and losing energy in favor of the tails of the .
In what follows, we compute the backreaction to Zel’dovich flow and derive the expressions for and in each domain. The calculation for the corrected motion is rather straightforward but needs several steps. Here, we give a brief sketch of the calculation, deferring details to Appendix C. The final results are given in Eqs. (58) and (60) for the outer part, Eqs. (62) and (64) for the inner part, together with the coefficients in Tables 1 and 2. Note that although our expansion is rigorously valid only at third-order in , higher-order contributions will appear in the corrections to provide a continuous solution, up to fifth- and seventh-order for the velocity and the position, respectively.
3.3.1 Velocity and position in the outer part:
In the outer part of the multi-valued region, (ii), the correction to the velocity can be expressed as
[TABLE]
In the above, while the first integral is performed with the help of formulae in Appendix E, the second integral is computed exactly. The whole derivation is provided in Appendix C.1. The resultant expression can be summarized as
[TABLE]
with the time-dependent coefficients given in Table 1. Note that the coefficient implicitly depends on the Lagrangian position through but is not Taylor expanded with respect to , for simplicity. A fully analytical theory, in particular to predict the power spectrum of the projected density field, would in principle require such a Taylor expansion. Performing it should not change significantly the performances of post-collapse PT as presented in this article.
Once is obtained, the expression for the correction is evaluated by further integrating Eq. (58) over time:
[TABLE]
The resultant expression becomes (see Appendix C.2 for derivation):
[TABLE]
with the time-dependent coefficients presented in Table 1. Note again the dependence on of parameter through and the fact that we did not Taylor expand it in polynomials of , as it would be required for a fully analytical theory.
3.3.2 Velocity and position in inner part:
Let us consider now the inner part of the multi-valued region (iii). For the correction to the velocity, , the expression to be evaluated is divided into three contributions:
[TABLE]
Evaluating each contribution in Appendix C.3, we obtain
[TABLE]
with the time-dependent coefficients given in Table 2. Note that the coefficient in the above equation is the same one as in Eq. (58).
Similarly, the correction to the position, , is divided into three pieces:
[TABLE]
where the function represents the correction of the velocity given in either Eq. (58) or Eq. (62), but without the contribution corresponding to the Zel’dovich motion, i.e., the last term in their expressions. The term is the same quantity as in Eq. (60) (see Table 1).
Based on the calculation presented in Appendix C.4, we obtain the expression for :
[TABLE]
with the time dependent coefficients given in Table 2.
3.4 Higher-order corrections
The perturbative description of post-collapse dynamics introduces two kinds of approximations around the singularities (shell-crossing points). Firstly, we assumed the position and velocity to be described by third-order polynomials in in the multi-stream regions. To preserve continuity between the multi-stream and the single-stream regime, we added higher-order terms in to the solutions for and , but strictly speaking, these solutions are valid locally only at third-order in around the singularities. Secondly, these solutions are correct only very shortly after collapse time up to second-order in . Our post-collapse dynamics can thus be improved at a twofold level:
(a)
Improvement at the spatial level: we could try, starting from initial conditions, to use a description of the shape of Fig. 1 at higher-order in , e.g. fourth- or fifth- instead of third-order. In this latter case, deriving the solution of equation comes down to find the roots of a third or fourth-order polynomial instead of a second-order one, which makes, in addition to the symmetry breaking of the shape, the calculations much more cumbersome, even though a series expansion of the roots at fourth or fifth-order in remains tractable and might bring interesting improvements on the description of the tails of the shape.
(b)
Higher-order in time: we can introduce higher-order corrections in time to follow more accurately the evolution of the central part of the singularity after collapse. Note that after integration of the equations of motion, such corrections bring out higher-order terms in . Strictly speaking, the right way to proceed would consist of performing “post-post-collapse” dynamics, i.e. of computing a counter-term to the Zel’dovich+leading-order post-collapse solution. This procedure is cumbersome, and seems to make sense only if performed simultaneously with a description at higher-order in , in particular because the tails of the structure in Fig. 1 also influence the dynamics of the central part of the system (but this remains to be proved).
Going beyond leading-order in post-collapse dynamics thus seems rather challenging even in the simple one-dimensional case considered in this article.
However, we can still try to investigate some potential improvements of the theory while maintaining the complexity of the calculations at an acceptable level. For this purpose, we consider here mainly two approaches. The first one just consists in incorporating next-to-leading order terms into the expressions for and (and and ). These quantities, which determine the boundary between the single-stream and the multi-stream regime, are indeed crucial because they correspond to the main scale/timescale of the system. The second approach consists, in addition, of improving the evaluation of the time integrals in Eqs. (49) and (50). In Sec. 3.3, when evaluating these time integrals, a part of the integrands is Taylor-expanded and the results are presented at leading-order in time. It is thus possible to include higher-order corrections from the Taylor-expanded integrands.
Again, one has to be aware that the two corrections we propose might not improve the results, since a correct higher-order time approach, as discussed in point (b) above, would actually require to compute counter terms to the standard post-collapse result, which we do not do here. In fact, our higher-order corrections might even worsen the results since the corrections brought by such additional terms might have opposite sign compared to the iterated counter-term approach.
The derivation of the higher-order corrections we propose here and summarize just below is rather straightforward compared to its rigorous alternatives (a) and (b) above, but still cumbersome. The resultant expressions now include higher-order polynomials of as well as new contributions, which are given in Appendix D. To be more specific, we consider the following cases:
Higher-order for the expressions of critical time and position (hc): based on the basic results in Sec. 3.3, we include the contributions arising from the next-to-leading order expressions for and (as well as and ). The calculation to derive higher-order terms is almost the same as presented in Sec. 3.3, but we replace the leading-order expressions for and in Eqs. (35) and (37) with those including the next-to-leading order, given in Eqs. (120) and (122). The higher-order results, denoted by and , correspond to the expressions summarized in Appendix D, but where the terms involving the quantities or are all set to zero.
Higher-order in time (ho): in addition to the above treatment, we also include higher-order corrections to the time-integrals, based on the formulas given in Appendix E. The resultant expressions, and , are summarized in Appendix D.
Finally, we also consider a semi-analytic approach in order to smoothly connect the inner part of the multi-stream region to the Zel’dovich solution:
Third-order spline (spl): basic post-collapse PT is applied to the inner part, , but we take the higher-order expression for . In the outer part, , on the other hand, a third-order spline interpolation is used to smoothly connect the inner post-collapse region to the Zel’dovich solution, with again, the higher-order expression for .
4 Improved treatment with adaptive smoothing
Post-collapse PT as described in § 3 provides a way to describe the dynamics around an initial density peak shortly beyond collapse time. In practice, the quality of the description is expected to be good until next-crossing time, at least this is what is known in the non-cosmological case (Colombi, 2015). Using the leading-order expression given in Eq. (64), the duration of time spent between collapse and next shell-crossing can be estimated from Eq. (55),
[TABLE]
which is derived by solving the equation .
Beyond time , post-collapse PT is expected to fail, similarly as Zel’dovich approximation fails beyond . The only way to fix this problem at the fine grain level would be to use an iterative procedure similar to what was proposed by Colombi (2015) to reconstruct the spiral structure built up during the course of the dynamics. In addition, post-collapse PT assumes the initial density distribution to be sufficiently smooth. While this is certainly expected at very small scales in the standard Cold Dark Matter scenario, the issue of mergers arises. Indeed, post-collapse PT is valid only if the forming halo is isolated. If there are two adjacent collapsing regions, they might, at some point, overlap in Lagrangian space, hence merge with each other in Eulerian space. After the merger, a separate description of each halo becomes inadequate. In this case, the prediction of post-collapse PT, which is in principle able to handle multiple but non-overlapping structures in Lagrangian space, will break down much earlier than expected.
One simple way to fix both limitation in time and mergers issues, already applied in the past to the Zel’dovich approximation (see, e.g. Coles et al., 1993; Sahni & Coles, 1995, and references therein), consists in trying to describe the dynamics at the coarse-grained level by introducing a smoothing procedure in Lagrangian space, so that post-collapse PT remains applicable, in particular that the condition remains valid for the system with coarse-grained initial conditions. This means that single or composite halos as well as mergers are locally summarized by a simplified shape (Fig. 1) corresponding to a single peak of the coarse-grained initial density field and such that , while retaining the correct large-scale dynamics. The hope is that this simplified structure still captures the most important features of the nonlinear system, e.g. its mass and its typical size.
Of course, the coarsening scale should depend on environment, hence on Lagrangian position , so for best results, adaptive smoothing of initial conditions has to be introduced. In this sense our method is very similar to the “peak-patch” treatment by Bond & Myers (1996) (see also Monaco et al., 2002), but we implemented it in the post-collapse PT framework. To be precise, the practical set up of post-collapse dynamics with adaptive smoothing is summarized as follows:
The initial density field is smoothed at various scales, employing the sharp- filter function in Fourier space:
[TABLE]
with being the Heaviside step function. The cutoff wavenumber ranges from the fundamental mode to a certain value , where is an integer corresponding to the effective resolution of the simulation, as discussed more in details in § 5.1. In our investigations, we covered all the possible values of , i.e. with , but this is costly and probably not required, in practice.
Zel’dovich solutions are constructed from the smoothed density fields. Then we identify the critical points in Lagrangian space that have undergone shell-crossing at time , where is the time of interest for the analyses. For each critical point, various quantities needed for perturbative calculations, such as and , are also computed.
Starting from the largest smoothing scale (i.e. ), we compute the next-crossing time associated to each critical point. For the points satisfying with a constant parameter, the post-collapse correction to the Zel’dovich flow is obtained in the interval , which is labeled.
The procedure (iii) is repeated, but with a smaller smoothing scale, i.e., larger , and only for peaks lying in the unlabeled regions, until the cutoff wavenumber reaches . Note that if two Lagrangian regions overlap, e.g. and , corresponding to two different values of the smoothing scale, , the interval corresponding to the smallest smoothing scale, i.e. largest wavenumber , prevails. This is the way we solve the cloud-in-cloud problem, but not the only possible choice.
At the end of the iterative process, post-collapse dynamics is applied as well for the rest of the critical points corresponding to the highest resolution scale that still require treatment, i.e. for those with . The procedure is the same as described in point (iv) above, i.e. only critical points which do not fall in a labeled region are treated, and if it is the case, then their Lagrangian interval prevails for applying post-collapse dynamics. For the rest of the unlabeled regions, standard Zel’dovich dynamics is applied and the mapping uses the initial density corresponding, of course, to the highest resolution scale.
Note that the approximation of the phase-space density distribution constructed this way is generally discontinuous, because we collect perturbative solutions for the displacement field in phase space coming from different coarse-graining scales without imposing smoothness at the transition between these solutions. In this respect, the procedure given above is not entirely consistent and remains to be improved in the future. Nevertheless, the discontinuities only affect the results at small scales and in regions where the projected density is small. Indeed, the discontinuities correspond to transitions between tails of the shape representing a halo and a non collapsed region or transitions between tails of two shapes in the case of a merger. As illustrated below through numerical examples, these regions of phase space have small density contrast and thus contribute weakly to second- or higher-order statistics, hence to the power spectrum which is the statistics under focus in this work. With an appropriate choice of , we shall show that the adaptive smoothing procedure provides a substantial improvement on the post-collapse prediction for the overall phase-space structure of the system and the power spectrum of the projected density.
5 Comparison with controlled -body experiments
We are now in position to compare the predictions of post-collapse PT with simulation results. In § 5.1, we briefly describe the one-dimensional cosmological -body code we wrote to perform the numerical analyses of this article and discuss various important parameters settings. We then present the results of the comparison between theory and numerical experiments in various cases with cold initial conditions: formation of an isolated single halo (§ 5.2), merging of two neighboring halos (§ 5.3), as well as random initial conditions with Cold Dark Matter-like (CDM) and power-law power spectra (§ 5.4 and 5.5). In the later case, we shall discuss more in details highly nonlinear dynamics, in particular predictions from stable clustering, the expected outcome on the power spectrum from halos with power-law density profiles, as well as the statistical properties of the phase-space sheet resulting from a caustic distribution.
5.1 Simulations
To solve numerically Vlasov-Poisson equations in one dimension, we wrote a Fortran90 particle-mesh code, Vlafroid. In this code, publicly available through the following web page, www.vlasix.org, Poisson equation is solved on a periodic grid of fixed resolution using Fast Fourier Transform, after projecting the particles on the grid using simple Cloud-in-Cell interpolation (see, e.g. Hockney & Eastwood, 1988). The particles positions and velocities are updated at each time step using a standard predictor-corrector integration scheme. After the predictor step,
[TABLE]
Poisson equation is solved and the force is estimated on the grid by direct derivation in Fourier space. The acceleration of each particle is computed by linear interpolation of the force. This is followed by the corrector step:
[TABLE]
The value of the time step is slowly varying with time and is bounded by 3 constraints:
- (i)
a limit on the relative variations of the expansion factor during a time step to have a good description of the linear growing modes at early stages of the simulation,
[TABLE] 2. (ii)
a dynamical condition which states that the time step should be a small fraction of the smallest harmonic oscillator time scale (see, e.g. Colombi & Touma, 2014),
[TABLE]
where is the maximum value of the projected density (normalized to be unity in average), 3. (iii)
the standard Courant-Friedrichs-Lewy (CFL) condition
[TABLE]
with the maximum magnitude of the velocity, stating that a particle should not travel more than a fraction of a grid cell size during a time step.
For all the runs performed in this article, our practical and safe choice of the time step parameters was:
[TABLE]
The evolution of single or composite objects as studied in § 5.2 and 5.3 does not pose any problem from the numerical point of view. Indeed, we consider here only a small number of dynamical times, so we do not expect significant defects due to the discrete nature of the particle distribution (for instance, such as described in Melott et al., 1997). The case of a random field, as studied in § 5.4 for CDM and in § 5.5 for power-law power spectra, is less trivial. Indeed, one needs accurate numerical modeling of the dynamics of each single structure, which has to be sampled with many particles and be resolved with a sufficient number of grid elements. In particular, to avoid contamination by the discrete nature of the particle distribution, we require, in average, many particles per grid element:
[TABLE]
To avoid the formation of small structures that would be insufficiently resolved by the computational grid, there must be a coarse-graining scale below which the initial conditions are smooth. This also facilitates, of course, the calculation of (post-collapse) PT predictions. Here we impose a cut-off on the power spectrum of initial conditions using the sharp- filter introduced in Eq. (66), with
[TABLE]
In practice, a choice which works well and adopted in this paper is
[TABLE]
which means that is five time smaller than the Nyquist frequency of the grid. Since the cut-off scale on the power spectrum can affect the dynamics compared to the true, un-smooth system, we shall perform some convergence studies to determine the available dynamic range in which the results are not influenced by .
5.2 Formation of a single structure
Let us first present the simplest case of the formation of a single structure seeded with the following sinusoidal density perturbation:
[TABLE]
We adopt the Einstein-de Sitter cosmology (, with ) and the box size of the simulation is arbitrarily set to . The simulation was started at expansion factor with initial amplitude and was run using . PT analyses were performed on the fly as part of routines of the code, when needed.
Fig. 2 shows snapshots of the phase-space structure (upper inserts) and projected density profiles (lower inserts) obtained from the simulation (red) and post-collapse PT (blue) based on the basic prescription described in Sec. 3.3. Also, for reference, the Zel’dovich solution is plotted as a green dashed line on each panel. Although the three curves are identical at early time, differences show up, as expected, after collapse, which occurs at . A closer look at the central part of the system reveals however that post-collapse PT captures the trends of the -body simulation quite well and succeeds to reproduce the structure of the inner part of the system until next-crossing time, which occurs at . This is in marked contrast with Zel’dovich solution, which, without the back reaction contribution present in the post-collapse correction, overshoots the displacement of the phase-space sheet after first crossing time, in particular of the caustics. After next-crossing time, , the phase-space structure given by post-collapse PT goes away from the -body result. Nevertheless, the predicted projected density profile still provides a reasonably good approximation of the evolved density structure, which is clearly not the case for the Zel’dovich solution.
Fig. 3 presents results of post-collapse PT when including the higher-order corrections discussed in § 3.4. While the three prescriptions examined here, namely hc, ho and spl, account for higher-order terms in different manners, they are all based on the force expression in Eq. (40), which is valid only shortly after collapse. Indeed, as discussed in § 3.4, our implementation of higher-order corrections is not performed in a fully rigorous way, that would require to compute a back reaction to post-collapse PT, similarly as post-collapse PT introduces a back reaction correction to Zel’dovich solution. It also does not propose an iteration procedure that would allow us, similarly as in Colombi (2015), to follow the system behind next-crossing time. Hence, all of the treatments break down at the end of the simulation. Nevertheless, taking the higher-order corrections to and into account improves the prediction of the outer boundary of the multi-valued region, and the phase-space description becomes visually better than that of the basic PT prescription in Fig. 2. Even after next-crossing time (), the predicted phase-space structure still remains consistent, although the size of the halo is prone to be over-estimated, leading to an extended density profile, the best result being obtained with the spl prescription. Still, one can notice that if only the projected density is under consideration, basic post-collapse PT does better than any of its higher-order counter-part, except for spl which, not surprisingly, provides nearly equivalent results, since it is exactly the same as post-collapse PT in the center of the system. As far as statistics depending on the projected density are concerned, such as the power spectrum that will be studied later, we thus expect from the analyses of this single halo that the simplest prescription of post-collapse PT is going to perform best.
5.3 Merger
We examine now the second simplest case of two merging structures in Einstein-de Sitter cosmology (, and ). The initial density contrast is given by
[TABLE]
The constant is determined so as to satisfy the condition . With , and , we have . The overall amplitude is set to and the simulation, of box size unity (), was started at and was run using .
The results are shown in Fig. 5 for basic post-collapse PT and Fig. 5 for the higher-order variants of post-collapse PT. Here, we focus on the merger phase which takes place after collapse of each individual structure. Figs. 5 and 5 demonstrate how the adaptive smoothing introduced in Sec. 4 improves PT prescriptions to describe the overall phase-space structure after the merger. In each figure, the upper and lower panels show the results without and with adaptive smoothing, respectively. The free parameter controlling adaptive smoothing, , is set here to for post-collapse PT and for the Zel’dovich solution.
Without adaptive smoothing, as soon as the two clusters cross each other (), both post-collapse PT and Zel’dovich solution start failing to describe the real dynamics in the -body simulation. Indeed, while post-collapse PT tries to provide a local perturbative correction for each single cluster, it is clearly unable to account for nonlinear couplings taking place at larger scales, where it simply follows Zel’dovich motion. As a result, the location of multi-valued regions predicted by post-collapse PT, as well as their shape, largely deviates from the correct one.
On the other hand, implementing adaptive smoothing in the PT calculations significantly improves the predicted phase-space structure, starting already from the beginning of the merger phase (). After merging, dynamics of the central part of the system is now described by the displacement field derived from the smoothed initial density field. By construction, it summarizes the composite structure made of two sub-halos into a single cluster. While this is a rough approximation of the real dynamics, the description of the outer part of the system is substantially improved as well as the prediction for projected density profiles. Introducing both adaptive smoothing and higher-order corrections to post-collapse PT further improves the results from the visual point of view (Fig. 5), although the predicted size of the halo tends to be slightly overestimated. As already noticed in § 5.2, the best results are obtained for spl, which also gives, after adaptive smoothing, comparable projected density profiles to basic post-collapse PT.
The above results demonstrate that post-collapse PT with adaptive smoothing is effective in capturing the main features of phase-space structures. Strictly speaking, it does not give an accurate prescription for the fine structure of high-density regions, but it provides a way of regularizing or mitigating the impact of highly nonlinear dynamics, keeping the location and size of halos reasonably accurate. As we will see next, adaptive smoothing brings a drastic improvement on the prediction of power spectra for systems with random initial conditions. Furthermore, the introduction of adaptive smoothing makes the analytic calculations much less sensitive to the choice of the small-scale cutoff in the initial conditions so that PT predictions become more robust. In these respects, the determination of the extension of multi-streaming regions for a given smoothing scale of initial conditions, as performed in step (iii) of § 4, is essential, making the choice of crucial. Concerning the particular merger experiment studied here, and seem to be the best choices for tuning post-collapse PT and Zel’dovich solutions, respectively, and we shall adopt these values in subsequent section. We shall however discuss again about the choice of later.
5.4 Random initial conditions: CDM-like spectrum
Let us now consider random CDM initial conditions. Although there is no realistic setup in 1D, a way of mimicking as well as possible the 3D case may be to consider initial conditions given by a random Gaussian field with the following power spectrum (see, e.g. McQuinn & White, 2016):
[TABLE]
with being the 3D matter power spectrum of the initial density fluctuations. Here we use for the transfer function of Eisenstein & Hu (1998) who wrote the corresponding Fortran module that we included in Vlafroid. We set the cosmological parameters to those of the base CDM model determined by Planck (Planck Collaboration et al., 2016): total matter density parameter , cosmological constant , baryonic matter density parameter , Hubble constant km s*-1* Mpc*-1*, effective slope and root mean square deviation .111To be precise, the power spectrum amplitude is normalized by through
\displaystyle=\int\frac{{\rm d}^{3}\mbox{\boldmathk}}{(2\pi)^{3}}\,P_{\rm 3D}(k)\{W_{\rm th}(kR_{8})\}^{2}
(80)
(81)
with Mpc, and . The simulations were performed with a box size Mpc and started at initial redshift . Convergence of simulations results was tested by varying the number of particles , the resolution of the grid used to compute the force and the cutoff wavenumber of the initial power spectrum, , introduced in § 5.1.222Namely, we performed three sets of 50 simulations for , and . While the highest resolution runs with only allowed us to test the -body results and Zel’dovich solution without adaptive smoothing, lower resolution runs included full PT predictions with and without adaptive smoothing. Here, we present the results of the analyses for , , and Mpc*-1*. Our convergence study shows, unless specified otherwise below, that the results can be trusted for Mpc*-1*, i.e. in the available dynamic range shown in left panels of Figs. 7 and 7. To have sufficient statistics for measuring power spectra, we ran simulations with different random seeds.
Figs. 7 and 7 show the measured power spectra at different times (left panel of each figure) as well as snapshots of phase space for a given realization (right panel). Right panel of Fig. 7 is supplemented with Fig. 8, which displays in approximately the same scale interval, Eulerian positions and velocities as functions of Lagrangian coordinates. One can first notice that, in contrast to the 3D case, the amplitude of the power spectrum at small scales is not strongly enhanced in 1D, but damped. For the CDM like initial conditions considered here, function becomes nearly flat in the highly nonlinear regime, which is probably a consequence of stable clustering as discussed further in details in next section, where the case of scale-free initial conditions is considered.
In Fig. 7, the simulations results are compared to Basic PT predictions (solid blue) and Zel’dovich solution (dotted green), while, in Fig. 7, they are compared to various higher-order variants of post-collapse PT. Note, that except for the solid line that we explain below, all power spectra predicted by post-collapse PT or Zel’dovich solution are computed directly from the particle distribution evolved according to the displacement field predicted by each theoretical model under scrutiny. So it is important to mention that we did not actually compute analytically power spectra in these cases, although it would be in fact possible. The thin black solid lines, on the other hand, provide the analytic results for the Zel’dovich solution without adaptive smoothing, (e.g., Bond & Couchman, 1988; Schneider & Bartelmann, 1995; McQuinn & White, 2016):
[TABLE]
where is the un-truncated power spectrum, i.e. without suppression of power at . The black solid curves therefore differ slightly from the green dotted ones on left panel of Fig. 7, because of the effect of the cut-off at (and not because of discreteness effects or finite grid resolution). While the cut-off wavenumber has, as expected, a noticeable effect on the un-smoothed Zel’dovich solution as well as the un-smoothed post-collapse PT prediction and its variants (a property that we do not show but that we verified), its influence becomes negligible in the range of values of under consideration, both for the simulations and the theoretical models with adaptive smoothing.
As expected, without adaptive smoothing, both post-collapse PT prescriptions and Zel’dovich solution deviate from the simulations as soon as the nonlinear scale becomes larger than the cut-off scale , which is already the case in upper insert of left panel of Figs. 7 and 7, i.e. at . Still, at this redshift, post-collapse PT (or its variants) without adaptive smoothing is able to capture the main features of phase space better than Zel’dovich solution, although this is not obvious to decipher on upper left insert of right panel of Figs. 7 and 7. At lower redshift, halos relax to a highly nonlinear state and can merge together, which cannot be described well by Zel’dovich solution that gives elongated structures in phase space (middle-left and lower-left inserts of right panels of Fig. 7. As already argued in previous section, post-collapse PT provides only local corrections to Zel’dovich flow and therefore cannot account for nonlinear dynamics at large scales, hence presents the same defect. This elongation effect obviously leads to a strongly underestimated power spectrum. Note however that this depends also on the choice of the cut-off scale : if this scale would be kept, for the theoretical predictions, variable so that it remains close to the typical transition scale to the nonlinear regime, the results would improve greatly both for Zel’dovich solution and for post-collapse PT, as it is the case for upper insert of left panels of Figs. 7 and 7. But here, instead of varying the effective smoothing scale for the PT predictions, we now focus on adaptive smoothing, which should be even better because it provides locally the “optimal” softening of initial conditions.
When adaptive smoothing is employed, agreement between PT predictions and simulation improves strikingly both from the visual and statistical points of view, even for Zel’dovich solution which performs comparably well to post-collapse PT at , reproducing simulations results in the whole available dynamic range, including the highly nonlinear regime, with a slight underestimation of the power spectrum, though. Of course, superiority of post-collapse PT shows up at but Zel’dovich solution remains quite good. Keep in mind, however, that Zel’dovich solution with adaptive smoothing as presented here uses the calculation of next-crossing time to estimate the local smoothing scale, hence relies as well partly on post-collapse PT.
These results confirm the conclusions of § 5.3: adaptive smoothing allows one to summarize composite structures into a single S shape like halo with the right size so that two-point statistics matches the true one even in the highly nonlinear regime. In fact, PT predictions reproduce, at least partly, the plateau that can be observed at small scales in middle and lower insert of left panels of Figs. 7 and 7. In this regime that we will discuss more in details in next section, the system is most probably following stable clustering. Of course, because of the way we implement our adaptive smoothing, continuity of the phase-space sheet is not preserved anymore, as illustrated by right panel of Figs 7 and 7, but the overall description of phase-space structures is nevertheless improved tremendously. This is illustrated even more clearly by Fig. 8, which also shows, as expected, that positions predicted by post-collapse PT or Zel’dovich solution match the simulation significantly better than predicted velocities. This is a trivial consequence of the fact that positions correspond to velocities integrated over time: by construction, velocities are better approximated by theoretical models at early than at late times, hence their integrated counterpart, which corresponds to some averaged behavior, will compare better to the simulations. Note thus on right panel of Fig. 8 the excellent agreement with the simulation for the predicted position from basic post-collapse PT, even when significant mergers take place: one just needs to compare bottom right insert of right panel of Fig. 7 to bottom insert of right panel of Fig. 8 to be convinced by this state of fact. As a consequence, we obtain an excellent match between basic post-collapse PT with adaptive smoothing and simulations measurements for the power spectrum, which is a two-point statistics not depending directly on the quality of representation of the velocity. These arguments also apply, to a lesser extent, obviously, to the Zel’dovich solution.
When focusing again on the power spectrum (left panel of Figs. 7 and 7), a detailed comparison between various variant of post-collapse PT models reveals that the basic prescription remains the best among all of them, although spl performs nearly as well. This confirms the results of § 5.2 and 5.3, where we noticed that higher-order prescriptions for post-collapse PT tended to over-predict the size of halos and mergers (although this is not obvious at first sight when examining phase-space diagrams of right panels of Figs. 7 and 7), leading to a smaller power spectrum than the basic post-collapse PT prediction at small scales.
Note that, in agreement with intuition, the results of PT with adaptive smoothing are rather insensitive to the choice of cutoff wavenumber (as long as it is kept large enough compared to the wavenumber corresponding to transition towards nonlinearity), but the quality of the agreement with simulations depends, on the other hand, on the value of the parameter introduced in Sec. 4. Here, we adopt the intuitive choices for post-collapse PT and for Zel’dovich solution, but, as illustrated in Appendix F, it is possible to optimize the value of . For instance at , we find for the analyses performed in this section that and would provide better agreement with simulation measurements. However we also noticed that the optimal value of depends slightly on redshift as well as on initial conditions, this is why we decided for the analyses presented in this article to stick to our fiducial, slightly suboptimal choice.
5.5 Random initial conditions: power-law power spectrum
We now turn to random initial conditions with a power-law power spectrum:
[TABLE]
This allows us to test post-collapse PT for different kinds of initial conditions and to understand more in details small-scale dynamics. Assuming Einstein-de Sitter universe (, and ), we performed simulations for and , with the same resolution as our main CDM runs, i.e. , and Mpc*-1*. Although the initial conditions are scale free, we set an actual size for the simulation box, Mpc, which allows us to normalize the power spectrum using (see footnote 1). With such a normalization, the level of nonlinearity of at the end of the simulations is roughly the same as in the CDM runs studied in previous section. To measure power spectra with sufficient statistics, we performed sets of 50 simulations with different initial random seed for each value of .
Fig. 10 displays the power spectra (left panel) and phase-space diagrams (right panel) measured at in the scale-free simulations. The numerical results are compared to basic post-collapse PT (blue solid) and to Zel’dovich solution (green dotted) with adaptive smoothing. For reference, linear theory (black dashed) and the analytic prediction given by Eq. (82) for the Zel’dovich solution without smoothing (black solid) are also plotted in left panel.333While computing in Eq. (82), we introduced the actual cut-offs in the initial conditions of the simulations, i.e. below and above Mpc*-1*. In addition, predictions from stable clustering (Joyce & Sicard, 2011; Benhaiem et al., 2013) are displayed as a thin solid line and will be discussed more in details below. To complete visual inspection of phase space, Fig. 10 plots velocities and positions as functions of Lagrangian coordinate in the same region of phase space as in right panel of Fig. 10. Note that the simulations considered in right panel of Fig. 10 and in Fig. 10 have the same random initial seed as for right panel of Fig. 7 and Fig. 8, so one can compare directly phase-space features between scale-free and CDM-like simulations.
Examination of Figs. 10 and 10 demonstrates again the success of post-collapse PT with adaptive smoothing. While the results are nearly as good for as for the CDM-like case, post-collapse PT does slightly less well for . The same trend can be seen for Zel’dovich solution with adaptive smoothing. This result is not surprising, for two reasons. Firstly, the level of nonlinearity in the simulations is slightly larger than for and the CDM-like initial conditions. Indeed, deviation from linear theory or pure Zel’dovich solution happens at smaller values of for than for other cases. Secondly and more importantly, adaptive smoothing using a sharp low-pass filter is less efficient, from the dynamical point of view, when considering a power spectrum with a high slope. The results would look even worse in the case that we simulated but do not show here. To explain this, one can just take the example of the single mode system, e.g. initial conditions with a single sine wave as shown in Fig. 2. In this case, with a sharp- filter, the fluctuation at is either kept, hence post-collapse PT becomes invalid after the system reaches next-crossing time, or it is completely suppressed, which means that the procedure with adaptive smoothing is at the end equivalent to no smoothing at all. Increasing the power spectrum slope basically makes us approaching the single mode case with , since we always introduce a cut-off at some value to have a sufficiently smooth representation of the phase-space sheet at very small scales. Of course, this situation is unrealistic in practice, because no such problem occurs when considering CDM-like power spectra, but one has to stay aware of the fact that adaptive smoothing does not always work. A possible way to improve it could consist in changing the filtering window so that all the wavenumbers are affected by the softening procedure, even if it is performed predominantly at a given scale: to this respect, adaptive smoothing with a Gaussian window might represent an interesting alternative to our sharp- filter.
Another interesting result of our measurements is the excellent agreement between the simulations and the stable clustering prediction from Joyce & Sicard (2011) and Benhaiem et al. (2013) in the nonlinear regime444We also note here that the asymptotic slope seen in the simulations is rather different from the one predicted by the so-called adhesion model based on Burgers’ equation, which gives at very high-, irrespective of the initial power-law slope (Valageas, 2009a, b; Valageas & Bernardeau, 2011)., who predict that the slope of the power spectrum should be equal to
[TABLE]
In fact, even in the CDM-like case, which roughly corresponds, in one dimension, to an asymptotic value of equal to at large , stable clustering prediction works, since it predicts a plateau, , which is indeed what we observe on left panel of Fig. 7. Note that what Joyce & Sicard (2011) and Benhaiem et al. (2013) call “stable clustering” is very specific. In the 1D case considered here, it is equivalent to assuming that relaxed objects become of constant size in the coordinate instead the physical coordinate as normally supposed in the three-dimensional case (see, e.g. Davis & Peebles, 1977; Peebles, 1980). It is important to notice here that post-collapse PT seems to be able to reproduce correctly the measured power spectrum in the beginning of the stable clustering regime. This suggests that stable clustering regime is established only in a few dynamical times, since post-collapse PT is only able to follow the dynamics of a single halo until next-crossing time.
On the other hand, if one considers a “halo model” approach (see, e.g. Ma & Fry, 2000; Peacock & Smith, 2000; Scoccimarro et al., 2001; Cooray & Sheth, 2002, and references therein) and assumes that at large the power spectrum is dominated by the “one-halo term”, interesting properties are expected, at least when initial density fluctuations are sufficiently smooth at small scale. Indeed, one dimensional simulations of single cold halos suggest that their central density profile nearly behaves like a power-law
[TABLE]
with . The actual value of is controversial. While Binney (2004) postulated it to be equal to , detailed -body simulations suggest instead (Schulz et al., 2013). On the other hand, pure Vlasov simulations results using the waterbag method suggest (Colombi & Touma, 2014). Other values to keep in mind are from the singularity occurring exactly at crossing times and from the caustics (see, e.g. Shandarin & Zeldovich, 1989). These singularities imply, as discussed below, some expected behavior for at large (Gouda & Nakamura, 1989).
A slope for the density profile corresponds, after some trivial algebraic calculations, to a power-law slope
[TABLE]
for the power spectrum. For instance, one expects the Zel’dovich approximation to be dominated by the caustics at very small scales, i.e. constant at large , which is indeed the case on left panels of Figs. 7 and 10 in both analytic (black) and numerical (green) cases, irrespective of adaptive smoothing. Turning to post-collapse PT with adaptive smoothing, since, by using , we push the (smoothed) system exactly to next-crossing time, we are dominated by the singularity with . Hence, in this case one expects close to when becomes very large, as can be observed on all the figures. With a smaller value of , the system is dominated, similarly as for Zel’dovich solution, by the caustics at very small scales and then reaches a plateau at large as can be seen in Fig. 11 of appendix F. Without adaptive smoothing, one also intuitively expect post-collapse PT power spectrum to be dominated by caustics at very small scales, however the objects have a more complex structure than Zel’dovich solution as shown by lower right panel of Fig. 2, hence the plateau at large is not clearly visible on left panel of Figs. 7, 7 and 10 probably because not reached as quickly as for the pure Zel’dovich solution when increasing .
Turning to the actual nonlinear evolution of single objects, assuming as found in the literature and for the density profile of a halo would respectively give and [math] for the “one halo” power spectrum. This of course supposes that statistical averaging over many halos of different masses does not affect the slope of the power spectrum predicted by the single halo. Again, we have to be aware of the fact that at very large , one must be dominated by the caustics, since we introduce a cut-off at large . Such a cut-off enforces some level of smoothness of the curve representing the system in phase space and thus actually imposes the existence of well-defined caustic structures in the nonlinear regime,555Things would be more difficult to apprehend if the curve representing the phase-space sheet would be un-smooth. hence a plateau for at very large , even for the simulations. However this plateau is not visible in the simulations, because it lies outside the range of shown in the figures. On the other hand, there must be, according to our “halo model” prescription, a regime where might decrease () or present another plateau () or at least some kind of inflection at some intermediate values of . This is what can indeed be observed e.g. on lower insert of left panel of Fig. 7 for Mpc*-1* for the CDM-like case and on lower insert of left panel of Fig. 10 for Mpc*-1* for the simulation. There is no such signature obviously visible in the simulation, although the measurements are too noisy to conclude in this case. Interestingly, when the inflection is visible, its extension in space is correctly predicted by basic post-collapse PT with adaptive smoothing. This just shows that basic post-collapse PT is able to predict approximately the correct size for the halos, which we already knew.
Note that here, we did not try to position our single halo term discussion in the context of stable clustering predictions. These latter can be considered as the outcome of some statistical averaging after weighting the one-halo term by the mass function. They are also expected to be valid only in a finite range of values of if one assumes a cut-off at : what we mean here is that our “one-halo term” discussion is in fact not incompatible with the concept of stable clustering although the link between the halo model and the stable clustering predictions remain to be performed in details in the present case. Clearly, however, the small inflection discussed above marks a small deviation from stable clustering.
6 Conclusion
In this paper, we have developed a Lagrangian perturbation theory method for solving, beyond crossing time, Vlasov-Poisson equations in the expanding universe for initially cold systems in the infinite parallel planes geometric set up. The proposed approach captures post-collapse dynamics by computing at leading-order a counter term to the Zel’dovich solution just after first crossing time. It extends earlier work of Colombi (2015) to the cosmological case and, thanks to adaptive smoothing, to random initial conditions instead of a single halo. By performing a local Taylor expansion of the velocity and position of each mass element as functions of Lagrangian coordinate around each initial density peak, we are able to compute the force field in the multi-valued region. A correction to the single-stream flow motion, described by the Zel’dovich solution, is systematically obtained. The results are given as polynomial functions of Lagrangian position and the coefficient of each power is expressed in terms of quantities associated to density peaks. A simple analysis of the formation of a single cluster shows that post-collapse PT reproduces well the phase-space structure up to next-crossing time and can describe the projected density profile quite well even at later time.
To deal with the bulk properties of phase-space structures, we propose an adaptive filtering scheme. The idea is very close to the peak-patch treatment by Bond & Myers (1996). In particular, the local smoothing scale of initial conditions is determined by the calculation of a dynamical time, but in post-collapse PT this dynamical time is given by next-crossing time instead of collapse time. This allows us to produce a realistic coarse-grained version of the phase-space structure, even in the presence of mergers. After demonstrating that adaptive smoothing indeed works well in a simple setup with coalescing clusters, we applied it to more general, random initial conditions that include CDM-like as well as scale-free power spectra. The results of our analyses show that post-collapse PT predictions with adaptive smoothing reproduce the power spectrum of our numerical simulations remarkably well even at scales sufficiently nonlinear to probe the beginning of the stable clustering regime.
On important thing to keep in mind for further work is that the critical quantity in post-collapse PT is the time of next crossing, which can in fact be used to also create a spectacularly good pure Zel’dovich prescription with adaptive smoothing. While the calculations of post-collapse PT can be cumbersome, the expression (65) of this time is relatively simple and this can be very interesting when trying to generalize results of post-collapse PT to the three-dimensional case. When analyzing equation (65), one indeed notes that depends both on collapse time and on the height of the initial density peak. Since depends solely on , we can see that basically the dynamical time we consider for defining the local smoothing scale is a single function of , larger than . In fact, things are very simple: when performing adaptive smoothing one enforces to be equal to present time, or in other words, adaptive smoothing of initial conditions is performed so that the local peak height given by linear theory at present time, , is a number larger than unity (instead of exactly unity for the collapse time). For instance, in the Einstein-de Sitter case, and , we have
[TABLE]
with chosen in this article to be equal to unity for post-collapse PT and one-half for the Zel’dovich solution. Note that the optimal value of is smaller than these intuitive choices (see Appendix F).
The present post-collapse PT treatment together with the adaptive smoothing scheme seems a very promising tool. A successful extension to the 3D case should lead to a major breakthrough in the precision PT calculations of large-scale structure statistics beyond the single-stream approximation. Our 1D study is however only a first step. Toward a practical application to the 3D case, there are several issues to be addressed. One is the analytic calculation of statistical quantities such as the power spectrum. Indeed, in the present work, we measured directly the power spectrum from a particle distribution following the dynamics prescribed by the various PT prescriptions under scrutiny. In fact, analytic calculations of the power spectrum are tractable, at least in the presence of fixed smoothing. Adaptive smoothing itself can in principle be implemented in the analytic framework by inspiring for instance from the theory of excursion sets discussed in e.g. Bond et al. (1991), although we can foresee that actual calculations will probably be very cumbersome.
Another issue, when generalizing to the 3D case, is that the Zel’dovich solution no longer provides a sufficiently accurate description of the dynamics before collapse. At least third-order Lagrangian PT is required to estimate correctly collapse times as well as the power spectrum beyond tree-level approximation, which clearly adds some significant level of complexity. Furthermore, the variety and complexity of topological configurations of singularities appearing in the 3D case might in fact represent a nearly impassable obstacle in three dimensions (see, e.g. Hidding et al., 2014). We will tackle these issues in future work.
Finally, while we have focused in this article on the impact of post-collapse dynamics in 1D on the phase-space structure and on the power spectrum of the matter distribution, it would be interesting to see quantitatively how the small-scale modes responsible for the multi-stream flows couple to large-scale fluctuations described by the single-stream flows. This can be presumably clarified by measuring the response of the evolved power spectrum to a small initial perturbation of high- modes, through the comparison between a pair of particle distributions modeled with post-collapse PT that slightly differ in their initial conditions (Nishimichi et al., 2016). In the 3D case, a strong damping is found for the mode transfer from small to large scales, but the physical origin of it is still unclear. A detailed analysis in the 1D cosmological setup would provide important clues on the nature of mode transfers.
Acknowledgments
AT thanks Patrick Valageas for references and useful comments on the adhesion model. SC greatly acknowledges hospitality of Yukawa Institute for Theoretical Physics (YITP), where this work was initiated. We also acknowledge the support of YITP in organizing the workshop “Vlasov-Poisson: towards numerical methods without particles” in Kyoto, funded by grant YITP-T-15-02, ANR grant ANR-13-MONU-0003 and by Institut Lagrange de Paris (ANR-10-LABX-63 and ANR-11-IDEX-0004-02). This work was also supported in part by MEXT/JSPS KAKENHI Grant Numbers JP15H05889 and JP16H03977 (AT) as well as ANR grant ANR-13-MONU-0003 (SC).
Appendix A Derivation of integral expression for the potential
In this appendix, we derive the expression in Eq. (19). To do this, we first consider the 1D Green function which satisfies the following equation:
[TABLE]
With this Green function, the solution of the Poisson equation is formally expressed as
[TABLE]
To solve Eq. (88) under periodic boundaries condition, we move to harmonic space and express the Green function as:
[TABLE]
Recalling that the Fourier transform of Dirac’s delta function reads ,666Here we assume . substitution of Eq. (90) into Eq. (88) gives
[TABLE]
Thus, we obtain
[TABLE]
Note that there is a subtlety in the above equation. In deriving it, we ignore the contribution from mode in , while the mode does not vanish in Dirac’s delta function. The result remains however correct as long as the density field does not contain a mode, or equivalently, .
Keeping this point in mind, we can further simplify the expression. We use the following formula (see, e.g. Chap.1.44 of Gradshteyn & Ryzhik, 1980; Moriguchi et al., 1987):
[TABLE]
The final expression for the Green function is then given by
[TABLE]
Substituting this into Eq. (89), we obtain the expression in Eq. (19).
Appendix B Derivation of the expression for the force
In this Appendix, we derive the expression for the force in the multi-stream region, given in Eqs. (40), (47) and (48).
To start with, we divide the integral in Eq. (20) into several pieces according to Eq. (39):
[TABLE]
We shall below evaluate these integrals one by one.
Consider the first integral in right-hand side of Eq. (95). We evaluate it at position in Fig. 1. Using Eq. (16), the integral is rewritten
[TABLE]
Exploiting the geometric properties of the configuration in Fig. 1, substitution of the explicit form of (Eq. 28) into the above leads to
[TABLE]
Here we used the fact that . Recalling that and are related with each other through Eqs. (32) and (33), Eq. (97) can be recast as a single function of either or , which respectively provides the general expression valid for and . Omitting the subscript or in either cases, we finally obtain
[TABLE]
Next, we consider the second integral in Eq. (95). Rewriting the integrand in terms of the variable with the help of Eqs. (16) and (28), it is straightforward to obtain
[TABLE]
The remaining integrals to be computed, i.e., the third and fourth terms in Eq. (95), are performed over intervals outside the multi-valued region, which allows us to use Zel’dovich solution to evaluate the integrands. For the third integral in Eq. (95), we obtain
[TABLE]
Although this expression is exact, our perturbative treatment implicitly assumes that the multi-stream region is small and the quantities inside this region are written in polynomial forms of . We may expand the last two terms of Eq. (104) as
[TABLE]
where we used the fact that (see Eq. 27). Additionally, periodic boundaries imply . Therefore, Eq. (104) can be recast as
[TABLE]
In similar manner, we evaluate the fourth integral in Eq. (95) as follows:
[TABLE]
The first of these two integrals is performed exactly to give
[TABLE]
where we exploited properties of periodic boundaries, and , as well as momentum conservation, which implies . Turning to the second integral in Eq. (107), we use again the fact that the domain of integration is assumed to be narrow to Taylor expand the integrand using the formula
[TABLE]
to obtain the following expression, valid at third-order in :
[TABLE]
Here, Eqs. (29)–(31) are used to rewrite the expression in the last line. Summing up the results in Eqs. (108) and (109), we obtain, at the end,
[TABLE]
Plugging Eqs. (102), (103), (106) and (110) into Eq. (95), we finally obtain the expression of the force in the multi-valued region (Eqs. 40, 47 and 48).
Appendix C Derivation of basic post-collapse PT results
In this Appendix, we detail the calculations of basic post-collapse PT predictions summarized in Eqs. (58), (60), (62) and (64), together with the coefficients in Table 1 and 2.
C.1 Derivation of
Our starting point is Eq. (57), which involves two integrals. To evaluate the first integral, the approximate formula given by Eq. (164) is applied. Recalling that we are considering the outer part of the multi-valued region, , we have
[TABLE]
with . Here, we used the leading-order expression for critical time , , and ignore in Eq. (164) the higher-order terms involving the quantities or .
On the other hand, the second integral in Eq. (57) is analytically performed. Using the evolution equation in Eq. (22), we obtain
[TABLE]
Summing up the above two contributions, we obtain the expression for summarized in Eq. (58).
C.2 Derivation of
To derive the expression for , we just need to integrate over time Eq. (58), as written in Eq. (59):
[TABLE]
Making use of the expressions summarized in Table 1, we evaluate the first integral. A straightforward calculation leads to
[TABLE]
Note again that we used the leading-order expression for (see Eq. 37).
The second integral in Eq. (113) is analytically performed to give
[TABLE]
Collecting Eqs. (114) and (115), we obtain Eq. (60).
C.3 Derivation of
For , the expression to be evaluated is given by Eq. (61), which involves three integrals. The first integral is obtained from Eq. (111) by setting (Eq. 38). Recalling the fact that (Eq. 34), we have
[TABLE]
For the second integral in Eq. (61), we use the approximate formula given by Eq. (164) in the regime . Ignoring higher-order terms involving the quantities or , we obtain
[TABLE]
Finally, the last integral in Eq. (61) is evaluated using the formula given in Eq. (112). Collecting the three contributions above, we obtain the expression for (Eq. 62).
C.4 Derivation of
The starting expression for is given in Eq. (63), which involves two integrals. To evaluate the first integral in Eq. (63), we use Eq. (114). Setting , an expression at seventh-order in is obtained:
[TABLE]
To compute the second integral in Eq. (63), we substitute Eq. (62) into the integrand, but without the term . Using the explicit expression for the coefficients in Table 2, a straightforward calculation gives
[TABLE]
Combining the above two results, we obtain Eq. (64), which gives the final expression for along with the coefficients in Table 2.
Appendix D Higher-order corrections to post-collapse PT
In this Appendix, we detail calculations of the various higher-order corrections to post-collapse PT we propose in § 3.4.
One possible way to improve on basic post-collapse PT is to Taylor expand at next-to-leading order the expression for (as well as ) that determines the boundary of the post-collapse region (Eqs. 34 and 35):
[TABLE]
with
[TABLE]
Inverting the relation [], we obtain the corresponding critical times at next-to-leading order,
[TABLE]
and
[TABLE]
Using these expressions, we repeat below the same calculations as in Sec. 3.3. Higher-order expressions for the correction to Zel’dovich flow are derived on top of basic post-collapse PT results and are denoted by and . In Sec. 3.4, a semi-analytic treatment designed by spl is also discussed. In this case, we adopt the basic post-collapse PT results in the inner part of the multi-stream region, but the extension of this latter is calculated using the higher-order expression (120) for . Then, a straightforward third-order spline interpolation, that we do not detail here, is used to connect the inner part to the Zel’dovich solution, assuming again that the Lagrangian boundary of the multi-stream region is given by the higher-order expansion (120).
Another improvement may come from the calculation of the integrals in Eqs. (49) and (50). In Sec. 3.3, a part of the integrands is Taylor-expanded in time and the integration is performed for the leading-order terms. The resultant expressions include time-dependent terms up to for the velocity and for the position. Here, using the approximate integral formulae for the Taylor-expanded integrands in Appendix E, we keep terms up to for the velocity, for the position, and derive the corresponding higher-order corrections to the Zel’dovich flow, and . These corrections assume, of course, higher-order expansions for critical times and .
Below, we present the expressions for and . The expressions for and are obtained by simply setting and to zero.
D.1 Velocity and position in the outer part of the multi-stream region:
The expression for the higher-order correction to the velocity, , is
[TABLE]
with
[TABLE]
Here, we have introduced the following function:
[TABLE]
where is the value of the Hubble parameter at time . For the higher-order correction to the position, , we have
[TABLE]
with
[TABLE]
Note the term proportional to , which is necessary to enforce continuity of the solution.
D.2 velocity and position in the inner part of the multi-stream region:
The expression for the higher-order corrections to the velocity, , becomes
[TABLE]
with
[TABLE]
Finally, the higher-order corrections to the position, , are given by
[TABLE]
with
[TABLE]
Appendix E Approximate formulae for integrals
In this appendix, we provide useful formulae for primitive integrals over time, which are used to derive the post-collapse PT results in Appendices C and D.
Our aim here is to evaluate the time-integral of the force in the multi-valued region, given in Eq. (40). The time-integral of the second term, involving function , is performed exactly (see Eq. 112). On the other hand, the integral of the first term needs an approximate treatment which will be valid as long as the multi-valued region is sufficiently small, or equivalently, shortly after collapse time . Below, we evaluate the integral by Taylor-expanding the integrand. We first give the results term by term:
[TABLE]
where is given by Eq. (129). For the term with fractional power of time, we obtain
[TABLE]
where we used the leading-order approximation for , i.e., .
Summing up the above formulas, we obtain the primitive of the time integrals of the force involving function . The result including polynomials of time up to is
[TABLE]
Appendix F On the choice of parameter in adaptive smoothing
In this Appendix, we study the effect of varying the parameter on the performances of PT predictions with adaptive smoothing (see Sec. 4 for definition). In the main text, we adopt the intuitive setup for post-collapse PT and for Zel’dovich solution. However, PT predictions are prone, after shell-crossing, to get worse over time, so we can naively expect that choosing a smaller value of will provide a better agreement with simulations.
Fig. 11 shows the results obtained in the case of the CDM-like power spectrum studied in § 5.4 when varying the value of . Here, the calculations are performed using , and Mpc*-1* for a boxsize Mpc, which is a relatively low resolution set up, but this will not have any consequence on the discussion that follows. Left panel shows the power spectrum at and the resultant predictions for different values of are plotted as solid lines with different colors. As anticipated, decreasing results in an enhancement of small-scale power and a better agreement between post-collapse PT and the simulation is obtained. But a too small value of overshoots the simulation and an optimal choice of roughly corresponds to . The same trend can be also seen for the Zel’dovich solution, but the improvement of the power spectrum amplitude is rather mild: while seems to provide the best overall behavior for the Zel’dovich solution, the result is clearly still not as good as what post-collapse PT can provide, in particular in the regime Mpc*-1*.
Right panel of Fig. 11 shows the phase-space portraits and compares the simulation results (red) to predictions obtained with our “standard” value of for post-collapse PT (upper insert) and for Zel’dovich solution (lower insert) and to predictions obtained with respective nearly optimal values of and . It seems difficult to judge by eye whether the prediction with a smaller value of really improves the description at small scales, although more structures seems to be captured for post-collapse PT with , which contradict the results obtained for the power spectrum. However the simple examination of this figure corresponding to a single realization of the random initial conditions is not conclusive, obviously: it does not preclude the fact that a smaller value of can give a better result for the power spectrum obtained from averaging over many realizations.
Finally, to conclude this section, although we do not show the results here, we noticed as well that the best choice of could vary slightly according to redshift or initial condition: for example, for the CDM cosmology considered here, a better choice of at redshift is and respectively for post-collapse PT and Zel’dovich solution. This finally explains why we decided to keep, for simplicity, the generic values and respectively for post-collapse PT and Zel’dovich solution, even if they are sub-optimal.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Alard (2013) Alard C., 2013, MNRAS , 428, 340 · doi ↗
- 2Baldauf et al. (2015) Baldauf T., Mercolli L., Zaldarriaga M., 2015, Phys. Rev. D , 92, 123007 · doi ↗
- 3Baldauf et al. (2016) Baldauf T., Codis S., Desjacques V., Pichon C., 2016, MNRAS , 456, 3985 · doi ↗
- 4Baumann et al. (2012) Baumann D., Nicolis A., Senatore L., Zaldarriaga M., 2012, J. Cosmology Astropart. Phys. , 7, 51 · doi ↗
- 5Benhaiem et al. (2013) Benhaiem D., Joyce M., Sicard F., 2013, MNRAS , 429, 3423 · doi ↗
- 6Bernardeau & Valageas (2010) Bernardeau F., Valageas P., 2010, Phys. Rev. D , 81, 043516 · doi ↗
- 7Bernardeau et al. (2002) Bernardeau F., Colombi S., Gaztanaga E., Scoccimarro R., 2002, Phys. Rept. , 367, 1 · doi ↗
- 8Bernardeau et al. (2008) Bernardeau F., Crocce M., Scoccimarro R., 2008, Phys. Rev. D , 78, 103521 · doi ↗