Renormalization group and UV completion of cosmological perturbations:
Gravitational collapse as a critical phenomenon

Current and upcoming surveys of the cosmic large-scale structure are expected to test the cosmological concordance model at a precision of about one percent Beutler et al. (2011); T. Abott et al. (2016) (Dark Energy Survey Collaboration); H. Aihara et al. (2018) (HSC collaboration); P. A. Abell et al. (2009) (LSST Science Collaboration); R. Laureijs et al. (2011) (Euclid collaboration). Therefore, having fast and accurate theoretical predictions for the large-scale structure is of fundamental importance. For example, this is relevant in the context of field-level forward models for reconstruction of tracers of the matter distributions, such as retrieved from the galaxy density field Jasche and Wandelt (2013); Kitaura (2013); Wang et al. (2013); Ata et al. (2021) or from quasar spectra obtained from the Lyman-$\alpha$ forest data Kitaura et al. (2012); Metcalf et al. (2018); Porqueres et al. (2020); Ravoux et al. (2020). Another important application relates to pairing theoretical predictions with numerical simulations Tassev et al. (2013); Howlett et al. (2015); Feng et al. (2016); Chartier et al. (2021); Kokron et al. (2021); Zennaro et al. (2021); Aricò et al. (2022).

Cosmological perturbation theory (CPT; Peebles (1980); Fry (1984); Bernardeau et al. (2002); Blas et al. (2014); Rampf (2021)) provides highly accurate predictions on large cosmological scales, and in particular plays a crucial role in connecting early- with late-time cosmology, such as by providing initial conditions for numerical simulations Klypin and Shandarin (1983); Efstathiou et al. (1985); Scoccimarro (1998); Crocce et al. (2006); Michaux et al. (2021). However, CPT struggles to accurately predict the small-scale collapse to gravitationally bound structures, a process that is intimately tied to the shell-crossing singularity—the crossing of trajectories of collisionless matter, which comes with extreme matter densities.

To remedy the problem, many different approaches beyond standard CPT have been developed, for example by applying several renormalization or resummation techniques, as well as effective field-theory methods or path-integral formalisms (see e.g. Valageas (2004); Crocce and Scoccimarro (2006); McDonald (2007); Matarrese and Pietroni (2007, 2008); Pietroni (2008); Matsubara (2008); Carrasco et al. (2012); Carlson et al. (2013); Blas et al. (2013); Bartelmann et al. (2019)). However, these approaches are either of statistical nature and thus cannot be directly tested against deterministic solutions or, to our knowledge, such critical tests have not yet been performed (see however McQuinn and White (2016) for an exception). Other approaches circumvent the shortcomings of CPT, by combining CPT predictions with those of the spherical or ellipsoidal collapse model Kitaura and Hess (2013); Bernardeau (1994); Mohayaee et al. (2006); Monaco et al. (2002, 2013); Stein et al. (2019); Neyrinck (2016); Tosone et al. (2021) by performing a large-scale/small-scale split—see e.g. Ref. Monaco (2016) for a review, and Ref. Lippich et al. (2019) for performance tests of such hybrid methods for modeling galaxy clustering.

Still, these studies do not attempt to tackle the obvious problem, namely that CPT is highly inefficient for predicting the nonlinear collapse. We believe that this problem has been partially addressed by Refs. Tatekawa (2007); Yoshisato et al. (1998); Matsubara et al. (1998); Taruya et al. (2022), who tested Padé approximants or Shanks transforms to accelerate the convergence of CPT. As we will outline in this article, even a faster convergence acceleration can be achieved, provided one exploits the asymptotic structure of the underlying collapse problem.

In all these matters, the choice of coordinates can be crucial. Indeed, the first nontrivial shell-crossing solutions have been found using Lagrangian coordinates Rampf and Frisch (2017); Saga et al. (2018, 2022), specifically within the framework of Lagrangian perturbation theory (LPT). Reasons for LPT performing better in comparison to its Eulerian counterpart are various, here we just name two: First, Lagrangian approaches are very efficient in resolving advected motion, essentially since the involved nonlinearities are absorbed into the Lagrangian time derivative. Second, the Lagrangian map acts as a de-singularization transformation, meaning that the infinity in the matter density at shell crossing gets converted into a vanishing of the Jacobian determinant. Obviously, resolving the latter is technically easier than the former, especially within the context of perturbation theory.

Despite these advantages, however, even LPT converges extremely slowly for spherical or quasispherical collapse Munshi et al. (1994); Yoshisato et al. (1998); Rampf (2019); Saga et al. (2018, 2022). Another, somewhat surprising, limitation of LPT is related to the late-time evolution of voids, which shows clear signs of loss of convergence (e.g. Sahni and Shandarin (1996); Karakatsanis et al. (1997)). This divergent behavior has been analyzed by the concept of “mirror symmetry” Nadkarni-Ghosh and Chernoff (2011, 2013), which elucidates that the radius of convergence of the LPT series is identical for both over- and underdense regions, essentially since the radius of convergence is independent of the sign of the spatial curvature parameter.

In this article, we revisit in detail the convergence issues of LPT for spherical over- and underdensities, and provide a systematic analysis of two asymptotic techniques. We will see that the poor convergence of the LPT series can be remedied by employing a technique dubbed UV completion, and this twofold: First, the convergence of the UV-completed LPT series is vastly accelerated in overdense regions and, second, the underdense evolution does not display any divergent behavior anymore. At the core, the UV completion exploits the asymptotic structure of the LPT series at order infinity. The functional form of this asymptotic behavior appears to be quite generic Rampf and Hahn (2021), which raises the justified hope that this technique might be adaptable in the near future also to the case of gravitational collapse with cosmological random field initial conditions.

In addition to the UV completion, we also study an alternative approach by applying a renormalization group (RG) technique to the spherical-collapse problem. Historically, the RG approach was introduced in quantum electrodynamics to regularize infinities in quantum field theory. RG is also vital in statistical physics, for example when investigating critical phenomena in phase transitions Wilson and Fisher (1972), which is also an instance where singularities appear (e.g., in the heat capacity). RG techniques are also heavily employed in non-equilibrium physics to detect singularities (e.g. Goldenfeld (1992); Chen et al. (1994, 1996); DeVille et al. (2008); Kirkinis (2008)), which is particularly relevant for investigating critical phenomena in general fluids. Crucially, as we show here by means of the spherical-collapse problem, the mere knowledge of the underlying singular structure provides a superior theoretical model as obtained through conventional perturbative techniques.

This paper is organized as follows. In the following we provide the basic fluid equations, first for generic initial conditions (Sec. II.1) and then limited to spherical symmetry (Sec. II.2). Afterwards, in Sec. III, we briefly review LPT. In Sec. IV, we discuss a particularly simple implementation of a renormalization-group approach (see App. B for complementary derivations exploiting the RG-flow condition), including also the pairing of RG and Padé approximants (Sec. IV.4). Then, we provide an asymptotic analysis of the LPT series (Sec. V.1), whose findings are then implemented to achieve a UV completion of the LPT series (Sec. V.2). Subsequently, we apply our findings to the calculation of the nonlinear density (Sec. VI.1) and of its one-point probability distribution function (Sec. VI.2). Finally, we summarize our findings and provide concluding remarks in Sec. VII.

The central quantity in LPT is the displacement field

which in the case of spherical symmetry can be written in a simplified manner as $\bm{\psi}(\bm{q},t)=\bm{q}\,\psi(t)$, where the scalar displacement $\psi(t)$ depends only on time. The fastest-growing mode of this displacement is expanded as a perturbative series (e.g. Moutarde et al. (1991); Buchert and Ehlers (1993); Bouchet et al. (1995); Ehlers and Buchert (1997); Zheligovsky and Frisch (2014); Saga et al. (2018); Rampf (2019)),

where $\psi(t)=x(t)-1$, and the (scalar) comoving trajectory $x(t)$ is defined via $x_{i,j}=\delta_{ij}x$. Furthermore, $k$ is a free parameter fixed by the initial conditions which, physically, amounts to an effective curvature parametrizing the local departure from a spatially flat Universe; see e.g. Ref. Rampf (2019) for calculational details.

In LPT, the Ansatz (6) is used to solve Eq. (4) at subsequent perturbation orders $n$. The $\psi_{n}$’s occurring in Eq. (6) are easily determined by the following recursive relations ($n\geq 1$) Rampf (2019)

which in particular leads to the first few coefficients

See Sec. V for investigating the large-order asymptotic properties of the displacement series.

In Fig. 1 we show the comoving trajectory $x(a)=1+\psi(a)$ for several LPT truncation orders (various colored lines), and compare it against the exact parametric solution (black dotted line, based on App. A). Specifically, the positive $a$-time branch depicts the collapse of a spherically overdense region with curvature $k=3/10$, while the negative time branch reflects the evolution of an underdense region with $k\!=\!-3/10$ where, for convenience, we exploit the sign symmetry of the tuple $(a,-k)\leftrightarrow(-a,k)$ within the definition of $\psi(a)$. In the present case, LPT convergence is lost at the critical time $|a|=a_{\star}=(3\pi\sqrt{2})^{2/3}\simeq 5.622$, which coincides exactly with the time of spherical collapse. More generally, for given curvature $k$, the time of collapse is $a_{\star}=\delta_{\rm c}/k$, where $\delta_{\rm c}=(3/5)(3\pi/2)^{2/3}\simeq 1.686$ is the usual critical density at collapse time (see e.g. Peebles (1967); Tomita (1969); Gunn and Gott (1972); Bertschinger and Jain (1994)).

Moreover, since LPT is a Taylor series in powers of $a$, the temporal regime of convergence is spanned by $|a|<a_{\star}$, where the radius of convergence $a_{\star}$ is set by the nearest singularity(ies) in the complex-time plane around the expansion point $a=0$. In the case of spherical symmetry, the nearest singularity occurs at the real-valued collapse time $a_{\star}$ where, specifically, the velocity $v=\dot{x}$ blows up Rampf (2019, 2021). This explains the exact coincidence between the collapse time and the loss of convergence in the case of spherical symmetry.

The above argument also implies that the LPT series in the void case must have an identical radius of convergence of $a_{\star}$, which is also shown in Fig. 1 in the negative time branch (exploiting the aforementioned sign symmetry). Clearly, the lack of convergence in the void case beyond $a_{\star}$ is a mathematical artifact since, physically, the void underdensity should not be influenced by a singular velocity appearing in the spherical collapse. To our knowledge, the fact that the LPT series for over- and underdensities have an identical radius of convergence was first pointed out in Ref. Nadkarni-Ghosh and Chernoff (2011).

Multiplying Eq. (4) by $\dot{r}/r^{2}$ and integrating the resulting equation in time, one obtains

where $C$ is an integration constant. Now, changing from cosmic time to scale-factor time $a$ with $\partial_{t}=\dot{a}\partial_{a}$ and using the Friedmann equation $(\dot{a}/a)^{2}=8\pi G\bar{\rho}_{0}/(3a^{3})$ to get an expression for $\dot{a}$, we arrive after some straightforward calculations at our main evolution equation in the RG approach,

where a prime denotes a partial derivative w.r.t. $a$-time, and we have defined $\epsilon=3C/(4\pi G\bar{\rho}_{0})$. Actually, in the RG approach $\epsilon$ acts as a perturbative bookkeeping parameter; physically, it can be interpreted as a curvature perturbation locally induced through a spherical over- or underdense region in an otherwise spatially flat Universe (see also App. A). In the following, we seek solutions for (10) by means of a standard perturbative expansion, i.e.,

This perturbative Ansatz allows us to obtain a set of differential equations for $r_{0}$, $r_{1}$, etc. that are easily integrated in time. The resulting approximation for $r$ will have so-called secular terms that grow unboundedly for large times, which is an indication of ill-posed asymptotic behavior. To cure the perturbative results from these secular terms, we then apply suitable renormalization techniques.

In the following we provide a reduced RG approach that, for the present problem, is particularly simple in its implementation; see App. B for more traditional avenues exploiting the RG-flow equation, leading however to identical results as quoted here.

We have just seen above that RG techniques can be used to circumvent inherent limitations of the LPT series. Here we shall follow an orthogonal approach: Instead of avoiding LPT, we develop a framework that allows us to extract and exploit the asymptotic properties of the displacement series

thereby providing a theoretical description that we call UV completion [see Eq. (8) for the first few coefficients $\psi_{s}$].

The general idea of the UV completion is as follows [see Eq. (37) for the final result]. Suppose that there exists an explicit solution for the nonperturbative displacement field, valid deep in the ultraviolet (short wavelength) regime. The functional form of that nonperturbative displacement might not be known in general, but let us assume that we know at least some of its “intrinsic” properties, which are encapsulated in the singular term

Here, $\nu$ is neither zero nor a positive integer, and $\mathfrak{a}_{\star}$ denotes a certain time value when a singularity in the problem arises. For example, if $\nu<0$, then the displacement itself would “blow up” at $\mathfrak{a}=\mathfrak{a}_{\star}$ which of course would be unphysical. If instead $\nu>0$ but $\nu\notin\mathbb{N}$, then the displacement will remain bounded, but the $n$th derivatives with $n>\lfloor\nu\rfloor$ exhibit singular behavior at $\mathfrak{a}=\mathfrak{a}_{\star}$, which is an instance of non-analyticity ($\psi\notin C^{\infty}$, i.e., the solution is not infinitely differentiable, and therefore not globally representable by a Taylor series).

Having the above in mind, we conjecture that the UV-completed solution of the LPT displacement field (29) is then obtained by splitting off the non-analytic piece $\psi_{\infty}$ as follows,

Here, the first term on the r.h.s. is the truncated LPT series, while $\psi_{\text{$\infty$}}^{(\!n-1\!)}$ denotes the Taylor expansion of $\psi_{\text{$\infty$}}$ about $a=0$ to truncation order $n-1$. This last term is needed to avoid the double counting of certain low-order coefficients. We remark that such UV completing schemes are well known in the field of general fluid dynamics and asymptotic analysis; early related work can be found e.g. in Refs. Kuo (1953); Domb and Sykes (1957); van Dyke (1974).

Before going into the calculations, let us briefly summarize the steps needed to determine the unknowns $\nu$ and $\mathfrak{a}_{\star}$. First of all, the physical meaning of these unknowns can be obtained by considering the Taylor expansion of (30) about $a=0$, and comparing subsequent ratios of these Taylor coefficients with ratios of the LPT coefficients $\psi_{n}$. A straightforward analysis then reveals that $\mathfrak{a}_{\star}$ is nothing but the radius of convergence of the LPT series, if and only if the involved limit in d’Alembert’s ratio test

exists. Similar asymptotic considerations performed on the “graphical” level (Fig. 4) then lead to the conclusion that $\nu$ is related to the slope of the ratio of coefficients at order infinity.

In the present article, we have analyzed two asymptotic approaches that remedy these drawbacks of LPT. One of the avenues employs renormalization-group techniques, where the evolution equation is first recast such as to identify an effective expansion parameter—in the present case the rescaled spatial curvature parameter $\epsilon=10k/3$. The evolution equation is then solved with a naîve perturbation Ansatz in powers of $\epsilon$. Subsequently, the perturbation equations are renormalized in order to expose the secular terms, i.e., terms that would grow unboundedly for large times. After evaluating the RG-flow condition, which removes the arbitrariness of certain integration constants, one then obtains the final renormalized result. These steps have been effectively implemented in a particularly concise way in Sec. IV; we refer the interested reader to App. B for more elaborate technical details.

The analyzed RG technique predicts already in the first iteration a singular behavior of the solution with the correct critical exponent within the leading-order asymptotics (cf. App. A for complementary asymptotic considerations by means of the exact parametric result). Subsequent higher-order iterations within the RG method are particularly fruitful when paired with Padé approximants (Sec. IV.4). Overall, we find that the Padé approximant (2,2) of the 4RG prediction (dubbed 4RG+PD22) delivers the most accurate RG result. We find some evidence that Padé approximants are not just better approximations, but in fact are crucial to restore a symmetry property of the solution which, in the present case, would otherwise lead to unphysical collapse of void solutions.

The other asymptotic approach that we considered is the UV completion (Sec. V). For this we first analyzed the large-order asymptotic properties of the LPT series, for which we drew the so-called Domb–Sykes plot (Fig. 4). From this plot it is seen that subsequent ratios of LPT coefficients settle into a linear relationship at large orders, from which one can deduce that the radius of convergence of the LPT series is limited by a singularity at the (real-valued) collapse time. Even more, the graphical analysis reveals a measured critical exponent that closely resembles the one obtained from the RG method. With this input at hand, one can complete the LPT series in a highly efficient manner (eq. 37). Crucially, the UV-completed prediction does not come with the above mentioned pitfalls of LPT, since the remainder of the series is not a Taylor series anymore. We note however that the UV completion should be performed at sufficiently low truncation orders, otherwise unwanted features in the void evolution are “reactivated”. For both over- and underdense regions, we find that the UV completion with third-order truncation in LPT (dubbed 3UV) leads to the most convincing predictions.

Finally, we have tested our two asymptotic methods by means of predicting the nonlinear density evolution and corresponding one-point distribution (Sec. VI). Characteristically, we find that the UV method works marginally better for predicting the nonlinear density near the collapse (which comprises a huge challenge for LPT) in comparison to our flagship RG method, while the situation is opposite in voids; see Fig. 6. In contrast, when predicting the one-point probability distribution function of matter, we found only negligible differences between the tested UV and RG methods (Fig. 7).

The obvious challenge of the discussed asymptotic methods is their potential applicability for predicting the evolution of collisionless matter for cosmological initial conditions in three space dimensions. For the UV method, one first needs to find the precise nature of the convergence limiting singularity(ies) for the Lagrangian displacement field, which requires high-order LPT solutions that are however already available Zheligovsky and Frisch (2014); Matsubara (2015); Schmidt (2021). In this context, we note that in Ref. Rampf and Hahn (2021), we obtained some numerical evidence that the norm of the displacement contains a non-analytic term that is in structure formally identical to the one in the case of spherical collapse.

Regarding a UV completion within an Eulerian-coordinates formulation, it is likely that the leading-order asymptotic behavior of the fluid variables is characterized by a pair of complex-conjugated singularities in time; see Ref. Rampf et al. (2022) for highly related avenues applied to the inviscid Burgers equation. However, in such a case, the asymptotic structure of the fluid variables would be just $\propto(a_{\star}-a)^{\nu}+(\bar{a}_{\star}-a)^{\nu}$, where $a_{\star}$ is now complex and $\bar{a}_{\star}$ denotes its complex conjugate. Hence, the remainder of the respective series has a different structure as in the spherical case, but this can easily be handled by suitable alterations within the framework.

Similarly, the RG method needs to be suitably adapted when applied to cosmological initial conditions. First of all, the underlying fluid equations are partial differential equations in time and space. Therefore, as a preparatory step, the fluid equations could be first represented in a Fourier basis, which then leads to a spatially decoupled ODE in time for each Fourier mode of the fluid variable. The solutions to these ODEs could then be renormalized by the methods outlined in this paper, provided a suitable expansion parameter is identified. Regarding the latter, we remind the reader that our starting point for the RG method, Eq. (10), is obtained by time integrating the fluid equations for spherical collapse. This time integration then revealed as an integration constant the spatial curvature, which indeed acts as the expansion parameter in the present RG approach. Whether similar derivations hold for the Fourier coefficients of the fluid variables for random initial conditions remains to be investigated.

Future applications of the asymptotic methods include the accurate modeling of (non-spherical) void and overdense regions in deterministic or statistical contexts (e.g. excursion sets and data inference for field-level forward modelling), as well as for hybrid approaches where the methods are paired with numerical simulation- or machine learning techniques. Lastly, in this paper we did not consider post-shell-crossing effects Taruya and Colombi (2017); Pietroni (2018); Chen and Pietroni (2020); Rampf et al. (2021); Saga et al. (2022), but the RG and UV methods are in principle able to handle such critical behavior, which however requires further investigation. In the long term, asymptotic methods have the potential for reducing the gap between theoretical and numerical methods, while at the same time enhance the physical insight into highly nonlinear problems.