© The Authors 2025

1. Introduction

In the concordant model of cosmology, most of the mass of the Universe ∼84% (see, e.g., Planck Collaboration VI 2020, and references therein) is attributed to cold dark matter (CDM), a non-relativistic collisionless fluid with negligible initial velocity dispersion that exhibits gravitational interactions. The large-scale structure of the Universe is therefore dominated by the clustering of dark matter, with collapsed halos serving as the basic units of cosmological structures. The potential wells formed by these halos aid in the clustering of baryons, making them the birthing ground for stars and galaxies. The halo profile, crucial for modeling observables such as the annihilation signal of WIMPs (e.g., Slatyer 2022), is very sensitive to the dynamics and the nature of dark matter particles. Thus, the study of dark matter halos and associated dynamics has important implications for dark matter indirect detection, the nature of dark matter particles, cosmological structure as well as galaxy formation and evolution.

The complex and multifaceted dynamics of dark matter particles is challenging to capture into a single-encompassing theory. Several analytical theories, each with its own set of assumptions and limitations, have been suggested for each stage of dynamics, allowing us to piece together a sketch of CDM structure formation in our Universe. It begins with density perturbations set at the end of inflation that grow linearly, with particle motion dominated by Hubble expansion. After reaching a critical density, the particles decouple from the expansion, turn around and start collapsing. The single stream flow until the first shell-crossing (the intersection of particle trajectories) can be accurately captured by the Eulerian or Lagrangian perturbation theory (see, e.g., Bernardeau et al. 2002, and references therein). For simplified initial conditions composed of three crossed sine waves, Lagrangian perturbation theory (LPT) has also been shown to predict the crossing time and structure of phase space and caustics (the loci of points where particle trajectories cross) (see, e.g., Saga et al. 2022). After shell-crossing, the flow inside the splashback radius (the radius of the outermost caustic) turns multistream and LPT loses its validity. Nevertheless, the motion up to the next crossing time can still be followed by a post-collapse perturbative approach (Colombi 2015; Taruya & Colombi 2017; Rampf et al. 2021; Saga et al. 2023). Subsequent shell-crossings along transverse directions and violent relaxation in the multistream regime (Lynden-Bell 1967) lead to the formation of primordial halos. In N-body simulations, these primordial halos are found to have a cusp-like density profile at the center obeying the approximate power law: ρ∝r^−1.5 (e.g., Diemand et al. 2005; Ishiyama et al. 2010; Angulo et al. 2017; Colombi 2021; Delos & White 2022). Accretion and mergers further drive the evolution of these halos towards a dynamical attractor, the well-known universal NFW profile (Navarro et al. 1996, 1997) or its recent improvements (e.g., Einasto profile, see Einasto 1965; Navarro et al. 2004).

Despite multiple attempts, there exists no complete analytical theory capable of fully describing the multistream dynamics and accurately predicting the slope and evolution of the density profile of CDM halos. One still promising class of such analytical approaches is to invoke self-similarity (e.g., Fillmore & Goldreich 1984; Bertschinger 1985; White & Zaritsky 1992; Henriksen & Widrow 1995; Sikivie et al. 1997; Nusser 2001; Zukin & Bertschinger 2010a, b; Lithwick & Dalal 2011; Alard 2013; Sugiura et al. 2020, but this list is not exhaustive). For a monolithically growing self-gravitating collisionless CDM halo in a matter-dominated Ω_m = 1 Universe, the only characteristic scale is that set at turn-around, which is the instant of decoupling from cosmological expansion. Therefore, if position, time, and mass profile were to be scaled with respect to the turn-around radius, time, and mass inside the turn-around region, all the particles would trace the same trajectory regardless of their initial conditions; that is, the system would be “self-similar”. Fillmore & Goldreich (1984), hereby FG, developed self-similar solutions for planar, cylindrical and spherical symmetries. They were able to compute position-time and phase-space trajectories as well as mass and density profiles. For the cylindrically symmetric (2D) case, which will be the focus of this work, the predicted asymptotic density profile is ρ∝r⁻¹, independent of the model parameters and for the spherically symmetric (3D) case, it varies between ρ∝r^−2.25−r⁻² depending on the halo mass. Bertschinger (1985) independently developed a self-similar theory of secondary spherical infall around an already collapsed density perturbation, which predicts a profile ρ∝r^−2.25. There lies a clear contention between the density profiles predicted by these self-similar theories and those of prompt cusps and universal NFW profiles measured in CDM simulations. Many extensions to FG exist that include, for example, angular momentum (e.g., White & Zaritsky 1992; Sikivie et al. 1997; Nusser 2001), tidal torque (e.g., Zukin & Bertschinger 2010a, b), and ellipsoid dynamics (Lithwick & Dalal 2011). In this paper, we focus on the FG self-similar model, aware of the inherent limitations due to the assumed symmetry and purely radial motion.

CDM is modeled as a self-gravitating and collisionless fluid obeying Vlasov-Poisson equations:

$$\frac{\partial f}{\partial t}+u·{\nabla }_{r}f-{\nabla }_r\varphi ·{\nabla }_{u}f=0,$$(1)

$${\mathrm{\Delta }}_r\varphi =4\pi G\rho =4\pi G\int f(r,u,t){\mathrm{d}}^{3}u,$$(2)

where f(r, u, t) represents the phase-space density at position r, velocity u and time t and ϕ is the gravitational potential. These are traditionally resolved with an N-body approach, that is, representing the phase space by an ensemble of particles that follow the standard Lagrangian equations of motion in an expanding Universe:

$$\frac{{\mathrm{d}}^{2}x}{\mathrm{d}{t}^2}+2H\frac{\mathrm{d}x}{\mathrm{d}t}=-\frac{1}{a^2}{\nabla }_x\varphi (x);{\mathrm{\Delta }}_x\varphi (x)=4\pi G\overline{\rho }{a}^2\delta (x),$$(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.

In this work, we opt to study dark matter halo dynamics in 2D Vlasov simulations of simplified initial conditions composed of two crossed sine waves run using the ColDICE Vlasov solver of Sousbie & Colombi (2016). We choose these sine-wave initial conditions since they are symmetrical, smooth, and controllable configurations that still retain generality, as they can be assimilated to high peaks of a smooth random Gaussian field (Bardeen et al. 1986). Varying the amplitudes of the sine waves allows us to span a realistic range of initial configurations, from quasi-unidimensional to axis-symmetric. Additionally, the halo centers being predetermined and the absence of mergers and angular momentum make it convenient to test simplistic self-similar models like that of FG. For comparison to our 2D simulations, we use the cylindrically symmetric FG solutions that are known to lead to a parameter-independent asymptotic profile ρ∝r⁻¹, which can be easily verified with our simulations. Though the 2D dynamics are much simpler than in an actual 3D CDM cosmology, it makes the study of fine deviations from self-similar behavior and their causes much more feasible. FG self-similar solutions have already been used to compare against phase-space structures (Sugiura et al. 2020) and density profiles (Enomoto et al. 2023a, b) of late-time equilibrated CDM halos in the multistream regime. Our focus shall be on probing deeper at the level of individual particle trajectories starting from the first shell-crossing itself, in addition to investigating phase space, transverse motion, and profiles of mass, density, and anisotropy parameter as measured in our simulations. As an additional note, dynamics in 2D corresponds to the interaction between infinite lines in 3D, which can be viewed as an approximation of the dynamics of filaments in the context of 3D CDM cosmological dynamics.

The plan of the paper is as follows. In Section 2, we detail the methodology used to conduct numerical simulations and generate data for isolated halos. In Section 3, we revisit the FG self-similar model and show the semi-analytic solutions for trajectories, mass and density profiles. In Section 4, we outline our choice of observables and the procedure of fitting self-similar solutions to our numerical data to obtain bounds and best-fit parameters, additionally discussing the causes behind deviations from self-similarity. The resulting distribution of best-fit parameters and their comparison across simulations with different initial set-ups are discussed in Section 5. A summary of our key qualitative inferences about dark matter halo dynamics in 2D, possible theoretical extensions and its implications on actual 3D CDM cosmological simulations and observations are presented in Section 6.

2. Numerical simulations

The simulations under study were carried out by Saga et al. (2022) using the public Vlasov solver ColDICE (Sousbie & Colombi 2016). It models the CDM phase-space distribution f as a 2D (or 3D) sheet, assuming negligible initial velocity dispersion, in a 4D (or 6D) phase space:

$$f(r,u,t=t_{\mathrm{ini}})={\rho }_{\mathrm{ini}}(r){\delta }_{\mathrm{D}}(u-u_{\mathrm{ini}}),$$(4)

where ρ_ini, u_ini are the density and peculiar velocity fields initialized by Zeldovich flow. The phase-space sheet is adaptively tessellated with simplices (triangles in 2D and tetrahedra in 3D) whose vertices [x(t), v(t)] are then evolved according to the Lagrangian equations of motion (3). The matter is linearly distributed inside each simplex instead of being transported by the vertices, unlike the N-body approach. For more details on refinement and measurements, refer to Sousbie & Colombi (2016) and Saga et al. (2022).

We choose to study highly symmetric cases, where the displacement field Ψ is initialized by sine waves with amplitudes (ϵ_x, ϵ_y) along the x and y axes:

$$x(q,t)=q+\mathrm{\Psi }(q,t),$$(5)

$${\mathrm{\Psi }}_{\mathrm{i}}(q,t_{\mathrm{ini}})=\frac{L}{2\pi }{D}_{+}(t_{\mathrm{ini}}){ϵ}_{\mathrm{i}}sin\left(\frac{2\pi }{L}{q}_{\mathrm{i}}\right),$$(6)

where q is the comoving initial position, D₊ is the linear growth factor, L is the comoving size of the simulation box with periodic boundaries and i is an index for x, y. The three sets of (ϵ_x, ϵ_y) = (−18, −3), (−18, −12), (−18, −18) are used to set the initial conditions for the simulations, which we denote as quasi-1D (Q1D), anisotropic (ANI), and axial-symmetric (SYM), respectively. The parameters used for each simulation are tabulated in Table 1. The first shell-crossing (self-intersection of the phase-space sheet) occurs along the x-axis, followed by shell-crossing along the transverse y direction that leads to the formation of a single monolithically growing halo precisely at the center, which is the ideal scenario to test for self-similar particle trajectories unaffected by mergers. The simulations are in increasing order of resemblance to circular symmetry that is assumed in FG self-similar solutions. Our comparative study across the three simulations allows us to investigate how the halo properties vary with initial conditions.

Figure 1 shows the density projections of the halos seeded by the different initial conditions in each simulation at gradually increasing times. The structure of caustics, defined by the points corresponding to the folds in the phase-space sheet with singular density, is clearly visible. The most notable feature about their evolution is that despite the stark symmetry difference in the shape of the caustics at shell-crossing, they grow roughly similar and close to circular at late times.

Fig. 1. Density color maps of our numerically simulated halos. In the left column, we have Q1D simulation snapshots for the expansion factors a = 0.055, 0.26, 0.7. In the middle, we have ANI simulation snapshots for a = 0.05, 0.125, 0.42. In the right column, we have SYM simulation snapshots for a = 0.045, 0.1, 0.32. The top row consists of the closest available snapshots after the first shell crossing in each of the simulations and the bottom row consists of the last available snapshots.

3. Theory

We briefly recap the FG self-similar solutions. A starting expansion factor a_i in the matter domination era wherein the flow is dominated by Hubble expansion and a scale-free (power-law) initial perturbation δ are assumed:

$$\delta \equiv \frac{\delta M_{\mathrm{i}}}{M_{\mathrm{i}}}={\left(\frac{M_{\mathrm{i}}}{M_0}\right)}^{-ϵ};M_{\mathrm{i}}\equiv M(r,t_{\mathrm{i}})={\displaystyle {\int }_0^r}{r}^{\prime n-1}\rho (r\prime ,t_{\mathrm{i}})\hspace{0.5em}\mathrm{d}r\prime ,$$(7)

where M_i is the initial profile of cumulative mass in planar (n = 1), cylindrical (n = 2) or spherical (n = 3) shells and M₀, ϵ are the model parameters. ϵ is related to the halo mass and accretion rate and M₀ is a reference mass. We select the units $\overline{\rho }(t_{\mathrm{i}})(a_{\mathrm{i}}L{)}^n$ for M_i, M₀ so that their values are independent of the length and mass scales used and physically represent the fraction of the total mass. For the first output snapshots of our simulations with a_i = 0.005, the assumption of the initial motion being dominated by the Hubble flow remains valid, but not that of the initial perturbation δ obeying a power law. However, it is a plausible assumption that the trajectories with small δ do approach a self-similar solution gradually.

Thus, we start in the single-stream regime, where the particles (analogous to points in the phase-space sheet in the Vlasov simulations) continue to expand away from the center. Once the local gravitational pull dominates over the Hubble expansion, the particles turn around. Since the average density decreases with distance from the center, particles initially farther away turn around later, as a result of which the turnaround radius grows with time. Subsequently, particles undergo first shell-crossing symmetrically along all directions, leading to the formation of the first caustics. The flow inside the splashback radius (radius of the outermost caustic) turns multi-stream. The particles in this regime continue to oscillate about the center with an asymptotically converging amplitude. As particles farther away from the center continue to turn around and infall, a power-law density profile builds up in the multistream regime.

As discussed in Section 1, the system possesses only one characteristic scale, which is defined by the turnaround. The variables in the system – position r and time elapsed t for particles (labeled by the initial mass M_i enclosed in concentric shells on which they are located) and mass profile M(r, t) – are normalized w.r.t. the values at turnaround:

$$\tau =\frac{t}{t_{\mathrm{ta}}(M_{\mathrm{i}})}=C_t^{-3/2}\left(\frac{t}{t_{\mathrm{i}}}\right){\left(\frac{M_0}{M_{\mathrm{i}}}\right)}^{3ϵ/2},$$(8)

$$\lambda =\frac{r}{r_{\mathrm{ta}}(M_{\mathrm{i}})}=C_r^{-1}\left(\frac{r}{r_{\mathrm{i}}}\right){\left(\frac{M_0}{M_{\mathrm{i}}}\right)}^{ϵ},$$(9)

$$\mathcal{M}\left(\frac{r}{r_{\mathrm{ta}}(t)}\right)=\frac{M(r,t)}{M_{\mathrm{ta}}(t)};M_{\mathrm{ta}}(t)=\frac{M_0}{C_t^{1/ϵ}}{\left(\frac{t}{t_{\mathrm{i}}}\right)}^{2/3ϵ-2(3-n)/3},$$(10)

where the indices ta and i refer to turnaround and initial. Our case of interest is that of cylindrical symmetry n = 2, for which the coefficients C_r, C_t relating turnaround radius r_ta and time t_ta to the initial density perturbation δ are 0.74 and 1.39, respectively. It is important to recognize that r_ta(M_i), the turnaround radius of the shell enclosing initial mass M_i, is different from r_ta(t), the radius of the shell turning around at the instant t. The coupled differential equations in terms of the scaled variables are as follows:

$$\frac{{\mathrm{d}}^2\lambda }{\mathrm{d}{\tau }^2}=\frac{2(3-n)}{9n}\frac{\lambda }{{\tau }^2}-\frac{2}{3n}{\left(\frac{C_t}{C_r}\right)}^n{\tau }^{2/3ϵ-2(3-n)/3}\frac{\lambda }{|\lambda {|}^n}\mathcal{M}\left(\frac{\lambda }{\mathrm{\Lambda }}\right),$$(11)

$$\mathcal{M}\left(\frac{\lambda (\tau )}{\mathrm{\Lambda }(\tau )}\right)=\frac{2}{3ϵ}{\displaystyle {\int }_1^{\infty }}\frac{\mathrm{d}\tau \prime }{{\tau }^{\prime 1+2/3ϵ}}H\left[\left|\frac{\lambda (\tau )}{\mathrm{\Lambda }(\tau )}\right|-\left|\frac{\lambda (\tau \prime )}{\mathrm{\Lambda }(\tau \prime )}\right|\right],$$(12)

where Λ(τ) = τ^2(1+1/3ϵ)/3 and H is the Heaviside step function. The boundary conditions are λ(τ = 1) = 1 and dλ/dτ(τ = 1) = 0. The equations have been slightly modified by taking the absolute magnitude of λ to change the sign of the particle's position upon crossing through the center λ = 0. This allows for smooth integration of the equations about the center.

The fact that Eqs. (11) and (12) as well as the boundary conditions are independent of initial time t_i, position r_i or enclosed mass M_i implies that all the particle trajectories trace the same curve in the scale-free space, that is, they are “self-similar”. We need to solve the equations only once and then appropriately scale the solution to obtain the trajectory of any particle in our simulation.

Before turnaround (τ<1), we integrate backwards over a grid of τ∈[0, 1] with dτ = 0.01. We substitute ℳ=τ^−2/3ϵ since the initial mass enclosed remains conserved M(r(t), t) = M_i for single stream flow. After turnaround (τ≥1), the equations are solved iteratively up to the desired convergence starting with an initial guess for ℳ(λ/Λ), which is a monotonic function between ℳ(0) = 0 and ℳ(1) = 1. Eq. (11) is integrated forwards over a grid of τ∈[1, τ_f] with dτ = 0.01 and the solution λ(τ) is linearly interpolated. τ_f is chosen such that the amplitude of oscillations decreases below λ/Λ = 10⁻³, allowing us to precisely resolve the mass profile down to r/r_ta = 10⁻³. The mass profile ℳ is computed over a grid of 500 points ∈[10⁻³, 1] and extrapolated down to ℳ(0) = 0 before being used for the next iteration.

The equation for computing the mass profile (12) can be further simplified by breaking the integral into intervals of τ_k's that satisfy λ(τ)/Λ(τ) = λ(τ_k)/Λ(τ_k) (Bertschinger 1985):

$$\mathcal{M}(\lambda /\mathrm{\Lambda })={\mathrm{\Sigma }}_{k=1}^N(-1{)}^{k+1}{\tau }_k^{-2/3ϵ}.$$(13)

Numerically integrated solutions to the coupled Eqs. (11) and (13) corresponding to the parameter ϵ = 0.4, 0.45, 0.5, 0.55, 0.6 are presented in Figure 2. The parameter M₀ enters only while scaling the self-similar solutions back to the appropriate particle. From the plots of self-similar trajectories, we note that lower ϵ corresponds to lower amplitude and greater frequency of oscillations as well as lower velocity at center crossing. A higher frequency of oscillations implies a larger number of caustics, seen as spikes in mass and density profiles. The density profile is given by:

$$\frac{\rho (r,t)}{\overline{\rho }(t)}=\frac{1}{n}{\left(\frac{C_t}{C_r}\right)}^n{\left(\frac{r_{\mathrm{ta}}}{r}\right)}^{n-1}\frac{\mathrm{d}\mathcal{M}(r/r_{\mathrm{ta}}(t))}{\mathrm{d}(r/r_{\mathrm{ta}})}·$$(14)

The asymptotic logarithmic slope for the density profile is −1 for cylindrical symmetry n = 2 and is independent of ϵ (refer Equation 39 of FG), which agrees with the numerical solution.

Fig. 2. FG self-similar solutions assuming cylindrical symmetry (n = 2) for ϵ∈[0.4, 0.6]. The solutions are depicted in scale-free spaces(variables normalized w.r.t. the turnaround). Top left: position – time trajectories. Top right: phase-space trajectories. Bottom left: mass profiles. Bottom right: density profiles.

4. Data comparison

In this section, we test the FG self-similar solutions assuming cylindrical symmetry against data from the 2D Vlasov simulations detailed in Sect. 2. We particularly focus on particle trajectories, mass, and density profiles while also investigating the phase space, transverse motion, and anisotropy parameter. The key questions we want to address pertain to

1. The extent of self-similarity in the motion of particles in our simulations.

2. Verifying if the mass and density profiles are in agreement with self-similar predictions.

3. The deviations from self-similarity and their causes.

4. The distribution of best-fit parameters M₀ and ϵ, comparison between the three simulations, and verifying if the halos become more circular over time.

5. The implications on CDM halos seeded by Gaussian random fields in 3D cosmological simulations.

4.1. Particle trajectories

There are two ways to interpret the self-similar solution for λ(τ). Refer Eq. (8). If t is varied while M_i is fixed, τ acts as time elapsed for a particle with initial position q on a shell enclosing an initial mass M_i that oscillates about the center. In this case, the self-similar solution λ−τ can be scaled using Eqs. (8) and (9) to get the predicted r−t trajectory of the corresponding Lagrangian point q on the phase-space sheet in our simulation. The initial mass profile M_i(q) used to map to the Lagrangian point q is measured from the Vlasov simulation.

Figure 3 shows self-similar fits (blue) to the position-time trajectories (red) of selected Lagrangian points along the x-axis in SYM simulations. The fits are made using the standard least squares method. The data points are then cubically interpolated and compared with the fits to identify the points beyond which the relative differences between the data curves and the fits exceed 10%. The time intervals within these bounding points are shown in green. In each case, the trajectory follows the fit only after 1−2 oscillations after shell-crossing. A possible explanation could be that at t_i, there is no collapsed material and the initial perturbation is not a power-law as assumed in the self-similar model (see Eq. (7)). It takes a brief period of relaxation for a power-law profile to build up before the trajectories approach a self-similar solution. The trajectories deviate again after tracing the fits for 3−4 oscillations, the reason for which could not be verified due to the sparsity of available snapshots at late times.

Fig. 3. Position-time trajectories of three Lagrangian points along the x-axis: q/L = 0.072, 0.087, 0.103 (left to right) in the SYM simulation. The data points (red, labeled “Vl” as in Vlasov solver used for the simulations) have been cubically interpolated. The self-similar fits, made using the least-squares method with ϵ = 0.5, are shown in blue. The time intervals (green) within which the trajectories follow the fits were computed assuming a threshold: Δr/rdata≤10%. The spatial resolution of the tessellation grid is ∼0.0005 comoving box units.

In the second interpretation, if M_i is varied keeping t fixed, τ acts as a label running over particles in the snapshot t. The self-similar solution λ−τ can be scaled using Eqs. (8) and (9) to get the predicted locus of positions r of Lagrangian points q in the snapshot t. The initial mass profile M_i(q) measured from the Vlasov simulation is used to map τ to the Lagrangian points q. This interpretation has an advantage over the former – the Lagrangian points q can be sampled down to the grid resolution whereas the number of snapshots we have are limited and sparse.

Figure 4 shows self-similar fits (blue) to the r−q curves (red) of Lagrangian slices along the x-axis for Q1D, ANI, and SYM simulations in increasing order of expansion factor. The fits were made keeping ϵ fixed, allowing only M₀ to vary and then repeated for several values of ϵ. For the purpose of presentation, ϵ = 0.8, 0.7, and 0.5 for Q1D, ANI, and SYM, respectively, were chosen as representative values. Further details regarding the distribution of best-fit parameters are in Section 5. The subsets of particles whose positions differ by ≤10% relative to the fits are shown in green. Similar to what we observed in the former approach, in each of the three simulations, the particles initially deviate for ∼1−2 oscillations after shell-crossing, then follow the fit for 3−4 oscillations in Q1D, 5−6 oscillations in ANI and 8−9 oscillations in SYM simulations and eventually, deviate again. Put differently, the subset of particles that follow the FG self-similar fit (green), increases in the order Q1D < ANI < SYM and it is clearly due to the lower bound decreasing in the same order. This suggests that the reason for the eventual deviation is correlated to the extent of non-radial dynamics arising from symmetry and not an artifact of numerical simulations. As expected, the FG solutions work best for SYM simulations where the trajectories are highly radial. The ANI and Q1D cases have increasingly elliptical collapses and transverse motion, which lead to greater deviations. From the last available snapshots (bottom row), we may concur that the initial conditions do leave their imprint on the particle dynamics for 10−11 shell crossings at the very least.

Fig. 4. r−q curves of Lagrangian slices along the x-axis for the three simulations – Q1D (left), ANI (middle), SYM (right) in increasing order of expansion factor from top to bottom. The data points (red, labeled “Vl” as in Vlasov solver used for the simulations) have been cubically interpolated. The self-similar fits (blue) are made using the least-squares method with the parameter ϵ fixed at 0.8, 0.7, and 0.5 for Q1D, ANI, and SYM simulations, respectively. The subsets of particles (green) that follow the self-similar fits were computed assuming a threshold: Δr/rdata≤10%. The axes have been scaled to make the features at lower values of q more prominent. The spatial resolution, ∼0.0005 comoving box units, is shown using dashed black curves.

Also, we note the increase in frequency of oscillations from Q1D to SYM case. Consider the Lagrangian point q/L = 0.1 as an example. It has undergone 3 oscillations in Q1D, 5 oscillations in ANI and 6 oscillations in SYM case by a = 0.32. This is the reason why the best-fit self-similar solutions have correspondingly lower ϵ across the 3 simulations, even though we have the same mass in all of them.

Each halo ought to be characterized by a single set of (M₀, ϵ) in the FG model. What we find in our tests is rather a distribution of best-fit parameters at different times and spatial regions of each halo, which we examine in Section 5.

4.2. Deviations

In this subsection, we study the deviations of particle trajectories from the FG self-similar fits in greater detail. In particular, we look at where the deviations occur, and how they change with time as well as across the three simulations to try and identify the causes behind them.

The subset of particles bounded in green denotes the portion of the phase-space sheet whose Eulerian position in the simulation is matched by the self-similar fits within a 10% error. It implies that an entire halo does not “turn self-similar” as a whole. By repeating this fitting exercise for Lagrangian slices along different directions θ for a given snapshot of a simulation, we find the subset of particles to be confined between two concentric rings in Lagrangian space q_x−q_y whose radii are equal to the computed bounds. Figure 5 depicts the bounds in Lagrangian space of the SYM simulation at one of the snapshots as an example.

Fig. 5. Bounds in Lagrangian space of the SYM simulation at the snapshot a = 0.2. The space has been color-mapped to corresponding Lagrangian density $\sim {\left|\partial x/\partial q\right|}^{-1}$. The contours of high density represent the caustics.

Fig. 6 shows the evolution of the bounds computed for Lagrangian slices along the x-axis in each of the three simulations. We note that the first snapshot wherein a self-similar fit could converge is always after shell-crossing along the y-axis and at increasingly late times a_r = 0.075, 0.125, 0.26 going from SYM to ANI to Q1D simulations (refer Table 1). As discussed earlier, a power-law profile needs to build up before the trajectories can approach a self-similar solution. The lower amplitude of initial displacement ϵ_y delays the shell-crossing along the y-axis and as a consequence, relaxation and profile build-up are also delayed.

Fig. 6. Evolution of bounds in Lagrangian space computed using r−q curves of Lagrangian slices along the x-axis for three cases: Q1D (left), ANI (middle), SYM (right). The parameter ϵ was fixed at 0.8, 0.7, and 0.5, respectively. The multicolored patches labeled with percentages depict the subset of particles in Lagrangian space for which the relative difference in their positions between the data curves and the self-similar fits are within thresholds of 5%, 7.5%, 10%, 12.5%, and 15%. The spaces excluded by bounds arising from relaxation, periodic boundaries, and anisotropy have been shown in blue, yellow, and red patches, respectively.

Looking at the subset of particles bound in green, the upper bound grows initially (clearly observed in ANI and SYM cases), but then stagnates, even reducing afterward. To track the initial growth, we hypothesize that after shell-crossing, the particles take 1−2 oscillations to relax to the self-similar fits. Relaxation (Lynden-Bell 1967) refers to the process of phase-space mixing following which particle motion can be described using acceleration computed from a smooth potential field. The subset of such particles grows as more particles starting farther away from the center (at higher q) gradually relax and turn self-similar. This argument is tricky for the Q1D and ANI cases, as they have different shell-crossing times along the x and y axes, whereas the FG model just has one owing to circular symmetry. In the simulations, we observed that the particles begin to relax to self-similar fits only after the last shell-crossing has occurred (along the y-axis in our cases). Since we cannot track the two shell-crossings separately using the FG model, we simply use the shell-crossing from theory (which might not correspond to the actual shell-crossing along the y-axis) to build our argument. Referring to the top left panel of Fig. 2, 1−2 oscillations after shell-crossing roughly corresponds to τ = 7−20 for ϵ = 0.8, τ = 6−18 for ϵ = 0.7, τ = 5−13 for ϵ = 0.5 in Q1D, ANI, and SYM cases, respectively. Using Eq. (8), the M_i's corresponding to these values of τ were computed, which were then mapped to the bounding q's in Lagrangian space for each snapshot. The subsets of particles that lie above these bounds are shown in blue and they seem to coincide well, at least for the SYM case, with the particles that initially deviate.

As all particles would eventually relax, the upper bound should have kept increasing until q/L = 0.5, however, it stagnates for Q1D and ANI cases and decreases for the SYM case after it reaches q∼0.3. Taking a look at the r−q curves in Figure 4 at the last available snapshots (bottom row), we notice that the particles with q/L roughly ≥0.3 have not completed as many oscillations as they ought to if they had traced the self-similar fits. This is clearly due to a difference in the force between FG theory, which assumes a single isolated halo, and our simulations, which have periodic boundaries. As a hand-waving argument, consider the motion of a particle with the same initial position q in the two setups, one with a single isolated halo at the center and the other having an additional halo at x/L = 1. Initially, their relative difference in positions is negligible but grows gradually. Thus, the initial deviation from self-similarity is dominantly due to the time taken for relaxation, after which the deviation caused by periodic boundaries takes over once the relative difference in positions crosses the % threshold we set to compute the bounds in q. The relative difference between the force fields in the two setups scales as x/(L−x); x/L∈[0.0, 0.5], which implies that if the initial position q was closer to the boundary, the particle would suffer from a greater erroneous force arising from the halo image. Hence, it would accumulate error in its trajectory faster and deviate from self-similar fit earlier. This explains why the upper bound stagnates and moves to lower q later in time, most prominently observed for the SYM case. The subsets of particles that experience a relative error in force ≥40% are shown in yellow in Fig. 6.

The lower bound remains roughly constant over time in each simulation, but its magnitude decreases going from Q1D to ANI to SYM. Since the upper bound stagnates at q/L∼0.3 for each simulation, it implies that the subset of particles in agreement with FG solutions increases in the same order. We hypothesize this to be correlated to the extent of transverse motion across the simulations. Figure 9 shows that the degree of transverse motion is indeed the highest for particles in Q1D simulation. Since FG solutions assume fully radial orbits, the fits for r−q curves of Lagrangian slices along non-axial directions barely converge in Q1D and ANI simulations. In Figure 10, the transition of the anisotropy parameter β from 1 to 0 as we move closer to the halo center roughly demarcates the halo into two regions with different dynamics – the outer region, dominated by radial infall and the inner region, where the velocity distribution approaches isotropy after violent relaxation due to transverse motions. For Lagrangian slices along axial directions, we could expect to see deviations from FG solutions once the amplitude of oscillations falls within the radius of the inner region r_trans where the transverse motion is no longer negligible. From Fig. 4, we note that the asymptotic amplitude of oscillations can be parameterized by a power law: r/aL=A(a)q^α, where α = 1.85, 1.8, 1.7 for Q1D, ANI, and SYM cases, respectively, and A(a)∼0.3−0.5. The transition radius r_trans(a) was computed assuming β(r_trans) = 0.5 using sigmoid fits to the anisotropy profile at each snapshot. The subset of particles whose amplitudes of oscillations are less than r_trans at a given snapshot, that is, A(a)q^α≤r_trans(a), are shown in red in Fig. 6.

Summing up the deviations seen in particle trajectories:

• –

  Particles typically take 1−2 oscillations after the last shell-crossing (along the y-axis) to relax to self-similar motion, resulting in the initial upper bounds. Therefore, this deviation is due to physical reasons.

• –

  On the other hand, the upper bounds in the later stages are due to the periodic boundary condition that exerts artificial forces on the particles. The difference between FG trajectories and simulated trajectories of particles gradually increases, faster for the particles that start closer to the boundaries. This deviation is an artifact of our simulations.

• –

  Finally, the lower bound is associated with the fact that we are comparing FG self-similar fits assuming fully radial motion to simulated trajectories in halos that have significant transverse motion in their central regions. Once the amplitude of oscillations of the particles decreases down to the radius inside which transverse motion turns significant, they start to deviate. It is to be expected that the lower bound decreases across the simulations in increasing order of circular symmetry, which means that the SYM simulations are the most consistent with the FG model of self-similarity. It would therefore be interesting to consider more general models of self-similarity that take transverse motion as well as elliptical collapse into account.

4.3. Mass and density profiles

After analyzing the motion of particles in our simulations, we turn to mass and density profiles. The self-similar prediction for mass profile inside the turnaround region ℳ(r/r_ta) is obtained by iteratively solving Eqs. (11) and (13), from which we can derive the density profile using Eq. (14). These are, again, characterized by M₀ and ϵ. From the simulations, we measure the cumulative mass and circularly averaged density in radial bins. Using least squares, we obtain the best-fit M₀ for each snapshot, keeping ϵ fixed at 0.8, 0.7, and 0.5 for Q1D, ANI, and SYM, respectively. Figure 7 shows the self-similar solution overplotted with the mass and density profiles from the three simulations normalized w.r.t. the mass M_ta and radius r_ta of the turnaround region using the best-fit M₀(ϵ) at multiple snapshots, color coded from blue to red in increasing order of expansion factor. We note that the first snapshot for which self-similar fit could converge to the mass (or density) profile is at increasingly late times a_M = 0.055, 0.1, 0.2 going from SYM to Q1D, while satisfying a_{SC, y}≤a_M≤a_r in each case (refer Table 1). Non-convergence of self-similar fit means that a power-law profile has not been fully built up yet. This verifies that lower ϵ_y delays the shell-crossing along y, which delays relaxation and profile build-up, which in turn delays the approach of trajectories to self-similar solutions.

Fig. 7. Mass (top) and density (bottom) profiles from the three simulations: Q1D (left), ANI (middle) and SYM (right) normalized w.r.t. to the mass Mta and radius rta of the turnaround region using the best-fit M0(ϵ) at different snapshots, color coded from blue to red in increasing order of expansion factor. The dashed black curve denotes the self-similar predictions for the mass and density profiles corresponding to ϵ = 0.8, 0.7, 0.5 for Q1D, ANI, and SYM, respectively.

The magnitude and slopes of the mass and density profiles are in good agreement. The number and location of caustics (spikes in mass and density profiles) differ initially but eventually conform to the self-similar fits. This can be understood from the r−q curves of Lagrangian slices in Fig. 4. Caustics correspond to the extrema of oscillations. As the self-similar fits to the number and amplitude of oscillations improve with time, so does the number and location of caustics. At late times, however, we note a dip in the slope of mass and density profiles at outer region of the halos. This stems from the fact that the FG model assumes indefinite infall of mass onto the halo whereas in our simulations, once the mass M_ta (or radius r_ta) of turnaround region grows beyond the mass (or box size) of our simulation, there is a deficit of infalling matter. The periodic boundaries holding back the particles closer to the boundaries also aggravate this issue.

The time-averaged slopes of the mass profiles in the range 0.001≤r/r_ta≤0.01 for Q1D, ANI, and SYM cases are 1.05, 1.10, and 1.00, respectively. These are in good agreement with an asymptotic slope of 1.0 predicted in FG for cylindrical symmetry. Since all the simulations produce nearly the same slope despite differing in their values of best-fit ϵ, it verifies that the slope is indeed independent of ϵ (in turn, the halo mass) for the 2D case.

4.4. Scale-free space

Self-similarity implies that trajectories in position and phase space basically have the same shape. Upon scaling w.r.t. the characteristic length and time – the turnaround radius and time in the FG model, the trajectories should overlap. As a consistency check, after obtaining the best-fit parameter M₀(ϵ) using fits to particle trajectories in Sect. 4.1, we scale the trajectories w.r.t. to turnaround using Eqs. (8) and (9), and overplot all the snapshots together with the self-similar solutions for ϵ = 0.8, 0.7, 0.5 for Q1D, ANI, SYM cases, respectively, in Figure 8.

Fig. 8. Superposition of position-time (top) and phase-space (bottom) trajectories of particles of Lagrangian slices along the x-axis from all the snapshots (where a fit could converge) after being scaled w.r.t. turnaround computed using Equations (8) and (9) using the best fit M0(ϵ) obtained from fitting r−q curves (refer Fig. 4) in Sect. 4.1. Zoomed panels of the central regions in phase space have also been added. Self-similar solutions for ϵ = 0.8, 0.7, 0.5 corresponding to Q1D (left), ANI (middle) and SYM (right) are shown in red.

The superposition of scale-free position-time trajectories from simulations, shown in the top row of panels, is strikingly good. The curves corresponding to later snapshots follow the self-similar curves longer in each simulation. As expected, the SYM simulations are the most consistent, with most of the curves being in agreement for roughly 7−8 oscillations. We also note that the oscillation frequency in each curve decreases with τ, that is, the more the number of oscillations a particle completes, the less its oscillation frequency. The r−q curves for snapshots at later times would thus be better fit by self-similar solutions corresponding to higher values of ϵ (refer to the top-left panel of Fig. 2). This is indeed what we found upon redoing the fits for different values of ϵ, discussed in detail in Section 5. The parameter ϵ is inversely related to the mass accretion rate, which means that the mass accretion rate in our simulations is less than that expected from the FG model at later snapshots. This is clearly due to the deficit of infalling matter in our simulations that also leads to the dips in mass and density profiles at larger radii at later snapshots.

In the phase space curves shown in the bottom row of panels, the self-similar spiral pattern is remarkably consistent across the snapshots for each simulation. It is crucial to note that even though the phase-space spirals in Q1D and ANI show deviations from FG self-similar solutions, the simulated curves superpose quite well within themselves. This insinuates the existence of a homothetic transform (or self-similar solution) more general than the FG solution alone, one that incorporates transverse motion. In the zoomed panels of the central regions, we can clearly see the resonant modes compromising the mean-field limit, especially in Q1D. This is a numerical defect. The radius below which the motion of particles is contaminated by these resonant modes is approximately r/aL∼1−2×10⁻³, which is greater than the grid resolution ∼5×10⁻⁴. However, this numerical bound is still less than the physical lower bound arising from transverse motion below the transition radius, the least of which is for the SYM case: r_trans/aL≥5×10⁻³ (refer Fig. 10 and Section 4.2).

4.5. Transverse motion and anisotropy parameter

To corroborate our hypothesis that non-radial dynamics cause particle trajectories in the interior of halos to deviate from FG solutions, we probe the extent of transverse motion and anisotropy in our simulations.

Fig. 9 shows the x−y trajectories of particles with initial positions along different directions in the three simulations. For the SYM case, only the particles starting along the x = 0, x=y, y = 0 directions have completely radial trajectories. The other trajectories, while being mostly radial, do show deviations. For ANI and Q1D cases, only the particles starting along axial directions are radial. Q1D being the most asymmetrical has non-axial trajectories exhibiting the greatest extent of transverse motion.

Fig. 9. x−y trajectories of five Lagrangian points (qx/L, qy/L) = (0.094, 0.0), (0.094, 0.008), (0.094, 0.094), (0.023, 0.094), (0.0, 0.094) for the 3 cases – Q1D (left), ANI (middle), SYM (right). The axes have been logarithmically scaled to feature the motion close to the center prominently.

To determine the region of the halos where transverse motion starts to be significant, we look into the anisotropy parameter, which is defined as $\beta (r)=1-{\sigma }_{\perp }^2(r)/{\sigma }_r^2(r)$, where σ_⊥, σ_r are the transverse and radial velocity dispersions, respectively. For radial orbits, ${\sigma }_{\perp }^2=0\Rightarrow \beta =1$. For virialized orbits, ${\sigma }_{\perp }^2={\sigma }_r^2\Rightarrow \beta =0$. The higher the value of β, the more radial the particle trajectories and hence, we would expect better fits to FG solutions. Fig. 10 shows the radial profile of the anisotropy parameter for the three simulations at several snapshots. The decrease in β from 1 to 0 as we move inwards suggests that in the outer region, the dynamics is dominated by radially infalling matter, whereas in the inner region, due to isotropization in velocity space, the particle trajectories have a non-negligible transverse component. To estimate the radius r_trans inside which transverse motion starts to be significant, sigmoid fits to the radial profile of anisotropy parameter were made for each snapshot and r_trans was determined using β(r_trans) = 0.5. The fits were made in the range 0.0005≥r/aL≥0.2. The gray regions denote the ranges of values of r_trans over all the snapshots in each simulation. The key observation is that the transition from β = 1 to β = 0 happens at correspondingly smaller radii the closer the simulation is to circular symmetry – r_trans/aL∈[0.08, 0.1], [0.01, 0.04], [0.005, 0.02] for Q1D, ANI, and SYM, respectively. We also note the unusually high value of β at r<r_trans for Q1D. Since Q1D deviates the most from circular symmetry, one would expect greater isotropization in velocity space close to the center and hence, lower β, which is not what we observe. One hypothesis could be that the Q1D particles indeed show a greater extent of transverse motion until they undergo the first crossing along the x-axis, after which they move more or less along the y-axis and their collapse into the center of the halo (shell-crossing along the y-axis) is more radial than that of the ANI particles, leading to a higher value of β in the central region. It does not seem to be a numerical artifact like the resonant modes since its extent r_trans/aL∼10⁻¹ far exceeds the extent of the resonant modes r/aL∼1−2×10⁻³.

Fig. 10. Radial profiles of anisotropy parameter $\beta =1-{\sigma }_{\perp }^2/{\sigma }_r^2$ at several snapshots for the 3 cases – Q1D (left), ANI (middle), SYM (right). The grid resolution = 0.0005 comoving box units is marked in dashed gray lines. The values of rtrans corresponding to β = 0.5 over different snapshots are marked as the gray patch. Illustrative sigmoid fits to the profiles are shown in dashed black lines.

In actual CDM cosmologies seeded by Gaussian random fields, the collapse of initially overdense perturbations leading to the formation of halos is better modeled by ellipsoidal instead of spherical collapse and the particles exhibit significant non-radial motion as well. Such halos would be closer in resemblance to the Q1D than the SYM case. Therefore, a self-similar analysis of 3D CDM halos would entail the inclusion of ellipsoidal collapse and transverse motion to build a more generic model of self-similarity.

5. Parameter distribution

The self-similar fits to the data from simulations were characterized by two parameters: M₀ and ϵ, defined in Eq. (7). Recall that ϵ characterizes the mass accretion rate and frequency of oscillations (refer Fig. 2), whereas M₀ sets the scale of the initial perturbation that seeds the halo. As the convergence was slow for a two-parameter fit, we kept ϵ fixed and obtained the best fits for M₀ using the least-squares method, then repeated the process for several values of ϵ. Through the procedure outlined in Sections 4.1 and 4.3, we have two sets of best-fit parameters corresponding to – (i) r−q curves of Lagrangian slices (ii) radial profiles of mass and density – for the three simulations at the available snapshots. We now analyze the trend and spread in their distributions and their variation across the three simulations.

Fig. 11 shows the trends in the best-fit M₀ for ϵ fixed at 0.8, 0.7, and 0.5 for Q1D, ANI, and SYM simulations, respectively. The best-fit M₀ for mass-density profiles is typically higher (for ANI and SYM cases at least) than that for r−q curves at earlier snapshots. Looking at the fits for r−q curves in Figure 4, we note that the amplitudes of the first 1−2 oscillations in the simulations are higher than that of the self-similar fits, which means that the radii of splashback (outermost caustic) and subsequent few caustics in the simulations are greater than that of their fits. This implies that the initial perturbation δ in the simulations falls off slower, in other words, is broader at larger radii than what the self-similar fits to r−q curves would suggest. In the fits to mass and density profiles (refer Fig. 7) however, the positions of the outer caustics are predominantly matched at earlier snapshots resulting in the best-fit M₀ being higher. At later times, when more particles have relaxed to lower amplitudes and more caustics have formed, the weight on matching the first caustics is less and the best-fit M₀ decreases.

Fig. 11. Trends in best fit M0 for r−q curves of Lagrangian slices along different angles θ from the x-axis and mass and density profiles for the three cases: Q1D (left), ANI (middle), SYM (right). The parameter ϵ was fixed at 0.8, 0.7, and 0.5, respectively.

We expected to see the best-fit M₀ for r−q curves of Lagrangian slices along the orthogonal x, y axes to converge for ANI and Q1D cases, suggesting that the dynamics inside the halo turns circularly symmetric with time. But, this is not what we observe. The r−q curves along the y-axis are better fit by lower M₀ across all the snapshots. It thus seems that the initial conditions do leave their imprints on halo dynamics for 10−11 oscillations about the center at the very least. But since the caustic structure seemed to grow roughly similar at later times (refer Fig. 1), comparing directly the axial density profiles along orthogonal directions might provide a better demonstration of halos gradually turning circular.

Fig. 12 shows the histograms of best-fit M₀ for r−q curves of Lagrangian slices, mass and density profiles with ϵ fixed for each simulation. Their average and standard deviations (over snapshots at different times) are tabulated in Table 2. The best-fit M₀'s are of the order M₀∈[10^−8.0, 10^−6.9] in units of total simulation mass. The relative spread in the initial perturbation δ (refer Eq. (7)) assumed in the FG model corresponding to the spread in best-fit M₀ from r−q curves is 12%, 20%, and 13% for Q1D, ANI, and SYM, respectively, which is not unreasonably high given the differences between our simulation setups and the ideal theoretical setup.

Fig. 12. Distribution of best fit M0 for r−q curves of Lagrangian slices along different angles, mass and density profiles for the three cases: Q1D (left), ANI (middle), SYM (right). The parameter ϵ was fixed at 0.8, 0.7, and 0.5, respectively.

To get an idea of the trends in best-fit ϵ in each simulation, we compare between different values of ϵ, the sum of squared residuals $\left(=\frac{1}{N}\hspace{0.5em}\mathrm{\Sigma }\mathrm{\Delta }{r}^2\right)$ corresponding to the fits of r−q curves of Lagrangian slices along the x-axis at each snapshot in Fig. 13. N is the number of data points in each r−q curve and Δr is the difference between the position of a Lagrangian point measured from simulation and that expected from the self-similar fit to the curve. The residuals indeed decrease over time implying better fits, which means that the system behaves increasingly in accordance with FG's model of self-similarity. The ϵ corresponding to the least residue does vary over the span of the simulations: 0.4−0.6 for Q1D, 0.8−1.0 for ANI and 0.45−0.55 for SYM. To demonstrate our fitting procedure and results, we choose the self-similar fits corresponding to ϵ = 0.8, 0.7, and 0.5 for Q1D, ANI, and SYM simulations as representatives. The variation in residuals between different ϵ's for the Q1D and ANI cases is not as drastic as that for the SYM case. The residuals corresponding to ϵ = 0.5 remain roughly intermediate throughout the SYM simulation. A better choice of ϵ's, the ones with lower residues, does improve the bounds in Fig. 6 slightly, but the deviations due to relaxation, periodic boundaries and transverse motion are still observable and lead to similar conclusions. The crucial thing to note, at least in Q1D and SYM cases, is that the value of ϵ that gives the least residue increases with time. This is due to the deficit of infalling matter resulting in a lowering of oscillation frequency with time, as observed in Section 4.4.

Fig. 13. Trends in residuals $\left(=\frac{1}{N}\mathrm{\Sigma }(\mathrm{\Delta }r{)}^2\right)$ computed using r−q curves of Lagrangian slices along the x-axis with varying ϵ for the three cases: Q1D (left), ANI (middle), SYM (right).

6. Conclusions

In this work, we have aimed to demonstrate the applicability of self-similarity, particularly the FG model, to dark matter halo dynamics using 2D Vlasov simulations of monolithically growing halos seeded by sine wave initial conditions of three different symmetries: quasi-1D (Q1D), anisotropic (ANI), and symmetric (SYM). By testing the FG model against particle trajectories and mass and density profiles, we have been able to determine not only the subset of particles that followed the model's predictions at a given instant but also the causes behind the deviations for others. Lastly, the trends and spread of the resulting distribution of best-fit parameters have been studied and compared across the three simulations. Below, we summarize our important findings together with the key inferences regarding 3D CDM halos seeded by Gaussian random fields:

1. The initial perturbation δ (refer Eq. (7)) in the simulations certainly does not follow a power law as assumed in FG model. Shell crossing is followed by relaxation that leads to the build-up of a power-law density profile, after which the self-similar fits begin to converge. Particles typically relax after 1−2 oscillations after the last shell-crossing has occurred (which is along the y-axis in our cases) and then move in agreement with self-similar fits (refer Fig. 4). The subset of such particles grows (refer Fig. 6) as the particles starting farther out gradually relax and the model predictions improve. This feature is common to all three simulations. Therefore, the dynamics of relaxation and formation of prompt cusp in 3D CDM halo simulations cannot be explained by the FG model of self-similarity.

2. The particles closer to the boundaries, however, experience additional forces due to the periodic boundaries, which slows down their in-fall, causing them to eventually deviate from the self-similar fits. This is purely an artifact of our simulations and does not concern the CDM halos in 3D cosmological simulations.

3. Even after relaxing to self-similar fits, the particle trajectories trace them only for a limited number of oscillations and deviate again. The number of such oscillations is the highest for the SYM simulation. We believe this to be correlated to the varying degrees of transverse motion in our three simulations. Indeed, the SYM case has particle trajectories with the least transverse motion in the inner region of the halo (refer Fig. 9). Since the FG solutions assume purely radial particle trajectories, it is, but obvious that the SYM simulations appear to be the most “self-similar” under the FG model. We posit that the particle trajectories deviate once their amplitudes decrease below the radius where the anisotropy parameter transitions from β∼1 to β∼0, where transverse motion starts being significant (refer Fig. 10). The deviation of particle trajectories embedded deep inside the halo interior is therefore a consequence of our choice of a self-similar model. Since CDM halos formed from Gaussian initial conditions are closer to ANI and Q1D cases, we need to consider a more general model that accounts for ellipsoidal collapse and transverse motion of particles while studying CDM halo dynamics in actual cosmologies. Looking at the superposition of particle trajectories at different snapshots in scale-free position-time and phase spaces in Fig. 8, we reach to a similar conclusion. Though there exist deviations from the FG model in Q1D and ANI cases, the simulated curves are strikingly self-consistent, which indicates a self-similar pattern more general than the purely radial FG solution alone.

4. The normalization and positions of caustics in mass and density profiles are well captured by the self-similar fits (refer Fig. 7). The measured slopes of the central density profiles – −0.95, −0.90, and −1.00 in Q1D, ANI, and SYM simulations, respectively – are also in agreement with FG's prediction of −1 for the asymptotic slope of the density profile, which is independent of the model parameters. However, at late times, the slope at outer radii dips signifying a deficit of infalling mass in the simulations. Though an artifact of our simulations, we speculate that it might actually be relevant for CDM cosmologies with Gaussian initial conditions. Since most of the mass eventually accretes onto some or the other halo, a similar situation of a lack of infalling matter will arise that will cause the self-similar power-law profile to dip towards the outer region of the halo and lead to a running power-law profile, which might explain the universal attractor – NFW profile.

5. Instead of a single set of best-fit parameters (M₀, ϵ), we obtained a distribution across the snapshots for each simulation (refer Figures 12 and 13), which was to be expected since our setups were not as ideal as that assumed in theory. Nevertheless, the spreads in the distributions using particle trajectories were narrow (refer Table 2) in the sense that the corresponding error implied on the initial perturbation δ (refer Eq. (7)), which defines the model parameters, is within reasonable limits. The fact that we could trace 30−60% of the particle trajectories measured from our simulations within 10% error for periods spanning over several oscillations after shell-crossing as well as the mass and density profiles using FG self-similar solutions corresponding to a narrow range of parameters justifies its relevance in the study of halo dynamics, at least in the 2D case examined in the present work.

6. Since the best-fit parameters for particle trajectories of ANI and Q1D along orthogonal directions do not converge within the simulation time (refer Fig. 11), we could not conclude if the halos gradually turn circular. This suggests that particle trajectories carry the imprints of initial conditions for 11−12 oscillations at the very least. The CDM halos turning increasingly circular is pretty evident from the caustics’ structure in the simulations, see Fig. 1. Therefore, looking for the convergence of axial density profiles along orthogonal directions instead might provide better insights in this regard.

7. Within a particular simulation, the frequency of oscillations of the particles gradually decreases, which we verify by showing that the parameter ϵ corresponding to the self-similar fit to the particle trajectories with the least residue gradually increases (refer Fig. 13). Since ϵ relates inversely to the mass accretion rate, its gradual increase corroborates our claim of a deficit of infalling matter at late times in our simulations.

In conclusion, by performing this exercise of studying 2D Vlasov simulations using the FG self-similar approach, one of the basic models of self-similarity, we have demonstrated self-similar patterns in individual particle trajectories. Though there are deviations, they can be accounted for. The averaged observables, like the mass and density profiles, are still very well explained by the model. It has helped us gain deeper insights into the signatures of self-similarity to look for in 3D CDM simulations with Gaussian fields. Future work would involve performing similar CDM numerical simulations with sine waves as well as Gaussian initial conditions but in 3D. It might also be interesting to study non-spherical power-law initial conditions in 3D because the 3D case is more complex – the asymptotic slope of the radial density profile varies with the FG model parameters as opposed to the 2D case – a prediction that could be tested using numerical simulations. On the analytical front, we would investigate more general self-similar models that include ellipsoidal collapse (e.g., Lithwick & Dalal 2011) and transverse motion of particles and then analyze the particle trajectories, phase space, and mass-density profiles using the procedure established in this work. Therefore, self-similarity could potentially be one of the keys to decoding the complex dynamics in dark matter halos.

Acknowledgments

This work was supported in part by the French Doctoral school ED127, Astronomy and Astrophysics Ile de France (AP), as well as Programme National Cosmology et Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES (SC & AP). This work was also financially supported by the “PHC Sakura” program (project number: 51203TL, grant number: JPJSBP120243208), implemented by the French Ministry for Europe and Foreign Affairs, the French Ministry of Higher Education and Research and the Japan Society for Promotion of Science (JSPS). This work is additionally supported by the JSPS KAKENHI Grant No. JP23K19050 and No. JP24K17043 (SS); and JP20H05861 and JP21H01081 (AT). Numerical computation with ColDICE was carried out using the HPC resources of CINES (Occigen supercomputer) under the GENCI allocations c2016047568, 2017-A0040407568 and 2018-A0040407568. Post-treatment of ColDICE data was performed on the INFINITY cluster of Institut d’Astrophysique de Paris.

References

All Tables

All Figures

Fig. 1. Density color maps of our numerically simulated halos. In the left column, we have Q1D simulation snapshots for the expansion factors a = 0.055, 0.26, 0.7. In the middle, we have ANI simulation snapshots for a = 0.05, 0.125, 0.42. In the right column, we have SYM simulation snapshots for a = 0.045, 0.1, 0.32. The top row consists of the closest available snapshots after the first shell crossing in each of the simulations and the bottom row consists of the last available snapshots.

In the text

	Fig. 2. FG self-similar solutions assuming cylindrical symmetry (n = 2) for ϵ∈[0.4, 0.6]. The solutions are depicted in scale-free spaces(variables normalized w.r.t. the turnaround). Top left: position – time trajectories. Top right: phase-space trajectories. Bottom left: mass profiles. Bottom right: density profiles.
In the text

Fig. 3. Position-time trajectories of three Lagrangian points along the x-axis: q/L = 0.072, 0.087, 0.103 (left to right) in the SYM simulation. The data points (red, labeled “Vl” as in Vlasov solver used for the simulations) have been cubically interpolated. The self-similar fits, made using the least-squares method with ϵ = 0.5, are shown in blue. The time intervals (green) within which the trajectories follow the fits were computed assuming a threshold: Δr/rdata≤10%. The spatial resolution of the tessellation grid is ∼0.0005 comoving box units.

In the text

Fig. 4. r−q curves of Lagrangian slices along the x-axis for the three simulations – Q1D (left), ANI (middle), SYM (right) in increasing order of expansion factor from top to bottom. The data points (red, labeled “Vl” as in Vlasov solver used for the simulations) have been cubically interpolated. The self-similar fits (blue) are made using the least-squares method with the parameter ϵ fixed at 0.8, 0.7, and 0.5 for Q1D, ANI, and SYM simulations, respectively. The subsets of particles (green) that follow the self-similar fits were computed assuming a threshold: Δr/rdata≤10%. The axes have been scaled to make the features at lower values of q more prominent. The spatial resolution, ∼0.0005 comoving box units, is shown using dashed black curves.

In the text

	Fig. 5. Bounds in Lagrangian space of the SYM simulation at the snapshot a = 0.2. The space has been color-mapped to corresponding Lagrangian density $\sim {\left|\partial x/\partial q\right|}^{-1}$. The contours of high density represent the caustics.
In the text

Fig. 6. Evolution of bounds in Lagrangian space computed using r−q curves of Lagrangian slices along the x-axis for three cases: Q1D (left), ANI (middle), SYM (right). The parameter ϵ was fixed at 0.8, 0.7, and 0.5, respectively. The multicolored patches labeled with percentages depict the subset of particles in Lagrangian space for which the relative difference in their positions between the data curves and the self-similar fits are within thresholds of 5%, 7.5%, 10%, 12.5%, and 15%. The spaces excluded by bounds arising from relaxation, periodic boundaries, and anisotropy have been shown in blue, yellow, and red patches, respectively.

In the text

Fig. 7. Mass (top) and density (bottom) profiles from the three simulations: Q1D (left), ANI (middle) and SYM (right) normalized w.r.t. to the mass Mta and radius rta of the turnaround region using the best-fit M0(ϵ) at different snapshots, color coded from blue to red in increasing order of expansion factor. The dashed black curve denotes the self-similar predictions for the mass and density profiles corresponding to ϵ = 0.8, 0.7, 0.5 for Q1D, ANI, and SYM, respectively.

In the text

Fig. 8. Superposition of position-time (top) and phase-space (bottom) trajectories of particles of Lagrangian slices along the x-axis from all the snapshots (where a fit could converge) after being scaled w.r.t. turnaround computed using Equations (8) and (9) using the best fit M0(ϵ) obtained from fitting r−q curves (refer Fig. 4) in Sect. 4.1. Zoomed panels of the central regions in phase space have also been added. Self-similar solutions for ϵ = 0.8, 0.7, 0.5 corresponding to Q1D (left), ANI (middle) and SYM (right) are shown in red.

In the text

	Fig. 9. x−y trajectories of five Lagrangian points (qx/L, qy/L) = (0.094, 0.0), (0.094, 0.008), (0.094, 0.094), (0.023, 0.094), (0.0, 0.094) for the 3 cases – Q1D (left), ANI (middle), SYM (right). The axes have been logarithmically scaled to feature the motion close to the center prominently.
In the text

Fig. 10. Radial profiles of anisotropy parameter $\beta =1-{\sigma }_{\perp }^2/{\sigma }_r^2$ at several snapshots for the 3 cases – Q1D (left), ANI (middle), SYM (right). The grid resolution = 0.0005 comoving box units is marked in dashed gray lines. The values of rtrans corresponding to β = 0.5 over different snapshots are marked as the gray patch. Illustrative sigmoid fits to the profiles are shown in dashed black lines.

In the text

	Fig. 11. Trends in best fit M0 for r−q curves of Lagrangian slices along different angles θ from the x-axis and mass and density profiles for the three cases: Q1D (left), ANI (middle), SYM (right). The parameter ϵ was fixed at 0.8, 0.7, and 0.5, respectively.
In the text

	Fig. 12. Distribution of best fit M0 for r−q curves of Lagrangian slices along different angles, mass and density profiles for the three cases: Q1D (left), ANI (middle), SYM (right). The parameter ϵ was fixed at 0.8, 0.7, and 0.5, respectively.
In the text

	Fig. 13. Trends in residuals $\left(=\frac{1}{N}\mathrm{\Sigma }(\mathrm{\Delta }r{)}^2\right)$ computed using r−q curves of Lagrangian slices along the x-axis with varying ϵ for the three cases: Q1D (left), ANI (middle), SYM (right).
In the text