Skip to main content
Sandbox Physics
From ripples to the cosmic web

From ripples to the cosmic web · L2

How Mesh Gravity Follows Folding Matter

Sandbox Physics ·

A particle-mesh method deposits the current mass on nodes, solves for its gravitational field and samples that force back at the particles. Repeating these operations lets the changed density influence motion after crossing.

Distribute mass between the nodes

Assigning all mass to the nearest node would make the deposited field jump when a particle crossed the halfway point. Cloud-in-cell, or CIC, shares its mass between neighbouring nodes: two in one dimension and eight corners in three. The particle stays where it is; these are weights for computing a field.

Let ff be the fractional distance between two nodes. The left and right weights are 1−f1-f and ff. A unit mass one-quarter of the way across deposits 0.750.75 and 0.250.25 mass units.

w0=1−f,w1=f,w0+w1=1w_0=1-f,\qquad w_1=f,\qquad w_0+w_1=1
The weights are nonnegative and sum to one. Three-dimensional corner weights multiply the three directional factors, so their sum is still one.
Cloud-in-cell weights for a particle at quarter and four-fifths fractional positions
Original assignment diagram. The white dot is one particle; coloured circles are node masses. Arrows show mass weights for field calculation, not particle motion. At fractional positions 0.250.25 and 0.80.8, the two weights sum to one. Open the full figure

Reading: Springel (2005), The cosmological simulation code GADGET-2, astro-ph/0505010v1 · Original finite-mesh calculations and independent small-wave Python replay

Check mass before removing the mean

With NgN_g nodes per axis and NpN_p equal-mass particles, let wgpw_{gp} denote the weight deposited by particle pp on node gg. Normalising by a grid volume gives density relative to the mean at the same epoch.

Adding all node masses must recover the original total. Check that first, then form the contrast δ=ρ/ρˉ−1\delta=\rho/\bar\rho-1. Forcing a measured mean back to unity could conceal missing periodic weights.

ρgρˉ=Ng3Np∑pwgp,∑gwgp=1\frac{\rho_g}{\bar\rho}=\frac{N_g^3}{N_p}\sum_p w_{gp},\qquad\sum_gw_{gp}=1
The expansion factors in physical cell volume and mean density cancel in the ratio. Total mass remains fixed.

Reading: Original finite-mesh calculations and independent small-wave Python replay · Springel (2005), The cosmological simulation code GADGET-2, astro-ph/0505010v1

Solve the potential wave by wave

Each periodic Fourier wave has an integer vector n\boldsymbol n. Its dimensionless wavenumber is ky=2πn\boldsymbol k_y=2\pi\boldsymbol n, while the dimensional comoving wavenumber is k=2πn/L\boldsymbol k=2\pi\boldsymbol n/L. The constant, zero mode belongs to the mean background that sets expansion.

The scaled potential obeys ∇y2U=3δ/2\nabla_y^2U=3\delta/2. A second spatial derivative of a Fourier wave contributes minus its squared wavenumber, letting us solve each mode separately. CIC deposition smooths the density, and interpolating the force applies another smoothing. The declared operator explicitly corrects those two average windows.

Define sinc⁡x=sin⁡x/x\operatorname{sinc}x=\sin x/x, with continuous value one at zero. The window and retained modes appear below. FFT is a fast way to perform the discrete transform; the mesh and cutoff still determine its physical resolution.

Wn=∏j=13sinc⁡2 ⁣(πnjNg),Un=−3δn2∣ky∣2Wn2W_{\boldsymbol n}=\prod_{j=1}^3\operatorname{sinc}^2\!\left(\frac{\pi n_j}{N_g}\right),\qquad U_{\boldsymbol n}=-\frac{3\delta_{\boldsymbol n}}{2|\boldsymbol k_y|^2W_{\boldsymbol n}^2}
Use this for nonzero modes satisfying ∣n∣≤Ng/3|\boldsymbol n|\leq N_g/3. Other modes and the zero potential mode vanish. The low-pass force retains finite-grid errors.

Reading: Springel (2005), The cosmological simulation code GADGET-2, astro-ph/0505010v1 · Original finite-mesh calculations and independent small-wave Python replay

Bring the potential gradient back to particles

After the inverse transform, force points down the potential gradient. We use a four-point central derivative along each axis, with dimensionless grid spacing h=1/Ngh=1/N_g, then interpolate with the same CIC weights used for deposition.

For Ng=16N_g=16 and L=20 MpcL=20\,\mathrm{Mpc}, nodes are 1.25 Mpc1.25\,\mathrm{Mpc} apart in comoving units, and the Fourier cutoff further limits short waves. Particles can crowd much more closely than this even though the force remains coarse. Internal details of dense knots need a separate resolution study.

(∂jU)i≃8(Ui+1−Ui−1)−(Ui+2−Ui−2)12h(\partial_jU)_i\simeq\frac{8(U_{i+1}-U_{i-1})-(U_{i+2}-U_{i-2})}{12h}
The Fourier potential solve and this derivative together define the actual operator. The reduced PM calculation does not add GADGET-2’s short-range tree force.

Reading: Springel (2005), The cosmological simulation code GADGET-2, astro-ph/0505010v1 · Original finite-mesh calculations and independent small-wave Python replay

Kick, move, then kick again

Force changes momentum and momentum advances position. Symmetric kick–drift–kick stepping first updates part of the momentum using the old force, moves the positions, solves for their new force and completes the momentum update.

For a step from a0a_0 to a1a_1, take am=a0a1a_m=\sqrt{a_0a_1}, consistent with logarithmic spacing. The EdS background lets us integrate the expansion factors exactly. Write the drift and kick factors as D\mathcal D and K\mathcal K to distinguish them from a density growth factor.

D(a0,a1)=2(a0−1/2−a1−1/2),K(a0,a1)=2(a1−a0)\mathcal D(a_0,a_1)=2(a_0^{-1/2}-a_1^{-1/2}),\qquad\mathcal K(a_0,a_1)=2(\sqrt{a_1}-\sqrt{a_0})
These are the dimensionless integrals appearing in the position and momentum equations.
uh=u0+K(a0,am)F(y0),y1=y0+D(a0,a1)uh\boldsymbol u_h=\boldsymbol u_0+\mathcal K(a_0,a_m)\boldsymbol F(\boldsymbol y_0),\qquad\boldsymbol y_1=\boldsymbol y_0+\mathcal D(a_0,a_1)\boldsymbol u_h
Recompute force at the new position, then finish with u1=uh+K(am,a1)F(y1)\boldsymbol u_1=\boldsymbol u_h+\mathcal K(a_m,a_1)\boldsymbol F(\boldsymbol y_1).

Reading: Springel (2005), The cosmological simulation code GADGET-2, astro-ph/0505010v1 · Original finite-mesh calculations and independent small-wave Python replay

Refine time and space separately

Check the factors on an arithmetic example: a0=0.25a_0=0.25 to a1=0.36a_1=0.36 gives D=2(2−1/0.6)=2/3\mathcal D=2(2-1/0.6)=2/3 and K=2(0.6−0.5)=0.2\mathcal K=2(0.6-0.5)=0.2. Actual trajectories divide their interval into many small steps and update density and force throughout.

At fixed mesh, doubling the step count tests temporal error in the same operator. Keeping the continuous initial condition and refining particles and mesh tests spatial and sampling errors. An independent adaptive integrator checks the former; planar collapse and linear growth also supply analytic or independent-force limits.

Mass must be conserved. Kinetic plus potential energy is generally not constant in an expanding background, and symmetric stepping alone does not prove exact energy conservation. A useful readout must still be checked under the refinement and reference appropriate to it.

Reading: Original finite-mesh calculations and independent small-wave Python replay · Springel (2005), The cosmological simulation code GADGET-2, astro-ph/0505010v1

Exercises and explained answers

  1. How does a unit mass at fractional coordinate f=0.8f=0.8 split between its two nodes?

    Read the explanation

    The left receives 0.20.2 and the right 0.80.8, summing to one. Motion changes these weights continuously.

  2. With 1616 cells per axis and L=20 MpcL=20\,\mathrm{Mpc}, does saving 257257 frames instead of 129129 improve force resolution?

    Read the explanation

    It changes how often motion is recorded. Node spacing remains 1.25 Mpc1.25\,\mathrm{Mpc} and the cutoff is unchanged. Force resolution needs an appropriate mesh and particle refinement.

  3. Why check raw deposited mass before subtracting the mean background?

    Read the explanation

    It reveals lost particles or periodic weights. Removing the mean isolates gravity relative to the background; rescaling the deposit first could hide a mass error behind a correct-looking mean.

Try the prediction

Compare both evolutions at the same epoch

  1. Choose An oblique web and predict how an updated force might change motion near a fold.
  2. Overlay the displacement reference and pause late in the evolution. Compare positions and long-wave power at the same initial condition and expansion factor.
  3. Export the full native record, inspect its saved epochs and explain the difference between denser output times and finer integration steps.

Sources and calculation scope

  1. Springel (2005), The cosmological simulation code GADGET-2, astro-ph/0505010v1

    Focused reading: Section 2.1, equations 1–3; Section 3.2/3.2.1, equations 20–23; Section 4.1, equations 26–31, with original pages 3 and 11 viewed. Supplies cosmological momentum, PM windows and time-integration context. The teaching calculation specifies its own cutoff and finite gradient and adds no short-range tree force.

  2. Original finite-mesh calculations and independent small-wave Python replay

    Original accepted kernel/linear, nonlinear and five-preset F64 calculations support the teaching figures and calibration. This small public ZIP contains only an exact current-state small-wave fixture, an independent NumPy/SciPy trajectory replay and expected output. The fixture was generated through the frozen model and canonical serializer in Node, not downloaded in a browser. Large raw preset archives remain separately retained; no Engine or source paper is included.

FeedbackSomething broken, confusing, or missing on this page?Tell us
What kind of feedback?

Sent with this page’s address, its Lab ID, the site language and your device type. Please leave out personal details. Feedback is kept for up to 12 months.Privacy policy