{'fraction_non_alphanumeric': 0.037225663048404155, 'fraction_numerical': 0.012257493687822352, 'mean_word_length': 4.557687694890861, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dark matter clustering: a simple renormalization group approach

I compute a renormalization group (RG) improvement to the standard beyond-linear-order Eulerian perturbation theory (PT) calculation of the power spectrum of large-scale density fluctuations in the Universe. At z=0, for a power spectrum matching current observations, lowest order RGPT appears to be as accurate as one can test using existing numerical simulation-calibrated fitting formulas out to at least k~=0.3 h/Mpc; although inaccuracy is guaranteed at some level by approximations in the calculation (which can be improved in the future). In contrast, standard PT breaks down virtually as soon as beyond-linear corrections become non-negligible, on scales even larger than k=0.1 h/Mpc. This extension in range of validity could substantially enhance the usefulness of PT for interpreting baryonic acoustic oscillation surveys aimed at probing dark energy, for example. I show that the predicted power spectrum converges at high k to a power law with index given by the fixed-point solution of the RG equation. I discuss many possible future directions for this line of work. The basic calculation of this paper should be easily understandable without any prior knowledge of RG methods, while a rich background of mathematical physics literature exists for the interested reader.


I. INTRODUCTION
Statistics of the large-scale density fluctuations in the Universe provide an invaluable source of information on the very early Universe, where the initial density perturbations that grow by gravitational instability were created, and on the material content and gravitational laws of the Universe, which determine the growth of structure and how it looks to us from a distance [1]. The standard practical toolbox we use for interpreting observations consists of linear theory on very large scales or at very early times [2], and everywhere else numerical simulations which compute the fully non-linear evolution of a realization of the density field in a patch of Universe. Fitting formulas, motivated by physical intuition but ultimately heuristic, are often used to more efficiently interpolate between simulation results [3,4].
Beyond linear order perturbation theory (PT) has been considered for a long time [5,6,7,8,9,10,11], but does not currently play an important role in our interpretation of precision observations (see [12] for a thorough review of largescale structure perturbation theory). The desire to measure dark energy properties using observations of baryonic acoustic oscillation features [13,14,15,16,17,18,19,20,21,22,23,24,25,26,27], which fall more or less exactly in the range of scales where higher order perturbation theory could be most useful, gives this line of work much more pressing practical relevance than it has had in the past [28]. More traditional uses of galaxy clustering to measure cosmological parameters, long reliant on the assumption that galaxy power is simply proportional to linear mass power, have also reached the point where weakly non-linear effects are a substantial limitation [29,30]. Additionally, other probes of large-scale structure, like the Lyα forest [31,32], weak lensing [33], galaxy cluster/Sunyaev-Zel'dovich effect (SZ) measurements [34], and possibly future 21cm surveys [35], could all benefit from improved computational techniques.
The problem with higher order perturbation theory is that it does not work very well in exactly the regime where it could currently be most useful: at low redshift and wavenumber k ∼ 0.1 h Mpc −1 , where galaxy clustering is significantly but still only weakly non-linear [36]. Recently, [37,38] proposed a renormalization technique aimed at improving PT calculations. They re-sum an infinite series of perturbation theory terms (Feynman diagrams) which determine the memory of perturbations to initial conditions as a function of scale (see also [39]). This method is apparently computationally difficult because they have not yet computed a power spectrum using it. The renormalization group (RG) method I present here is different, but surely related. There has been some work toward using RG techniques to model large-scale structure [40,41,42]. The methods used are different than mine, and have not yet seen much practical application for comparison with observations. My implementation of the RG is inspired more by its use for dealing with secular divergences in perturbative solutions of differential equations [43,44] than by the probably better known applications in high energy and statistical physics [45,46,47,48], although all of these uses are related. The basic idea is to identify a quantity in the problem that is not directly observable and modify it to remove breakdowns in the perturbation theory, leading to a wellbehaved renormalized observable. In this paper I renormalize the initial (linear theory) power spectrum to remove the relative (not actually infinite for realistic power spectra) divergence of the lowest order correction to linear theory. More detailed mathematical clarification and explanation of the meaning of the method of [43,44], in terms of the theory of envelopes, can be found in [49,50]. Even more detailed treatments can be found in [51,52]. These papers should be useful to the reader who is unsatisfied with my presentation, or wants to extend it. During the last several years, [53,54,55,56,57,58,59] have extended and generalized this dynamical renormalization group method to compute the real time evolution of expectation values of quantum fields.
The rest of the paper is as follows: In §II I derive a simple formula for the RG-improved power spectrum of mass density fluctuations. In §III I show numerical results and compare them to standard PT and fitting formulas calibrated using simulations. Finally, in §IV I discuss possible future directions for this kind of work.

II. CALCULATION
In this section I present the basic calculation. In §II A I discuss the time evolution equations for cold dark matter, which are the starting point for perturbation theory. In §II B I review standard Eulerian perturbation theory. Finally, in §II C I apply the renormalization group improvement.

A. Evolution equations
We are interested in the statistics of the mass density field, , where x is the comoving position and τ = dt/a is the conformal time, with a = 1/(1 + z) the expansion factor. In particular, I will compute the power spectrum, P (k, τ ), defined by For now I will assume only cold dark matter, i.e., collisionless particles with insignificant initial velocity dispersion. The exact evolution is is described by the Vlasov equation [5] ∂f ∂τ with where f (x, p, τ ) is the particle density at phase-space position (x, p), m is the particle mass (which plays no role in the final results), and p = a mv (v here is a particle's peculiar velocity, not to be confused with the mean peculiar velocity of all particles at some point in space, used everywhere else). Except when otherwise indicated ∇ = ∂/∂x. The density field is obtained by averaging the distribution function over momentum: and the bulk (mean) velocity and higher moments of the velocity distribution e.g., the dispersion of particle velocities around their bulk velocity, can be similarly obtained by multiplying the distribution function by any number of p's (e.g., one to obtain bulk velocity) before integrating over p. As discussed by [5], taking moments of the Vlasov equation with respect to momentum leads to a hierarchy of evolution equations for these quantities. Dropping the higher moments of the velocity distribution gives the usual hydrodynamic equations for density and bulk velocity: and where H = d ln a/dτ = aH with H the usual Hubble parameter.

B. Standard Eulerian perturbation theory
Perturbation theory consists of writing the density and velocity fields as a series of terms of at least formally increasing order of smallness, i.e., δ = δ 1 + δ 2 + δ 3 + .... The evolution equations are solved order-by-order, with lower order solutions appearing as sources in the higher order equations so that δ n is of order δ n 1 [12]. The power spectrum for Gaussian initial conditions is given by where no terms 3rd order in δ 1 appear because the expectation value of any term cubic in a Gaussian field is zero. At this point I assume an Einstein-de Sitter Universe for simplicity, and define where D(τ ) = δ 1 (τ )/δ initial is the linear theory growth factor. The Einstein-de Sitter assumption is needed to avoid more complicated time dependence of P 13 and P 22 , but the real Universe is of course not Einstein-de Sitter. Fortunately, this is not a significant problem because, to percent level accuracy [28,60,61,62,63,64,65,66], the effect of changing the background model can be included by simply using the correct linear growth factor in Eq. (8). I will sometimes refer to P 13 and P 22 as 2nd order, meaning in the initial power spectrum amplitude, not to be confused with the fact that they are 4th order in δ 1 and require calculating the evolution of δ to 3rd order, i.e., δ 3 .
[67] derived the following useful form of the equations for P 13 (k) and P 22 (k): and Note that, due to many cancellations, these terms are not as divergent as they might appear at first glance, their sum being convergent at high k for power law P 11 (k) with n = d ln P/d ln k < −1 and at low k when n > −3 [67]. Their sum is zero for n ≃ −1.4 [12], a fact that will have interesting consequences for our RG calculation.

C. Renormalization group improvement
The problem with the standard calculation outlined in §II B, which leads us to renormalization, is that the 2nd term in Eq. (8) diverges relative to the first (although it does not literally become infinite for realistic power spectra), at increasingly large scales (small k) as time progresses. This divergence is not exactly unphysical -we expect nonlinearities to become important during gravitational collapse -however, the basic premise of perturbation theory is violated when higher order terms become large. There is no a priori reason to expect the results to be accurate. I employ a renormalization group calculation to cure this divergence. I start in §II C 1 with a quick derivation motivated by the method of [43]. Then I discuss the mathematical background, physical interpretation, approximations, and limitations of the method in §II C 2.

Simple calculation
To simplify the presentation of the calculation, I rewrite Eq. (8) in a more compact, schematic form, where A ∝ D 2 (τ ),P (k, τ ) ≡ P (k, τ )/A, P L (k) is the initial condition power, and [P 2 L ](k) ≡ P 13 (k) + P 22 (k) is the higher order correction term, which is quadratic in P L (k), as given by Eqs. (9) and (10). I now follow the method described in [43] to deal with the problem of the 2nd order term becoming large (deferring a detailed explanation of the meaning of this method to §II C 2). I rewrite A = A − A ⋆ + A ⋆ , where A ⋆ is just an arbitrary constant, so that I then absorb the constant A ⋆ P 2 L (k) into the initial conditions, P L (k), defining the renormalized initial conditions P ⋆ (k, A ⋆ ), so thatP where I can freely replace P L with P ⋆ in the 2nd term because the change is formally 3rd order (in the initial power spectrum amplitude). The observable power spectrumP (k) obviously must be independent of the completely artificial constant A ⋆ , so I impose dP (k)/dA ⋆ = 0 and find where I have dropped the derivative of P 2 ⋆ (k) because it is higher order (this is seen to be self-consistent from Eq. 14 itself). This equation is easy enough to solve numerically, given P ⋆ (k, A ⋆ ) at some value of A ⋆ . Since A ⋆ is completely arbitrary, I will in the end choose A ⋆ = A, removing the 2nd term in Eq. (13) entirely. The result forP (k, τ ) is then simply P ⋆ (k, A ⋆ = A), obtained by solving Eq. (14). To reproduce the linear theory result at very early times, the initial conditions for the solution of Eq. (14) must be P ⋆ (k, A ⋆ = 0) = P L (k).

Explanation
The method of [43] and related papers was, by the authors own admission, somewhat lacking in rigorous motivation. This situation was clarified by others [49,50,51,52]; however, their explanations may be opaque. A very simple example should serve to explain and clarify the method. Consider the differential equatioṅ The exact solution is easy to obtain, but suppose one did not know that and solved Eq. (15) using perturbation theory for small initial δ, as in the cosmological situation. The solution to the linearized equation is δ 1 = g 1 e t . The equation for the first non-linear correction, δ 2 , is theṅ with solution The approach generally followed in the cosmological case is to set the term that is the solution to the homogeneous equation, g 2 e t , to zero, i.e., to assume that at asymptotically early times δ → δ 1 ; however, this is a choice, not necessary. The key observation here is that solving two differential equations has produced two parameters in the solution, g 1 and g 2 , but we only need one to satisfy the physical boundary conditions. The RG method exploits and fixes the otherwise arbitrary extra parameter of the perturbative solution (note that this approach gives exactly the same result as the approach of "splitting" the linear solution into two pieces, one of which is considered to be second order and used to cancel the non-linear correction). For the standard choice g 2 = 0, δ 2 inevitably becomes larger than δ 1 , invalidating the perturbation theory. Note, however, that using g 2 we can actually guarantee that at any one chosen time, t ⋆ , δ 2 (t ⋆ ) = 0, and very near t ⋆ , δ 2 ≪ δ 1 . This requires g 2 = −g 2 1 e t⋆ , which leads to full solution Clearly the value of g 1 must depend on the choice of t ⋆ , i.e., g 1 = g 1 (t ⋆ ), for a fixed physical situation (fixed full solution δ(t)); however, nothing would be accomplished by simply fixing g 1 (t ⋆ ) to match the initial conditions, because this generally requires using the perturbative solution outside its range of validity. The RG method is to determine g 1 (t ⋆ ) by enforcing the idea that the physical solution should be independent of the arbitrary time t ⋆ , while using the perturbative solution only for t infinitesimally close to t ⋆ , where it should be valid, i.e., or The exact solution to this equation is g 1 (t) = c/(1 − ce t ) where c is a constant which will allow matching of the physical initial conditions. Finally, referring back to Eq. (18), evaluated at t ⋆ = t, gives the final solution While such good fortune would generally not be expected, it happens that this is the exact solution to the original differential equation, Eq. (15). The intuitive picture one should have of this calculation is of stepping forward in time slightly using the perturbative solution centered at the present time, then taking the result and using it as the initial condition for a perturbative solution around the new time, followed by another step, and so on. [49] showed that this method gives the envelope equation for the set of solutions obtained by applying the perturbation theory near various times t ⋆ (the envelope is the curve which is tangent to each of the different solutions).
To make the analogy with original case of the power spectrum clear, note that changing variables to A = exp (t), A ⋆ = exp (t ⋆ ), andδ = δ/A in Eq. (18) produces an equation with the same form as Eq. (13). The rest of the derivation follows identically until, unfortunately, Eq. (14) is not analytically solvable because P (k) represents an infinite number of variables rather than one.
While the resulting equations look similar, the skeptical reader may still question the relation between the derivation of Eq. (14) for the power spectrum and Eq. (20) for a single variable governed by a differential equation, e.g., δ would be most naturally be identified with the cosmological density field, not the power spectrum. To complete the connection one must return to the original derivation of cosmological perturbation theory. At linear order δ 1 (k) = D g 1 (k) where g 1 (k) is the normalization and D = a is the growth factor (assuming an EdS Universe as usual), and I have dropped the decaying mode solution as usual (I will discuss this further below). At 2nd order where the D g 2 (k) term is set to zero in the standard calculation (see [68] for explicit forms of J below). Again, this equation is missing terms which grow like D −1/2 or slower (naively one would expect them to be unimportant). δ 3 is needed to compute the power spectrum, but computing the bispectrum at this point is informative: The last term here is the usual perturbative result. The first term is normally zero, for Gaussian initial conditions. The middle term is normally zero because g 2 (k) is chosen to be zero. One can look at Eq. (23) two ways: The easy way is to say that there is no breakdown in perturbation theory, where a higher order term becomes larger than a lower order term. The last term is the leading term, i.e., the bispectrum is fundamentally O δ 4 1 , and there is no need for the first two terms to enter this discussion at all. The second approach is to argue that all possible terms should be included in the calculation -the fact that the first term happens to be zero initially is a peculiarity of the initial conditions, not fundamental to the calculation. In this approach, the new middle term, since it is arbitrary, could be set in a way that would cancel the last term at some chosen D ⋆ . The first term becomes a function of D ⋆ , and requiring that the solution be independent of D ⋆ would produce an RG equation for g 1 (k 1 )g 1 (k 2 )g 1 (k 3 ) (D ⋆ ), just as above in the simple example. The second approach is clearly more thorough, however, in this paper I will take the easy route and declare the usual bispectrum term to be the leading order, with no need for renormalization, i.e., I will set g 2 (k) ≡ 0 as usual (we will see that this term is not needed to renormalize the power spectrum). While my approach may not give the best possible result, it is not fundamentally incorrect, at least not any more than standard perturbation theory is fundamentally incorrect. Renormalization is a tool, not an imperative, and renormalizing part of the calculation but not all of it should be better than nothing.
Moving on to δ 3 , and again including the homogeneous solution, gives L ](k) above, just a quick way of writing the usual O δ 4 1 perturbative correction terms (note that if I had promoted the leading order bispectrum to O(δ 3 1 ) there would be more new terms in the power spectrum). Following the simple example above, we know just what to do here: eliminate the O δ 4 1 term at D ⋆ by a specific choice of g 3 (k), i.e., which finally leads to where P 11 (k) must be a function of D ⋆ . After the same renamings, this is exactly the starting point for the original RG calculation above, Eq. (13). While the approach described here apparently takes care of the relative divergence of the corrections to the power spectrum of density fluctuations, there is a hidden approximation that guarantees some imperfection. In the usual calculation the initial conditions are specified by simply specifying the density-density power spectrum P δδ (k), because for growing mode initial conditions the power spectrum of velocity divergence fluctuations θ ≡ −H −1 ∇ · v, P θθ (k), and the cross-correlation between δ and θ, P δθ , are equal to P δδ (k). However, the non-linear corrections to these three spectra are not equal, i.e., decaying modes are created. To simultaneously cancel the corrections to all three spectra, decaying modes would need to be included in the solutions for δ, e.g., δ 1 = g 1 (k)D + d 1 (k)D −3/2 . They do not decay instantly, so the non-linear terms they produce can not be safely dropped. It is hard to estimate how large the net effect of these changes would be without actually doing the full calculation, although [39] showed that the decaying mode contribution is a modest, although not negligible, fraction of the non-linear evolution (in a somewhat different context). Extension of the calculation in this paper to properly account for decaying modes should not be too difficult. One simply needs to write down solutions for P δδ (k), P θθ (k), and P δθ (k), including the possibility of decaying modes in the initial conditions. This will produce terms similar in form to P 13 (k) and P 22 (k), with various changes in numerical factors. It will then be possible to simultaneously cancel all of the corrections to P δδ (k), P θθ (k), and P δθ (k), and obtain RG equations for each of them. If the basic principle behind this method is useful, doing the calculation more carefully should result in even better accuracy.
While the potential usefulness of this RG method appears to be clear, and versions of it have been successful in many previous applications [43,44,49,50,51,52,53,54,55,56,57,58,59,69,70], the exact conditions in which it will succeed or fail are not completely understood. However, the presence of an attractive fixed point (or manifold) is clearly a positive feature [49,50]. Eq. (14) has such a fixed point at n = −1.4, which appears to be attractive (I haven't proved this but it works for every power spectrum I've tried). In the asymptotic (long time) limit the arbitrarily many degrees of freedom in a general power spectrum are reduced to a single amplitude of this power law. The system is naturally driven towards a regime where the perturbation theory approximation works best, i.e., the 2nd order correction to the linear theory power spectrum vanishes. Furthermore, the correction to the bispectrum also vanishes for n = −1.4, so we are driven to a place where the Gaussianity of the field is well-preserved. On the other hand, in the simple example of Eq. (15) we obtained an exact solution despite the fact that there is no fixed point -in fact, the relative size of the non-linear term diverges with time. Ultimately, the best way to see whether this rearrangement of the perturbation series is useful or not is probably to simply compute higher order terms to see if the result converges better than standard PT or not.

III. RESULTS
In Figs. 1-4 I show the results obtained by solving Eq. (14). As discussed above, I am assuming an Einstein-de Sitter Universe for the background time evolution, to avoid even small time dependence of the P 2 L (k) term. I compare to the simulation-calibrated fitting formulas of [3], including both their implementation of the formula from Peacock and Dodds (hereafter PD96) [71] (note that they did not redetermine the parameters using their expanded set of simulations) and their newer HALOFIT formula. These formulas were not calibrated for models with baryonic acoustic oscillations, so I compare using their standard transfer function from [72]. I choose parameters to produce a power spectrum similar to recent best fits to observations [1]: σ 8 = 0.85, n = 0.96, and Γ = 0.15. I reiterate that this Γ is chosen to give the closest possible match to the best fit power spectrum from [1], i.e., it has nothing to do with the Einstein-de Sitter background evolution, and is less than Ω m h in the concordance model because of broad-band baryonic power suppression which is not otherwise included in the transfer function of [72]. Figure 1(a) shows essentially perfect (better than percent level) agreement at z = 0 between RGPT and the Peacock and Dodds fitting formula [71] for k 0.3 h Mpc −1 . This does not mean that RGPT is perfect -it is likely that the agreement with PD96 is partially a coincidence, because we have no reason to expect their formula itself to be accurate to this level. In fact, because of the approximation discussed above related to decaying modes, RGPT must inevitably have some error as well. I note here that [3]'s HALOFIT formula is not suitable for testing perturbation theory at low-k because the fit was actually constrained by PT calculations for k < 0.1 h Mpc −1 (R. E. Smith, private communication). The fitting formula is a smooth function, so it is not surprising that it does not agree perfectly with standard PT at k < 0.1 h Mpc −1 , i.e., if the PT power is generally greater than the simulation power, the region near k = 0.1 h Mpc −1 will be some average of the two, because a step function is not allowed. It is not surprising that this imperfect fit would be missed by [3], because the deviation is < 5%, i.e., smaller than their claimed accuracy.
Additionally, the k < 0.1 h Mpc −1 power was not constrained exclusively by PT, because the criterion for using scale free simulation results was different (R. E. Smith, private communication). In any case it is clear that RGPT is an improvement over standard PT. Figure 1(b) shows the comparison at z = 0 again, with rescaled axes, this time focusing on higher k, including deep in the non-linear regime. I note optimistically that, while the prediction is not precise at high k, it continues to be a dramatic improvement over linear theory, i.e., even when the non-linear power is a factor of ∼ 10 greater than the linear power, the perturbation theory result correctly accounts for most of the difference (admittedly, this was also true of standard PT). If similar improvements are obtained by removing approximations or simply with each additional order of perturbation theory, one could imagine converging quickly to a precise result. However, it is possible, or even likely, that missing higher moments of the velocity distribution will undermine the calculation at high k.
One of the notable features of the results is their convergence at high k to the d ln P/d ln k ≃ −1.4 fixed point of Eq. (14), where the 2nd order PT term vanishes for a pure power law power spectrum [11]. At z = 0 this limit is reached for k 0.3 h Mpc −1 , while at earlier times it is only reached at higher k (e.g., k 3 h Mpc −1 at z = 3). Qualitatively, this seems like a positive development. Whatever small-scale structure exists initially, it appears to be effectively taken out of play. This may alleviate one common worry about this kind of PT, that one is integrating over high-k modes that are highly non-linear, violating the premise of the perturbation theory. This fixed-point behavior provides another reason to hope that higher order calculations can significantly extend the range of accuracy of the calculation: it may be that the high k behavior of P (k) can be thought of as a slow rolling with scale, as higher order terms become important, of this kind of fixed-point power law. The fixed point will also be modified by the full inclusion of decaying modes. Figure 2 shows the comparison at z = 1, a typical redshift for planned galaxy redshift surveys aimed at measuring the baryonic acoustic oscillation feature [22]. As at z = 0, RGPT agrees well with the fitting formulas in the mildly non-linear regime, although now the agreement is better with HALOFIT than PD96 at the high k end of this regime. An optimistic interpretation of this figure would be that PD96 is more accurate at low k, because it was not constrained by standard PT, but HALOFIT becomes more accurate at higher k, as the authors [3] claim it should be. Dedicated numerical simulations and/or higher order PT calculations testing convergence are needed to determine the true power. Figure 2b again shows that the RGPT prediction traces the fitting formulas reasonably well deep into the highly non-linear regime. In this case, the agreement is actually as good as the agreement between the two fitting formulas, out to k ∼ 5 h Mpc −1 .
Moving on to z = 3, relevant to the Lyα forest [73,74,75], and again potentially future galaxy redshift surveys [22], we see clear deviation between the PT predictions and the fitting formulas, even in the weakly non-linear regime. This is especially puzzling considering that the RG and standard version of PT actually give very similar predictions in the weakly non-linear regime. This agreement means that the higher order terms that are being effectively re-summed by the RG calculation are not important. Considering the substantial disagreement between the two fitting formulas at very low k, I am not ready to conclude that PT is wrong. It may be that the fitting formulas were not well-calibrated in this regime. Focused simulations are needed to test this. Deeper in the non-linear regime at z = 3, Fig. 3(b) shows that the effect of renormalization becomes substantial, and that RGPT agrees quite well with HALOFIT, while both disagree significantly with PD96. The optimistic reading is that RGPT is very successful, but considering my willingness to dismiss HALOFIT in the weakly non-linear regime, it is probably best to wait for more simulations and/or higher order PT calculations before becoming too pleased with RGPT.
Finally, we consider z = 6. In the weakly non-linear regime shown in Fig. 4(a), RG and standard PT are essentially identical, and agree quite well with PD96. HALOFIT differs substantially, so again the results are ambiguous. The same can be said for the non-linear regime in Fig. 4(b), where the fitting formulas disagree substantially.

IV. DISCUSSION AND SPECULATION
The central result of this paper is Eq. (14), the RG equation for the renormalized power spectrum, demonstrated in Fig. 1(a) to give good results for the power spectrum in the quasi-linear regime at low redshift, where standard PT performs poorly. It is clear from Figs. 1-4 that the RG improvement to standard PT is generally helpful, in both the weakly and strongly non-linear regimes. Except in the most challenging case of z = 0 and high k, the simulationcalibrated fitting formula predictions of Peacock and Dodds [71] and HALOFIT [3] are sufficiently contradictory that it is not completely obvious that they are more reliable than RGPT.
It is possible that this method will turn out to be a computational curiosity, maybe leading to improved a posteriori understanding of what we see in simulations, but unable to achieve sufficient accuracy for practical problems, and thus leading to little fundamental change in how we carry out precision cosmology measurements. (It may appear that the method is guaranteed to be useful for describing weakly non-linear galaxy clustering, but this could still be foiled by bias and redshift-space distortions, which I discuss below.) On the other hand, I hope that this line of work is only in its infancy, and will lead to a diversification of computational methods, beyond standard N-body simulations. I am distinguishing here between systematic approximation methods and heuristic ones like the halo model, which rely heavily on calibration by simulations. One of the implicit goals of this work has been to avoid any approximations where the path to improved accuracy is unclear. A good quality for an approximation to have is a clear limit in which it can be expected to return to the exact result. Note that standard simulations are themselves exactly the kind of systematic approximation that I would like to see more of. Simulating a finite volume with a finite number of particles and limited force resolution are all approximations, but each can be tested for convergence as their controlling parameters are changed.
With the hope of stimulating further work, I now discuss a variety of possible directions that could be pursued, although I can't guarantee that all of these are good ideas: The first priority should be to include the velocity-velocity and velocity-density power spectra as independent functions to be renormalized, as I discussed extensively in §II C 2. The main requirement for doing this properly is a straightforward generalization of the standard PT results to allow for general initial conditions, including decaying modes. The n ≃ −1.4 fixed-point power spectrum will be modified, hopefully to a more interesting, accurate fixedpoint shape.
After that, brute force computation to higher order in perturbation theory should produce more accurate results. The fact that the current calculation is a substantial improvement over linear theory at all k (as opposed to diverging wildly at some point), well beyond the scale where it is no longer very precise, makes me hopeful that each additional order could extend the effective scale significantly. The possibility of estimating the errors in the calculation by comparing results at different orders is an equally strong motivation for going to higher order. [In the same sense, the present calculation should show where linear theory breaks down much more accurately than the common practice of guessing based on the value of ∆ 2 (k) = k 3 P L (k)/2π 2 .] These calculations give clear motivation and targets for high accuracy simulation work. RGPT, being a first principles method containing the possibility of internal error control, may in places lead the comparison rather than simply being tested itself. If a higher order of PT is computed, and the results agree between different orders to some level of precision, it is reasonable to expect that the result is correct, independent of simulations. (Missing higher moments of the particle velocity distribution function may be a loophole in this argument, as I discuss below.) With the fixed-point concept in mind, it may be interesting to look more carefully at the evolution of the logarithmic slope of the power spectrum with time and scale in simulations.
Abandoning the first principles philosophy, RGPT results may provide a useful template for a new simulation fitting formula. Because the errors in RGPT are not too large, and should be slowly varying functions of k and cosmological parameters, the new formula could be just a very simple parameterization of corrections to RGPT computed using simulations. One can obviously apply the techniques of this paper to higher order statistics such as the bispectrum [76,77,78], trispectrum, etc. [79]. As discussed in §II C 2, one can consider promoting the leading order for these statistics to the more generally possible O(δ 3 1 ) and O(δ 4 1 ) (for the bispectrum and trispectrum, respectively). Bispectrum results can be compared to the simulation fitting formula of [80]. If the trispectrum computation is sufficiently accurate, RGPT should be useful for computing expected statistical errors on the power spectrum, which is notoriously painful to do with simulations [81]. Similarly, errors on higher order statistics may be computable.
Mostly for convenience, I have assumed that the only effect of deviations from Einstein-de Sitter background evolution is to change the linear growth factor. This is probably sufficient in the weakly non-linear regime, but if RGPT is to achieve much accuracy at higher k, this assumption inevitably must be relaxed. This is clear, even though the corrections in standard PT are very small, because the dependence of the power spectrum on the equation of state of dark energy, w, at fixed observation-time linear theory power [82], can not be reproduced under the approximation used here. There is simply no way the dependence quantified in [82] can enter the current calculation. Fortunately, even small changes in the right hand side of Eq. (14) can lead to significant changes in the solution of this non-linear equation. It will also be necessary to move beyond the very simple time dependence in Eq. (8) if one wants to include any non-gravitational effects that break this form (e.g., gas pressure).
To describe direct tracers of gas, like the Lyα forest or SZ, gas dynamics could be included in the calculation. Even for weak lensing, which traces mass directly, gas physics has been shown to matter at a level problematic for high precision experiments [83]. Very rudimentary pressure approximations should be relatively straightforward to include [84]. While it may be that very sophisticated implementations of heating and cooling are hopeless, one should remember that there is some non-trivial inclusion of non-linearity in this calculation, so this is at least worth considering.
To describe clustering of galaxies, or other unavoidably inexact tracers of density, some concept of bias will need to be included in the calculation. I think this may be the area where RGPT is most likely to provide a real fundamental improvement over simulations, or at least complement to them. With a large but probably achievable amount of computer power, one can simulate the power spectrum of dark matter and its dependence on a relatively small number of cosmological parameters more or less perfectly. This will not be possible in the foreseeable future for galaxies, so robust, high precision cosmological measurements using them will rely on developing a comprehensive understanding of how relatively microscopic galaxy formation details translate into large-scale clustering patterns. While the separation of scales probably is not clear enough to make the analogy perfect, this is reminiscent of the classic uses of the renormalization group to understand how complex microscopic theories lead to simple macroscopic behavior [46]. The simplest way to introduce bias is to follow [68,85] in writing the galaxy density fluctuation, δ g , as a Taylor series in the mass density, δ g = i b i δ i /i!, which leads to an expression for the PT power spectrum of galaxies involving a few more terms than the mass power spectrum, but not fundamentally more complicated. This path has been followed by [86]. Much more ambitiously, one could write local formation rate equations for a galaxy density field and apply perturbation theory to them (something like this was discussed using linear theory in [87]). For example, ρ g = ρ f (ρ), withρ = −ρ g to represent depletion of gas, and with both components obeying the usual gravitational evolution equations. The mean galaxy density as a function of time would be an output of the calculation as well as fluctuations. The goal here would not be to write down a realistic local model and compute exact predictions as much as to see what kind of relations between observables can be generically expected. An intriguing possibility is that one could show that, regardless of the small-scale galaxy formation details, bias can only take certain universal forms, described by a small number of free parameters. Note that the validity of the scale-independent bias assumption has never really been proven even in the perfectly linear regime (see [88] for one attempt). This entire discussion of bias of galaxies applies equally well to many other traces of density, e.g., the Lyα forest [89].
For many applications, redshift space distortions will need to be included. Like bias, these should be relatively straightforward to include in the present RG formulation using the extension of [90]'s linear theory calculation described by [68]. [36] discusses the imperfections of the [90] approach, and proposes alternative ideas which may be useful.
It may be interesting, as a computational curiosity, to see if RGPT can give sensible results for n ≥ −1 power law initial conditions. In this case the integrals giving the first correction in standard PT truly diverge, in the sense of being infinite rather than just large [11]. The integrals can be cut off at wavenumber k c , but the results then depend on k c , which is arbitrary. One may be able to show that, using the RG as described in this paper, the results will converge as k c is taken to infinity. The small-A ⋆ integration of Eq. (14) will need to be done increasingly carefully as k c is increased, but hopefully convergence to the n = −1.4 fixed point will quickly erase all trace of the cutoff. Note that k c could affect the normalization of the asymptotic power law, which would spoil this idea (e.g., it may only work over some limited range of n).
One potential limitation in RGPT as presented here is the absence of higher moments of the velocity distribution function, e.g., the dispersion in the velocities of particles around their mean velocity at a point. I wrote the Vlasov equation (2) instead of cutting straight to the hydrodynamic equations to point out this problem, shared by standard PT (subsequently discussed by [91]). Dropping higher moments of the velocity distribution, also known as the singlestream approximation, is not really an added approximation in standard Eulerian perturbation theory. This is surely known to experts (e.g., [92]), but may not be universally recognized. One can easily include these moments as variables in the calculation, for example multiplying the Vlasov equation by p i p j and integrating over p gives (after some substitution to remove bulk velocity terms) where ∂ k = ∂/∂x k , σ ij (x, τ ) ≡ δv i δv j , q ijk (x, τ ) ≡ δv i δv j δv k , with δv i the deviation of a particle's velocity from the local mean velocity (note that a σ ij term would appear in Eq. 6 if I had not dropped it). Unfortunately, the lowest order evolution equations contain only the Hubble drag term, so any initial velocity dispersion will selfconsistently disappear from the perturbation theory. Furthermore, higher order equations contain no source terms, meaning that if the velocity dispersion is initially zero it can never become non-zero. In fact, [93] showed that even perturbation theory using the distribution function directly leads to the same result. This is a vexing problem because stream crossing is clearly ubiquitous in the real Universe; however, [39] estimated that it only becomes significant at k 1 h Mpc −1 , so it may not be fatal for intermediate scales. It seems likely that, if these variables really matter, there should be some way to recover them using a RG method. The basic idea is that, no matter how small the effect is in the bare perturbation theory, if it is dynamically relevant it will grow through a RG equation like Eq. (14). Note that even WIMP cold dark matter has a very small seed velocity dispersion [94]. The solution is probably to simply retain the dispersion variable in the calculation, even though it is decaying, in the same way I have argued above that one should retain the decaying combination of δ and θ -because it can be rapidly fed by higher order terms. Vorticity, ∇ × v, falls out of the standard PT calculation for the same reason velocity dispersion does [12]. The standard calculation makes use only of the velocity divergence, ∇ · v. A method that successfully reintroduces velocity dispersion may also work for vorticity.
I started this line of work with the intention of using a significantly different type of renormalization group method, based on integrating out (averaging over) small-scale (high-k) modes while modifying the equations describing evolution of the large-scale modes in a way that preserves their statistics [47]. This method usually starts with a path integral formulation of the problem. [93] presented a path integral approach to large-scale structure using the full distribution function, and it should not be too hard to write something similar for the usual moments of the distribution function. Velocity dispersion should arise inevitably in this approach, because smoothing a dispersionless field obviously produces dispersion. It is less inevitable that this dispersion will actually fundamentally change the outcome of the calculation, because the effect of these velocities is already present in the standard PT calculation. Among other possibly interesting features, the noise discussed in previous RG approaches to large scale structure [40,41] may be generated naturally using this approach (i.e., if high-k modes are eliminated, the time evolution of the remaining modes will no longer be perfectly deterministic).
Further work could elucidate the connection between this work and the re-summation method of [37,38]. Related renormalization group procedures, starting from general relativistic perturbation theory, are used by [95] and [96] to compute the back-reaction of inhomogeneities on the global evolution of the Universe.
Finally, the success of the RGPT calculation in this paper, which leads to what looks like a time evolution equation for the power spectrum, suggests that it might be useful to consider the exact time evolution equation for the power spectrum,Ṗ (k) ∝ 2 Re δ k δ −k , whereδ k is obtained from the usual evolution equations. The usual argument against this would be that δ k δ −k involves higher order terms like the bispectrum, which then require their own time evolution equations, and so on, leading to an infinite hierarchy of equations. Firstly, this could be a useful way to look at perturbation theory, especially if we are going to be renormalizing statistics like the power spectrum instead of the field itself. Second, one may be able to make progress by simply solving these equations numerically, in the same sense that simulations are usually used. One could, for example, truncate the hierarchy by replacing the connected spectra at some level by the appropriate functions of lower order spectra given by perturbation theory, and look for convergence in the results as higher order terms are added. The advantage of working with statistics of the density field rather than the density field itself is enormous. One obtains the desired results directly rather than averaging over fluctuations in simulations, where there is usually an unavoidable conflict between the need for large boxes to limit finite volume effects like sample variance, and the need for high force and mass resolution to properly compute small-scale structure. The fact that we are working with continuous, relatively slowing varying functions, rather than wildly fluctuating fields, allows many fewer points to be computed, e.g., we can compute a power spectrum over many decades of dynamical range by interpolating between a few tens of logarithmically spaced computation points when millions or billions of particles would be needed to cover the same range with a simulation. Tasks that at first glance seem arduous, e.g., parameterizing the trispectrum or even higher spectra and setting up their evolution equations, should be considered in contrast to the decades of person-power expended on simulations.