Simulating Galaxies in a Cluster Collapse

In this paper, we introduce a new approach to cluster simulation with a large enough dynamic range to resolve galaxies within a cluster and examine galaxy merging in a cosmological context. We simply assume that disk galaxies form instantly in the centre of galactic-mass dark halos early on in the evolution of the dark matter cluster in a cosmological N-body simulation. At an early time, we replace galactic dark halos with equilibrium, N-body galaxy models scaled to the appropriate mass and dimensions and resolved with 10-100 times as many particles. We assume that the first galaxies are disks embedded in dark halos with flat rotation curves similar to the Milky Way and other nearby galaxies (Kuijken & Dubinski 1995)[34]. To guarantee that the chosen galaxies end up in the final cluster we use the following procedure:

1. A dark matter simulation of the cluster is run until the current epoch, z=0.
2. The particles in the final virialized cluster are labelled and identified at the earlier epoch, z=2.
3. From this subset of particles, all galactic dark halos are identified using the friends of friends linking method (Davis et al. 1985)[15] and replaced by randomly-oriented disk galaxy models of the same mass scale placed on the same trajectories as the halos.
4. The simulation of the cluster is then rerun with the higher resolution galaxy population until the present epoch.

A similar method has been used in the study of galaxy harassment in clusters (Moore et al. 1996a)[41]. The idea is to enhance the resolution of the simulation selectively by using galaxy models as in phenomenological studies while retaining the cosmological character of the mass distribution and large-scale kinematics.

The experimental cluster is chosen from a cold dark matter (CDM) simulation of periodic cubic volume with L=32 Mpc on a side, assuming km/s/Mpc and normalized to . With this normalization, the CDM model is a good description of the clustering of galaxies for the scale we are examining, although it is known to lack sufficient power on larger scales. The initial conditions were generated by applying the Zel'dovich approximation to a random realization of a CDM density field generated with a Fourier transform. The simulation is first run at resolution to identify the site of cluster formation. A spherical volume of comoving radius R = 11 Mpc, associated with the virialized cluster, is identified in the initial conditions and then resampled at the full resolution. The tidal boundary of the collapsing cluster is adequately handled using two concentric shells in the radial ranges of 11<R<16 Mpc and 16<R/L<27 Mpc, sampled at and resolution respectively. The simulation therefore contains a total of 4.3 million particles of which about 1 million end up in the virialized cluster halo. All of the simulations are run with a parallel N-body treecode adapted for both periodic and vacuum boundaries (Dubinski 1996)[20].

The cluster has a virial radius of 1.2 Mpc, a mass of M within this radius and a spherically averaged, central line-of-sight velocity dispersion of 550 km/s. This cluster would be classified as a poor cluster or a large group by observers. At z=2, the 100 most massive dark halos associated with the cluster with masses in the range of to M are identified and replaced with N-body galaxy models. The model used for each galaxy (with different scaling) is composed of an exponential disk, a truncated King model bulge and King model dark halo with the potential and orbital distribution derived from a self-consistent distribution function (Model B of Kuijken & Dubinski 1995)[34]. By design, the model has a flat rotation curve out to 10 exponential scale-lengths and declines beyond that distance.

The 20 most massive halos are replaced with ``high'' resolution galaxy models including 50000 disk particles, 10000 bulge particles and 40000 dark halo particles with a softening length of 0.32 kpc for the stars and 0.64 kpc for the dark halo particles. The remaining 80 halos are sampled at ``low'' resolution with one tenth as many particles as above in the same ratios and a softening length twice as large. The remaining cluster dark matter is retained with a softening length of 3.2 kpc.

The models are scaled according to the value of the circular velocity and mass of the halos in the dark matter simulation at about 1/2 the virial radius. The scale-lengths, h, of the disks are determined by the measured mass and velocity () and fall in the range of observed disks. The 100 galaxies roughly follow a Tully-Fisher relation (Tully & Fisher 1977)[58] in their circular velocity vs. mass profiles with (Figure 1).

  


Figure 1: Relation between mass and circular velocity for the initial galaxy population. The relation is similar to the Tully-Fisher relation.

The disk scale-lengths also vary according to observed laws with again cf. with (Freeman 1970)[25] (Figure 2).

  


Figure 2: Relation between mass and exponential scale-length for the initial galaxy population. The relation follows the prediction of Freeman's (1970) law for exponential disks.

The mass function also has a Schechter form with and M (stellar mass) similar to the observed local galaxy luminosity functions but perhaps somewhat steeper (Loveday et al. 1992)[36]. The simulation is run for 10.5 Gyr with a single leapfrog timestep Myr for a total of 4700 steps. This timestep allows the resolution of structures down to about 0.5 kpc.

One problematic feature of the distribution is that the 3 most massive galaxies have unusually large circular velocities for normal disk galaxies, with km/s. These galaxies are probably ellipticals rather than disks at the time they are selected. Their exact morphology probably makes little difference to the final outcome since they quickly merge at the outset. The simulation should eventually be rerun with elliptical galaxies to check for possible discrepancies. However, the remaining galaxies follow a realistic mass distribution having properties similar to the observed high surface brightness disks.