© The Authors 2026

1. Introduction

The large-scale structure of the Universe exhibits a striking and complex pattern known as the cosmic web. This is a vast network of halos, walls, filaments, and voids that emerged from tiny primordial density perturbations (Zel’dovich 1970; Peebles 1980; Bond et al. 1996; Springel et al. 2005). In the standard ΛCDM cosmological framework, approximately 84% of the matter content in the Universe (see, e.g., Planck Collaboration VI 2020) is constituted of cold dark matter (CDM). This is a nonrelativistic collisionless component that drives structure formation. As the Universe expanded and evolved under gravity, overdense regions grew, giving rise to a hierarchical structure in which CDM clustered into halos that are connected by elongated filaments and walls and are surrounded by underdense voids. The baryons having decoupled from the early hot plasma condensed into collapsed halos, which became the primary sites of galaxy formation. This web-like arrangement of galaxies has been revealed in a series of pioneering redshift surveys, from the early work of Jõeveer et al. (1978) and de Lapparent et al. (1986) to large-scale surveys such as 2dF (Colless et al. 2003), Sloan Digital Sky Survey (SDSS) (Gott et al. 2005), and Two Micron All Sky Survey (2MASS) (Huchra et al. 2012). Beyond its role in structuring the Universe, the cosmic web serves as a sensitive probe of fundamental physics: its morphology and connectivity carry imprints of the nature of dark matter, the dynamics of structure formation, and the statistical properties of the primordial density field, including possible non-Gaussianities (Libeskind et al. 2018; Codis et al. 2012; Cautun et al. 2014; Kitaura et al. 2019). The study of the cosmic web therefore is of growing importance in the era of precision cosmology and large-scale surveys such as Dark Energy Spectroscopic Instrument (DESI) (DESI Collaboration 2016) and Euclid (Euclid Collaboration: Mellier et al. 2025).

As CDM is the dominant driver of structure formation, analyses of the cosmic web are naturally framed in terms of dark matter dynamics. CDM is modeled as a self-gravitating and collisionless fluid obeying the Vlasov–Poisson equations,

$$\begin{array}{cc}& \frac{\partial f}{\partial t}+\mathbf{u}·{\mathrm{\nabla }}_{\mathbf{r}}f-{\mathrm{\nabla }}_{\mathbf{r}}\mathrm{\Phi }·{\mathrm{\nabla }}_{\mathbf{u}}f=0,\hfill \end{array}$$(1)

$$\begin{array}{cc}& {\mathrm{\Delta }}_{\mathbf{r}}\mathrm{\Phi }=4\pi G\rho =4\pi G{\displaystyle \int f}(\mathbf{r},\mathbf{u},t){\mathrm{d}}^3\mathbf{u},\hfill \end{array}$$(2)

where f(r, u, t) represents the phase-space density at position r, velocity u, and time t; and Φ is the gravitational potential. Owing to negligible velocity dispersion, the CDM distribution can be thought of as a 2D (or 3D) sheet (Abel et al. 2012; Shandarin et al. 2012) in 4D (or 6D) phase space. The system of equations is usually numerically resolved with an N-body approach, that is, representing the CDM sheet by an ensemble of particles that follow the standard Lagrangian equations of motion in an expanding Universe,

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2\mathbf{x}}{\mathrm{d}{t}^2}+2H\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}=-\frac{1}{a^2}{\mathrm{\nabla }}_{\mathbf{x}}\mathrm{\Phi }(\mathbf{x});{\mathrm{\Delta }}_{\mathbf{x}}\mathrm{\Phi }(\mathbf{x})=4\pi G\overline{\rho }{a}^2\delta (\mathbf{x}),\end{array}$$(3)

where x is the comoving position, a is the cosmological expansion factor corresponding to time t, and δ is the matter density contrast, defined as $\delta =\rho /\overline{\rho }-1$, with $\overline{\rho }$ being the background matter density.

On the analytical front, the single-stream evolution from the primordial Gaussian random field to the onset of shell crossing can be accurately described using Lagrangian perturbation theory (LPT; see, e.g., Bernardeau et al. 2002, and references therein). The pioneering Zel’dovich approximation (ZA; also 1-LPT; Zel’dovich 1970) was the first to capture the anisotropic nature of collapse and the emergence of pancake-like structures that evolve into a hierarchy of halos, filaments, and walls in a cosmological setting. Building on this, Arnold et al. (1982) framed the emerging structures in terms of the topology of singularities in the CDM manifold in 1D and 2D. Arnold (1982, 1983) presented the normal forms of all generic singularities in 3D, laying the foundation of catastrophe theory (Arnold 1984). Under this formalism, the structural hierarchy (halos, filaments, walls, and voids) as well as the singularities that seed them are characterized through the eigenvalues and eigenvectors of the deformation tensor. On the statistical side, Doroshkevich (1970), Doroshkevich & Shandarin (1978) derived the distribution of these eigenvalues and applied it to make statistical predictions about the pancake area and curvature, while Bardeen et al. (1986) developed the statistics of peaks in Gaussian random fields. Later, Bond et al. (1996) coined the term “cosmic web”, linking the statistics of eigenvalues to predictions about the filamentary network connecting halos.

Considerable effort has been devoted to tracing the cosmic web by reducing it to the underlying skeleton or spine of the large-scale structure. Early descriptions include the adhesion model (Gurbatov et al. 1989), while nonlocal ridge-tracing techniques using the density field in 2D and 3D were developed by Novikov et al. (2006) and Sousbie et al. (2008), respectively. Incorporating the information from the deformation tensor, more recent approaches (Hidding et al. 2014; Feldbrugge et al. 2018) used catastrophe theory to define the web as the locus of singularities. Extending this framework by combining Gaussian statistics, Feldbrugge & van de Weygaert (2023), Feldbrugge et al. (2023), Feldbrugge & van de Weygaert (2025) probed the statistical properties of the progenitors of the hierarchical structures in 2D, such as the correlation function of singular points, filament formation times, and multistream properties.

We focus on the 2D case, which is analytically more tractable, but retains the essential features of the hierarchical structures present in 3D. In Sect. 2 we draw on catastrophe theory (Arnold et al. 1982) to derive the 2D normal form (valid up to O(q⁴), where q is the Lagrangian coordinate) for the subset of singularities that seed pancakes, which we then use to characterize their shape and related properties. Section 3 extends the distribution of deformation-tensor eigenvalues (Doroshkevich 1970) to include O(q⁴) derivatives of the Gaussian field, treated as parameters in the normal form. Adopting the peak statistics formalism of Bardeen et al. (1986), we compute their conditional distribution at singular points where pancakes form. Together, these ingredients yield the distributions of several observable pancake features. Our framework is in line with recent advances (Feldbrugge & van de Weygaert 2025) and their precursors, but places particular emphasis on pancake geometry, especially curvature and higher-order corrections that produce C- and S-type morphologies. Finally, Sect. 4 summarizes the main results of our 2D analysis, outlines extensions to 3D, and discusses the observational ramifications. We have relegated the extended formulae appearing in the derivation of our analytical results in Sects. 2 and 3 to Appendices A and B, respectively. In Appendix C we corroborate our analytical and statistical predictions with measurements from numerical simulations of Zeldovich flow initialized with Gaussian random fields.

2. Analytical model using catastrophe theory

2.1. Lagrangian formalism and shell crossing

In Lagrangian formalism, we track the motion of particles x(q, t) in terms of their displacement Ψ(q, t) from initial positions q,

$$\begin{array}{c}\hfill \mathbf{x}(\mathit{q},t)=\mathit{q}+\mathrm{\Psi }(\mathit{q},t).\end{array}$$(4)

In Lagrangian q-space, the density is uniform and equal to the background density in our expanding universe, that is the Eulerian x-space. The Eulerian density is given by the Jacobian of the transformation between Lagrangian and Eulerian spaces,

$$\begin{array}{cc}\hfill J_{\mathit{ij}}(\mathit{q},t)& =\frac{\partial x_i}{\partial q_j}(\mathit{q},t)={\delta }_{\mathit{ij}}^D+\frac{\partial {\mathrm{\Psi }}_i}{\partial q_j}(\mathit{q},t)\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill \rho (\mathbf{x},t)& =\overline{\rho }(t)(1+\delta (\mathbf{x},t))=\frac{\overline{\rho }(t)}{J(\mathit{q},t)},\hfill \end{array}$$(6)

where J is the determinant of the Jacobian, and δ^D is the Dirac-delta function. We define the deformation matrix d_ij(q, t) = δ_ij^D − J_ij(q, t) =  − ∂Ψ_i / ∂q_j(q, t) with eigenvalues α > β and corresponding eigenvectors ${\widehat{n}}^{\alpha },{\widehat{n}}^{\beta }$, which together encode the complete dynamical information. In terms of eigenvalues,

$$\begin{array}{c}\hfill \rho (\mathbf{x},t)\propto J{(\mathit{q},t)}^{-1}={[(1-\alpha (\mathit{q},t))(1-\beta (\mathit{q},t))]}^{-1}.\end{array}$$(7)

Within the framework of catastrophe theory, the particle motion can be viewed as the evolving shape of the CDM sheet in x − q space (see Fig. 1). The topology of the sheet is determined by the set of critical points: J = ∥∂x_i/∂q_j∥ = 0. The tangent to the sheet at critical points is perpendicular to the Eulerian space and hence, the density at these points projected onto the Eulerian space is singular, as implied by Eq. (6). From Eq. (7), the density being singular implies

$$\begin{array}{c}\hfill \alpha (\mathit{q},t)=1\end{array}$$(8)

Fig. 1. 2D shell crossing in Eulerian x − y space (top row) along with a qx axis slice through the shell-crossing point qc in qx − x space (bottom row) to help visualize the folding of the CDM sheet. Left column: Snapshot during the single-stream flow prior to shell crossing (t < tc). Middle column: Moment of shell crossing (t = tc). A catastrophe event marked by the appearance of an A3+ singularity (red dot) seeds the 2D pancake-like structure, shown in the right column (t > tc). The boundary of the pancake encompassing the multi-streaming region is defined by A2 singularities (blue), which grows along the A3+ spine (green).

for points undergoing shell crossing for the first time (see Sect. 5(a) in Doroshkevich & Shandarin 1978 and Eq. (10) in Arnold et al. 1982). This condition identifies the A₂ class of singularities in 2D.

When shell crossing occurs for the first time in a local neighborhood, the sheet starts folding onto itself, turning perpendicular to the Eulerian space and forming a cusp. The first critical point (or singularity) appears, indicating a change in topology (or catastrophe). We refer to this point as the shell-crossing point (q_c, t_c). It is the first point to satisfy the condition α(q, t) = 1 in its local neighborhood and is therefore a local maximum of the eigenvalue field, that is ∇_qα(q, t) = 0. In 2D, it is a singularity of class A₃⁺, also called cusp singularity. shell crossing typically happens along one axis; it is nearly impossible to have shell crossings along multiple axes simultaneously (see Eq. (12) in Doroshkevich 1970). Therefore, as the CDM sheet folds further, its projection onto the Eulerian space is shaped like a pancake, bulging towards the center. In 2D, it is shaped like an eye. The cusps at either end of the fold in CDM sheet branch out in opposite directions (${\widehat{n}}_{\beta }$) transverse to the shell-crossing direction (${\widehat{n}}_{\alpha }$) and follow the locus of A₃⁺ singularities (see Sect. 5(b) in Doroshkevich & Shandarin 1978 and Eq. (13) in Arnold et al. 1982), given by

$$\begin{array}{c}\hfill {\overline{\mathrm{\nabla }}}_{\mathit{q}}\alpha (\mathit{q},t)·{\widehat{n}}^{\alpha }(\mathit{q},t)=0.\end{array}$$(9)

Thus, the A₃⁺ locus serves as the spine along which the pancake grows, with its boundary (lateral edges of the fold in CDM sheet) being defined by the A₂ locus Eq. (8). The A₂ boundary divides the multi-stream and single-stream regions in Eulerian space. As more particles enter the multi-stream region, the size of the pancake (A₂ boundary) grows, with the two A₃⁺ singularities at either tips branching out farther.

In fact, the locus of A₃⁺ singularities passes through all the critical – local maxima, minima, and saddle - points of the eigenvalue field, since ${\overline{\mathrm{\nabla }}}_{\mathit{q}}\alpha (\mathit{q},t)=0$ automatically satisfies Eq. (9). These are the points around which collapses occur and structures emerge. Therefore, the A₃⁺ locus traces not only the spine of pancakes, but also the entire the skeletal network of structures in our 2D setup, akin to the cosmic web of the universe in 3D. Put differently, the study of the cosmic web can be reduced to analyzing the critical points of the eigenvalue field α(q, t) within the framework of Morse-Smale theory (Smale 1960).

Growing pancakes form only around points of local maxima with positive eigenvalue α. These points, which we refer to as shell-crossing points, are the focus of our formalism. If there is no shell crossing along the transverse direction(s), the pancake grows into a filament in 2D (wall in 3D). If there is shell crossing along the secondary transverse direction, the pancake turns into a cluster in 2D (filament in 3D). Further, in 3D, shell crossing along the tertiary transverse direction turns it into a cluster. Clusters, followed by violent relaxation in the multistream regime (Lynden-Bell 1967), lead to the formation of primordial halos.

On the other hand, saddle points with positive eigenvalue α undergo shell crossing to form tube-shaped structures, whereas points of local minima correspond to no physical structures. The saddle points, however, are not the first points to undergo shell crossing in their local neighborhood, which is already deep in the multi-stream regime by the time these points turn singular. Although the same formalism, with minor tweaks discussed in the next subsection, holds true for these points asymptotically, its applicability in studying the large-scale structure around these points is limited.

Equations (3)–(6) form a complete set of equations that can be solved for a given set of initial conditions. For cosmological setups, the initial conditions are generally provided in terms of the matter power-spectrum P_δ(k) of a Gaussian random field for the overdensity δ(q), defined as $⟨\delta (\overline{k}){\delta }^{\ast }(\overline{k}\prime )⟩={(2\pi )}^2{\delta }_D^2(\overline{k}-\overline{k}\prime )P_{\delta }(k)$, where k is the Fourier mode and δ_D is the Dirac-delta function.

If there is no vorticity in Lagrangian space in our setup, then the displacement field derives from a potential ϕ

$$\begin{array}{c}\hfill \mathrm{\Psi }(\mathit{q},t)=-\overline{\mathrm{\nabla }}\varphi (\mathit{q},t),\end{array}$$(10)

which is related to the overdensity δ(q) in the linear regime: Δ_qϕ(q) = δ(q). The displacement potential ϕ(q) is therefore Gaussian distributed as well, with the power spectrum

$$\begin{array}{c}\hfill P_{\varphi }(k)=k^{-4}{P}_{\delta }(k).\end{array}$$(11)

Eq. (6) is valid as long as the flow is single-stream and one can obtain recursive solutions using LPT. After shell crossing, the flow is multistream inside the boundary of the pancake defined by A₂ singularities, and LPT breaks down. Equation (6) must be modified into

$$\begin{array}{c}\hfill \rho (\mathbf{x},t)=\overline{\rho }(t){\displaystyle \sum _{\mathit{q}:\mathbf{x}(\mathit{q},t)=\mathbf{x}}}J{(\mathit{q},t)}^{-1}\end{array}$$(12)

to take into account the contribution from the overlapping folds of the CDM sheet. This is the basis of post-collapse perturbation theory (PCPT; Colombi 2014; Taruya & Colombi 2017; Rampf et al. 2021; Saga et al. 2023) in which the corrections to the force field from the overlapping folds are computed that allows the post-collapse motion to be tracked deeper into the multi-stream regime.

2.2. Normal form of the shell-crossing singularity

For our purpose of studying the structure of an emerging pancake in 2D, we derived the normal form of the first singularity that appears in a local neighborhood at shell crossing. We start by expanding the motion around the shell-crossing location and time (q_c, t_c) into a Taylor series up to O(q³, t),

$$\begin{array}{c}\hfill \mathbf{x}(\mathit{q},t)={\displaystyle \sum _{0\le k\le 1}\sum _{0\le i+j\le 3}}\frac{{\partial }^{i+j+k}\mathbf{x}}{\partial q_x^i\partial q_y^j\partial t^k}{|}_{{\mathit{q}}_c,t_c}{(q_x-q_{x,c})}^i{(q_y-q_{y,c})}^j{(t-t_c)}^k.\end{array}$$(13)

The first assumption made is that of ballistic approximation, that is, expansion in x up to linear order in time. Higher-order expansion in time would require computation of the force field as is done in PCPT. We shall show in a following subsection (see Eq. (29) and the subsequent discussion) that in order to be self-consistent under this assumption, the expansion in velocity must be linear order in space: dx/dt ∼ O(q). The second assumption is that the motion of particles is vorticity-free in Lagrangian space,

$$\begin{array}{c}\hfill \mathbf{x}(\mathit{q},t)={\overline{\mathrm{\nabla }}}_q(q^2/2-\varphi (\mathit{q},t)),\end{array}$$(14)

which is obtained by substituting Eq. (10) in Eq. (4). Thus, the coefficients of our expansion comprise derivatives of ϕ up to fourth order in q and first order in t. This assumption is justifiable for displacement fields derived from 1- or 2-LPT motion. For more generic cases, such as N-body simulations or 3-LPT and beyond, the effect of vorticity in Lagrangian space is non-negligible, but is left for future work (see Sect. 4). We use the following notation for the derivatives of ϕ:

$$\begin{array}{c}\hfill {\varphi }^{i,j,k}(\mathit{q},t)=\frac{{\partial }^{i+j+k}\varphi }{\partial q_x^i\partial q_y^j\partial t^k}(\mathit{q},t).\end{array}$$(15)

In equations where the function input (q, t) has been dropped to make them concise, it is to be assumed that the function is evaluated at the shell-crossing point (q_c, t_c). These coefficients govern the dynamics in the vicinity of the shell-crossing point (q_c, t_c) and therefore determine the morphology of the pancake that forms. The first-order derivatives describe the displacement of the shell-crossing point from the origin. The second-order derivatives correspond to the curvatures of the displacement potential ϕ, i.e., the eigenvalues; their signs indicate whether the resulting structure is a halo or a filament, while their magnitudes control the growth rate of caustics along each axis. The third-order derivatives represent gradients of the eigenvalues and set the relative positions of their peaks. They introduce asymmetries, such as unequal extension along the y-axis or the emergence of a C shape. The fourth-order derivatives describe the curvatures of the eigenvalue field, determining the ellipticity of the emerging pancake and producing features such as an S shape. Based on the geometry of the CDM sheet in the neighborhood of the shell-crossing point (q_c, t_c), we impose constraints on the set of coefficients of our Taylor expansion Eq. (13).

• O(q⁰) derivatives: ϕ^0, 0, 0, ϕ^0, 0, 1 These neither appear in Eq. (13) nor affect the dynamics and hence are immaterial.

• O(q¹) derivatives: ϕ^1, 0, 0, ϕ^0, 1, 0, ϕ^1, 0, 1, ϕ^0, 1, 1 Upon translating on to the shell-crossing point in Lagrangian space q_c = 0 and in Eulerian space x(q_c, t_c) = 0, the O(q¹, t⁰) derivatives vanish: ϕ^1, 0, 0, ϕ^0, 1, 0 = 0. Similarly, upon subtracting the velocity at the shell-crossing point globally from the velocity field such that dx/dt (q_c, t_c) = 0, the O(q¹, t¹) derivatives vanish as well: ϕ^1, 0, 1, ϕ^0, 1, 1 = 0.

• O(q²) derivatives: ϕ^2, 0, 0, ϕ^1, 1, 0, ϕ^0, 2, 0, ϕ^2, 0, 1, ϕ^1, 1, 1, ϕ^0, 2, 1 With the assumption of a vorticity-free displacement field in Lagrangian space, the deformation matrix d_ij = δ_ij^D − J_ij can be rewritten by substituting Eq. (14) in Eq. (5),

  $$\begin{array}{c}\hfill d_{\mathit{ij}}(\mathit{q},t)=\left[\begin{array}{cc}{\varphi }^{2,0}(\mathit{q},t)& {\varphi }^{1,1}(\mathit{q},t)\\ {\varphi }^{1,1}(\mathit{q},t)& {\varphi }^{0,2}(\mathit{q},t)\end{array}\right],\end{array}$$(16)

  which is simply the Hessian of the potential ϕ and has the eigenvalues α > β. Since d_ij is symmetrical, its eigenvectors ${\widehat{n}}^{\alpha },{\widehat{n}}^{\beta }$ are orthogonal (see Eqs. (A.1) and (A.2) for the full expressions of the eigenvalues and eigenvectors). By performing rotation, both in Lagrangian and Eulerian spaces, to align the reference axis ${\widehat{q}}_x$ along the shell-crossing axis ${\widehat{n}}^{\alpha }({\mathit{q}}_c,t_c)$, we diagonalize d_ij(q_c, t_c), which implies ϕ^1, 1, 0 = 0 and ϕ^2, 0, 0 > ϕ^0, 2, 0. The eigenvalues and eigenvectors evaluated at (q_c, t_c) are

  $$\begin{array}{cc}& \alpha ({\mathit{q}}_c,t_c)={\varphi }^{2,0,0}\beta ({\mathit{q}}_c,t_c)={\varphi }^{0,2,0},\hfill \\ \hfill & {\mathit{n}}^{\alpha }({\mathit{q}}_c,t_c)=\left[\begin{array}{c}2({\varphi }^{2,0,0}-{\varphi }^{0,2,0})\\ 0\end{array}\right]{\mathit{n}}^{\beta }({\mathit{q}}_c,t_c)=\left[\begin{array}{c}0\\ 2({\varphi }^{0,2,0}-{\varphi }^{2,0,0})\end{array}\right].\hfill \end{array}$$(17)

  The density ρ(x_c, t_c)∝[(1 − α(q_c, t_c))(1 − β(q_c, t_c))]⁻¹, Eq. (7). The shell-crossing point corresponds to a singularity in the density field, which implies α(q_c, t_c) = ϕ^2, 0, 0 = 1. If we assume Zeldovich flow prior to shell crossing, then ϕ(q, t) = D₊(t)ϕ(q) implies ϕ^i, j, 1 ∝ ϕ^i, j, 0. This suggests that the constraints on the signs of ϕ^i, j, 1 are the same as that of ϕ^i, j, 0, which are ϕ^2, 0, 1 > 0, ϕ^1, 1, 1 = 0, ϕ^2, 0, 1 > ϕ^0, 2, 1. Even for higher-order LPT motion, the constraints on the signs should remain the same, as the magnitude of the corrections would typically be less than that of the leading order. Therefore, the constraints on the signs of ϕ^i, j, 1 are not just limited to Zeldovich flow.

• O(q³) derivatives: ϕ^3, 0, 0, ϕ^2, 1, 0, ϕ^1, 2, 0, ϕ^0, 3, 0 The expression for ${\overline{\mathrm{\nabla }}}_{\mathit{q}}\alpha (\mathit{q},t)$ is provided in the appendix, Eq. (A.3). Evaluating it at (q_c, t_c) using the previously obtained constraints ϕ^1, 1, 0 = 0, ϕ^2, 0, 0 = 1, ϕ^2, 0, 0 > ϕ^0, 2, 0 (see Eq. (109) in Feldbrugge et al. 2023),

  $$\begin{array}{c}\hfill {\overline{\mathrm{\nabla }}}_{\mathit{q}}\alpha ({\mathit{q}}_c,t_c)=\left[\begin{array}{c}{\varphi }^{3,0,0}\\ {\varphi }^{2,1,0}\end{array}\right].\end{array}$$(18)

  Since the shell-crossing point is a local maximum of the eigenvalue α(q, t), ${\overline{\mathrm{\nabla }}}_{\mathit{q}}\alpha ({\mathit{q}}_c,t_c)$ must vanish, which implies ϕ^3, 0, 0 = ϕ^2, 1, 0 = 0. The other two derivatives ϕ^1, 2, 0, ϕ^0, 3, 0 are unconstrained. Since we restrict dx/dt up to O(q) (see Eq. (29) and the subsequent discussion), the O(q^≥3, t¹) derivatives of ϕ are not considered.

• O(q⁴) derivatives: ϕ^4, 0, 0, ϕ^3, 1, 0, ϕ^2, 2, 0, ϕ^1, 3, 0, ϕ^0, 4, 0 The expression for the Hessian of the higher eigenvalue field, defined as H_ij(α(q, t)) = ∂²α(q, t)/∂q_i∂q_j, is provided in the appendix, Eq. (A.4). Evaluating it at (q_c, t_c) using the previously obtained constraints: ϕ^1, 1, 0 = 0, ϕ^2, 0, 0 = 1, ϕ^2, 0, 0 > ϕ^0, 2, 0, ϕ^3, 0, 0 = ϕ^2, 1, 0 = 0 (see Eq. (6.5) in Feldbrugge & van de Weygaert 2023),

  $$\begin{array}{c}\hfill H_{\mathit{ij}}(\alpha ({\mathit{q}}_c,t_c))=\left[\begin{array}{cc}{\varphi }^{4,0,0}& {\varphi }^{3,1,0}\\ {\varphi }^{3,1,0}& {\varphi }^{2,2,0}+2{({\varphi }^{1,2,0})}^2/(1-{\varphi }^{0,2,0})\end{array}\right].\end{array}$$(19)

  It is symmetric but not diagonal, unless ϕ^3, 1, 0 = 0. Consequently, its eigenvectors ${\widehat{n}}_H^{\alpha }({\mathit{q}}_c,t_c)$ and ${\widehat{n}}_H^{\beta }({\mathit{q}}_c,t_c)$ are orthogonal, though not necessarily aligned with the eigenvectors of d_ij(q_c, t_c), namely ${\widehat{n}}^{\alpha }({\mathit{q}}_c,t_c)$ and ${\widehat{n}}^{\beta }({\mathit{q}}_c,t_c)$, which serve as our reference axes. The expressions for the eigenvalues α_H > β_H and corresponding eigenvectors at (q_c, t_c) are provided in the appendix, Eqs. (A.5) and (A.6). For the shell-crossing point to be a non-degenerate local maximum of α(q, t), both the eigenvalues of its Hessian, which are also the curvatures at that point, must be negative. This is possible only if the trace and the determinant of the Hessian are negative and positive,

  $$\begin{array}{cc}& {\varphi }^{4,0,0}+{\varphi }^{2,2,0}+\frac{2{({\varphi }^{1,2,0})}^2}{1-{\varphi }^{0,2,0}}<0,\hfill \\ \hfill & {\varphi }^{4,0,0}({\varphi }^{2,2,0}+\frac{2{({\varphi }^{1,2,0})}^2}{1-{\varphi }^{0,2,0}})>{({\varphi }^{3,1,0})}^2,\hfill \end{array}$$(20)

  which together imply that ϕ^4, 0, 0 < 0 and ϕ^2, 2, 0 < 0. The remaining derivatives ϕ^3, 1, 0, ϕ^1, 3, 0, ϕ^0, 4, 0 are unconstrained. Furthermore, the shell-crossing point (q_c, t_c) automatically satisfies the condition for an A₃⁺ singularity: ${\overline{\mathrm{\nabla }}}_{\mathit{q}}\alpha (\mathit{q},t_c)·{\widehat{n}}^{\alpha }(\mathit{q},t_c)=0$, as well as that for a ridge point: ${\overline{\mathrm{\nabla }}}_{\mathit{q}}\alpha (\mathit{q},t_c)·{\widehat{n}}_H^{\beta }(\mathit{q},t_c)=0$, by virtue of being a local maximum of α(q, t_c) (see Fig. 6). The tangent to the locus of A₃⁺ points at (q_c, t_c) is $-{\varphi }^{3,1,0}{\widehat{n}}^{\alpha }({\mathit{q}}_c,t_c)+{\varphi }^{4,0,0}{\widehat{n}}^{\beta }({\mathit{q}}_c,t_c)$, Eq. (30), which is aligned along neither the reference axis ${\widehat{n}}^{\beta }$ nor the ridge ${\widehat{n}}_H^{\alpha }$ in Lagrangian space unless ϕ^3, 1, 0 = 0. However, it is aligned with the reference axis in Eulerian space. This is because ${\widehat{n}}^{\alpha }({\mathit{q}}_c,t_c)$ gives the direction of shell crossing of particles, or the folding of the CDM sheet, which is quasi-1D for an infinitesimally small neighborhood and time interval. Therefore, the only portion of the CDM sheet in the immediate neighborhood that can continue folding lies along the transverse direction, which is ${\widehat{n}}^{\beta }({\mathit{q}}_c,t_c)$. It is thus more straightforward to characterize the Eulerian shape of the shell-crossing structure in the frame of $({\widehat{n}}^{\alpha }({\mathit{q}}_c,t_c),{\widehat{n}}^{\beta }({\mathit{q}}_c,t_c))$, and therein lies the reason behind our choice of the reference axes.

The final coefficient matrices ϕ^i, j, k after compiling all the constraints we discussed are

$$\begin{array}{cc}& {\varphi }^{i,j,0}=\left[\begin{array}{ccccc}0& 0& (-\infty ,1)& (-\infty ,\infty )& (-\infty ,\infty )\\ 0& 0& (-\infty ,\infty )& (-\infty ,\infty )& \_\\ 1& 0& (-\infty ,0)& \_& \_\\ 0& (-\infty ,\infty )& \_& \_& \_\\ (-\infty ,0)& \_& \_& \_& \_\\ \end{array}\right],\hfill \\ \hfill & {\varphi }^{i,j,1}=\left[\begin{array}{ccc}0& 0& (-\infty ,\infty )\\ 0& 0& \_\\ (0,\infty )& \_& \_\\ \end{array}\right],\hfill \end{array}$$(21)

where the dashes correspond to the higher order derivatives we do not consider in our expansion.

To reiterate, these constraints are derived for shell-crossing points where pancakes emerge, that is, the local maxima of the eigenvalue field α(q, t). The pancakes can be further classified into halos and filaments based on the sign of the lower eigenvalue β(q_c, t_c) = ϕ^0, 2, 0. If ϕ^0, 2, 0 ∈ (0, 1) and ϕ^0, 2, 1 ∈ (0, ∞), shell crossing occurs along the transverse axis ${\widehat{n}}^{\beta }({\mathit{q}}_c,t_c)$, and the pancake collapses into a 2D halo. Conversely, if ϕ^0, 2, 0, ϕ^0, 2, 1 ∈ ( − ∞, 0], the particles continue to drift apart along the transverse axis without undergoing a second shell crossing, resulting in a 2D filament.

Extension to local minima and saddle points of the eigenvalue field α(q, t) is obtained simply by changing the sign constraints on the eigenvalues α_H > β_H of the Hessian H_ij(α(q_c, t_c)). For saddle points, α_H(q_c, t_c) > 0 and β_H(q_c, t_c) < 0, which is possible only when the determinant of the Hessian is negative. If ϕ^0, 2, 0 > 0, the saddle point undergoes shell crossing along the transverse axis by merging into a halo. If ϕ^0, 2, 0 ≤ 0, the saddle point does not undergo further shell crossing and becomes part of a filament. For local minima, α_H(q_c, t_c), β_H(q_c, t_c) > 0, which is possible only when both the trace and the determinant of the Hessian are positive.

After plugging in the constraints, Eq. (21), into our expansion Eq. (13), we obtain the following reduced expression with a total of ten remaining coefficients:

$$\begin{array}{cc}\hfill x(\mathit{q},t)& =-{\varphi }^{2,0,1}{q}_x(t-t_c)-\frac{1}{2}{\varphi }^{1,2,0}{q}_y^2-\frac{1}{6}{\varphi }^{4,0,0}{q}_x^3\hfill \\ \hfill & -\frac{1}{2}{\varphi }^{3,1,0}{q}_x^2{q}_y-\frac{1}{2}{\varphi }^{2,2,0}{q}_x{q}_y^2-\frac{1}{6}{\varphi }^{1,3,0}{q}_y^3,\hfill \\ \hfill y(\mathit{q},t)& =(1-{\varphi }^{0,2,0})q_y-{\varphi }^{0,2,1}{q}_y(t-t_c)-\frac{1}{2}{\varphi }^{0,3,0}{q}_y^2\hfill \\ \hfill & -{\varphi }^{1,2,0}{q}_x{q}_y-\frac{1}{6}{\varphi }^{0,4,0}{q}_y^3-\frac{1}{2}{\varphi }^{1,3,0}{q}_x{q}_y^2\hfill \\ \hfill & -\frac{1}{2}{\varphi }^{2,2,0}{q}_x^2{q}_y-\frac{1}{6}{\varphi }^{3,1,0}{q}_x^3.\hfill \end{array}$$(22)

In terms of catastrophe theory, these are the normal form derivatives corresponding to the subset of A₃⁺ singularities that appear at the local maxima of the eigenvalue α(q, t).

2.3. Pancake structure

Eq. (22) governs the motion of particles in the neighborhood of the shell-crossing location q_c = 0 and time t_c, from which we can derive the loci of A₃⁺ and A₂ singularities that define the structure of the emerging pancake as shown in Fig. 1.

The deformation matrix d_ij = δ_ij^D − J_ij can be rewritten substituting our expansion form, Eq. (22), into Eq. (5). Its elements, up to O(q², t), are

$$\begin{array}{c}\hfill \left[\begin{array}{cc}1+{\varphi }^{2,0,1}(t-t_c)+\frac{1}{2}{\varphi }^{4,0,0}{q}_x^2& {\varphi }^{1,2,0}{q}_y+\frac{1}{2}{\varphi }^{3,1,0}{q}_x^2\\ +{\varphi }^{3,1,0}{q}_x{q}_y+\frac{1}{2}{\varphi }^{2,2,0}{q}_y^2& +{\varphi }^{2,2,0}{q}_x{q}_y+\frac{1}{2}{\varphi }^{1,3,0}{q}_y^2\\ & \\ {\varphi }^{1,2,0}{q}_y+\frac{1}{2}{\varphi }^{3,1,0}{q}_x^2& {\varphi }^{0,2,0}+{\varphi }^{0,2,1}(t-t_c)\\ +{\varphi }^{2,2,0}{q}_x{q}_y+\frac{1}{2}{\varphi }^{1,3,0}{q}_y^2& +{\varphi }^{0,3,0}{q}_y+{\varphi }^{1,2,0}{q}_x+{\varphi }^{1,3,0}{q}_x{q}_y\\ & +\frac{1}{2}{\varphi }^{0,4,0}{q}_y^2+\frac{1}{2}{\varphi }^{2,2,0}{q}_x^2\end{array}\right].\end{array}$$(23)

Its eigenvalues α > β and corresponding eigenvectors ${\widehat{n}}^{\alpha },{\widehat{n}}^{\beta }$ are expanded up to O(q², t) to be consistent,

$$\begin{array}{cc}\hfill \alpha (\mathit{q},t)& =1+{\varphi }^{2,0,1}(t-t_c)+\frac{1}{2}{\varphi }^{4,0,0}{q}_x^2+{\varphi }^{3,1,0}{q}_x{q}_y\hfill \\ \hfill & +\frac{1}{2}({\varphi }^{2,2,0}+\frac{2{({\varphi }^{1,2,0})}^2}{(1-{\varphi }^{0,2,0})})q_y^2,\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill {\widehat{n}}^{\alpha }(\mathit{q},t)& =\left[\begin{array}{c}1-{\varphi }^{0,2,0}+({\varphi }^{2,0,1}-{\varphi }^{0,2,1})(t-t_c)-{\varphi }^{0,3,0}{q}_y\\ -{\varphi }^{1,2,0}{q}_x+({\varphi }^{3,1,0}-{\varphi }^{1,3,0})q_x{q}_y-\frac{1}{2}({\varphi }^{2,2,0}-{\varphi }^{4,0,0})q_x^2\\ -\frac{1}{2}({\varphi }^{0,4,0}-{\varphi }^{2,2,0}-\frac{2{({\varphi }^{1,2,0})}^2}{(1-{\varphi }^{0,2,0})})q_y^2\\ \\ \\ {\varphi }^{1,2,0}{q}_y+\frac{1}{2}{\varphi }^{3,1,0}{q}_x^2+{\varphi }^{2,2,0}{q}_x{q}_y+\frac{1}{2}{\varphi }^{1,3,0}{q}_y^2\end{array}\right].\hfill \end{array}$$(25)

The higher the value of 1 − ϕ^0, 2, 0, the faster the convergence of these expansions. Thus, our formalism works better for collapse scenarios that have a greater difference between the eigenvalues, and hence, are closer to being quasi-1D. Substituting Eq. (24) in Eq. (8), we obtain the expression for the locus of A₂ singularities at a given instant t > t_c up to O(q², t) (see Eq. (17) in Arnold et al. 1982):

$$\begin{array}{c}\hfill {\varphi }^{4,0,0}{q}_x^2+2{\varphi }^{3,1,0}{q}_x{q}_y+({\varphi }^{2,2,0}+\frac{2{({\varphi }^{1,2,0})}^2}{(1-{\varphi }^{0,2,0})})q_y^2=-2{\varphi }^{2,0,1}(t-t_c),\end{array}$$(26)

which is of the form Aq_x² + Bq_xq_y + Cq_y² + F = 0. It represents an ellipse if 4AC − B² > 0. Since F > 0 for t > t_c, the ellipse is real only if A, C < 0. These conditions are equivalent to Eq. (20) for the shell-crossing point to be a local maximum of the eigenvalue field α(q, t). For A, C < 0,

$$\begin{array}{cc}& 2{\theta }_{A2}=arctan\left(\frac{B}{|C-A|}\right)=arctan\left(\frac{2{\varphi }^{3,1,0}}{|{\varphi }^{4,0,0}-{\varphi }^{2,2,0}-\frac{2{({\varphi }^{1,2,0})}^2}{(1-{\varphi }^{0,2,0})}|}\right),\hfill \\ \hfill & a={\left(\frac{2F(-(A+C)+\sqrt{{(C-A)}^2+B^2})}{4AC-B^2}\right)}^{1/2}\propto {(t-t_c)}^{1/2}\hfill \end{array}$$(27)

$$\begin{array}{cc}& b={\left(\frac{2F(-(A+C)-\sqrt{{(C-A)}^2+B^2})}{4AC-B^2}\right)}^{1/2}\propto {(t-t_c)}^{1/2},\hfill \end{array}$$(28)

where θ_A2 is the angle by which the axes of the A₂ ellipse $({\widehat{n}}_H^{\alpha },{\widehat{n}}_H^{\beta })$ are rotated with respect to the reference frame $({\widehat{n}}^{\alpha },{\widehat{n}}^{\beta })$ at the shell-crossing point (q_c, t_c). The semi-major and semi-minor axes, of lengths a and b, are oriented along the directions of less and more negative curvature of α at (q_c, t_c), denoted by ${\widehat{n}}_H^{\alpha }$ and ${\widehat{n}}_H^{\beta }$, respectively, and satisfy 0 < b/a ≤ 1. Consequently, the α-ridge aligns with the semi-major axis of the A₂ ellipse. For the case where ϕ^3, 1, 0 = B = θ_A2 = 0, that is, the axes of the A₂ ellipse being aligned with the reference frame, the shell-crossing direction ${\widehat{n}}^{\alpha }$ coincides with the minor axis (more negative curvature) if A < C < 0. Interestingly, if C < A < 0, it coincides with the major axis (less negative curvature).

From Eq. (28), all points on the A₂ ellipse satisfy |q_A2 − q_c|∝(t − t_c)^1/2. Using the parameterization h = t − t_c and q_A2 ∝ h^1/2 in Eqs. (22) to obtain the Eulerian shape of the A₂ caustic,

$$\begin{array}{cc}\hfill x_{A2}(h)& =Ph+Q{h}^{3/2}\hfill \\ \hfill y_{A2}(h)& =R{h}^{1/2}+Sh+T{h}^{3/2}.\hfill \end{array}$$(29)

The A₂ caustic demarcates the multi-streaming region in Eulerian space, within which our expansion form of the motion around the shell-crossing point, Eq. (22), remains valid. Therefore, Eq. (22) is consistent up to O(h^3/2). Terms of O(q²) in dx/dt would have introduced some, but not all the terms of O(h²) in x_A2 leading to inconsistency. Hence, they were neglected.

For saddle points of the eigenvalue field α(q, t), the determinant of the Hessian Eq. (19) is negative, which implies 4AC − B² < 0. Equation (26) then represents a hyperbola. For points of minima, the determinant is positive ⟹4AC − B² > 0. However, the trace A + C being positive implies A, C > 0 for which Eq. (26) does not have real solutions after shell crossing t > t_c.

Returning to our shell-crossing point, we substitute Eqs. (24) and (25) in Eq. (9) and keep only the terms up to O(q, t) to obtain (see Eq. (5.12) in Feldbrugge & van de Weygaert 2023),

$$\begin{array}{cc}& \left[\begin{array}{c}{\varphi }^{4,0,0}{q}_x+{\varphi }^{3,1,0}{q}_y\\ {\varphi }^{3,1,0}{q}_x+({\varphi }^{2,2,0}+\frac{2{({\varphi }^{1,2,0})}^2}{(1-{\varphi }^{0,2,0})})q_y\end{array}\right]·{\widehat{n}}^{\alpha }=0\hfill \\ \hfill & \Rightarrow {\varphi }^{4,0,0}{q}_x+{\varphi }^{3,1,0}{q}_y=0.\hfill \end{array}$$(30)

This is the locus of A₃⁺ points in Lagrangian space, regardless of the time when they turn singular. It is a straight line at an angle θ_A3 with respect to the reference axis ${\widehat{q}}_y\equiv {\widehat{n}}^{\beta }({\mathit{q}}_c,t_c)$

$$\begin{array}{c}\hfill {\theta }_{A3}=arctan(-{\varphi }^{3,1,0}/{\varphi }^{4,0,0}),\end{array}$$(31)

and is time-independent for our expansion form, Eq. (22), of the order O(q³, t). Higher order expansions would indeed introduce O(q^≥2, t^≥1) corrections (see Fig. 6 to discern the alignments between A₂ ellipse, A₃⁺ line, and the reference axes).

To obtain the Eulerian shape of the A₃⁺ line, we substitute Eq. (30) in our expansion form Eq. (22). The expressions for x_A3(q, t) and y_A3(q, t) in parametric form are provided in the appendix, Eq. (A.7). It is to be noted that Eq. (30) is of O(q) and Eq. (22) is of O(q³). So, upon substitution, Eq. (A.7) for the A₃⁺ spine is strictly correct only up to O(q). To have the complete expression of Eq. (A.7) correct up to O(q³), we need to start with expansion form, Eq. (22), of the order O(q⁵) that will result in Eq. (30) of the order O(q³), that is, quadratic and cubic corrections to the A₃⁺ line in Lagrangian space. An O(q⁵) expansion in space necessitates O(t²) in time so that our expansion form is consistent up to O(h^5/2) (see Eq. (29)). However, O(t²) expansion requires PCPT for force field corrections, which is beyond the scope of this work. Nevertheless, only at shell crossing (t = t_c), the O(q) term in expansion of x(q, t), Eq. (22), vanishes and hence, the A₃⁺ spine trajectory in Eulerian space,

$$\begin{array}{cc}\hfill x_{A3}(y_{A3},t_c)=& -\left[\frac{{\varphi }^{1,2,0}}{2{(1-{\varphi }^{0,2,0})}^2}\right]y_{A3}^2\hfill \\ \hfill & -[\frac{{\varphi }^{1,3,0}+{\varphi }^{3,1,0}(2{\left(\frac{{\varphi }^{3,1,0}}{{\varphi }^{4,0,0}}\right)}^2-3\frac{{\varphi }^{2,2,0}}{{\varphi }^{4,0,0}})}{6{(1-{\varphi }^{0,2,0})}^3}\hfill \\ \hfill & +\frac{{\varphi }^{1,2,0}}{{(1-{\varphi }^{0,2,0})}^4}(\frac{1}{2}{\varphi }^{0,3,0}-\frac{{\varphi }^{1,2,0}{\varphi }^{3,1,0}}{{\varphi }^{4,0,0}})]y_{A3}^3,\hfill \end{array}$$(32)

is indeed correct up to O(y²). The only missing term is at the order O(y³) and involves O(q⁵) derivative of the potential ϕ. On physical grounds, we can argue that if the field ϕ is Gaussian distributed, its O(q⁵) derivatives and their effects on the pancake shape are typically suppressed compared to those of the lower order derivatives for a relevant range of power spectra (see Fig. B.1 and the discussion therein). Equation (32) is therefore an accurate enough approximation at shell crossing, which is where we focus the derivation of expressions for properties related to the pancake shape.

We can transform the intersection of Eqs. (26) and (30) through Eq. (A.7) to obtain the Eulerian location of the tips of the growing pancake at instants shortly after shell crossing t > t_c. From Eq. (A.7), the leading order term O(q_y) in x_A3(q, t) is time dependent suggesting that the spine gradually rotates. At shell crossing t = t_c, however, it vanishes, which means the spine is aligned with our reference axis ${\widehat{q}}_y\equiv {\widehat{n}}^{\beta }({\mathit{q}}_c,t_c)$. The order O(q_y²) term depends on ϕ^1, 2, 0 and contributes to the C shape of the spine. Last, the order O(q_y³) term depends on ϕ^1, 3, 0 and ϕ^3, 1, 0, and contributes to the S shape of the spine.

From Eq. (32), we compute a few properties related to the shape of the pancake at shell crossing:

1. The curvature of the A₃⁺ spine,

   $$\begin{array}{c}\hfill \frac{d^2{x}_{A3}}{d{y}_{A3}^2}{|}_{q_y=0}(t_c)=\frac{{\varphi }^{1,2,0}}{{(1-{\varphi }^{0,2,0})}^2}·\end{array}$$(33)

   Doroshkevich & Shandarin (1978, in Sect. 5.3) approximated the mean curvature of the lateral surfaces (equivalent to the A₂ boundary in 2D) of the pancake in 3D Eulerian space, assuming it to be infinitely thin and symmetric. It is different from the A₃⁺ spine’s curvature Eq. (33), which is related to the leading order correction to the presumed symmetrical shape of the pancake.

2. From Eq. (32), we note that close to the shell-crossing point, the A₃⁺ spine is C-shaped given the O(y²) term. At larger scales, it turns S-shaped as the O(y³) term starts to dominate. The transition scale provides an estimate of the extent up to which the pancake can be considered C-shaped, which we estimate using the point y at which the spine starts to bend from C to S in Eulerian space at shell crossing, that is dx_A3/dy_A3 (t_c) = 0, which implies

   $$\begin{array}{cc}\hfill y_{C\to S}(t_c)& =-2(1-{\varphi }^{0,2,0}){\varphi }^{1,2,0}\times \hfill \\ \hfill & [{\varphi }^{1,3,0}+{\varphi }^{3,1,0}(2{\left(\frac{{\varphi }^{3,1,0}}{{\varphi }^{4,0,0}}\right)}^2-3\frac{{\varphi }^{2,2,0}}{{\varphi }^{4,0,0}})\hfill \\ \hfill & +\frac{6{\varphi }^{1,2,0}}{(1-{\varphi }^{0,2,0})}(\frac{1}{2}{\varphi }^{0,3,0}-\frac{{\varphi }^{1,2,0}{\varphi }^{3,1,0}}{{\varphi }^{4,0,0}}){]}^{-1}.\hfill \end{array}$$(34)

   For halos, ϕ^0, 2, 0, ϕ^0, 2, 1 > 0, which means the transverse axis contracts. Thus, the curvature increases and the y_C → S decreases with time. On the other hand, for filaments, ϕ^0, 2, 0, ϕ^0, 2, 1 ≤ 0, which means the transverse axis stretches, leading to decreasing curvature and increasing y_C → S.

2.4. Parameters and their effects

The shape of the shell-crossing structure in Lagrangian space is captured by the loci of A₂ and A₃⁺ points, which are governed by Eqs. (26) and (30), respectively. Their transformation to Eulerian space is governed by Eq. (22). These three sets of equations constitute our analytical model and contain ten parameters in total. Constraints on the signs of some of these parameters are listed in Eq. (21). We further isolate and investigate the effect of each of these parameters on the shape of the pancake. The values of parameters used in Figs. 2, 4, 3, and 5 are only for the purpose of presentation.

• ϕ^2,0,1: velocity gradient along shell-crossing axis. It determines how quickly the particles cross trajectories and turn singular, and hence, the rate at which the A₂ caustic grows in Lagrangian, Eq. (26), and Eulerian spaces.

• ϕ^2,0,0, ϕ^0,2,1: displacement and velocity gradient along the transverse axis. They determine if the pancake stretches or shrinks longitudinally in Eulerian space, Eq. (22).

• ϕ^4,0,0, ϕ^2,2,0: they are related to the two eigenvalues (the principal curvatures) of H_ij(α(q_c, t_c)), Eq. (19). They determine the axis ratio of the A₂ caustic, Eq. (28), which in turn affects how flat the pancake is in Eulerian space.

• ϕ^1,2,0, ϕ^1,3,0: they affect the shape of the A₃ spine Eq. (32) in Eulerian space. The former is responsible for the C shape, and the latter is responsible for the S shape.

• ϕ^0,3,0, ϕ^0,4,0: the former adds asymmetric and the latter adds symmetric higher order corrections to the longitudinal stretch of the pancake in Eulerian space, Eq. (22).

• ϕ^3,1,0: leads to rotation of the A₃⁺ spine at leading order and contributes to its S shape at order O(y³) in Eulerian space, Eq. (32).

Fig. 2. Effect of varying ϕ2, 0, 1 (left) and ϕ0, 2, 0 (right) on the shape of A2 caustic in Eulerian space without changing the other parameters and the time.

Fig. 3. Effect of varying the ratio ϕ2, 2, 0/ϕ4, 0, 0 on the shape of A2 caustic in Lagrangian (left) and Eulerian (right) spaces without changing the other parameters and the time.

Fig. 4. Effect of varying ϕ1, 2, 0 (left) and ϕ1, 3, 0 (right) on the shape of A2 caustic and A3 spine in Eulerian space without changing the other parameters and the time.

Fig. 5. Effect of varying ϕ0, 3, 0 (left) and ϕ0, 4, 0 (right) on the shape of A2 caustic in Eulerian space without changing the other parameters and the time.

Fig. 6. Effect of ϕ3, 1, 0 on the shape of the A2 caustic (solid) and the A3+ spine (dashed) in Lagrangian (blue) and Eulerian (red) spaces. The angles θA2, Eq. (27), and θA3, Eq. (31), are indicated in green and blue, respectively. The vector fields of the eigenvector ${\widehat{n}}^{\alpha }$ and the gradient ∇qα are shown as light and dark gray arrows, respectively.

3. Gaussian statistics

Using arguments from catastrophe theory, we modeled how the derivatives of the displacement potential ϕ at the shell-crossing point (q_c, t_c) relate to the shape of the emerging pancake. There are ten such derivatives, or, as we prefer to call them, the parameters in our model that appear in our expansion form, Eq. (22). In this section, starting from a Gaussian random field for the displacement potential ϕ(q) and its derivatives, we compute the conditional and marginal probabilities of our model parameters subject to the constraints of a shell-crossing point, and obtain the distribution and expectation values of the features defining the shape of the pancakes.

We first construct our random variables. As discussed earlier, the shape of the pancake is captured by the A₂ caustic and the A₃ spine. Properties like the curvature, Eq. (33), and the transition scale from C to S shape depend on O(q²),O(q³) and O(q⁴) derivatives of the displacement potential ϕ at shell crossing. For the purpose of computing statistics, we assume Zeldovich flow (ZA) before shell crossing. Coles et al. (1993) and Melott et al. (1994) showed that the ZA is a good approximation down to the nonlinearity scale, by which stage most pancakes have already formed yet remain relatively thin. So, ZA can be safely used to capture the time-evolution of the potential ϕ up to shell crossing: ϕ(q, t) = D₊(t)ϕ(q). With this approximation, the time derivatives ϕ^i, j, 1s are proportional to ϕ^i, j, 0s, which makes their constraints redundant in the computation of conditional probabilities at the shell-crossing point. For our random variables, it is thus sufficient to consider the set of O(q²),O(q³) and O(q⁴) spatial derivatives for computing the properties we are interested in. The Zeldovich approximation (ZA) also allows for a rescaling of these derivatives, thereby eliminating the explicit time dependence,

$$\begin{array}{c}\hfill {\varphi }^{i,j}(\mathit{q})=D_{+}^{-1}(t){\varphi }^{i,j,0}(\mathit{q},t).\end{array}$$(35)

We define two random variables: $\stackrel{\sim }{Y}$, which comprises the derivatives at a point in the generic frame, and Y, which comprises the derivatives in the reference frame (same as used in Sect. 2) rotated to align with the eigenvectors ${\widehat{n}}^{\alpha },{\widehat{n}}^{\beta }$ at that point, Eq. (A.2),

$$\begin{array}{cc}& \tilde{Y}=[{\tilde{\varphi }}^{2,0},{\tilde{\varphi }}^{1,1},{\tilde{\varphi }}^{0,2},{\tilde{\varphi }}^{3,0},{\tilde{\varphi }}^{2,1},{\tilde{\varphi }}^{1,2},{\tilde{\varphi }}^{0,3},{\tilde{\varphi }}^{4,0},{\tilde{\varphi }}^{3,1},{\tilde{\varphi }}^{2,2},{\tilde{\varphi }}^{1,3},{\tilde{\varphi }}^{0,4}],\hfill \\ \hfill & Y=[{\varphi }^{2,0},\theta ,{\varphi }^{0,2},{\varphi }^{3,0},{\varphi }^{2,1},{\varphi }^{1,2},{\varphi }^{0,3},{\varphi }^{4,0},{\varphi }^{3,1},{\varphi }^{2,2},{\varphi }^{1,3},{\varphi }^{0,4}],\hfill \end{array}$$(36)

where θ is the rotation angle between the two frames, and it diagonalizes the deformation matrix:

$$\begin{array}{c}\hfill \left[\begin{array}{cc}{\tilde{\varphi }}^{2,0}& {\tilde{\varphi }}^{1,1}\\ {\tilde{\varphi }}^{1,1}& {\tilde{\varphi }}^{0,2}\end{array}\right]=\left[\begin{array}{cc}cos\theta & sin\theta \\ -sin\theta & cos\theta \end{array}\right]\left[\begin{array}{cc}{\varphi }^{2,0}& 0\\ 0& {\varphi }^{0,2}\end{array}\right]\left[\begin{array}{cc}cos\theta & -sin\theta \\ sin\theta & cos\theta \end{array}\right].\end{array}$$(37)

The transformations of the O(q³) and O(q⁴) derivatives under rotation, as well as the Jacobian $d{\stackrel{\sim }{Y}}_i/d{Y}_j$, are given in the appendix (Eqs. (B.1) and (B.2); see also Eqs. (E.3)–(E.5) in Feldbrugge & van de Weygaert 2023). The determinant of the Jacobian is

$$\begin{array}{c}\hfill \left||\frac{d\tilde{Y}}{\mathit{dY}}\right||=|{\varphi }^{2,0}-{\varphi }^{0,2}|.\end{array}$$(38)

We know that $\stackrel{\sim }{\varphi }(\mathit{q})$ is Gaussian distributed with the power spectrum $P_{\stackrel{\sim }{\varphi }}(k)$, which is related to the matter power spectrum P_δ(k) through Eq. (11). Therefore, $\stackrel{\sim }{Y}$ comprising derivatives of ϕ(q) is Gaussian distributed as well. Its probability density is

$$\begin{array}{c}\hfill P(\tilde{Y})=\frac{1}{\sqrt{{(2\pi )}^N|M|}}exp(-\frac{1}{2}{\tilde{Y}}^T{M}^{-1}\tilde{Y}),\end{array}$$(39)

where N is the number of elements in $\stackrel{\sim }{Y}$ and M is the covariance matrix between $\stackrel{\sim }{Y}$ and Y, given by (see also Eqs. (57)–(65) in Feldbrugge & van de Weygaert 2025),

$$\begin{array}{cc}\hfill M_{2,2}& =\frac{{\sigma }_2^2}{8}\left[\begin{array}{ccc}3& 0& 1\\ 0& 1& 0\\ 1& 0& 3\end{array}\right],\hfill \\ \hfill M_{2,4}& =-\frac{{\sigma }_3^2}{16}\left[\begin{array}{ccccc}5& 0& 1& 0& 1\\ 0& 1& 0& 1& 0\\ 1& 0& 1& 0& 5\end{array}\right],\hfill \\ \hfill M_{3,3}& =\frac{{\sigma }_3^2}{16}\left[\begin{array}{cccc}5& 0& 1& 0\\ 0& 1& 0& 1\\ 1& 0& 1& 0\\ 0& 1& 0& 5\end{array}\right],\hfill \\ \hfill M_{4,4}& =\frac{{\sigma }_4^2}{128}\left[\begin{array}{ccccc}35& 0& 5& 0& 3\\ 0& 5& 0& 3& 0\\ 5& 0& 3& 0& 5\\ 0& 3& 0& 5& 0\\ 3& 0& 5& 0& 35\end{array}\right],\hfill \end{array}$$(40)

where M_i, j denotes the covariance matrix between O(qⁱ) and O(q^j) derivatives of the potential ϕ. The O(q³) derivatives have no correlation with O(q²) or O(q⁴) derivatives. σ_i refers to the ith moment of the power spectrum $P_{\stackrel{\sim }{\varphi }}(k)$,

$$\begin{array}{cc}\hfill {\sigma }_i^2& =\frac{1}{2\pi }{\displaystyle {\int }_0^{\infty }}{k}^{2i+1}{P}_{\tilde{\varphi }}(k)dk\hfill \\ \hfill \propto & \frac{1}{R_c^{n+2(i-1)}}\{\begin{array}{c}\mathrm{\Gamma }(n+2(i-1));\mathrm{Exponential}\mathrm{filter}{P}_{\tilde{\varphi }}\propto k^{n-4}{e}^{-k{R}_c}\hfill \\ \mathrm{\Gamma }(n/2+i-1);\mathrm{Gaussian}\mathrm{filter}{P}_{\tilde{\varphi }}\propto k^{n-4}{e}^{-k^2{R}_c^2},\hfill \end{array}\hfill \end{array}$$(41)

where n and R_c are the spectral index and smoothing scale, respectively. These σ_i encapsulate the dependence on the power spectrum parameters n, R_c as well as the smoothing filter used, thereby capturing all the statistical information about the initial potential field ϕ(q) that we require for our computations. Their dependence on the spectral index n is shown in Fig. B.1.

To obtain P(Y), we use Eqs. (39) and (38),

$$\begin{array}{cc}& P(Y)dY=P(\tilde{Y})d\tilde{Y}\hfill \\ \hfill \Rightarrow & P(Y)=\parallel \frac{d\tilde{Y}}{\mathit{dY}}\parallel P(\tilde{Y})\hfill \\ \hfill \Rightarrow & P(Y)\propto |{\varphi }^{2,0}-{\varphi }^{0,2}|exp(-\frac{1}{2}{\tilde{Y}}^T(Y)M^{-1}\tilde{Y}(Y)).\hfill \end{array}$$(42)

The full expression after substituting $\stackrel{\sim }{Y}$ in terms of Y, Eq. (B.1), and the covariance matrix M, Eq. (40) is provided in the appendix, Eq. (B.3). This is the probability density of Y at a generic spatial point q. We are interested in the conditional probability of Y given that it is a shell-crossing point q_c, that is, a peak in the eigenvalue field α(q) which seeds a pancake. For this, we take an approach similar to the computation of the number density of peaks in Bardeen et al. (1986) (see the constraints, Eq. (21)).

• At shell crossing, the higher eigenvalue α(q_c, t_c) = ϕ^2, 0, 0 is equal to 1, which implies ϕ^2, 0 = D₊⁻¹(t_c) > 0. For pancakes, referring to the total set of halos and filaments, it is sufficient to have α(q_c, t_c) > β(q_c, t_c)⟹ϕ^2, 0 > ϕ^0, 2. For the subset of halos, ϕ^0, 2 must be further restricted to ϕ^0, 2 > 0 to have shell crossing along the second axis as well. On the other hand, for the subset of filaments, ϕ^0, 2 ≤ 0 so that the second axis does not collapse. All these constraints are introduced as Heaviside functions.

• The shell-crossing point is a local maximum in the eigenvalue field α(q, t). This implies constraints on its gradient, ϕ^3, 0, ϕ^2, 1 = 0, which appear as Dirac-delta functions, and on the eigenvalues of its Hessian, Eq. (20), which appear as Heaviside functions in the conditional probability.

The conditional probability of Y can be expressed as

$$\begin{array}{cc}\hfill P(Y|& \mathrm{Pancake})\propto P(Y){\delta }_D({\varphi }^{3,0}){\delta }_D({\varphi }^{2,1})\hfill \\ \hfill & H[{\varphi }^{2,0}]H[{\varphi }^{2,0}-{\varphi }^{0,2}]H[-{\varphi }^{4,0}-{\varphi }^{2,2}-\frac{2{({\varphi }^{1,2})}^2}{{\varphi }^{2,0}-{\varphi }^{0,2}}]\hfill \\ \hfill & H[{\varphi }^{4,0}({\varphi }^{2,2}+\frac{2{({\varphi }^{1,2})}^2}{{\varphi }^{2,0}-{\varphi }^{0,2}})-{({\varphi }^{3,1})}^2],\hfill \end{array}$$(43)

where H refers to the Heaviside function. For halo and filament populations, we need to further multiply H[ϕ^0, 2] and H[−ϕ^0, 2], respectively.

3.1. Marginal probabilities

The conditional probability Eq. (43) is completely parameterized by the moments of power spectrum σ₂, σ₃, and σ₄ through Eq. (B.3). The choice of the power spectrum parameters n, R_c, and the smoothing filter affects the probability distribution only through the values of σ. If ϕ^i, j is uncorrelated with other derivatives, the expectation value ⟨|ϕ^i, j|⟩ is proportional to σ_i + j. Even for the correlated ones, ϕ^i, j / σ_i + j is only weakly dependent on the spectral index n (see Fig. B.2 and the discussion therein). It is thus useful to present the distributions in terms of dimensionless quantities that have been scaled by their corresponding σ_i + j, so that they are effectively independent of the length scale R_c, the power spectrum index n, and the smoothing filter.

From Eq. (43), we note that the marginal distribution of θ is simply constant since P(Y) Eq. (B.3) is independent of θ, and those of ϕ^2, 1 and ϕ^3, 0 are Dirac-delta functions centered at zero. So, we integrate them over to obtain

$$\begin{array}{cc}& P({\varphi }^{2,0},{\varphi }^{0,2},{\varphi }^{1,2},{\varphi }^{0,3},{\varphi }^{4,0},{\varphi }^{3,1},{\varphi }^{2,2},{\varphi }^{1,3},{\varphi }^{0,4}|\mathrm{Pancake})\propto \hfill \\ \hfill & \mathrm{Exp}[-10{\left(\frac{{\varphi }^{1,2}}{{\sigma }_3}\right)}^2-2{\left(\frac{{\varphi }^{0,3}}{{\sigma }_3}\right)}^2-\frac{1}{2}({\displaystyle \sum _{\mu ,\nu }}{f}_{\mu \nu }({\sigma }_2,{\sigma }_3,{\sigma }_4){\varphi }_{\mu }{\varphi }_{\nu })]\hfill \\ \hfill & |{\varphi }^{2,0}-{\varphi }^{0,2}|H[{\varphi }^{2,0}]H[{\varphi }^{2,0}-{\varphi }^{0,2}]H[-{\varphi }^{4,0}-{\varphi }^{2,2}-\frac{2{({\varphi }^{1,2})}^2}{{\varphi }^{2,0}-{\varphi }^{0,2}}]\hfill \\ \hfill & H[{\varphi }^{4,0}({\varphi }^{2,2}+\frac{2{({\varphi }^{1,2})}^2}{{\varphi }^{2,0}-{\varphi }^{0,2}})-{({\varphi }^{3,1})}^2],\hfill \end{array}$$(44)