© The Authors 2026

1. Introduction

Cosmological N-body simulations are fundamental tools in modern cosmology, providing a numerical framework for studying the non-linear evolution of structure from the linear Gaussian early Universe to the highly non-linear cosmic web observed today. They allow us to test cosmological models, interpret large-scale surveys, and explore scenarios beyond the reach of analytic theory. The first attempts to model gravitational clustering in a cosmological context date back to the works of Peebles (1970) and Press & Schechter (1974), who studied the formation of galaxy clusters and introduced a statistical framework for halo formation. These early works laid the foundation for the development of numerical simulations that could capture the complexity of non-linear structure formation.

In the 1970s, the first cosmological N-body simulations were performed with relatively few particles and often employed free (open) boundary conditions in spherical regions, as in White (1976) and Aarseth et al. (1979). At that time, the focus was on modelling galaxy clusters, where edge effects were less critical. By the early 1980s, however, periodic boundary conditions became the standard choice. This shift was driven by the need to mimic an infinite, statistically homogeneous Universe and to avoid artificial edge effects, while also enabling efficient particle-mesh (PM) and hybrid TreePM algorithms that rely on Fourier methods (Klypin & Shandarin 1983; Xu 1995). The three-torus (T³ =S¹ × S¹ × S¹) topological manifold is homeomorphic to the Cartesian product of three circles, and it is often described as a periodic box. Since the Lagrangian of a self-gravitating system in T³ is invariant to translation, the linear momentum is conserved in such a system. In the late 20th century, the notion that our Universe has this sort of topological manifold was considered a real possibility (Zeldovich & Starobinskij 1984; Broadhurst et al. 1990; Einasto et al. 1997), but modern large-volume galaxy surveys and cosmic microwave background (CMB) observations did not find any signs of periodicity (Hawkins et al. 2002; Planck Collaboration XVIII 2016) on the observed scales. Periodic boundary conditions are widely used in such codes as GADGET4 (Springel et al. 2022), AREPO (Weinberger et al. 2020), and PKDGRAV3 (Potter et al. 2017) and remain the dominant approach in large-scale cosmological simulations such as Millennium (Springel et al. 2005), Illustris (Genel et al. 2014), and Euclid flagship (Euclid Collaboration: Castander et al. 2025).

Despite their practicality, the three-torus topology is not an ideal representation of the Universe with trivial ℝ³ topological manifold. While the cubic T³ topology preserves statistical homogeneity, it breaks statistical isotropy, since the cube is invariant only under the octahedral group O_h, not the full continuous rotation group SO(3). According to Noether’s theorem, the angular momentum is generally not conserved in a periodic simulation box, since the Lagrangian of a particle system in this topology is not invariant under continuous rotation. This symmetry-breaking also introduces subtle anisotropies in the gravitational field and can bias clustering statistics, especially when the simulation box is not much larger than the scales of interest. In practice, these effects are mitigated by choosing sufficiently large volumes (Rácz et al. 2021).

To address these limitations of the most used T³ topology, alternative topological manifolds must be explored. Rácz et al. (2018) proposed the STEreographically Projected cosmological Simulations (StePS) algorithm, which uses a spherically symmetric ℝ³ topological manifold with spatial compactification. This approach preserves full rotational invariance and implements a radially varying resolution: high in the simulation centre and gradually decreasing towards the outskirts, reaching the scale of homogeneity to avoid boundary effects. In this space, the full angular momentum vector of a self-gravitating system is conserved. In this paper, we introduce a new hybrid simulation topology that combines the fully periodic and spherical StePS approaches: the cylindrical S¹ × ℝ² topological manifold. In this configuration, the simulation domain is periodic along one axis (the cylinder’s axis) and open in the radial direction. This design enables a zoom-in strategy similar to the spherical StePS approach, with high resolution near the central axis and decreasing resolution outward, while retaining periodicity in one dimension. The Lagrangian of a self-gravitating particle system in this space is invariant under translation along the cylinder’s axis and rotations in the ℝ² plane, which are described by the SO(2) group. As a consequence of Noether’s theorem, the linear momentum along the cylinder’s axis and the angular momentum along the same axis are both conserved. Such a topology is particularly well-suited for studying naturally anisotropic environments such as filamentary structures or anisotropic cosmological models with the same symmetries.

This paper is organised as follows: in Sect. 2, we describe the base cosmological simulation algorithm. In Sect. 3, we present the initial condition generation method. Then, in Sect. 4, we introduce the first lambda cold dark matter (ΛCDM) simulations in cylindrical S¹ × ℝ² topological manifold. Finally, in Sect. 5, we summarise our results.

2. Algorithm

2.1. Equations of motion

The Newtonian equations of motion in comoving coordinates of a system composed of N particles can be written as

$$\begin{array}{c}\hfill m_i\ddot{{\mathit{x}}_i}={\displaystyle \sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^N}\frac{m_i{m}_j\mathit{F}({\mathit{x}}_i-{\mathit{x}}_j,h_i+h_j)}{a{(t)}^3}-2m_{i}H(t){\dot{\mathit{x}}}_i,\end{array}$$(1)

where x_i, m_i, and h_i are the coordinate vector, mass, and softening length of the particle with i index, a(t) is the scale factor determined by the background cosmological model,

$$\begin{array}{c}\hfill H(t)=\frac{\dot{a}(t)}{a(t)}\end{array}$$(2)

is the Hubble parameter, and F(x_i − x_j, h_i + h_j) is the smoothed gravitational force between particles i and j. This force depends on the ℱ(r, h) force softening and the background topology. In the StePS code, the softening,

$$\begin{array}{c}\hfill \mathcal{F}(r,h)=\{\begin{array}{cc}\frac{32r^4}{h^6}-\frac{38.4r^3}{h^5}+\frac{32r}{3h^3}\hfill & \text{if}r<\frac{h}{2},\hfill \\ \&\hfill \\ -\frac{32r^4}{3h^6}+\frac{38.4r^3}{h^5}-\hfill \\ frac48r^2{h}^4+\frac{64r}{3h^3}-\frac{1}{15r^2}\hfill & \text{if}\frac{h}{2}<r<h,\hfill \\ \&\hfill \\ \frac{1}{r^2}\hfill & \text{if}r>h,\hfill \end{array}\end{array}$$(3)

is used (Rácz et al. 2019), which corresponds to the cubic spline kernel of Monaghan & Lattanzio (1985). The shape of the F(r, h) force function is determined by the background topology.

For a trivial, flat Euclidean ℝ³ space, the smoothed gravitational force law can be written as

$$\begin{array}{c}\hfill {\mathit{F}}_{{\mathbb{R}}^3}(\mathit{r},h)=-G\mathcal{F}(\left|\mathit{r}\right|,h)\frac{\mathit{r}}{\left|\mathit{r}\right|},\end{array}$$(4)

where G is the gravitational constant. This formula is equivalent to Newton’s law of universal gravitation if r > h. The standard flat ΛCDM cosmological model assumes this isotropic gravity, and it is possible to use this directly in the StePS simulation code. In the case of comoving simulations in this topology, the right hand side of Eq. (1) has to be extended with a simple radial term,

$$\begin{array}{c}\hfill {\mathit{F}}_r({\mathit{x}}_i)=+m_i\frac{4\pi G}{3}{\rho }_0{\mathit{x}}_i,\end{array}$$(5)

that prevents the system from collapsing, where ρ₀ is the mean matter density in the simulation volume. For the detailed calculations, we refer to Rácz et al. (2018).

While the cosmological principle demands isotropic gravity, most cosmological simulations use non-trivial T³ topology for numerical convenience. In a periodic cubic box with a linear size, L_box, multiple periodic images of the particles have to be taken into account with

$$\begin{array}{c}\hfill {\mathit{F}}_{{\mathrm{T}}^3}(\mathit{r},h)={\displaystyle \sum _{\mathit{n}}}-G\mathcal{F}(|\mathit{r}-\mathit{n}{L}_{\mathrm{box}}|,h)\frac{\mathit{r}-\mathit{n}{L}_{\mathrm{box}}}{|\mathit{r}-\mathit{n}{L}_{\mathrm{box}}|},\end{array}$$(6)

summation formula, where n = {n₁, n₂, n₃}∈ℤ³, |r − L_box ⋅ n|< N_cut ⋅ L_box, and N_cut is the cut-off distance of the summation in linear box size units. In practice, this formula converges notoriously slowly and the number of periodic images scales with the third power of N_cut. To overcome this limitation, the Ewald (1921) summation method is used, in which the force is split into rapidly convergent real-space and reciprocal-space parts. While a straightforward implementation of Ewald summation provides sufficient accuracy, evaluating these sums requires hundreds or even thousands of floating-point operations for every interaction. This can be dramatically reduced by pre-computing the difference between the direct (nearest periodic image) forces and the fully periodic Ewald forces on a cubic grid (Hernquist et al. 1991). By constructing a lookup table from these values, one can use interpolation to efficiently add the Ewald correction term to the quasi-periodic force as

$$\begin{array}{c}\hfill {\mathit{F}}_{{\mathrm{T}}^3}(\mathit{r},h)={\mathit{F}}_{{\mathbb{R}}^3}(\mathit{r},h)+{\mathit{D}}_{\mathrm{Ewald}}(\mathit{r}),\end{array}$$(7)

where r is the separation vector between the nearest image in the T³ manifold and D_Ewald(r) is the Ewald correction term. In practice, the long-range forces can be calculated with Fast Fourier Transform (FFT) (Cooley & Tukey 1965; Hockney & Eastwood 1988), which corresponds to N_cut → ∞, but this method is not implemented in StePS, since this topology is not the main focus of our code. In the T³ topological manifold, the linear momentum is conserved due to translational invariance. However, angular momentum is not conserved because the system’s Lagrangian lacks symmetry under continuous rotation. This is reflected in the force law: the force vector between two particles depends not only on their relative positions, but also on their orientation relative to the axes of the periodic box. For the detailed description of T³ periodic effects in cosmological simulations, we refer to Rácz et al. (2021).

In this work, we propose another non-trivial topology for cosmological simulations, the cylindrical S¹ × ℝ² topological manifold. This space is periodic only along one axis with a periodicity of L_z. For simplicity, we chose this as the ‘z’-axis and

$$\begin{array}{c}\hfill {\mathit{F}}_{{\mathrm{S}}^1\times {\mathbb{R}}^2}(\mathit{r},h)={\displaystyle \sum _{\begin{array}{c}n=-N_{\mathrm{cut}}\\ |r_z-n{L}_z|<N_{\mathrm{cut}}{L}_z\end{array}}^{N_{\mathrm{cut}}}}-G\mathcal{F}(|\mathit{r}-n{L}_z{\mathit{e}}_z|,h)\frac{\mathit{r}-n{L}_z{\mathit{e}}_z}{|\mathit{r}-n{L}_z{\mathit{e}}_z|}\end{array}$$(8)

summation describes the gravitational force field of a particle in this manifold, where e_z = {0, 0, 1} is the unit vector in the z direction, r_z is the z component of the r vector, and N_cut is the cut-off distance of the summation in linear cylinder size units. While this summation formula still converges slowly, the contribution of distant periodic images now can be directly summed, since the number of images scales only with N_cut. Nevertheless, the Ewald summation method still can be used in S¹ × ℝ² for a more accurate force calculation. A comparison between the gravitational forces in ℝ³, T³, and S¹ × ℝ² manifolds are shown in Fig. 1. It is important to note that while the Newtonian gravitational forces in a T³ space are always weaker than in ℝ³, the S¹ × ℝ² forces are weaker only in axial directions, whereas they are always stronger in radial directions.

Fig. 1. Comparison of gravitational force laws between the periodic cubic (T3) and the periodic cylindrical (S1 × ℝ2) geometry. Each plot demonstrates the deviation from the Newton’s law of universal gravitation due to the periodic topology. In both cases, a test particle was placed to the centre of the simulation volume, and the gravitational force field was calculated both in the corresponding and in ℝ3 topology. The heat maps are showing the |F|/|Fℝ3| ratio of the force field magnitudes. The vector plots are demonstrating the F − Fℝ3 difference of the gravitational force fields around the particle. The length of the arrows is proportional to the magnitude of the difference between the forces.

In comoving cosmological simulations, we only follow the evolution of cosmic structures in comoving coordinates. Since Newtonian gravity is an attractive force, we have to introduce an external force that keeps the simulation from collapsing into the central axis. The net force in every point should be exactly zero for homogeneous matter distribution in comoving coordinates. This can be achieved, if we compensate for the gravitational field, F_cylinder(r), of the homogeneous cylinder in every point. Due to the construction of the summation expressed in Eq. (8), all particles experience the gravitational pull of a cylinder of a height of 2 ⋅ N_cut ⋅ L_z, as if they were in the central plane of this cylinder. This gravitational force is calculated by Chang (1988) and takes the form of

$$\begin{array}{cc}& F_h(\varrho ,R,L)\hfill \\ \hfill & =2G{\rho }_0[-2{\displaystyle {\int }_0^R}ln\left(\frac{\sqrt{L^2+R^2+{\varrho }^2+2\varrho \sqrt{R^2-Y^2}}+L}{\sqrt{L^2+R^2+{\varrho }^2-2\varrho \sqrt{R^2-Y^2}}+L}\right)dY\hfill \\ \hfill & +{\displaystyle {\int }_0^R}ln\left(\frac{\sqrt{R^2+{\varrho }^2+2\varrho \sqrt{R^2-Y^2}}}{\sqrt{R^2+{\varrho }^2-2\varrho \sqrt{R^2-Y^2}}}\right)dY],\hfill \end{array}$$(9)

for particles within the simulation volume, where $\varrho =\sqrt{x^2+y^2}$ is the radial coordinate measured from the central axis, L = N_cut ⋅ L_z, while ρ₀ is the mean matter density and R is the radius of the simulation volume. If N_cut → ∞, then the formula above can be simplified to

$$\begin{array}{c}\hfill F_h(\varrho )=-2\pi G{\rho }_0\varrho .\end{array}$$(10)

Using these formulae, we can write the equations of motion of a self-gravitating N-body system in S¹ × ℝ² topology in comoving coordinates as

$$\begin{array}{cc}\hfill m_i\ddot{{\mathit{x}}_i}& ={\displaystyle \sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^N}\frac{m_i{m}_j{\mathit{F}}_{{\mathrm{S}}^1\times {\mathbb{R}}^2}({\mathit{x}}_i-{\mathit{x}}_j,h_i+h_j)}{a{(t)}^3}-2m_{i}H(t)·{\dot{\mathit{x}}}_i\hfill \\ & -m_i{F}_h(\varrho ,R,N_{\mathrm{cut}}{L}_z).\hfill \end{array}$$(11)

This coupled, second-order differential equation system is solved in our StePS implementation using a kick-drift-kick (KDK) leapfrog integrator (Klypin & Shandarin 1983), with cosmic time, t, as the integration parameter.

2.2. Ewald summation in S¹ × ℝ² topological manifold

To precisely calculate the periodic forces in cylindrical S¹ × ℝ² topological manifold, one must employ an Ewald summation along the one periodic direction. This problem has been explored in molecular dynamics (Grzybowski & Bródka 2002) and electrostatics (Tornberg 2015; Saffar Shamshirgar & Tornberg 2017), where these boundary conditions are useful in describing charged particles in strings, pores or channels. In our case, the gravitational force between masses m_i and m_j separated by r is split into real- and reciprocal-space parts as

$$\begin{array}{c}\hfill {\mathit{F}}_{{\mathrm{S}}^1\times {\mathbb{R}}^2}(\mathit{r})={\mathit{F}}_{{\mathrm{S}}^1\times {\mathbb{R}}^2}^{\mathrm{real}}(\mathit{r})+{\mathit{F}}_{{\mathrm{S}}^1\times {\mathbb{R}}^2}^{\mathrm{Fourier}}(\mathit{r})+{\mathit{F}}_{{\mathrm{S}}^1\times {\mathbb{R}}^2}^{m=0}.\end{array}$$(12)

The real-space part of the Ewald split,

$$\begin{array}{c}\hfill {\mathit{F}}_{{\mathrm{S}}^1\times {\mathbb{R}}^2}^{\mathrm{real}}(\mathit{r})=-G{\displaystyle \sum _{n=-n_{\max }}^{n_{\max }}}\frac{{\mathit{r}}_n}{|{\mathit{r}}_n{|}^3}[\text{erfc}(\alpha |{\mathit{r}}_n|)+\frac{2\alpha }{\sqrt{\pi }}|{\mathit{r}}_n|e^{-{\alpha }^2{|{\mathit{r}}_n|}^2}],\end{array}$$(13)

is identical in form to the T³ case, but the image set now runs over one integer along the z-axis, where r_n = r + nL_ze_z and α is the Ewald-split parameter. In Fourier space, the z component of the wavevector $k_z=\frac{2\pi m}{L_z}$ is discrete due to periodicity in z, but continuous in the other components. Carrying out the two-dimensional Fourier integrals, the force contribution at a 𝜚 radial and z axial separation becomes

$$\begin{array}{c}\hfill {\mathit{F}}_{{\mathrm{S}}^1\times {\mathbb{R}}^2}^{\mathrm{Fourier}}(\varrho ,z)=-\frac{2G}{L_z}{\displaystyle \sum _{m=1}^{m_{\max }}}\mathbf{\nabla }[cos(k_{m}z)\mathcal{B}(k_m,\varrho )],\end{array}$$(14)

where the kernel,

$$\begin{array}{c}\hfill \mathcal{B}(k_m,\varrho )=e^{k_m\varrho }\text{erfc}(\frac{k_m}{2\alpha }+\alpha \varrho )+e^{-k_m\varrho }\text{erfc}(\frac{k_m}{2\alpha }-\alpha \varrho ),\end{array}$$(15)

avoids the use of modified Bessel functions by utilising the analytical Gaussian filtering instead. The m = 0 mode represents the background contribution of the periodic mass line, yielding a radial force of

$$\begin{array}{c}\hfill F_{\varrho }^{m=0}=-\frac{2G}{L_z\varrho }(1-e^{-{\alpha }^2{\varrho }^2}).\end{array}$$(16)

While this Ewald summation converges faster than the direct summation of the periodic images, it is still too slow to calculate directly between every particle-particle interaction. To overcome this issue in the StePS implementation, we build a 2D lookup table,

$$\begin{array}{c}\hfill \mathit{D}(\rho ,z)={\mathit{F}}_{{\mathrm{S}}^1\times {\mathbb{R}}^2}(\rho ,z)-{\mathit{F}}_{{\mathbb{R}}^3}(\rho ,z),\end{array}$$(17)

on an equally spaced grid in 𝜚 and z cylindrical coordinates, where F_ℝ³(𝜚, z) is the nearest-image Newtonian force. This correction vector is smoothly interpolated during the simulation with cloud-in-cell (CIC) or triangular shaped cloud (TSC) method and added to the nearest image interactions.

Since the direct summation of the periodic images in S¹ × ℝ² is possible, we used N_cut = 10⁷ real-space periodic images in both direction to construct a brute-force reference lookup-table. We used this to find the optimal α Ewald-split parameters for fixed n_max and m_max values. This can be seen in Fig. 2.

Fig. 2. Accuracy of S1 × ℝ2 periodic summation methods. For all tests, the correction lookup table D(ρ, z) was computed on a 144 × 128 grid, varying the truncation parameter Ncut or the Ewald parameters α, nmax, and mmax. Using these tables, we reconstructed the full periodic force field FS1 × ℝ2(ρ, z) = Fℝ3(ρ, z)+D(ρ, z) and compared it against a reference field. This reference field was generated via direct summation using Ncut = 107 periodic images, requiring approximately 104 s of computation time. Left: Mean force errors of direct periodic real-space summation as a function of Ncut. Although summation occurs only along a single direction (the z-axis), convergence remains exceptionally slow. Center: S1 × ℝ2 Ewald summation force accuracy as a function of the Ewald splitting parameter α with fixed nmax = 5 and mmax = 12 values. Using only 5 periodic images in real space and 12 modes in reciprocal space is sufficient to achieve single-precision (32-bit) floating-point accuracy when the optimal splitting parameter is applied. Right: Same as the center panel, but with fixed nmax = 7 and mmax = 14. Increasing both nmax and mmax by two improves the force accuracy by nearly three orders of magnitude (a factor of ∼600).

2.3. Background evolution

The background expansion of the Universe is described by the scale factor, a(t), and the Hubble parameter, H(t), and used in Eq. (11). These can be calculated for any given time by solving the Friedmann equation (Peacock 1999; Frieman et al. 2008)

$$\begin{array}{c}\hfill {\left(\frac{\dot{a}}{a}\right)}^2=H_0^2[{\mathrm{\Omega }}_{r,0}{a}^{-4}+{\mathrm{\Omega }}_{m,0}{a}^{-3}+{\mathrm{\Omega }}_{k,0}{a}^{-2}+{\mathrm{\Omega }}_{\mathit{DE}}(a)],\end{array}$$(18)

where Ω_r, 0 is the radiation,Ω_m, 0 is the non-relativistic matter, Ω_k, 0 is the curvature, and

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathit{DE}}(a)={\mathrm{\Omega }}_{DE,0}exp[-3{\displaystyle {\int }_1^a}\frac{1+w_{\mathit{DE}}(a\prime )}{a\prime }da\prime ]\end{array}$$(19)

is the dark energy density parameter. In a Chevallier-Polarski-Linder (CPL) parametrisation (Chevallier & Polarski 2001; Linder 2003), the w_DE(a) dark energy equation of state parameter is written as

$$\begin{array}{c}\hfill w_{\mathit{DE}}(a)=w_0+w_a(1-a),\end{array}$$(20)

where w₀ and w_a are free parameters of the model. The w_a = 0 corresponds to the wCDM model, while the cosmological constant (ΛCDM) is recovered when w₀ = −1 and w_a = 0 in this parametrisation. For the ΛCDM, wCDM, and general CPL cosmologies, the scale factor, a(t), is determined by solving Eq. (18) with a fourth-order Runge-Kutta method in our StePS implementation. This Friedmann solver uses the same timestep length as the N-body integrator. Our current implementation also supports the use of tabulated expansion histories, which can be loaded from an ASCII file. When using these tabulated histories, the code interpolates between the provided values. The interpolation scheme can be configured to use linear, quadratic, or cubic algorithms.

2.4. Force calculation

To solve Eq. (11), we have to calculate the acceleration of each particle at every timestep. The calculation of the Hubble drag and the boundary term scales with 𝒪(N), where N is the total number of particles in the simulation. The gravitational interaction between the particles is considerably more complex, since a naive approach requires evaluating all pairwise forces, leading to a computational cost of 𝒪(N²).

2.4.1. Direct

The first public version of the StePS code supported only direct force calculation, and this option remains available. This was feasible because the original ℝ³ StePS algorithm employs a radially decreasing resolution, allowing a large dynamical range to be covered with only a few million particles. In such a configuration, the computational cost of the 𝒪(N²) pairwise summation is manageable, especially when combined with efficient GPU parallelization. This approach leverages the zoom-in nature of the method: high resolution is concentrated near the centre, while outer regions are represented by fewer, more massive particles, significantly reducing the total N compared to uniform-resolution simulations. A direct force calculation was implemented in StePS for both CPUs and GPUs using MPI-OpenMP and MPI-OpenMP-CUDA parallelisation, respectively. In the newly implemented S¹ × ℝ² topological manifold, which applies a similar zoom-in strategy in cylindrical coordinates, direct summation remains practical only for smaller cylinders. Because the resolution along the periodic axis is constant, simulations of larger volumes become less efficient with this method. Nevertheless, direct summation is retained in the code as a benchmark due to its simplicity and minimal force-computation errors.

2.4.2. Octree

For larger particle numbers, the 𝒪(N²) force calculation method becomes inefficient even on GPUs. To address this, the new StePS version implements a Barnes & Hut (1986) octree force calculation method for both the original ℝ³ and the new cylindrical S¹ × ℝ² topological manifold to reduce the complexity to 𝒪(N log N). In our implementation, the octree domains are always cubic and the size of the root node is always twice that of the largest linear dimension of the simulation volume: either the cylinder’s height or diameter in the S¹ × ℝ² topology, or the simulation diameter in the ℝ³ mode.

To minimize systematic errors arising from spatial correlations of the octree geometry, the code introduces random shifts and rotations of the domain at every timestep. In cylindrical topology, the domain center is shifted randomly along the cylinder axis and rotated around the same axis by a random angle. In ℝ³ mode, the domain center is shifted randomly within a sphere of a radius R_sim/2 from the simulation center, while the octree is rotated around a random axis by a random angle. These transformations help decorrelate force errors and improve the statistical isotropy of the simulation.

In contrast to direct summation, the octree method offers a dramatic improvement in performance for large particle numbers, but at the cost of introducing small force errors that depend on tree parameters and particle distribution. This is particularly important for maintaining radial balance in both the S¹ × ℝ² and ℝ³ cases, where long-range force accuracy governs the stability of the system. While a direct summation automatically satisfies this requirement, the octree-based calculations can lead to spurious radial collapse or expansion if uncorrected. Such a radial instability would change the background density over time during a simulation, which could alter the linear and non-linear power spectrum in a cosmological case. The current StePS implementation performs a calibration step before the simulation starts: the initial particle load (glass) is used to measure the radial force error under the same octree parameters as the simulation, applying multiple random shifts and rotations. From these measurements, a radial lookup table can be constructed, which is then used during the run to correct octree forces by subtracting the interpolated radial error at every timestep. This correction ensures that the large-scale radial equilibrium is preserved in comoving coordinates without sacrificing the efficiency of the octree. We demonstrate the implemented octree force calculation accuracy with and without a radial correction in a cosmological simulation at different redshifts in Fig. 3.

Fig. 3. Force accuracy as a function of the octree opening angle, Θ, in a cylindrical cosmological volume at different redshifts. We evolved 6 ⋅ 104 particles in standard ΛCDM cosmology with direct force calculation in a periodic cylinder with Lz = 200.0 Mpc height and Rsim = 100.0 Mpc radius for this accuracy test. After the simulation was finished, we calculated the FBH octree forces with our StePS simulation code. The FD direct forces were used as a ground truth in the comparison. We calculated the octree forces with and without the radial force correction. Left: Relative total ⟨|FBH − FD|/|FD|⟩ force error. The magnitude of the force error is behaving similarly as in the well tested cubic T3 case. The new radial force correction calculated from the initial particle load increases the accuracy in all redshifts. Center: Mean radial ⟨FBH, 𝜚/FD, 𝜚⟩ force error. While the raw radial octree force have a non-zero expected value at high Θ values, the implemented radial correction fixes this in all redshifts for Θ ≤ 0.5. It is critical to use this correction, since a systematic bias in radial direction will cause contraction or expansion of the simulated volume, which will alter the structure formation due to the changed background density. Right: Mean tangential ⟨FBH, ϑ/FD, ϑ⟩ force error. In contrast to the radial case, the measured tangential force errors are consistent with a zero mean. A non-zero systematic bias in this component would cause errors in the total angular momentum.

The octree implementation is available only with MPI-OpenMP parallelization, making it suitable for CPU clusters. The GPU accelerated octree force calculation will be implemented in the future.

3. Generating initial conditions

The first step in running a cosmological N-body simulation is to generate initial conditions that are consistent with the chosen cosmological parameters. This is typically done in two separate steps: (i) prepare a particle distribution that represents a homogeneous density field and (ii) imprint the density fluctuations predicted by linear theory for the adopted cosmological model. In this section, we outline these steps for the S¹ × ℝ² topological manifold.

3.1. Simulation geometry and initial particle load

The main simulation challenge in non-finite-volume topological manifolds is to handle the edge effects where the simulation transitions from the unresolved to the resolved volume. The original spherical StePS algorithm addressed this in ℝ³ by compactifying space via the inverse stereographic projection. In S¹ × ℝ² topological manifold, we applied the same technique in two dimensions to the non-compact (x, y) directions.

The inverse stereographic projection is a bijection between the compact n-sphere Sⁿ and ℝⁿ (excluding the projection point Q). Here, we compactify ℝ² to S² and work in cylindrical coordinates (𝜚,φ, z). Only the radial coordinate $\varrho =\sqrt{x^2+y^2}$ is transformed as

$$\begin{array}{c}\hfill \omega =2arctan\left(\frac{\varrho }{D_s}\right),\end{array}$$(21)

where D_s is the diameter of the S² sphere used in the mapping and ω ∈ [0, π) is the new compact coordinate. The point ω → π corresponds to 𝜚 → ∞ and is excluded from the finite computational domain. This corresponds to the last volume element around Q. Next, we tessellate the compactified domain uniformly in the coordinates (ω, φ, z). To avoid excessively steep variations in the non-compact resolution, we enforced a constant mass resolution inside a cylindrical core ω < ω_c, where ω_c is the maximum compact coordinate for which constant mass resolution is maintained. We then placed a single particle at a random location in each compact cell and assigned it a mass of $m_i=\overline{\rho }{V}_i$, where $\overline{\rho }$ is the mean comoving density and V_i is the corresponding cell volume in the non-compact (𝜚,φ, z) coordinates. This ensures that the compact-to-physical Jacobian of the inverse stereographic map is properly accounted for. Once the sample particles are prepared, we transform the compact radius back to the non-compact radius using a stereographic projection,

$$\begin{array}{c}\hfill \varrho =D_{s}tan\frac{\omega }{2}·\end{array}$$(22)

Following this transformation, we have a zoom-in geometry in the non-compact space where the spatial resolution smoothly decreases with the radial 𝜚 coordinate, until we reach the edge of the resolved region (R_sim). The distribution of the particle masses for an example realisation as a function of the radial 𝜚 coordinate can be seen in Fig. 4. Using Eq. (22), the non-compact radius of constant resolution 𝜚_c can be calculated.

Fig. 4. Mass resolution as a function of radial coordinate 𝜚 in an example cylindrical glass. The parameters of this glass are the following: DS = 70 Mpc, Rsim = 500.0 Mpc, ωc = 1.63973 RAD, Lz = 500 Mpc, and the total particle number is Npart = 1.2 ⋅ 106. The glass has constant resolution for 𝜚 < 𝜚c = DStan(ωc/2) = 75.0 Mpc to avoid unwanted mixing of particles with widely different masses in the simulations.

To use this distribution as an initial particle load during the initial condition generation, we required the net comoving gravitational acceleration on each particle to be approximately zero. The usual method for this is to use this particle distribution as an initial condition for glass-making (White 1994). The StePS simulation code is able to generate such glasses in S¹ × ℝ² topology by running a comoving cosmological N-body simulation with repulsive gravity. During the simulation, the particles are losing energy due to the Hubble-drag term in Eq. (11), and after 𝒪(10⁵ − 10⁶) expansion factors, they are settling down in a glass-like structure where the net gravitational acceleration of every particle is negligible.

3.2. Initial perturbations

Once the initial particle load is ready, we generated the initial perturbations using the stepsic2 code (Pál et al. in prep.)¹. This program supports T³, ℝ³, and S¹ × ℝ² topologies and can be applied in either the Zel’dovich approximation (1LPT) or second-order Lagrangian perturbation theory (2LPT). Given the cosmology, we can compute the linear matter power spectrum at the time of the initial scale factor (a_init) as

$$\begin{array}{c}\hfill P_{\mathrm{lin}}(k,a_{\mathrm{init}})=D^2(a_{\mathrm{init}})T^2(k)P_0(k),\end{array}$$(23)

where P₀(k) is the primordial power spectrum, T(k) is the transfer function, and D(a) is the linear growth function. We then realised a Gaussian random field, δ(k), following the standard Rayleigh–phase sampling (Sirko 2005). For each independent Fourier mode, we draw a random phase θ_k ∼ 𝒰(−π, π) and an amplitude, A_k, from a Rayleigh distribution with scale parameter, σ_k, chosen so that the real and imaginary parts of δ(k) end up independent Gaussian variates with zero mean and variance σ_k² = P_lin(|k|,a_init)/2. We then enforced the Hermitian symmetry δ(−k) = δ^*(k) to ensure a real configuration-space field and proceeded to compute the displacement and velocity fields.

The initial comoving positions, x_i, and velocities, ẋ_i, were obtained by displacing the Lagrangian particle coordinates, q_i, of the initial particle load according to the Zel’dovich expression,

$$\begin{array}{cc}\hfill {\mathit{x}}_i& ={\mathit{q}}_i+{\displaystyle \int \frac{i\mathit{k}}{k^2}}\delta (\mathit{k})e^{i{\mathit{q}}_i\mathit{k}},\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill {\dot{\mathit{x}}}_i& =\frac{\dot{D}(a_{\mathrm{init}})}{D(1)}{\displaystyle \int \frac{i\mathit{k}}{k^2}}\delta (\mathit{k})e^{i{\mathit{q}}_i\mathit{k}},\hfill \end{array}$$(25)

or, alternatively, with a 2LPT approximation (Buchert 1994; Jenkins 2010).

4. Results

4.1. ΛCDM simulation

To validate the new simulation method, we run a new cosmological N-body simulation with the updated version of our StePS code in S¹ × ℝ² topology. We used the best-fit Planck 2018 (Planck Collaboration VI 2020) ΛCDM cosmological parameters and a 2LPT method to generate the initial conditions of a periodic cylinder with a height of L_z = 1.0 Gpc and resolved radius of R_sim = 500 Mpc. As an initial particle load, we used a cylindrical glass with height 500 Mpc and radius = 500 Mpc, repeated along the z axis. The initial power spectrum was calculated by CAMB (Code for Anisotropies in the Microwave Background) (Lewis et al. 2000; Lewis & Challinor 2011) at z = 0, and backscaled to the z_init = 49 initial redshift with the linear growth function. The parameters of the simulation can be seen in Table 1. The simulation was run on the CSC Puhti² supercomputer using 640 CPU cores and 3085 timesteps. On average, one integration timestep took 60 seconds of wall-clock time. Overall 138 particle snapshots were stored between redshifts 31 and 0. The pySPHviewer (Benítez-Llambay 2017) visualisation of the dark matter distribution in this simulation in different slices and redshifts can be seen in Fig. 5. In our S¹ × ℝ² topology, a cosmic web forms from the small initial density fluctuations due to gravity, similarly to the T³ and ℝ³ simulations.

Fig. 5. Thin slices of the dark matter density field from the first S1 × ℝ2 StePS simulation at z = 11, z = 3, and z = 0 redshifts in comoving coordinates. Each slice has a thickness of 15.0 Mpc and each column shows a different slice orientation. The late-time non-linear cosmic web is visually similar to that found in T3 and ℝ3 simulations: voids, walls, filaments, and dark matter halos are clearly visible. Due to the zoom-in nature of the StePS S1 × ℝ2 simulation method, the spatial resolution is highest along the central axis, while only the largest structures are resolved near the Rsim simulation radius.

4.2. Power spectra

The isotropic P(k) power spectrum of the dark matter density field is a standard statistic used to quantify the clustering. This is defined as

$$\begin{array}{c}\hfill \langle \delta (\mathit{k}){\delta }^{\ast }(\mathit{k}\prime )\rangle ={(2\pi )}^3{\delta }_{\mathrm{D}}(\mathit{k}-\mathit{k}\prime )P(k),k\equiv |\mathit{k}|,\end{array}$$(26)

where δ_D is the Dirac delta function and the angle brackets denote an ensemble average and

$$\begin{array}{c}\hfill \delta (\mathit{k})={\displaystyle \int \delta }(\mathit{x})e^{-i\mathit{k}·\mathit{x}}{\mathrm{d}}^{3}x={\displaystyle \int \frac{{\rho }_{\mathrm{DM}}(\mathit{x})-{\overline{\rho }}_{\mathrm{DM}}}{{\overline{\rho }}_{\mathrm{DM}}}}{e}^{-i\mathit{k}·\mathit{x}}{\mathrm{d}}^{3}x\end{array}$$(27)

is the Fourier transform of the dark matter density contrast of the simulation. In a periodic T³ volume, P(k) can be estimated efficiently using FFTs on a uniform grid. In the S¹ × ℝ² topology, however, the StePS zoom-in particle load implies a strongly position-dependent effective sampling and a non-trivial window, so making a direct FFT-based estimate on a single uniform grid is not straightforward. We therefore estimated the power spectra from particle snapshots using the Feldman–Kaiser–Peacock (FKP) method (Feldman et al. 1994). The FKP estimator was originally developed to measure galaxy power spectra from redshift surveys with non-uniform selection functions and complex survey geometries. It constructs a weighted fluctuation field by subtracting a synthetic, unclustered ‘random’ catalogue that follows the same selection function and window as the data and then it applies a Fourier transform to this weighted field. This random catalogue in our case is the initial conditions that was used during glass making. The estimated power spectra of a few selected snapshots can be seen in Fig. 6. The estimated power spectra follow the back-scaled CAMB predictions on linear scales, while showing a good agreement with the non-linear Mead et al. (2021) halofit power spectrum emulator, which was calibrated on large cubic T³ simulations.

Fig. 6. Dark matter power spectra of four selected snapshots (solid lines) with the corresponding theoretical linear power spectra (dashed lines) and the non-linear Mead et al. (2021) halofit power spectra (dotted dashed lines). Our S1 × ℝ2 StePS simulation is capable of following both linear and non-linear structure formation.

4.3. Dark matter halos

The dark matter halos are gravitationally bound dark matter regions that have decoupled from the cosmic expansion and collapsed (Wechsler & Tinker 2018). While these halos are the prospective hosts of galaxies, the present runs are dark-matter-only, so individual galaxies are not explicitly resolved in our simulations. As a consequence of the S¹ × ℝ² zoom-in configuration, low-mass halos and the fine internal structure of halos can be captured with high fidelity near the axis, while more massive halos are resolvable within larger volumes, according to the M(𝜚) mass resolution profile in Fig. 4. As existing halo finders are not designed for S¹ × ℝ² zoom-in geometries, we implemented our own spherical-overdensity (SO) halo finder, StePS_HF (Rácz et al. in prep.).

StePS_HF follows the SO algorithm (Tormen et al. 1998; Despali et al. 2016) widely used in cosmological analyses, adapted to the S¹ periodicity along the z direction and the non-compact, radially varying the resolution in the transverse plane. After loading one snapshot, we reconstructed a local density for each particle using a k-nearest-neighbour estimator that is robust to inhomogeneous sampling. We then identified the highest-density unflagged particle as a provisional centre and grew a sphere around it, accumulating mass until the mean enclosed density first dropped below the chosen threshold, $\overline{\rho }(<R_{\mathrm{vir}})={\rho }_{\mathrm{vir}}(z)$, where ρ_vir is the redshift-dependent virial density from spherical collapse (Bryan & Norman 1998). Whenever the resulting candidate contained at least N_p, min particles, we refined its centre by replacing the seed position with the centre of mass of the innermost N_p, min particles and recomputing essential parameters (R_vir, M_vir), other mass definitions (M_200c, M_200b, M_500c, M_1000c, etc.), the corresponding radii, and other halo properties such as the peculiar velocity (V), Klypin scale radius (R_Klypin, Klypin et al. 2011), concentration (c), angular momentum (J), Bullock-spin (λ_B, Bullock et al. 2001), and total energy (E). All particles within R_vir were then flagged as members and the process was repeated on the remaining unflagged particles until no seeds above threshold remain.

We used N_p, min = 15 as the minimum threshold for halo particles and identified dark matter halos in late-time snapshots. While this minimal particle number is usually sufficient to identify halo positions, a higher threshold is needed for the reliable measurement of structural parameters (Knebe et al. 2011). At the final state (redshift 0), we found 1.1 ⋅ 10⁵ dark matter halos. The spatial distribution of these halos is visualised in Fig. 7 at three selected redshifts. Overall, 68 parameters were calculated for each identified halo.

Fig. 7. Spatial distribution of dark matter halos in the new S1 × ℝ2 StePS simulation at different redshifts, in comoving coordinates. Due to the zoom-in geometry, and the fact that low-mass halos form first, halos first appear close to the central axis, where the resolution is highest, and the smallest halos can be resolved. As structure formation proceeds towards higher masses, halos can later be found at larger radii, where the resolution is lower, and only more massive halos can be resolved. Towards larger $\varrho =\sqrt{x^2+y^2}$ only increasingly massive halos are resolved, so the apparent decline in number density at large 𝜚 is a selection effect.

4.4. Identifying signals of anisotropy

A recent analysis of the James Webb Space Telescope (JWST) Advanced Deep Extragalactic Survey (JADES) reported a statistically significant deviation from isotropy in galaxy spin directions (Shamir 2025). Motivated by these findings and by the fact that the gravity is inherently anisotropic in S¹ × ℝ² topology, we have aimed to detect corresponding late-time signatures of anisotropy in our simulation. In our setup, the compact direction is the periodic z-axis and the manifold retains an SO(2) rotational symmetry around this axis. Since we assumed statistical isotropy during the initial condition generation, any global anisotropy at the halo scales should appear as a deviation from uniformity in the distribution of the angle

$$\begin{array}{c}\hfill {\theta }_{\mathit{Jz}}={cos}^{-1}\left(\frac{J_z}{|\mathit{J}|}\right)\end{array}$$(28)

between the halo angular-momentum vector, J, and the z-axis. We measured J for every halo with N_part ≥ 100 particles to ensure accurate angular momenta, binned the halos by mass and tested the distribution against the isotropic expectation and plotted the results in Fig. 8. Isotropy was quantified by testing uniformity of cos θ_Jz with Kolmogorov–Smirnov tests in each mass bin. Our result shows that for L_z = 1.0 Gpc comoving periodicity, within our statistical precision and for all mass bins considered, the halo spin–axis distribution is consistent with statistical isotropy in our ΛCDM simulation. We emphasise that smaller compact lengths could in principle enhance anisotropy, but such values are strongly constrained by CMB observations (Planck Collaboration XVIII 2016) and by modern large-scale galaxy surveys (Roukema & Edge 1997; Hawkins et al. 2002; Luminet 2016). Our result shows that the S¹ × ℝ² topology, at least for the comoving L_z = 1.0 Gpc periodicity, does not provide a natural explanation for the JWST rotation-sense asymmetry reported by Shamir (2025).

Fig. 8. Distribution of the halo spin orientation angle, θJz, with respect to the simulation z-axis, shown for three mass bins at redshift z = 0.5 in our S1 × ℝ2ΛCDM simulation. We only used halos with Npart > 100 in this analysis to reliably measure the angular momentum vectors. The histogram uses equal-area binning on the unit sphere, so that an isotropic spin distribution implies equal expected counts per bin. Each panel is annotated with the one-sample Kolmogorov-Smirnov test on $u=\frac{1}{2}(1+cos({\theta }_{\mathit{Jz}}))\in [0,1]$ against the uniform null. The three mass bins show no departures from the isotropic expectation. We found similar distributions in all other output redshifts and mass bins with no clear signs of anisotropy.

5. Summary

We present a new capability of the StePS simulation code that enables compactified cosmological N-body simulations in S¹ × ℝ² topology. We summarise the large-scale modifications of the gravitational field induced by this slab periodic boundary condition, described the simulation algorithm, and outline the initial condition generation procedure. We present the first S¹ × ℝ²ΛCDM simulation run with our StePS code and show that it reproduces the expected linear and non-linear growth of structures. The ΛCDM simulation particle outputs, post-processed data, and example notebooks are freely available on our project webpage³.

The new simulation topology offers several practical advantages compared to the usual T³ setup:

• Since the angular momentum in this topology is conserved along the axis of periodicity, the S¹ × ℝ² StePS simulations are ideal for studying systems where this quantity plays an important role; for instance, in the context of cosmic filaments and their interactions with dark matter halos or anisotropic cosmological models, such as Bianchi type VII models.

• The zoom-in implementation improves the balance between small- and large-scale modes relative to a constant-resolution T³ run: small scales are resolved at high fidelity only within ω < ω_c around the central axis, while long-wavelength modes are sampled over substantially larger volumes.

• Similarly, the halo catalogues sample the smallest halos near the axis, while more massive and rarer halos are sampled from much larger volumes due to the zoom-in strategy in ℝ².

These benefits come with some limitations:

• Since almost all simulation post-processing pipelines assume T³ geometry, processing the outputs of StePS usually requires heavy modifications in existing pipelines.

• The zoom-in strategy introduces selection effects and sampling challenges that are absent in constant-resolution T³ simulations and these must be accounted for.

• The anisotropic force field induced by the topology can imprint artefacts, similarly to known issues in the T³ case (Rácz et al. 2021). Based on preliminary tests, such effects might become non-negligible at late times if the periodic length is small, L_z ≲ 200 Mpc. A detailed characterisation of these systematics will be presented in future work.

We do not expect that the new algorithm will replace the well-tested cubic T³ simulations; rather, the aim of this work is to provide an alternative for cases where cubic symmetry can limit the results. This new simulation method is now part of the version 2 release of StePS along the ℝ³ and T³ topological manifolds. This simulation code is freely available as an open-source project at our github⁴ repository under GNU General Public License v2.0.

Acknowledgments

GR and TS acknowledge support of the Research Council of Finland grant 354905 and support by the European Research Council via ERC Consolidator grant KETJU (no. 818930). This work has been supported by the Hungarian National Research, Development and Innovation Office (NKFIH, grant OTKA NN147550). IS acknowledges NASA ROSES grants 80NSSC24K1489 and 24-ADAP24-0074, contract number 80NM0018F0610 via a JPL sub-award. The authors wish to acknowledge CSC – IT Center for Science, Finland, for computational resources.

References

All Tables

All Figures

Fig. 1. Comparison of gravitational force laws between the periodic cubic (T3) and the periodic cylindrical (S1 × ℝ2) geometry. Each plot demonstrates the deviation from the Newton’s law of universal gravitation due to the periodic topology. In both cases, a test particle was placed to the centre of the simulation volume, and the gravitational force field was calculated both in the corresponding and in ℝ3 topology. The heat maps are showing the |F|/|Fℝ3| ratio of the force field magnitudes. The vector plots are demonstrating the F − Fℝ3 difference of the gravitational force fields around the particle. The length of the arrows is proportional to the magnitude of the difference between the forces.

In the text

Fig. 2. Accuracy of S1 × ℝ2 periodic summation methods. For all tests, the correction lookup table D(ρ, z) was computed on a 144 × 128 grid, varying the truncation parameter Ncut or the Ewald parameters α, nmax, and mmax. Using these tables, we reconstructed the full periodic force field FS1 × ℝ2(ρ, z) = Fℝ3(ρ, z)+D(ρ, z) and compared it against a reference field. This reference field was generated via direct summation using Ncut = 107 periodic images, requiring approximately 104 s of computation time. Left: Mean force errors of direct periodic real-space summation as a function of Ncut. Although summation occurs only along a single direction (the z-axis), convergence remains exceptionally slow. Center: S1 × ℝ2 Ewald summation force accuracy as a function of the Ewald splitting parameter α with fixed nmax = 5 and mmax = 12 values. Using only 5 periodic images in real space and 12 modes in reciprocal space is sufficient to achieve single-precision (32-bit) floating-point accuracy when the optimal splitting parameter is applied. Right: Same as the center panel, but with fixed nmax = 7 and mmax = 14. Increasing both nmax and mmax by two improves the force accuracy by nearly three orders of magnitude (a factor of ∼600).

In the text

Fig. 3. Force accuracy as a function of the octree opening angle, Θ, in a cylindrical cosmological volume at different redshifts. We evolved 6 ⋅ 104 particles in standard ΛCDM cosmology with direct force calculation in a periodic cylinder with Lz = 200.0 Mpc height and Rsim = 100.0 Mpc radius for this accuracy test. After the simulation was finished, we calculated the FBH octree forces with our StePS simulation code. The FD direct forces were used as a ground truth in the comparison. We calculated the octree forces with and without the radial force correction. Left: Relative total ⟨|FBH − FD|/|FD|⟩ force error. The magnitude of the force error is behaving similarly as in the well tested cubic T3 case. The new radial force correction calculated from the initial particle load increases the accuracy in all redshifts. Center: Mean radial ⟨FBH, 𝜚/FD, 𝜚⟩ force error. While the raw radial octree force have a non-zero expected value at high Θ values, the implemented radial correction fixes this in all redshifts for Θ ≤ 0.5. It is critical to use this correction, since a systematic bias in radial direction will cause contraction or expansion of the simulated volume, which will alter the structure formation due to the changed background density. Right: Mean tangential ⟨FBH, ϑ/FD, ϑ⟩ force error. In contrast to the radial case, the measured tangential force errors are consistent with a zero mean. A non-zero systematic bias in this component would cause errors in the total angular momentum.

In the text

Fig. 4. Mass resolution as a function of radial coordinate 𝜚 in an example cylindrical glass. The parameters of this glass are the following: DS = 70 Mpc, Rsim = 500.0 Mpc, ωc = 1.63973 RAD, Lz = 500 Mpc, and the total particle number is Npart = 1.2 ⋅ 106. The glass has constant resolution for 𝜚 < 𝜚c = DStan(ωc/2) = 75.0 Mpc to avoid unwanted mixing of particles with widely different masses in the simulations.

In the text

Fig. 5. Thin slices of the dark matter density field from the first S1 × ℝ2 StePS simulation at z = 11, z = 3, and z = 0 redshifts in comoving coordinates. Each slice has a thickness of 15.0 Mpc and each column shows a different slice orientation. The late-time non-linear cosmic web is visually similar to that found in T3 and ℝ3 simulations: voids, walls, filaments, and dark matter halos are clearly visible. Due to the zoom-in nature of the StePS S1 × ℝ2 simulation method, the spatial resolution is highest along the central axis, while only the largest structures are resolved near the Rsim simulation radius.

In the text

	Fig. 6. Dark matter power spectra of four selected snapshots (solid lines) with the corresponding theoretical linear power spectra (dashed lines) and the non-linear Mead et al. (2021) halofit power spectra (dotted dashed lines). Our S1 × ℝ2 StePS simulation is capable of following both linear and non-linear structure formation.
In the text

Fig. 7. Spatial distribution of dark matter halos in the new S1 × ℝ2 StePS simulation at different redshifts, in comoving coordinates. Due to the zoom-in geometry, and the fact that low-mass halos form first, halos first appear close to the central axis, where the resolution is highest, and the smallest halos can be resolved. As structure formation proceeds towards higher masses, halos can later be found at larger radii, where the resolution is lower, and only more massive halos can be resolved. Towards larger $\varrho =\sqrt{x^2+y^2}$ only increasingly massive halos are resolved, so the apparent decline in number density at large 𝜚 is a selection effect.

In the text

Fig. 8. Distribution of the halo spin orientation angle, θJz, with respect to the simulation z-axis, shown for three mass bins at redshift z = 0.5 in our S1 × ℝ2ΛCDM simulation. We only used halos with Npart > 100 in this analysis to reliably measure the angular momentum vectors. The histogram uses equal-area binning on the unit sphere, so that an isotropic spin distribution implies equal expected counts per bin. Each panel is annotated with the one-sample Kolmogorov-Smirnov test on $u=\frac{1}{2}(1+cos({\theta }_{\mathit{Jz}}))\in [0,1]$ against the uniform null. The three mass bins show no departures from the isotropic expectation. We found similar distributions in all other output redshifts and mass bins with no clear signs of anisotropy.

In the text