---
title: "How Mesh Gravity Follows Folding Matter"
description: "Follow mass assignment, Fourier gravity and symmetric time stepping to see which forces a finite mesh retains and which structures it smooths."
canonical_url: "https://sandboxphysics.com/stories/cosmic-web/how-mesh-gravity-follows-the-fold/"
language: "en"
page_type: "story"
alternate_languages:
  zh-CN: "https://sandboxphysics.com/zh/stories/cosmic-web/how-mesh-gravity-follows-the-fold/"
last_modified: "2026-11-12"
---

[From ripples to the cosmic web](https://sandboxphysics.com/stories/cosmic-web/)

From ripples to the cosmic web · L2

# How Mesh Gravity Follows Folding Matter

Sandbox Physics · 2026-11-12

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.

## Before you begin

Read the coordinate and shell-crossing chapters first. Fourier decomposition means expressing a periodic distribution as waves with different wavelengths, directions and phases.

## 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 $f$ be the fractional distance between two nodes. The left and right weights are $1-f$ and $f$. A unit mass one-quarter of the way across deposits $0.75$ and $0.25$ mass units.

$$
w_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](https://sandboxphysics.com/_astro/mesh.B9WW-ZDl.svg)

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.25$ and $0.8$, the two weights sum to one. [Open the full figure](https://sandboxphysics.com/_astro/mesh.B9WW-ZDl.svg)

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 $N_g$ nodes per axis and $N_p$ equal-mass particles, let $w_{gp}$ denote the weight deposited by particle $p$ on node $g$. 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 $\delta=\rho/\bar\rho-1$. Forcing a measured mean back to unity could conceal missing periodic weights.

$$
\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 $\boldsymbol n$. Its dimensionless wavenumber is $\boldsymbol k_y=2\pi\boldsymbol n$, while the dimensional comoving wavenumber is $\boldsymbol k=2\pi\boldsymbol n/L$. The constant, zero mode belongs to the mean background that sets expansion.

The scaled potential obeys $\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 $\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.

$$
W_{\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 $|\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/N_g$, then interpolate with the same CIC weights used for deposition.

For $N_g=16$ and $L=20\,\mathrm{Mpc}$, nodes are $1.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.

$$
(\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 $a_0$ to $a_1$, take $a_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 $\mathcal D$ and $\mathcal K$ to distinguish them from a density growth factor.

$$
\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.

$$
\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 $\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: $a_0=0.25$ to $a_1=0.36$ gives $\mathcal D=2(2-1/0.6)=2/3$ and $\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.8$ split between its two nodes?

   Read the explanation

   The left receives $0.2$ and the right $0.8$, summing to one. Motion changes these weights continuously.

2. With $16$ cells per axis and $L=20\,\mathrm{Mpc}$, does saving $257$ frames instead of $129$ improve force resolution?

   Read the explanation

   It changes how often motion is recorded. Node spacing remains $1.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](https://sandboxphysics.com/labs/cosmic-web-shell-crossing/)

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.

[Previous: Why One Position Can Have Several Velocities](https://sandboxphysics.com/stories/cosmic-web/when-one-position-has-several-velocities/) [Next: What a Power Spectrum Remembers](https://sandboxphysics.com/stories/cosmic-web/what-a-power-spectrum-remembers/)

## Sources and calculation scope

1. [Springel (2005), The cosmological simulation code GADGET-2, astro-ph/0505010v1](https://arxiv.org/abs/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](https://sandboxphysics.com/stories/cosmic-web/materials/cosmic-web-replay-v1.zip)

   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.