What is Gravitational softening length?

Question

I was reading the Aquarius simulation preprint (Springel et al. 2008) The Aquarius Project: the subhalos of galactic halos as a reference for my internship. I came across the term Gravitational softening length in the Initial conditions section and couldn't get a clear explanation when I searched online.

Can someone explain what it means in this context?

2.2 The Aquarius simulation suite

In Table 1, we provide an overview of the basic numerical parameters of our simulations. This includes a symbolic simulation name, the particle mass in the high-resolution region, the gravitational softening length, the total particle numbers in the high- and low-resolution regions, as well as various characteristic masses and radii for the final halos, and the corresponding particle numbers. Our naming convention is such that we use the tags “Aq-A” to “Aq-F” to refer to simulations of the six Aquarius halos. An additional suffix “1” to “5” denotes the resolution level. “Aq-A-1” is our highest resolution calculation with ∼ 1.5 billion halo particles. We have level 2 simulations of all 6 halos, corresponding to 160 to 224 million particles per halo.

We kept the gravitational softening length fixed in comoving coordinates throughout the evolution of all our halos. The dynamics is then governed by a Hamiltonian and the phase-space density of the discretized particle system should be strictly conserved as a function of time (Springel, 2005), modulo the noise introduced by finite force and time integration errors. Timestepping was carried out with a kick-driftkick leap-frog integrator where the timesteps were based on the local gravitational acceleration, together with a conservatively chosen maximum allowed timestep for all particles

Answer 1

Gravitational softening is a numerical trick introduced in particle simulations to avoid too close encounters that would otherwise result in unrealistic motion of particles.

The gravitational force between two particles of masses $m_i$ and $m_j$ (where in simulations, usually $m_i=m_2$), separated by a distance $r_{ij}$, is $$ F_\mathrm{real} = \frac{G m_i m_j}{r_{ij}^2}, $$ where $G$ is the gravitational constant. This expression diverges for $r\rightarrow0$. The particles in a simulation represent a large amount of mass, thousands or even millions of Solar masses, so in reality if two masses come close, other forces (e.g viscosity) will dominate, and masses will disperse. In a simulation, particles may come arbitrarily and unnaturally close, and numerical rounding may blow up, so instead we substitute the force by $$ F_\mathrm{soft} = \frac{G m_i m_j}{r_{ij}^2 + \varepsilon^2}, $$ where $\varepsilon$ is the gravitational softening length. For $r_{ij}\gg\varepsilon$, this approaches the regular force, whereas for $r\ll\varepsilon$, the force is constant.

The value of $\varepsilon$ sets a limit to how small structures you can resolve in your simulation, and choosing a proper value for $\varepsilon$ is an art; a too small value reduces the desired effect, while a too large value reduces the resolution. Often the value of $\varepsilon$ is set to a certain fraction of the mean interparticle distance, e.g. 1/10 or 1/100. The choice will also affect the density profile of structures, so in order to make sure that your result doesn't depend too heavily on your chosen value, you need to run convergence test with various values, varying number of particles, etc.

The paper that @James K links to, Zhang et al. (2019), has a good list of resources if you want to dig further into the effect of the softening length.

As an example, I've plotted the forces for particles of mass $m_{i,j}=10^5\,M_\odot$, and a gravitational softening length of $\varepsilon=100\,\mathrm{pc}$: