Dark matter particles

Dark matter particles are included in the simulation by setting

nyx.do_dm_particles = true

in the inputs file.

When dark matter particles are present, one has the option to run with or without the baryons; to build with no hydro capability set:

NO_HYDRO=TRUE

It is also possible to build with

NO_HYDRO=FALSE

but turn off hydro at run-time by setting

nyx.do_hydro = 0

in the inputs file.

Our default is to use double precision for the mesh data and single precision for the particles; this is set in the GNUMakefile by:

PRECISION = DOUBLE
USE_SINGLE_PRECISION_PARTICLES = TRUE

We do not recommend using single precision for the mesh data, as this is not tested and will potentially degrade the quality of the hydrodynamical solve. To use single precision for the mesh use:

PRECISION = FLOAT

Equations

If we define \({\mathbf x}_i\) and \({\bf u}_i\) as the location and velocity of particle \(i\), respectively, then we wish to solve

\[\begin{split}\begin{aligned} \frac{d {\mathbf x}_i}{d t} &=& \frac{1}{a} {\mathbf u}_i \\ \frac{d (a {\mathbf u}_i) }{d t} &=& {\mathbf g}_i\end{aligned}\end{split}\]

where \({\mathbf g}_i\) is the gravitational force evaluated at the location of particle \(i\), i.e., \({\mathbf g}_i = {\mathbf g}({\mathbf x}_i,t).\)

Particle Time Stepping: Move-Kick-Drift Algorithm

In each time step:

  • Solve for \({\mathbf g}^n\) (only if multilevel, otherwise use \({\mathbf g}^{n+1}\) from previous step)
  • \({\mathbf u}_i^{{n+\frac{1}{2}}} = \frac{1}{a^{{n+\frac{1}{2}}}} ( (a^n {\mathbf u}^n_i) + \frac{{\Delta t}}{2} \; {\mathbf g}^n_i )\)
  • \({\mathbf x}_i^{n+1 } = {\mathbf x}^n_i + \frac{{\Delta t}}{a^{{n+\frac{1}{2}}}} {\mathbf u}_i^{{n+\frac{1}{2}}}\)
  • Solve for \({\mathbf g}^{n+1}\) using \({\mathbf x}_i^{n+1}\)
  • \({\mathbf u}_i^{n+1} = \frac{1}{a^{n+1}} ( (a^{{n+\frac{1}{2}}} {\mathbf u}^{{n+\frac{1}{2}}}_i) + \frac{{\Delta t}}{2} \; {\mathbf g}^{n+1}_i )\)

Note that at the end of the timestep \({\bf x}_i^{n+1}\) is consistent with \({\bf g}^{n+1}\) becasue we have not advanced the positions after computing the new-time gravity. This has the benefit that we perform only one gravity solve per timestep (in a single-level calculation with no hydro) because the particles are only moved once.

Computing g

We solve for the gravitational vector as follows:

  • Assign the mass of the particles onto the grid in the form of density, \(\rho_{DM}\). The mass of each particle is assumed to be uniformly distributed over a cube of side \(\Delta x\), centered at what we call the position of the particle. We distribute the mass of each particle to the cells on the grid in proportion to the volume of the intersection of each cell with the particle’s cube. We then divide these cell values by \(\Delta x^3\) so that the right hand side of the Poisson solve will be in units of density rather than mass. Note that this is the comoving density.
  • Solve \(\nabla^2 \phi = \frac{4 \pi G}{a} \rho_{DM}\). We discretize with the standard 7-point Laplacian (5-point in 2D) and use multigrid with Gauss-Seidel red-black relaxation to solve the equation for \(\phi\) at cell centers.
  • Compute the normal component of \({\bf g} = -\nabla \phi\) at cell faces by differencing the adjacent values of \(\phi,\) e.g. if \({\bf g} = (g_x, g_y, g_z),\) then we define \(g_x\) on cell faces with a normal in the x-direction by computing \(g_{x,i-{\frac{1}{2}},j,k} = -(\phi_{i,j,k} - \phi_{i-1,j,k}) / \Delta x.\)
  • Interpolate each component of \({\bf g}\) from normal cell faces onto each particle position using linear interpolation in the normal direction.

Output Format

The particle positions and velocities are stored in a binary file in each checkpoint directory. This format is designed for being read by the code at restart rather than for diagnostics. We note that the value of \(a\) is also written in each checkpoint directory, in a separate ASCII file called comoving_a, containing only the single value.

The particle positions and velocities will be written in binary files in each plotfile directory. Dark matter particles will be in DM, active galactic nuclei particles will be in AGN, neutrino particles will be in NPC.

In addition, we can also visualize the particle locations as represented on the grid. There are multiple “derived quantities” which represent the particles. Including particle variables in the derived variables will make them be written as plotfile fields on the grid, i.e.:

amr.derive_plot_vars = particle_count particle_mass_density

in the inputs file will generate plotfiles with only two variables. particle_count represents the number of particles in a grid cell; particle_mass_density is the density on the grid resulting from the particles.

The same naming convention follows for particle velocities on the grid: particle_x_velocity, particle_y_velocity, particle_z_velocity

Derived variables with particle_ represent quantities from the Dark Matter Particle Container. Similar variables from the AGN particle container, and the Neutrino Particle Container are named agn_ and neutrino_. Note these are particle fields written to the grid, which are distinct from the density field in the plotfile, which is baryonic density on the grid.

We note that the value of \(a\) is also written in each plotfile directory, in a separate ASCII file called comoving_a, containing only the single value.