Spherical N-body halo disperses

[wolfram-schmidt]

I observed a discrepancy between Enzo-E and Enzo legacy (enzo-dev) when running simulations of gas disks (isolated galaxies) with spherical N-body halo. The Hernquist-profile halo is initialised with particle data produced with an ICICLE-like Python tool (NFW-profile is also supported). The gas disk doesn't really matter, it is just the setup I happened to use when I noticed this problem.

I produced the images below from the particles in the simulation output using yt. The colour map specifies the particle velocities:

1. In the case of the Enzo-E run, the halo disperses (top: t = 0, bottom: 500 Myr of evolution):
   
   
2. This is completely different from Enzo legacy, where halo remains stable with just a little bit of migration beyond the cutoff radius (top: t = 0, bottom: 500 Myr of evolution):
   
   

By postprocessing the t=0 data dumps, I confirmed that the initial halo has the expected (i.e. prescribed in the IC generator) density profile and the total mass and kinetic energy of the particles match the expected values (consquently, the halo should be nearly in virial equilibrium). So I think a bug in reading and converting the particle data in Enzo-E can be ruled out.

To run the simulations, I used my fork of Enzo-E, where I only made modifications in the routines for the gas disk (hydrostatic equilibrium following Wang et al. 2010) and particle initialization in enzo_EnzoInitialIsolatedGalaxy.cpp. Similarly, there is an Enzo fork that supports this setup.

I can provide parameter and input files for test runs.

[jobordner]

I played around a bit more with this. When velocities are initialized to zero it collapses, so I suspect a scaling/units bug somewhere, such that either gravity is too weak or velocities are too high. Most likely gravity is still the issue since there's more complexity in scaling the total densities for computing the gravitational potential. [Edit: actually it could very well be velocity scaling, considering your tests both run to 500 MYr] Rerunning Enzo-E and Enzo but with velocities initialized to zero might be a good next test, and comparing the collapse rates might provide some qualitative measurement of the error.

[wolfram-schmidt]

Here we go:

The two runs reported above are shown as orange/red crosses. In both cases, the total kinetic energy of the halo is close to the black dashed line. That's the expected kinetic energy if the halo is in virial equilibrium (computed from the analytic value of the gravitational energy for the initial density profile). Interestingly, the kinetic energy doesn't change much in the Enze-E run, although the halo disperses.

If the initial velocities are set to zero, the halo contracts, particles are accelerated and kinetic energy rises until the halos relaxes into equilibrium (not quite reached after 0.5 Gyr). However, the growth of kinetic energy is smaller by several orders of magnitude in the Enzo-E run. So that could indicate a bug in velocity scaling, but , if so, I find it hard to understand why the kinetic energy in the equilibrium Enzo-E run wouldn't change significantly.