© The Authors 2024

1. Introduction

We have long since known that the large-scale distribution of matter throughout the Universe is well-approximated by a filamentary structure dubbed the cosmic web (de Lapparent et al. 1986). From a theoretical point of view, the first building blocks for a description of the cosmic web date back to the seminal work of Zel’dovich (1970) and many others that followed (e.g., Arnold et al. 1982; Klypin & Shandarin 1983). Indeed, the so-called Zeldovich approximation, which describes the first-order ballistic trajectories of particles in a Lagrangian description – predicts the existence of ‘pancakes’ (i.e., sheet-like tenuous walls), filaments, and clusters due to the collapse of anisotropic primordial fluctuations through gravitational instabilities in our expanding Universe. In the 90s, Bond et al. (1996) and Klypin & Shandarin (1983) showed that this ensemble of objects forms a connected network called the cosmic web based on the correlations imprinted in the primordial fluctuations. Initial peaks led to the formation of clusters at the nodes of the web while initial correlation bridges in between later form filaments that lie within walls, which themselves surround nearly empty void regions. Such a web classification can be achieved through the eigenvalues of the linear deformation tensor, which states that if all eigenvalues are negative (or below a threshold that is taken here to be zero) then the region is expanding in 3D, thus describing a void region; if all values are positive then the region is contracting, thus describing a knot, with the other two configurations of eigenvalues leading to either walls or filaments. It is possible to classify the cosmic web by means of its tidal field (Hessian of the gravitational potential). Pioneering work in this direction was presented by Doroshkevich (1970), and later by van de Weygaert & Bertschinger (1996), Rossi (2012), Desjacques et al. (2018), Castorina et al. (2016), Feldbrugge et al. (2018).

Another classification method was later introduced by Hahn et al. (2007a), Aragón-Calvo et al. (2007a) and Forero-Romero et al. (2009) using tidal fields, whereby either the Hessian of the gravitational potential is evaluated theoretically in the linear regime, or the non-linear potential is estimated in dark-matter-only numerical simulations. In subsequent years, several authors developed different classification schemes, to improve the detection of these cosmic web structures in various types of data (continuous fields, simulated datasets, point-like galaxy surveys, etc). For instance, Aragón-Calvo et al. (2010a) used the SpineWeb topological framework to segment the density field, Shandarin (2011), Falck et al. (2012) and Abel et al. (2012) showed how to classify morphological structures using the Lagrangian phase space sheet to count for shell crossings, and various authors have applied techniques from (continuous and then discrete) Morse theory (Colombi et al. 2000) to identify topological structures in the cosmological density field (Sousbie et al. 2008; Aragón-Calvo et al. 2007a; Sousbie et al. 2009; Sousbie 2011). Alternatively, Hoffman et al. (2012) introduced the V-web classification scheme, this time based on the shear of the velocity field, showing that it is able to better resolve smaller structures, which in turn allows for the study of finer dark-matter halo properties. An extension of this idea to Lagrangian settings was later proposed by Fisher et al. (2016). We note that the velocity divergence and density fields are closely related (e.g. through the mass-conservation equation in the Vlasov-Poisson system) and their statistics are virtually equivalent in the linear regime of structure formation. We refer readers to Wang et al. (2014) and Wang & Szalay (2014) for a thorough investigation of the differences between the dynamical and kinematical classifications.

Beyond the strict motivation of wanting to describe the cosmic web from a mathematical point of view, the classification schemes allow us to explore many different environmental effects on the properties of dark-matter halos (Hahn et al. 2007a,b; Aragón-Calvo et al. 2007b, 2010b; Codis et al. 2012; Libeskind et al. 2013; Hellwing et al. 2021) and galaxies within them (Nuza et al. 2014; Metuki et al. 2015; Poudel et al. 2017; Kraljic et al. 2018; Codis et al. 2018a; Hasan et al. 2023). Cosmic web classification can also be used to discriminate between cosmological models, as shown for instance by Lee & Park (2009) or Biswas et al. (2010) using voids; Codis et al. (2013) with counts of cosmic web critical points; Feldbrugge et al. (2019) with Betti numbers; Codis et al. (2018b) using the connectivity of the filaments; Bonnaire et al. (2022) relying on the power spectrum of the various cosmic web environments; and Dome et al. (2023) with cosmic web abundances. Understanding how this cosmic web evolves with time and scale is therefore paramount for both cosmology and studies of the galaxy formation that takes place within this large-scale environment. This evolution has been studied theoretically in some contexts (e.g., the local skeleton; Pogosyan et al. 2009b) with semi-analytical approaches (e.g., Fard et al. 2019) and in various numerical works. One example of the latter is the recent work by Cui et al. (2017, 2019), who used simulations to investigate the abundances of the various environments and their time evolution, relying on the T- and V-web decompositions. However, to the best of our knowledge, no theoretical model in the quasi-linear regime, and based on first principles, has been explicitly derived so far for these cosmic web abundances. This objective is nevertheless within the reach of standard techniques used in large-scale structure theoretical studies.

In the present paper, we therefore propose to derive theoretical predictions for web classifications based on the T-web definition. We first focus on a Gaussian description of the cosmic density field before turning to mild non-Gaussian corrections.

To assess the validity regime of our theoretical formalism and predictions for the one-point statistics of all four cosmic environments, we compare our results with measurements made in a dark-matter-only N-body numerical simulation, the Quijote simulation (Villaescusa-Navarro et al. 2020). The simulated box used by these latter authors has a size of 1 Gpc h⁻¹, with 1024³ dark matter particles. We are using their high-resolution fiducial cosmology snapshots at five different redshifts: z = 0, 0.5, 1, 2, and 3. Their fiducial cosmology is {Ω_m,  Ω_b,  h,  n_s, σ₈} = {0.3175,  0.049,  0.6711,  0.9624,  0.834}. Finally, in order to maximise the statistical relevance of our analysis, we typically average our measurements over 11 realisations of the Quijote suite and use those realisations to estimate error bars.

The outline of the paper is as follows. We first describe the T-web and V-web classifications in Sect. 2. In Sect. 3, we then present the theoretical formalism and the predictions obtained for the abundance of the different environments (as a function of threshold, redshift, and smoothing scale) in the linear regime of structure formation, that is assuming a Gaussian matter density field. Section 4 introduces a Gram-Charlier expansion to include non-Gaussian (i.e., non-linear) corrections to the joint probability distribution function (joint-PDF) of the elements of the Hessian of the gravitational potential, and we finally discuss our results and draw conclusions in Sect. 5. In the remainder of this paper, we will use ‘abundance’ or ‘probability’ equivalently.

2. T-web classification of the cosmic web

Among many possibilities (Libeskind et al. 2017; Leclercq et al. 2016), one commonly used mathematical way to classify the different cosmological environments is based on the number of eigenvalues of a deformation tensor above a threshold. Several definitions can be used for this deformation tensor (strain tensor, mass tensor, etc.) but the classification can always follow the outline below:

• 0 eigenvalue above the threshold: void,

• 1 eigenvalue above the threshold: wall,

• 2 eigenvalues above the threshold: filament,

• 3 eigenvalues above the threshold: knot.

Notably, the T-web classification developed by Forero-Romero et al. (2009) is based on the tidal shear tensor T, which is the tensor of the second derivative of the gravitational potential

$$\begin{array}{c}\hfill T_{\mathit{ij}}=\frac{{\partial }^2\mathrm{\Phi }}{\partial r_i{r}_j},\end{array}$$(1)

where Φ is the gravitational potential and r_i with i = 1, 2, 3 are the spatial coordinates. Alternatively, the V-web classification developed by Hoffman et al. (2012) is based on the velocity shear tensor Σ:

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathit{ij}}=-\frac{1}{2}(\frac{\partial v_i}{\partial r_j}+\frac{\partial v_j}{\partial r_i})/H_0,\end{array}$$(2)

where H₀ is the Hubble constant and v the velocity. We note that both of the above classifications are similar to the one used by Hahn et al. (2007a). Given those definitions, in practice, one can start from field maps, compute the deformation tensor, diagonalise it, and obtain the number of eigenvalues above a given threshold Λ_th at each spatial position, thus dividing the cosmic web into four characteristic classes.

Usually, the value of this threshold is not fixed a priori in the literature and is often chosen in order to obtain satisfactory visual agreement between the log density field and the classification. Originally, Hahn et al. (2007a) took Λ_th = 0 to distinguish structures with inner or outer flows. Follow-up studies demonstrated that having a non-zero threshold leads to a better visual match (for example, Forero-Romero et al. 2009; Libeskind et al. 2017; Suárez-Pérez et al. 2021; Cui et al. 2017; Hoffman et al. 2012; Carlesi et al. 2014). The threshold is generally taken as a positive constant of between 0 and 1, and usually does not depend on the redshift or smoothing scale, but we note that a positive threshold will tend to underline the most extreme knots and filaments (which is usually what is meant by good visual agreement) while a negative threshold will then underline the most empty voids.

In the remainder of this paper, we focus on the T-web formalism, although we emphasise that the V-web description would yield identical results at first perturbative order in a Lagrangian framework (the so-called 1LPT or Zeldovitch approximation). This is indeed due to the usual expression for these ballistic trajectories where the velocity is given by the gradient of the gravitational potential up to a uniform time-dependent factor (Zel’dovich 1970). As an example of the accuracy of this classification scheme, Fig. 1 presents a comparison of the density contrast (first and third columns, with a continuous colour bar) at different redshifts and smoothing scales, with the classification obtained using the T-web (second and fourth column, with a discrete color bar) at the same redshifts and smoothing scales. In the second and fourth columns, voids are shown in dark blue, walls in blue, filaments in green, and nodes in red. As expected from the T-web classification, we obtain the environments of the simulation by computing the second derivative of the potential and looking at the number of eigenvalues above our chosen threshold. Here, we are using a threshold of Λ_th = 0.01 for every redshift and smoothing scale based on the value used in Cui et al. (2017), which was chosen in order to have good visual agreement. For instance, the knots and voids (respectively in red and dark blue) can easily be identified by eye in both visualisations and are at the same position as rare maxima and minima in the density contrast maps.

Fig. 1. Logarithm of the density contrast of one slice of 1 Gpc h−1 in height and width and 2 Mpc h−1 in thickness of the Quijote simulation (first and third column) compared with the classification obtained using the T-web (second and fourth column) at different redshifts and smoothing scales. For the T-web classification, the colour bar is the number of eigenvalues above our chosen threshold Λth = 0.01 chosen specifically for good visual agreement.

To build a theoretical description of the cosmic web as defined by the T-web description, let us first define the variance of the contrast of the density field, δ = (ρ − ⟨ρ⟩)/⟨ρ⟩,

$$\begin{array}{c}\hfill {\sigma }_0^2=⟨{\delta }^2⟩.\end{array}$$(3)

Following Pogosyan et al. (2009a), we choose to normalise the derivatives of the gravitational potential by their variance, such that

$$\begin{array}{c}\hfill \nu =\frac{1}{{\sigma }_0}\delta ,{\varphi }_{\mathit{ij}}=\frac{1}{{\sigma }_0}{\mathrm{\nabla }}_i{\mathrm{\nabla }}_j\mathrm{\Phi }.\end{array}$$(4)

Formally, given our classification method, the probability of each cosmic environment then depends on the joint probability distribution 𝒫 of the eigenvalues of the tidal tensor. Given our normalisation choice in Eq. (4), here we use the eigenvalues normalised by their variance (or equivalently the eigenvalues of the normalised tidal tensor), which we denote λ₁ ≤ λ₂ ≤ λ₃. The probabilities of the four environments can then be written as

$$\begin{array}{cc}& P_{\mathrm{void}}={\displaystyle \int \mathrm{d}}{\lambda }_1\mathrm{d}{\lambda }_2\mathrm{d}{\lambda }_3\mathcal{P}({\lambda }_1,{\lambda }_2,{\lambda }_3)Boole({\lambda }_1<{\lambda }_2<{\lambda }_3<{\lambda }_{\mathrm{th}})\hfill \end{array}$$(5)

$$\begin{array}{cc}& P_{\mathrm{wall}}={\displaystyle \int \mathrm{d}}{\lambda }_1\mathrm{d}{\lambda }_2\mathrm{d}{\lambda }_3\mathcal{P}({\lambda }_1,{\lambda }_2,{\lambda }_3)Boole({\lambda }_1<{\lambda }_2<{\lambda }_{\mathrm{th}}<{\lambda }_3),\hfill \end{array}$$(6)

$$\begin{array}{cc}& P_{\mathrm{filament}}={\displaystyle \int \mathrm{d}}{\lambda }_1\mathrm{d}{\lambda }_2\mathrm{d}{\lambda }_3\mathcal{P}({\lambda }_1,{\lambda }_2,{\lambda }_3)Boole({\lambda }_1<{\lambda }_{\mathrm{th}}<{\lambda }_2<{\lambda }_3),\hfill \end{array}$$(7)

$$\begin{array}{cc}& P_{\mathrm{knot}}={\displaystyle \int \mathrm{d}}{\lambda }_1\mathrm{d}{\lambda }_2\mathrm{d}{\lambda }_3\mathcal{P}({\lambda }_1,{\lambda }_2,{\lambda }_3)Boole({\lambda }_{\mathrm{th}}<{\lambda }_1<{\lambda }_2<{\lambda }_3),\hfill \end{array}$$(8)

where {λ_i}_{i = 1, 2, 3} are the orderly eigenvalues of the (ϕ_ij)_{1 ≤ i, j ≤ 3} matrix and Boole is a Boolean equal to 1 if the condition is satisfied, and 0 otherwise.

As we are now working with normalised eigenvalues, and keeping our earlier choice of a Λ_th = 0.01 value, our threshold will be given by λ_th = (Λ_th = 0.01)/σ(z), with

$$\begin{array}{c}\hfill {\sigma }^2(z)=4\pi {\displaystyle \int \mathrm{d}}k{k}^{2}P(k,z)W{(kR)}^2,\end{array}$$(9)

where W(kR) is the applied smoothing, which in this paper is Gaussian with a smoothing scale R such that

$$\begin{array}{c}\hfill W_{\mathrm{G}}(kR)=exp(-\frac{1}{2}{k}^2{R}^2),\end{array}$$(10)

and P(k, z) is the matter density power spectrum. We denote the linear variance ${\sigma }_0^2$ when using the linear power spectrum, which is computed using the Boltzmann code CAMB (Lewis et al. 2000), and ${\sigma }_{\text{NL}}^2$ the non-linear variance (e.g., measured in the simulation or predicted with emulators for example).

3. Cosmic web abundances in the linear regime

To obtain some theoretical information about the abundance of the different cosmic environments in the density field and their redshift evolution, let us first consider the density field at linear order in Eulerian perturbation theory, this is, let us assume that it is Gaussian. For a Gaussian density field, the associated gravitational potential and its successive derivatives are also Gaussian-distributed, which means that we can use the Doroshkevich formula (Doroshkevich 1970), which gives the joint probability distribution of eigenvalues of a Gaussian symmetric matrix:

$$\begin{array}{c}\hfill {\mathcal{P}}_D({\lambda }_1,{\lambda }_2,{\lambda }_3)=\frac{675\sqrt{5}{e}^{\frac{3}{4}{({\lambda }_1+{\lambda }_2+{\lambda }_3)}^2-\frac{15}{4}({\lambda }_1^2+{\lambda }_2^3+{\lambda }_3^2)}}{8\pi }({\lambda }_3-{\lambda }_2)({\lambda }_2-{\lambda }_1)({\lambda }_3-{\lambda }_1).\end{array}$$(11)

We note that the general expression for the multi-dimensional PDF of the tidal tensor in Cartesian coordinates X = {ϕ_ij} in the Gaussian case reads

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathrm{G}}(X)={(2\pi )}^{-N/2}{|C|}^{-1/2}exp(-\frac{1}{2}X{C}^{-1}X),\end{array}$$(12)

where C is its covariance matrix.

From Eq. (11), we can then numerically integrate Eqs. (5)–(8) to obtain the probabilities of the various environments. In practice, those 3D integrals can be reduced to 1D as two degrees of freedom can be analytically integrated out. For that purpose, one needs to express the three degrees of freedom of the tidal tensor that are rotationally invariant as polynomials of its Cartesian coordinates ϕ_ij (in contrast to eigenvalues, which are not rotationally invariant). Following for instance Pogosyan et al. (2009a), a natural option can be the three coefficients of its characteristic polynomial:

$$\begin{array}{cc}& I_1=Tr({\varphi }_{\mathit{ij}})={\varphi }_{11}+{\varphi }_{22}+{\varphi }_{33}={\lambda }_1+{\lambda }_2+{\lambda }_3=\nu ,\hfill \end{array}$$(13)

$$\begin{array}{cc}& I_2={\varphi }_{11}{\varphi }_{22}+{\varphi }_{22}{\varphi }_{33}+{\varphi }_{11}{\varphi }_{33}-{\varphi }_{12}^2-{\varphi }_{23}^2-{\varphi }_{13}^2={\lambda }_1{\lambda }_2+{\lambda }_2{\lambda }_3+{\lambda }_1{\lambda }_3,\hfill \end{array}$$(14)

$$\begin{array}{cc}& I_3=\det |{\varphi }_{\mathit{ij}}|={\varphi }_{11}{\varphi }_{22}{\varphi }_{33}+2{\varphi }_{12}{\varphi }_{23}{\varphi }_{13}-{\varphi }_{11}{\varphi }_{23}^2-{\varphi }_{22}{\varphi }_{13}^2-{\varphi }_{33}{\varphi }_{12}^2={\lambda }_1{\lambda }_2{\lambda }_3.\hfill \end{array}$$(15)

To get variables that are as uncorrelated as possible, a simplification proposed by Pogosyan et al. (2009a) is to use combinations of the {I_k}_{1 ≥ k ≥ 3} – denoted {J_k}_{1 ≥ k ≥ 3} – such that the first variable is again the trace of the tidal tensor but higher-order variables only depend on the traceless part of that tensor and are thus uncorrelated with the trace in the Gaussian case. After some algebra, one can easily find

$$\begin{array}{c}\hfill J_1=I_1,J_2=I_1^2-3I_2,J_3=I_1^3-\frac{9}{2}{I}_1{I}_2+\frac{27}{2}{I}_3.\end{array}$$(16)

For a Gaussian random field, one can easily show that the probability distribution function of the above-defined variables reads (Pogosyan et al. 2009a)

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathrm{G}}(J_1,J_2,J_3)=\frac{25\sqrt{5}}{12\pi }exp[-\frac{1}{2}{J}_1^2-\frac{5}{2}{J}_2].\end{array}$$(17)

where J₂ ≥ 0 and J₃ is uniformly distributed between its boundaries such that $-J_2^{3/2}\le J_3\le J_2^{3/2}$. Let us note that the probability distribution of Js can be easily mapped to the Doroskevitch formula for the distribution of the eigenvalues after taking into account the usual Vandermonde determinant.

To get the abundances for the different environments, we use a criterion on the number of eigenvalues superior to the threshold. As we are now working with variables (J₁, J₂, J₃), we have to translate the criterion on eigenvalues into this new set of variables. Following Pogosyan et al. (2009a) and Gay et al. (2012), we can obtain the integration limits for the four different environments as follows:

$$\begin{array}{cc}& P_{\mathrm{void}}={\displaystyle {\int }_0^{\infty }}\mathrm{d}{J}_2{\displaystyle {\int }_{-\infty }^{-2\sqrt{J_2}+3{\lambda }_{\mathrm{th}}}}\mathrm{d}{J}_1{\displaystyle {\int }_{-J_2^{3/2}}^{J_2^{3/2}}}\mathrm{d}{J}_3{\mathcal{P}}_{\mathrm{G}}(J_1,J_2,J_3)+{\displaystyle {\int }_0^{\infty }}\mathrm{d}{J}_2{\displaystyle {\int }_{-2\sqrt{J_2}+3{\lambda }_{\mathrm{th}}}^{-\sqrt{J_2}+3{\lambda }_{\mathrm{th}}}}\mathrm{d}{J}_1{\displaystyle {\int }_{-J_2^{3/2}}^{-\frac{1}{2}{(J_1-3{\lambda }_{\mathrm{th}})}^3+\frac{3}{2}(J_1-3{\lambda }_{\mathrm{th}})J_2}}\mathrm{d}{J}_3{\mathcal{P}}_{\mathrm{G}}(J_1,J_2,J_3),\hfill \end{array}$$(18)

$$\begin{array}{cc}& P_{\mathrm{wall}}={\displaystyle {\int }_0^{\infty }}\mathrm{d}{J}_2{\displaystyle {\int }_{-2\sqrt{J_2}+3{\lambda }_{\mathrm{th}}}^{\sqrt{J_2}+3{\lambda }_{\mathrm{th}}}}\mathrm{d}{J}_1{\displaystyle {\int }_{-\frac{1}{2}{(J_1-3{\lambda }_{\mathrm{th}})}^3+\frac{3}{2}(J_1-3{\lambda }_{\mathrm{th}})J_2}^{J_2^{3/2}}}\mathrm{d}{J}_3{\mathcal{P}}_{\mathrm{G}}(J_1,J_2,J_3),\hfill \end{array}$$(19)

$$\begin{array}{cc}& P_{\mathrm{filament}}={\displaystyle {\int }_0^{\infty }}\mathrm{d}{J}_2{\displaystyle {\int }_{-\sqrt{J_2}+3{\lambda }_{\mathrm{th}}}^{2\sqrt{J_2}+3{\lambda }_{\mathrm{th}}}}\mathrm{d}{J}_1{\displaystyle {\int }_{-J_2^{3/2}}^{-\frac{1}{2}{(J_1-3{\lambda }_{\mathrm{th}})}^3+\frac{3}{2}(J_1-3{\lambda }_{\mathrm{th}})J_2}}\mathrm{d}{J}_3{\mathcal{P}}_{\mathrm{G}}(J_1,J_2,J_3),\hfill \end{array}$$(20)

$$\begin{array}{cc}& P_{\mathrm{knot}}={\displaystyle {\int }_0^{\infty }}\mathrm{d}{J}_2{\displaystyle {\int }_{\sqrt{J_2}+3{\lambda }_{\mathrm{th}}}^{2\sqrt{J_2}+3{\lambda }_{\mathrm{th}}}}\mathrm{d}{J}_1{\displaystyle {\int }_{-\frac{1}{2}{(J_1-3{\lambda }_{\mathrm{th}})}^3+\frac{3}{2}(J_1-3{\lambda }_{\mathrm{th}})J_2}^{J_2^{3/2}}\mathrm{d}{J}_3{\mathcal{P}}_{\mathrm{G}}(J_1,J_2,J_3)+{\displaystyle {\int }_0^{\infty }}\mathrm{d}{J}_2{\int }_{2\sqrt{J_2}+3{\lambda }_{\mathrm{th}}}^{\infty }\mathrm{d}{J}_1{\int }_{-J_2^{3/2}}^{J_2^{3/2}}}\mathrm{d}{J}_3{\mathcal{P}}_{\mathrm{G}}(J_1,J_2,J_3).\hfill \end{array}$$(21)

The analytical integration over J₁ and J₃ can now be performed and we thus obtain the following probabilities of the four environments (written in terms of 1D integrals over J₂):

$$\begin{array}{cc}& P_{\mathrm{void}}={\displaystyle {\int }_0^{\infty }}\frac{25\sqrt{5}}{48\pi }{e}^{-\frac{5J_2}{2}}[-\sqrt{2\pi }(2J_2^{3/2}-9J_2{\lambda }_{\mathrm{th}}+9(3{{\lambda }_{\mathrm{th}}}^3+{\lambda }_{\mathrm{th}}))\mathrm{erf}\left(\frac{3{\lambda }_{\mathrm{th}}-2\sqrt{J_2}}{\sqrt{2}}\right)+\sqrt{2\pi }(2J_2^{3/2}-9J_2{\lambda }_{\mathrm{th}}+9(3{{\lambda }_{\mathrm{th}}}^3+{\lambda }_{\mathrm{th}}))\mathrm{erf}\left(\frac{3{\lambda }_{\mathrm{th}}-\sqrt{J_2}}{\sqrt{2}}\right)\hfill \\ \hfill & +4\sqrt{2\pi }{J}_2^{3/2}\mathrm{erfc}\left(\frac{2\sqrt{J_2}-3{\lambda }_{\mathrm{th}}}{\sqrt{2}}\right)-2e^{-\frac{1}{2}{(2\sqrt{J_2}-3{\lambda }_{\mathrm{th}})}^2}(6\sqrt{J_2}{\lambda }_{\mathrm{th}}+J_2+9{{\lambda }_{\mathrm{th}}}^2+2)+e^{-\frac{1}{2}{(\sqrt{J_2}-3{\lambda }_{\mathrm{th}})}^2}(6\sqrt{J_2}{\lambda }_{\mathrm{th}}-4J_2+18{{\lambda }_{\mathrm{th}}}^2+4)]\mathrm{d}{J}_2,\hfill \end{array}$$(22)

$$\begin{array}{cc}& P_{\mathrm{wall}}={\displaystyle {\int }_0^{\infty }}\frac{25\sqrt{5}}{48\pi }{e}^{-\frac{5J_2}{2}}[-\sqrt{2\pi }(2J_2^{3/2}+9J_2{\lambda }_{\mathrm{th}}-9(3{{\lambda }_{\mathrm{th}}}^3+{\lambda }_{\mathrm{th}}))\mathrm{erf}\left(\frac{3{\lambda }_{\mathrm{th}}-2\sqrt{J_2}}{\sqrt{2}}\right)+\sqrt{2\pi }(2J_2^{3/2}+9J_2{\lambda }_{\mathrm{th}}-9(3{{\lambda }_{\mathrm{th}}}^3+{\lambda }_{\mathrm{th}}))\mathrm{erf}\left(\frac{\sqrt{J_2}+3{\lambda }_{\mathrm{th}}}{\sqrt{2}}\right)\hfill \\ \hfill & +2e^{-\frac{1}{2}{(2\sqrt{J_2}-3{\lambda }_{\mathrm{th}})}^2}(6\sqrt{J_2}{\lambda }_{\mathrm{th}}+J_2+9{{\lambda }_{\mathrm{th}}}^2+2)+2e^{-\frac{1}{2}{(\sqrt{J_2}+3{\lambda }_{\mathrm{th}})}^2}(3\sqrt{J_2}{\lambda }_{\mathrm{th}}+2J_2-9{{\lambda }_{\mathrm{th}}}^2-2)]\mathrm{d}{J}_2,\hfill \end{array}$$(23)

$$\begin{array}{cc}& P_{\mathrm{filament}}={\displaystyle {\int }_0^{\infty }}-\frac{25\sqrt{5}}{48\pi }{e}^{-\frac{5J_2}{2}}[-\sqrt{2\pi }(2J_2^{3/2}-9J_2{\lambda }_{\mathrm{th}}+9(3{{\lambda }_{\mathrm{th}}}^3+{\lambda }_{\mathrm{th}}))\mathrm{erf}\left(\frac{2\sqrt{J_2}+3{\lambda }_{\mathrm{th}}}{\sqrt{2}}\right)+\sqrt{2\pi }(2J_2^{3/2}-9J_2{\lambda }_{\mathrm{th}}+9(3{{\lambda }_{\mathrm{th}}}^3+{\lambda }_{\mathrm{th}}))\mathrm{erf}\left(\frac{3{\lambda }_{\mathrm{th}}-\sqrt{J_2}}{\sqrt{2}}\right)\hfill \\ \hfill & -2e^{-\frac{1}{2}{(2\sqrt{J_2}+3{\lambda }_{\mathrm{th}})}^2}(-6\sqrt{J_2}{\lambda }_{\mathrm{th}}+J_2+9{{\lambda }_{\mathrm{th}}}^2+2)+e^{-\frac{1}{2}{(\sqrt{J_2}-3{\lambda }_{\mathrm{th}})}^2}(6\sqrt{J_2}{\lambda }_{\mathrm{th}}-4J_2+18{{\lambda }_{\mathrm{th}}}^2+4)]\mathrm{d}{J}_2,\hfill \end{array}$$(24)

$$\begin{array}{cc}& P_{\mathrm{knot}}={\displaystyle {\int }_0^{\infty }}\frac{25\sqrt{5}}{48\pi }{e}^{-\frac{5J_2}{2}}[-\sqrt{2\pi }(2J_2^{3/2}+9J_2{\lambda }_{\mathrm{th}}-9(3{{\lambda }_{\mathrm{th}}}^3+{\lambda }_{\mathrm{th}}))\mathrm{erf}\left(\frac{\sqrt{J_2}+3{\lambda }_{\mathrm{th}}}{\sqrt{2}}\right)+\sqrt{2\pi }(2J_2^{3/2}+9J_2{\lambda }_{\mathrm{th}}-9(3{{\lambda }_{\mathrm{th}}}^3+{\lambda }_{\mathrm{th}}))\mathrm{erf}\left(\frac{2\sqrt{J_2}+3{\lambda }_{\mathrm{th}}}{\sqrt{2}}\right)\hfill \\ \hfill & +4\sqrt{2\pi }{J}_2^{3/2}\mathrm{erfc}\left(\frac{2\sqrt{J_2}+3{\lambda }_{\mathrm{th}}}{\sqrt{2}}\right)-2e^{-\frac{1}{2}{(2\sqrt{J_2}+3{\lambda }_{\mathrm{th}})}^2}(-6\sqrt{J_2}{\lambda }_{\mathrm{th}}+J_2+9{{\lambda }_{\mathrm{th}}}^2+2)+e^{-\frac{1}{2}{(\sqrt{J_2}+3{\lambda }_{\mathrm{th}})}^2}(-6\sqrt{J_2}{\lambda }_{\mathrm{th}}-4J_2+18{{\lambda }_{\mathrm{th}}}^2+4)]\mathrm{d}{J}_2,\hfill \end{array}$$(25)

where erf(z) is the error function $\mathrm{erf}(z)=2{\displaystyle {\int }_0^z}{e}^{-t^2}\mathrm{d}t/\sqrt{\pi }$ and erfc(z) is the complementary error function erfc(z) = 1 − erf(z).

It is then possible to perform numerical integrations and obtain the abundance of each environment. This is illustrated in Fig. 2, where the probabilities of the different environments are shown as a function of redshift. Each panel corresponds to a different Gaussian smoothing scale and displays the probability of voids, walls, filaments, and knots in blue, green, yellow, and red, respectively. The dots are the measurements obtained from the simulation, the dashed lines represent the Gaussian prediction obtained with the formalism described in this section, and the solid lines – which can be ignored for now – are the predictions at next-to-leading order obtained with a Gram-Charlier expansion, which is described in Sect. 4, below. As expected, we observe that the higher the redshift and/or the larger the smoothing scale – and thus the closer the simulation gets to the linear regime –, the closer the Gaussian prediction to the measurements. At lower redshift and smaller smoothing scales, non-Gaussianities are more important and departures from the Gaussian prediction thus appear. Here, all probabilities are computed using a threshold λ_th = 0.01/σ_NL inspired from the literature (using simulations) for which the evolution of σ with redshift and smoothing scale is given in Table 1. We note that the variance used solely in the normalisation of the threshold is the non-linear variance measured from the simulation. This allows us to have a meaningful threshold – in terms of rarity – even though we describe the cosmic structures in linear theory. Given the good agreement with the simulation already, that is, simply with Gaussian theory for sufficiently large scale and redshift, this states that the probabilities of the different environments are roughly captured by the statistics of a Gaussian field, at least for redshifts and scales that typically correspond to typical variances of σ ≲ 0.1. Consequently, the redshift evolutions seen in previous works (Cui et al. 2017, 2019) with a fixed threshold can be roughly understood simply as the non-linear evolution of the amplitude of fluctuations for that threshold. This is because a fixed threshold in non-linear densities does not correspond to a fixed ‘rarity’ or ‘abundance’ threshold for which the cosmic evolution would be much less important. This interpretation is notably valid at sufficiently large scale and redshift. Beyond σ ∼ 0.1, such a Gaussian field approximation starts to break down. In the following section, we show how the accuracy of the theoretical model with the help of Gram-Charlier corrections.

Fig. 2. Probabilities of voids, walls, filaments, and knots as a function of the redshift. These probabilities are shown in blue, green, orange, and red, respectively. Each panel is obtained at a given smoothing scale. Dots are measurements from the simulation, dashed lines are the Gaussian prediction, and solid lines are the prediction obtained with the Gram-Charlier formalism at next-to-leading order. The error bars are errors on the mean but are too small to be distinguished.

Before turning to the case of non-Gaussian corrections, let us emphasise that our Gaussian theoretical formalism allows us to obtain the probability of voids, walls, filaments, and knots as a function of the threshold itself, as illustrated in Fig. 3 for redshift z = 0 and Gaussian smoothing R = 15 Mpc h⁻¹. The threshold used in the above analysis (Λ_th = 0.01) is the vertical black dashed line. For this choice of smoothing and redshift – which corresponds to a mildly non-linear regime – and for all the environments, we can draw the same conclusion as before: the qualitative picture is correctly captured but non-linear corrections are nonetheless necessary to improve the predictive power of our model. This figure could be helpful to guide the choice of threshold based on theoretical arguments. We note that an alternative approach could be to renounce defining a global threshold and choose it according to the studied environment(s). For example, one could determine the threshold that gives the 20% rarest knots and use the same threshold for the filaments, and obtain a threshold for the voids and walls in a similar fashion (which would be the opposite of the one used for knots and filaments in the simple Gaussian case). In this case, some spatial position may be in none of the environments and this will mean that this position is a transition between two environments. Another possibility would be to use the filling factor approach to fix a threshold in abundance (i.e., we fix the volume fraction occupied by the excursion above ν, Gott et al. 1987; Matsubara 2003): this would lead to a remapping of our environments.

Fig. 3. Probability of the environment as a function of the threshold for a Gaussian random field (dashed lines), Gram-Charlier (solid lines), and measurement in the simulations (dots) at R = 15 Mpc h−1 and z = 0. The vertical black line is the threshold Λth = 0.01 (thus λth = Λth/σNL) used in this study.

Hereafter, we adhere to the standard global-threshold strategy. The description and inclusion of the above-mentioned non-linear theoretical corrections are now the goals of the remainder of this paper.

4. Mild non-Gaussian corrections to cosmic web abundances

Relying on Gaussian random fields is valid only in the linear regime of structure formation. However, at low redshift and small scales, increasing numbers of non-Gaussianities appear in the density field and corrections to the Gaussian predictions need to be accounted for. To improve our previous Gaussian predictions, we now propose to work with a probability distribution function 𝒫(X) that is no more Gaussian but includes non-linearities in a perturbative manner. In practice, we use a Gram-Charlier expansion following previous works in the literature including Pogosyan et al. (2009a), Gay et al. (2012) and Codis et al. (2013).

4.1. Gram-Charlier expansion of the joint distribution

For a set X of random fields, the Gram-Charlier expansion of the joint PDF reads

$$\begin{array}{c}\hfill \mathcal{P}(X)={\mathcal{P}}_{\mathrm{G}}(X)[1+{\displaystyle \sum _{n=3}^{\infty }}\frac{1}{n!}Tr[{⟨X^n⟩}_{\mathrm{GC}}.h_n(X)]],\end{array}$$(26)

where 𝒫_G(X) is a Gaussian kernel as defined in Eq. (12), h_n(x) are the Hermite tensors $h_n(X)={(-1)}^n{\mathcal{P}}_{\mathrm{G}}^{-1}(X){\partial }^n{\mathcal{P}}_{\mathrm{G}}(X)/\partial X^n$ and ⟨Xⁿ⟩_GC = ⟨h_n(X)⟩ are the Gram Charlier coefficients.

For our rotation-invariant variables J₁, J₂, J₃, this translates into the following expression

$$\begin{array}{cc}\hfill \mathcal{P}(J_1,J_2,J_3)& ={\mathcal{P}}_{\mathrm{G}}(J_1,J_2,J_3)[1+{\displaystyle \sum _{n=3}^{\infty }}{\displaystyle \sum _{k,l}^{k+2l=n}}\frac{{(-1)}^l{5}^l\times 3}{k!(3+2l)!!}{⟨J_1^k{J}_2^l⟩}_{\mathrm{GC}}{H}_k(J_1)L_l^{(3/2)}\left(\frac{5}{2}{J}_2\right)\hfill \\ \hfill & +{\displaystyle \sum _{n=3}^{\infty }}{\displaystyle \sum _k^{k+3=n}}\frac{25}{k!\times 21}{⟨J_1^k{J}_3⟩}_{\mathrm{GC}}{H}_k(J_1)J_3+{\displaystyle \sum _{n=5}^{\infty }}{\displaystyle \sum _{k,l,m=1}^{k+2l+3m=n}}\frac{c_{\mathit{lm}}}{k!}{⟨J_1^k{J}_2^l{J}_3^m⟩}_{\mathrm{GC}}{H}_k(J_1)F_{\mathit{lm}}(J_2,J_3)],\hfill \end{array}$$(27)

where 𝒫_G(J₁, J₂, J₃) is the Gaussian case given by Eq. (17), H_n are the successive Hermite polynomials, $L_l^{(\alpha )}(x)$ are the generalised Laguerre polynomials, c_lm are the normalisation coefficients that we leave undetermined and the orthogonal polynomials F_lm associated with J₂ and J₃ are such that

$$\begin{array}{c}\hfill {\displaystyle {\int }_0^{\infty }}\mathrm{d}{J}_2{\displaystyle {\int }_{-J_2^{3/2}}^{J_2^{3/2}}}\mathrm{d}{J}_3{\mathcal{P}}_{\mathrm{G}}(J_1,J_2,J_3)F_{\mathit{lm}}(J_2,J_3)F_{l^{\prime }{m}^{\prime }}(J_2,J_3)={\delta }_l^{l^{\prime }}{\delta }_m^{m^{\prime }}.\end{array}$$

We have, for example, these special cases: $F_{l0}=\sqrt{3\times 2^l\times l!/(3+2l)!!}\times L_l^{3/2}(5J_2/2)$ and $F_{01}=5J_3/\sqrt{21}$.

Here, we focus on the first corrective term, n = 3, such that

$$\begin{array}{cc}\hfill \mathcal{P}(J_1,J_2,J_3)& ={\mathcal{P}}_{\mathrm{G}}(J_1,J_2,J_3)[1+\frac{1}{6}{⟨J_1^3⟩}_{\mathrm{GC}}{H}_3(J_1)-{⟨J_1{J}_2⟩}_{\mathrm{GC}}{H}_1(J_1)L_1^{(3/2)}\left(\frac{5}{2}{J}_2\right)+\frac{25}{21}{⟨J_3⟩}_{\mathrm{GC}}{J}_3]+o({\sigma }_0^2).\hfill \end{array}$$(28)

The Gram-Charlier cumulants read

$$\begin{array}{c}\hfill {⟨J_1^3⟩}_{\mathrm{GC}}=⟨H_3(J_1)⟩,{⟨J_1{J}_2⟩}_{\mathrm{GC}}=-\frac{2}{5}\langle H_1(J_1)L_1^{(3/2)}\left(\frac{5}{2}{J}_2\right)\rangle ,{⟨J_3⟩}_{\mathrm{GC}}=⟨J_3⟩,\end{array}$$(29)

such that we finally obtain the Gram-Charlier expression at first non-linear order (NLO):

$$\begin{array}{c}\hfill \mathcal{P}(J_1,J_2,J_3)={\mathcal{P}}_{\mathrm{G}}(J_1,J_2,J_3)[1+\frac{1}{6}⟨J_1^3⟩(J_1^3-3J_1)+\frac{5}{2}⟨J_1{J}_2⟩J_1(J_2-1)+\frac{25}{21}⟨J_3⟩J_3]+o({\sigma }_0^2).\end{array}$$(30)

Three cumulants appear in the next-to-leading-order Gram-Charlier probability distribution function (Eq. (30)): $\langle J_1^3\rangle$, ⟨J₁J₂⟩ and ⟨J₃⟩. In Eulerian perturbation theory, those cumulants are linear in σ₀ at tree order¹ and we therefore introduce the reduced cumulants $S_3=\langle J_1^3\rangle /{\sigma }_0$, U₃ = ⟨J₁J₂⟩/σ₀ and V₃ = ⟨J₃⟩/σ₀, which are constant in time at tree order. For example, S₃ is the usual cosmological skewness, whose analytical prediction at tree order is well known (for example, Peebles 1980; Juszkiewicz et al. 1993; Lokas et al. 1995; Colombi et al. 2000). The other two reduced cumulants can be computed in a similar manner, described in Appendix A. Figure 4 shows the resulting tree-order cumulants as a function of the smoothing scale in dashed black, compared to the measurements in the simulation at different redshifts (as shown with different colours). We see a good agreement at almost all smoothing scales and redshifts as the prediction is almost always within the error bars. For large smoothing scales, the error bars increase due to the finite volume of the Quijote simulation which thus misses large wave modes. As expected, at low redshift and small smoothing scales, departures from tree-order predictions are seen as the non-linearities increase. In this work, we focus on mildly non-linear scales (about 10 Mpc h⁻¹ and above) where a perturbative treatment is accurate, as illustrated by the values of the cumulants depicted in this figure.

Fig. 4. Cumulants S3, U3, and V3 as a function of the smoothing scale R. The dashed black line is the analytical prediction. The coloured lines are the measurements from the simulations at different redshifts.

4.2. Cosmic web abundances at next-to-leading order

We now turn to the explicit computation of the probability of different cosmic environments at next-to-leading order in the Gram-Charlier formalism. Let us first rewrite the Gram-Charlier expansion of the joint distribution of the rotational invariants of the tidal tensor as

$$\begin{array}{c}\hfill \mathcal{P}(J_1,J_2,J_3)={\mathcal{P}}_{\mathrm{G}}(J_1,J_2,J_3)+{\sigma }_0{\mathcal{P}}^{(3)}(J_1,J_2,J_3)+O({\sigma }_0^2),\end{array}$$(31)

where 𝒫_G is the Gaussian part and 𝒫⁽³⁾ is the first non-Gaussian corrective term

$$\begin{array}{c}\hfill {\mathcal{P}}^{(3)}(J_1,J_2,J_3)=\frac{25\sqrt{5}}{12\pi }exp[-\frac{J_1^2}{2}-\frac{5J_2}{2}](\frac{1}{6}(J_1^3-3J_1)S_3+\frac{5}{2}{J}_1(J_2-1)U_3+\frac{25}{21}{J}_3{V}_3).\end{array}$$(32)

Integrating this distribution function over the constraints given in Eqs. (18)–(21) will give us the expression of the probability of a given environment (E) as

$$\begin{array}{cc}\hfill P_{(\mathrm{E})}({\lambda }_{\mathrm{th}})& =P_{\mathrm{G},(\mathrm{E})}({\lambda }_{\mathrm{th}})+{\sigma }_0{P}_{(\mathrm{E})}^{(3)}({\lambda }_{\mathrm{th}})\hfill \\ \hfill & =P_{\mathrm{G},(\mathrm{E})}({\lambda }_{\mathrm{th}})+{\sigma }_0[s_{(\mathrm{E})}({\lambda }_{\mathrm{th}})S_3+u_{(\mathrm{E})}({\lambda }_{\mathrm{th},})U_3+v_{(\mathrm{E})}({\lambda }_{\mathrm{th}})V_3],\hfill \end{array}$$(33)

where s, u and v are threshold-dependent functions and where the Gaussian probability of the environments was derived in Eqs. (22)–(25).

Again, two degrees of freedom can be integrated out and one remains to be integrated numerically. The resulting expressions for s(λ_th), u(λ_th), and v(λ_th) in each environment are given below (written in terms of 1D integrals over J₂).

4.2.1. Void

For voids, we get

$$\begin{array}{cc}\hfill s_{\mathrm{void}}({\lambda }_{\mathrm{th}})& ={\displaystyle {\int }_0^{\infty }}\frac{25\sqrt{5}}{144\pi }{e}^{\frac{1}{2}(-9)(J_2+{\lambda }_{\mathrm{th}}^2)}[-3\sqrt{2\pi }{e}^{2J_2+\frac{9{\lambda }_{\mathrm{th}}^2}{2}}\mathrm{erf}\left(\frac{2\sqrt{J_2}-3{\lambda }_{\mathrm{th}}}{\sqrt{2}}\right)+3\sqrt{2\pi }{e}^{2J_2+\frac{9{\lambda }_{\mathrm{th}}^2}{2}}\mathrm{erf}\left(\frac{\sqrt{\mathrm{J}2}-3{\lambda }_{\mathrm{th}}}{\sqrt{2}}\right)\hfill \\ \hfill & +\sqrt{J_2}{e}^{3\sqrt{J_2}{\lambda }_{\mathrm{th}}}[e^{3\sqrt{J_2}{\lambda }_{\mathrm{th}}}(-48J_2^{3/2}{\lambda }_{\mathrm{th}}+16J_2^2+36J_2{\lambda }_{\mathrm{th}}^2-27\sqrt{J_2}{\lambda }_{\mathrm{th}}\hfill \\ \hfill & +14J_2+12)-6e^{\frac{3J_2}{2}}]]+\frac{25\sqrt{5}}{36\pi }{J}_2^{3/2}{e}^{6\sqrt{J_2}{\lambda }_{\mathrm{th}}-\frac{9J_2}{2}-\frac{9{\lambda }_{\mathrm{th}}^2}{2}}(1-{(2\sqrt{J_2}-3{\lambda }_{\mathrm{th}})}^2)\mathrm{d}{J}_2,\hfill \end{array}$$(34)

$$\begin{array}{cc}\hfill u_{\mathrm{void}}({\lambda }_{\mathrm{th}})& ={\displaystyle {\int }_0^{\infty }}-\frac{125\sqrt{5}}{96\pi }(J_2-1)e^{\frac{1}{2}(-9)(J_2+{\lambda }_{\mathrm{th}}^2)}[-3\sqrt{2\pi }{e}^{2J_2+\frac{9{\lambda }_{\mathrm{th}}^2}{2}}(J_2-9{\lambda }_{\mathrm{th}}^2-1)\times \mathrm{erf}\left(\frac{2\sqrt{J_2}-3{\lambda }_{\mathrm{th}}}{\sqrt{2}}\right)+3\sqrt{2\pi }{e}^{2J_2+\frac{9{\lambda }_{\mathrm{th}}^2}{2}}(J_2-9{\lambda }_{\mathrm{th}}^2-1)\mathrm{erf}\left(\frac{\sqrt{J_2}-3{\lambda }_{\mathrm{th}}}{\sqrt{2}}\right)\hfill \\ \hfill & -2e^{6\sqrt{J_2}{\lambda }_{\mathrm{th}}}(4J_2^{3/2}+6\sqrt{J_2}+9{\lambda }_{\mathrm{th}})+6e^{3\sqrt{J_2}{\lambda }_{\mathrm{th}}+\frac{3J_2}{2}}(\sqrt{J_2}+3{\lambda }_{\mathrm{th}})]-\frac{125\sqrt{5}}{12\pi }(J_2-1)J_2^{3/2}{e}^{6\sqrt{J_2}{\lambda }_{\mathrm{th}}-\frac{9J_2}{2}-\frac{9{\lambda }_{\mathrm{th}}^2}{2}}\mathrm{d}{J}_2,\hfill \end{array}$$(35)

$$\begin{array}{cc}\hfill v_{\mathrm{void}}({\lambda }_{\mathrm{th}})& ={\displaystyle {\int }_0^{\infty }}\frac{625\sqrt{5}}{4032\pi }{e}^{-\frac{5J_2}{2}}[-\sqrt{2\pi }[-4J_2^3+9J_2^2+243(5-2J_2){\lambda }_{\mathrm{th}}^4+81((J_2-4)J_2+5){\lambda }_{\mathrm{th}}^2-18J_2+729{\lambda }_{\mathrm{th}}^6+15]\hfill \\ \hfill & \times \mathrm{erf}\left(\frac{3{\lambda }_{\mathrm{th}}-2\sqrt{J_2}}{\sqrt{2}}\right)+\sqrt{2\pi }[-4J_2^3+9J_2^2+243(5-2J_2){\lambda }_{\mathrm{th}}^4+81((J_2-4)J_2+5){\lambda }_{\mathrm{th}}^2-18J_2+729{\lambda }_{\mathrm{th}}^6+15]\mathrm{erf}\left(\frac{3{\lambda }_{\mathrm{th}}-\sqrt{J_2}}{\sqrt{2}}\right)\hfill \\ \hfill & +2e^{-\frac{1}{2}{(\sqrt{J_2}-3{\lambda }_{\mathrm{th}})}^2}[J_2^{3/2}(-(45{\lambda }_{\mathrm{th}}^2+13))+12J_2^2{\lambda }_{\mathrm{th}}+4J_2^{5/2}-9J_2{\lambda }_{\mathrm{th}}(15{\lambda }_{\mathrm{th}}^2+7)+3\sqrt{J_2}(27{\lambda }_{\mathrm{th}}^4+36{\lambda }_{\mathrm{th}}^2+5)\hfill \\ \hfill & +9{\lambda }_{\mathrm{th}}(27{\lambda }_{\mathrm{th}}^4+42{\lambda }_{\mathrm{th}}^2+11)]-2e^{-\frac{1}{2}{(2\sqrt{J_2}-3{\lambda }_{\mathrm{th}})}^2}[J_2^{3/2}(4-36{\lambda }_{\mathrm{th}}^2)+3J_2^2{\lambda }_{\mathrm{th}}+2J_2^{5/2}+J_2(18{\lambda }_{\mathrm{th}}-54{\lambda }_{\mathrm{th}}^3)\hfill \\ \hfill & +6\sqrt{J_2}(27{\lambda }_{\mathrm{th}}^4+36{\lambda }_{\mathrm{th}}^2+5)+9{\lambda }_{\mathrm{th}}(27{\lambda }_{\mathrm{th}}^4+42{\lambda }_{\mathrm{th}}^2+11)]]\mathrm{d}{J}_2,\hfill \end{array}$$(36)

where erf(z) is, again, the error function. Once that analytical expression is obtained, we can perform the numerical integration over J₂. We display s(λ_th), u(λ_th), and v(λ_th) as a function of the threshold in Fig. B.1.

4.2.2. Wall

For walls, the same procedure gives

$$\begin{array}{cc}\hfill s_{\mathrm{wall}}({\lambda }_{\mathrm{th}})& ={\displaystyle {\int }_0^{\infty }}\frac{25\sqrt{5}}{48\pi }{e}^{-\frac{9J_2}{2}}[-\sqrt{2\pi }{e}^{2J_2}\mathrm{erf}\left(\frac{3{\lambda }_{\mathrm{th}}-2\sqrt{J_2}}{\sqrt{2}}\right)+\sqrt{2\pi }{e}^{2J_2}\mathrm{erf}\left(\frac{\sqrt{J_2}+3{\lambda }_{\mathrm{th}}}{\sqrt{2}}\right)+\sqrt{J_2}(-e^{\frac{1}{2}(-3){\lambda }_{\mathrm{th}}(2\sqrt{J_2}+3{\lambda }_{\mathrm{th}})})\hfill \\ \hfill & \times (e^{9\sqrt{J_2}{\lambda }_{\mathrm{th}}}(-9\sqrt{J_2}{\lambda }_{\mathrm{th}}+6J_2+4)+2e^{\frac{3J_2}{2}})]\mathrm{d}{J}_2,\hfill \end{array}$$(37)

$$\begin{array}{cc}\hfill u_{\mathrm{wall}}({\lambda }_{\mathrm{th}})& ={\displaystyle {\int }_0^{\infty }}\frac{125\sqrt{5}}{32\pi }{e}^{-\frac{9J_2}{2}}(J_2-1)[-\sqrt{2\pi }{e}^{2J_2}(J_2-9{\lambda }_{\mathrm{th}}^2-1)\mathrm{erf}\left(\frac{\sqrt{J_2}+3{\lambda }_{\mathrm{th}}}{\sqrt{2}}\right)+\sqrt{2\pi }{e}^{2J_2}(J_2-9{\lambda }_{\mathrm{th}}^2-1)\mathrm{erf}\left(\frac{3{\lambda }_{\mathrm{th}}-2\sqrt{J_2}}{\sqrt{2}}\right)\hfill \\ \hfill & -2e^{\frac{1}{2}(-3){\lambda }_{\mathrm{th}}(2\sqrt{J_2}+3{\lambda }_{\mathrm{th}})}[e^{9\sqrt{J_2}{\lambda }_{\mathrm{th}}}(2\sqrt{J_2}+3{\lambda }_{\mathrm{th}})+e^{\frac{3J_2}{2}}(\sqrt{J_2}-3{\lambda }_{\mathrm{th}})]]\mathrm{d}{J}_2,\hfill \end{array}$$(38)

$$\begin{array}{cc}\hfill v_{\mathrm{wall}}({\lambda }_{\mathrm{th}})& ={\displaystyle {\int }_0^{\infty }}-\frac{625\sqrt{5}}{4032\pi }{e}^{-\frac{5J_2}{2}}[-\sqrt{2\pi }[-4J_2^3+9J_2^2+243(5-2J_2){\lambda }_{\mathrm{th}}^4+81((J_2-4)J_2+5){\lambda }_{\mathrm{th}}^2-18J_2+729{\lambda }_{\mathrm{th}}^6+15]\mathrm{erf}\left(\frac{3{\lambda }_{\mathrm{th}}-2\sqrt{J_2}}{\sqrt{2}}\right)\hfill \\ \hfill & +\sqrt{2\pi }[-4J_2^3+9J_2^2+243(5-2J_2){\lambda }_{\mathrm{th}}^4+81((J_2-4)J_2+5){\lambda }_{\mathrm{th}}^2-18J_2+729{\lambda }_{\mathrm{th}}^6+15]\mathrm{erf}\left(\frac{\sqrt{J_2}+3{\lambda }_{\mathrm{th}}}{\sqrt{2}}\right)\hfill \\ \hfill & -2e^{-\frac{1}{2}{(2\sqrt{J_2}-3{\lambda }_{\mathrm{th}})}^2}[J_2^{3/2}(4-36{\lambda }_{\mathrm{th}}^2)+3J_2^2{\lambda }_{\mathrm{th}}+2J_2^{5/2}+J_2(18{\lambda }_{\mathrm{th}}-54{\lambda }_{\mathrm{th}}^3)+6\sqrt{J_2}(27{\lambda }_{\mathrm{th}}^4+36{\lambda }_{\mathrm{th}}^2+5)+9{\lambda }_{\mathrm{th}}(27{\lambda }_{\mathrm{th}}^4+42{\lambda }_{\mathrm{th}}^2+11)]\hfill \\ \hfill & +2e^{-\frac{1}{2}{(\sqrt{J_2}+3{\lambda }_{\mathrm{th}})}^2}[J_2^{3/2}(45{\lambda }_{\mathrm{th}}^2+13)+12J_2^2{\lambda }_{\mathrm{th}}-4J_2^{5/2}-9J_2{\lambda }_{\mathrm{th}}(15{\lambda }_{\mathrm{th}}^2+7)-3\sqrt{J_2}(27{\lambda }_{\mathrm{th}}^4+36{\lambda }_{\mathrm{th}}^2+5)+9{\lambda }_{\mathrm{th}}(27{\lambda }_{\mathrm{th}}^4+42{\lambda }_{\mathrm{th}}^2+11)]]\mathrm{d}{J}_2.\hfill \end{array}$$(39)

4.2.3. Filament

For filaments, we get

$$\begin{array}{cc}\hfill s_{\mathrm{filament}}({\lambda }_{\mathrm{th}})& ={\displaystyle {\int }_0^{\infty }}-\frac{25\sqrt{5}}{144\pi }{e}^{\frac{1}{2}(-3)(4\sqrt{J_2}{\lambda }_{\mathrm{th}}+3J_2+3{\lambda }_{\mathrm{th}}^2)}[3\sqrt{2\pi }{e}^{\frac{1}{2}{(2\sqrt{J_2}+3{\lambda }_{\mathrm{th}})}^2}\times \mathrm{erf}\left(\frac{2\sqrt{J_2}+3{\lambda }_{\mathrm{th}}}{\sqrt{2}}\right)+3\sqrt{2\pi }{e}^{\frac{1}{2}{(2\sqrt{J_2}+3{\lambda }_{\mathrm{th}})}^2}\mathrm{erf}\left(\frac{\sqrt{J_2}-3{\lambda }_{\mathrm{th}}}{\sqrt{2}}\right)\hfill \\ \hfill & -3\sqrt{J_2}(2e^{9\sqrt{J_2}{\lambda }_{\mathrm{th}}+\frac{3J_2}{2}}+9\sqrt{J_2}{\lambda }_{\mathrm{th}}+6J_2+4)]\mathrm{d}{J}_2,\hfill \end{array}$$(40)

$$\begin{array}{cc}\hfill u_{\mathrm{filament}}({\lambda }_{\mathrm{th}})& ={\displaystyle {\int }_0^{\infty }}-\frac{125\sqrt{5}}{48\pi }(J_2-1)[3\sqrt{\frac{\pi }{2}}{e}^{-\frac{5J_2}{2}}(J_2-9{\lambda }_{\mathrm{th}}^2-1)\times \mathrm{erf}\left(\frac{3{\lambda }_{\mathrm{th}}-\sqrt{J_2}}{\sqrt{2}}\right)-\frac{3}{2}{e}^{\frac{1}{2}(-3)(4\sqrt{J_2}{\lambda }_{\mathrm{th}}+3J_2+3{\lambda }_{\mathrm{th}}^2)}\hfill \\ \hfill & \times [\sqrt{2\pi }{e}^{\frac{1}{2}{(2\sqrt{J_2}+3{\lambda }_{\mathrm{th}})}^2}(J_2-9{\lambda }_{\mathrm{th}}^2-1)\times \mathrm{erf}\left(\frac{2\sqrt{J_2}+3{\lambda }_{\mathrm{th}}}{\sqrt{2}}\right)+6{\lambda }_{\mathrm{th}}(e^{9\sqrt{J_2}{\lambda }_{\mathrm{th}}+\frac{3J_2}{2}}-1)+2\sqrt{J_2}(e^{9\sqrt{J_2}{\lambda }_{\mathrm{th}}+\frac{3J_2}{2}}+2)]]\mathrm{d}{J}_2,\hfill \end{array}$$(41)

$$\begin{array}{cc}\hfill v_{\mathrm{filament}}({\lambda }_{\mathrm{th}})& ={\displaystyle {\int }_0^{\infty }}\frac{625\sqrt{5}}{4032\pi }{e}^{-\frac{5J_2}{2}}[-\sqrt{2\pi }[-4J_2^3+9J_2^2+243(5-2J_2){\lambda }_{\mathrm{th}}^4+81((J_2-4)J_2+5){\lambda }_{\mathrm{th}}^2-18J_2+729{\lambda }_{\mathrm{th}}^6+15]\mathrm{erf}\left(\frac{3{\lambda }_{\mathrm{th}}-\sqrt{J_2}}{\sqrt{2}}\right)\hfill \\ \hfill & +\sqrt{2\pi }[-4J_2^3+9J_2^2+243(5-2J_2){\lambda }_{\mathrm{th}}^4+81((J_2-4)J_2+5){\lambda }_{\mathrm{th}}^2-18J_2+729{\lambda }_{\mathrm{th}}^6+15]\mathrm{erf}\left(\frac{2\sqrt{J_2}+3{\lambda }_{\mathrm{th}}}{\sqrt{2}}\right)-2e^{-\frac{1}{2}{(\sqrt{J_2}-3{\lambda }_{\mathrm{th}})}^2}\hfill \\ \hfill & \times [J_2^{3/2}(-(45{\lambda }_{\mathrm{th}}^2+13))+12J_2^2{\lambda }_{\mathrm{th}}+4J_2^{5/2}-9J_2{\lambda }_{\mathrm{th}}(15{\lambda }_{\mathrm{th}}^2+7)+3\sqrt{J_2}(27{\lambda }_{\mathrm{th}}^4+36{\lambda }_{\mathrm{th}}^2+5)+9{\lambda }_{\mathrm{th}}(27{\lambda }_{\mathrm{th}}^4+42{\lambda }_{\mathrm{th}}^2+11)]\hfill \\ \hfill & +e^{-\frac{1}{2}{(2\sqrt{J_2}+3{\lambda }_{\mathrm{th}})}^2}[8J_2^{3/2}(9{\lambda }_{\mathrm{th}}^2-1)+6J_2^2{\lambda }_{\mathrm{th}}-4J_2^{5/2}-36J_2{\lambda }_{\mathrm{th}}(3{\lambda }_{\mathrm{th}}^2-1)-12\sqrt{J_2}(27{\lambda }_{\mathrm{th}}^4+36{\lambda }_{\mathrm{th}}^2+5)+18{\lambda }_{\mathrm{th}}(27{\lambda }_{\mathrm{th}}^4+42{\lambda }_{\mathrm{th}}^2+11)]]\mathrm{d}{J}_2.\hfill \end{array}$$(42)