text:
# Scaling up Tides in Numerical Models of Galaxy- and Halo-Formation

## 1 Introduction

Cold, sub-virial distributions of stars undergo a phase of gravitational focusing during which binding energy is redistributed between them (Lynden-Bell’s ‘violent relaxation’ process, see Binney & Tremaine 1987 \[$``$ BT+87\]). The equilibria established through this process show density profiles which, when averaged over spherical shells, approach a de Vaucouleurs law applicable to elliptical galaxies (van Albada 1982, McGlynn 1984). This has since motivated studies of galaxy and galactic halo formation by some degree of gravitational relaxation (e.g. Henriksen & Widrow 1997, Weinberg 2001). A hands-on approach to this problem, free of geometric constraints, consists in integrating the equations of motion with $`N`$-body numerical codes. A crucial step when applying results from $`N`$-body experiments to actual galaxies and haloes consists in bridging the gap between simulation particle numbers and the actual number of stars (or generally, mass elements) in galaxies, which still differ by five orders of magnitude or more in present-day simulations of collisionless dynamics (Athanassoula 2000). It is therefore essential to establish the scaling of $`N`$-body results with particle number.

The following example in spherical symmetry brings the problem to focus. A star at rest converges to the centre of gravity of a free-falling distribution of mass $`M`$ in an interval of time

$$t_{\mathrm{ff}}=\sqrt{\frac{3\pi }{32G<\rho (a,0)>}}$$
(1)
where $`<\rho (a,0)>=4\pi M/3a^3`$ is the mean density inside the star’s initial radius, $`a`$. (The result holds when stars accrete at the origin, so that shell-crossing is suppressed, see Lynden-Bell 1973.) When the density profile is flat initially, all stars converge to a point in a free-fall time. The time-dependent gravitational potential along a radial orbit is

$$\mathrm{\Phi }(r,t)=\frac{GM}{2R}\left(3\frac{r^2}{R^2}\right);^2\mathrm{\Phi }|_R=\frac{GM}{R^3}$$
(2)
with $`r<R`$, the system radius, and double-differencing with respect to $`r`$ at fixed time yields a measure of the tidal field at $`r`$. We note that the tide is unbound as collapse proceeds and $`R0`$. In general one would not expect a flat density profile on the scales of a galaxy but rather a heterogeneous or clumpy matter distribution. Furthermore fragmentation modes develop on all scales in homogeneous, cold distributions (Aarseth, Lin & Papaloizou 1988, hereafter ALP+88). Bound clumps would survive violent relaxation if their binding energy is high (van Albada 1982, Tsuchiya 1998). However the above argument suggests that remnant structures would be severely affected by the maximal potential depth experienced during collapse (here, when $`R`$ reaches a minimum). What constrains this maximum, and how does it scale with particle number? The survival of bound substructures, the mixing of orbits and energy exchanges between stars all relate to constraints set on (2).

We explore these questions both analytically, through an idealised model of collapsing self-gravitating pressureless fluids, and with a set of numerical $`N`$-body calculations to bring out discreteness effects. We find that while collapse simulations in spherical symmetry reproduces theoretical expectations of vigorous infall (ALP+88), small departures from sphericity lead to much gentler collapse. Specifically, we perform $`N`$-body simulations of collapsing spheroids covering four decades in particle number $`N`$ (up to 16 millions) using an FFT code to integrate the equations of motion. We show empirically that during the collapse of spheroidal distributions, the maximum gravitational energy achieved scales with particle number as $`N^{1/6}`$ and so becomes infinite as $`N\mathrm{}`$. Triaxial initial configurations offer no such scaling with particle number. We find from the pressure-free fluid calculations that the maximum gravitational energy achieved depends sensitively on the initial axis ratios. We survey the parameter space of axis ratios on a mesh of 836 points integrated with high resolution to quantify energy maxima. We then use this to constrain the increase in potential energy of self-gravitating finite-$`N`$ systems. Thus for instance, unless triaxial mass profiles are initially rounder than E1, their linear size on average will contract by no more than a factor $`O(20)`$.

We cover the mathematics in Section 2 before giving details of our simulations in the section that follows. The implications of our study to scenarios of galaxy formation are presented in the last section.

## 2 Analysis

Our overall objective is to find out the maximum potential depth haloes and galaxies may attain during violent relaxation. For constant initial density, equation (1) shows that the first caustic occurs at the unique free-fall time $`t_{\mathrm{ff}}`$. This would involve the whole of the galaxy, as opposed to a subset of stars. It is therefore the profile of choice for our study. We justify this choice in part below and later with numerical modelling (see Section 4).

### 2.1 Results in Spherical Symmetry

We begin by considering in more details the collapse of uniform density spheres. The results of this subsection are those of ALP+88.

The solution for radial infall of uniform spheres takes the parametric form

$$\mathrm{sin}2\eta +2\eta =\left(\frac{8GM}{r(0)^3}\right)^{1/2}\times (tt_o)$$
(3)
with $`r(t)=r(0)\mathrm{cos}^2\eta ,\eta (t=t_o)=0`$ and free-fall is complete when $`\eta =\pi /2`$. If we perturb the density profile so $`\rho :\rho _o+\delta \rho `$, $`\delta \rho >0`$ (or, = 0) for $`r<a`$ (or, $`>a`$), the free-fall time is now a function of position and no singularity forms. From (1) we write the new collapse time

$$t_{\mathrm{ff}}:t_{\mathrm{ff}}\times (1\delta \rho /\rho +O[(\delta \rho /\rho )^2])t_{\mathrm{ff}}\delta t_{\mathrm{ff}}.$$
(4)
Introducing

$$\eta (t_{\mathrm{ff}})=\frac{\pi }{2}ϵ$$

and inserting in (3), we find on truncating the Taylor expansion to second order in $`\delta t_{\mathrm{ff}}`$

$$ϵ^3=\left(\frac{8GM}{r(0)^3}\right)^{1/2}\times \delta t_{\mathrm{ff}}.$$
(5)
Since the motion is pressure-less, spherical mass shells re-expand radially once they have reached the centre. Those that originate outside $`r=a`$ meanwhile continue inwards. This spread in arrival time means that the linear dimension of the sphere reaches a positive minimum. This minimum is known in terms of the original system size and $`\delta t_{\mathrm{ff}}`$ :

$$𝒞^1\frac{W(0)}{W(t_{\mathrm{ff}})}=\frac{r(t_{\mathrm{ff}})}{r(0)}=\mathrm{cos}^2\eta (t_{\mathrm{ff}})=\mathrm{sin}^2ϵϵ^2\left(\delta t_{\mathrm{ff}}\right)^{2/3}.$$
(6)
where $`𝒞`$ is the collapse factor, and $`WGM^2/r`$ is the system’s gravitational energy. Note that the definition of $`𝒞`$ applies equally to non-spherical systems; the first equality in (6) is valid only for spheres.

So far we have not specified the form of $`\delta \rho `$ in the region $`0<r<a`$ of the initial configuration, only that it be positive. If furthermore $`\delta \rho `$ is a non-monotonic function of position, shell-crossing will occur before any shell has reached the centre. This would contribute to smooth out fluctuations by orbit mixing, but would not affect infall of the system as a whole. In the case of a point-mass realisation of a uniform density stellar system, discreteness introduces Poissonian noise so that $`\delta \rho /\rho 1/\sqrt{N}`$, with $`N`$ the total particle number. The ratio (6) now scales with the inverse one-third power of $`N`$, or

$$𝒞N^{1/3}.$$
(7)
The density perturbations $`\delta \rho >0`$ leading to (7) are Jeans-type fragmentation modes of instability: the enhanced gravity pulls in the matter which condenses faster at the origin, as in the classic Jeans condensation of star-formation studies. The scaling relation (7) recovers the solution for a cold Newtonian fluid in the limit $`N\mathrm{}`$, for which $`r(t_{\mathrm{ff}})0`$. Systems with small particle numbers experience relatively larger density fluctuations, which act as seeds for fragmentation modes contributing to halt radial infall. Infall stops once orbit crossing occurs in the centre, as encapsulated by the spread in radial collapse time (6).

The above analysis would not apply to initially cuspy profiles, since in that case shell-crossing takes place immediately at the centre. Analysis with shell-crossing is beyond our scope. However if we view a peaked profile as a perturbed, uniform-density distribution, where a large-amplitude perturbation is necessary to distinguish it from Poissonian noise, the spread in free-fall times, $`\delta t_{\mathrm{ff}}`$, would then be even larger, so presumably (6) is minimised for initially flat profiles.

### 2.2 Non-Spherical Collapse

#### 2.2.1 Small deviations

Pressure-less, self-gravitating oblate or prolate structures collapse first down the shortest axis (Lynden-Bell 1964, Lin, Mestel & Shu 1965). Hence perturbations breaking the spherical symmetry while preserving internal density homogeneity grow in time. Writing

$$𝒓(t)=r(0)\mathrm{cos}^2\eta \widehat{𝒓}+𝝃$$
(8)
where $`\eta (t)`$ is known from (3) and $`𝝃(r,\theta ,\varphi ,t)=`$ vector displacement of the spherical surface with respect to the radial direction $`\widehat{𝒓}`$. Lagrangian linear analysis shows growth rates for these modes (cf. eq.  of ALP+88)

$$\frac{\text{d}^2𝝃}{\text{d}t^2}=\frac{4\pi }{3}\frac{G\rho _o}{(\mathrm{cos}^2\eta )^3}\frac{n1}{2n+1}𝝃$$
(9)
where $`n>1`$ is the principal number of an angular decomposition of the displacement in spherical harmonics ($`n=1`$ corresponds to a homogeneous radial contraction). Thus the right-hand side in the above equation is positive and $`|𝝃|`$ becomes larger in time. Since (9) is linear in $`𝝃`$, angular and radial components of each mode (or, value of $`n`$) grow in time at the same rate. We may solve numerically for $`|𝝃|`$ as a function of time using (3). However an approximate solution is found immediately if we note that for collapse in spherical symmetry the time-averaged square cosine,

$$<\mathrm{cos}^2\eta >=\frac{1}{2}$$

gives a mean ratio $`r(t)/r(0)=1/2`$ averaged over $`t_{\mathrm{ff}}`$. Substituting this in (9) and writing $`|𝝃|=\xi (t)`$ we find on integrating

$$\xi (t)=\xi _o\mathrm{exp}\left(\sqrt{8\frac{4\pi }{3}G\rho _o\frac{n1}{2n+1}}[tt_o]\right)$$
(10)
Thus high-order $`n1`$ modes grow faster with increasing $`tt_o`$.

At the time when $`\xi (t)r(t)`$, the linear size of the system may be compared with (6) in order to determine which type of perturbations develop the fastest for a given particle number. For this purpose we truncate to third order a Taylor expansion with respect to $`tt_o`$ of (10) in the limit $`\eta \pi /2(\mathrm{i}.\mathrm{e}.,tt_ot_{\mathrm{ff}})`$. In the case of Poissonian fluctuations, we find from (6) and (7) that the fragmentation modes in spherical symmetry develop faster than surface modes (8) when

$$N^{1/3}\mathrm{}<\left(\frac{\pi }{2}\sqrt{\frac{n1}{2n+1}}+\frac{\pi ^2}{2}\frac{n1}{2n+1}+1\right)^1\left(\frac{\xi _o}{r[0]}\right)^1$$
(11)
For a particle realisation of a uniform-density sphere, the statistics will be Poissonian. The surface mode should initially rise above the noise level to be effective. We therefore set

$$\frac{\xi _o}{r(0)}=\frac{1}{\sqrt{N}},$$

from which we find a critical particle number,

$$N_c^{1/6}=\frac{\pi }{2}\sqrt{\frac{n1}{2n+1}}+\frac{\pi ^2}{2}\frac{n1}{2n+1}+1,$$
(12)
such that for $`N\mathrm{}>N_c`$ fragmentation modes (‘clumpiness’) halt radial infall before $`\sqrt{N}`$-seed surface modes have developed fully and led to pancaking. If we make $`n1`$ we compute the maximum value possible for $`N_c`$:

$$N_c4.6^69475.$$
(13)
Otherwise said, this $`N_c`$ is the largest possible particle number for which discreteness effects (noise) may significantly distort the flow of a collapsing sphere through a pancaking mode.

The results (12) and (13) apply to initially small amplitude deviations from spherical symmetry. For sufficiently large deviations from sphericity, linear analysis shows that small particle number simulations may yet reproduce the pancaking collapse of cold fluids (Lynden-Bell 1964, Lin, Mestel & Shu 1965). Consider for example an axisymmetric displacement of amplitude $`\xi _or(0)/4`$ mapping a sphere to a spheroid of aspect ratio 3:4. For this case we compute from (11) $`N\mathrm{}<4^3/(4.6)^31`$. Thus for axisymmetric cold distributions, of initial aspect ratio $`<3/4`$, any sensible simulation particle number will adequately reproduce the Lin-Mestel-Shu flow. The time-evolution of axis-ratio of collapsing triaxial systems with $`N=10^5`$ particles performed by Hozumi et al. (1996) shows growth of surface modes (pancaking) in agreement with (11). The initial axis ratios of their systems were $`\xi _o/r(0)0.01`$ and $`0.005`$, or three times the Poisson noise level for this number of particles.

We stress that uncorrelated Poisson noise of a uniform spherical distribution is not sufficient in itself to lead to appreciable flattening during infall, for any particle number. Thus the scaling (7) is well recovered from simulations with as few as $`N10^2`$ particles (see ALP+88, Boily et al. 1999), therefore (12) does not invalidate the interpretation of previous studies of small-$`N`$ collapse simulations in terms of one-dimensional radial motion (e.g. van Albada 1982, Aguilar & Merritt 1990, Cannizzo & Hollister 1992).

#### 2.2.2 Large deviations: ellipsoidal figures

An uniform-density ellipsoid collapses down the minor axis before major-axis collapse is complete, followed by re-expansion when the fluid is also pressure-less. Since all orbits are synchronous in a homogeneous distribution, phase-mixing is minimal. Minor-axis cyclic motion continues while the ellipsoid collapses down the major axis, until it too rebounds while the minor-axis assumes a finite value. For a perfect fluid, such cycles of collapse/expansion may repeat themselves without loss of cohesion. For a fluid made up of stars however, the stars exchange kinetic energy and phase along their orbits, causing damping. In the case where the system is initially oblate spheroidal, axes $`a_1=a_2>a_3`$, Boily et al. (1999) found in $`N`$-body simulations the time-evolution of the aspect ratio after plane crossing to be approximately given by

$$\frac{a_3(t)}{a_1(t)}=\frac{a_3(0)}{a_1(0)}\left(\frac{a_1(t)}{a_1(0)}\right)^{1/3}$$
(14)
i.e. the spheroid becomes rounder as collapse continues and the major axis $`a_1(t)0`$. In practice only a few minor-axis oscillations are detected before phase mixing and violent relaxation lead to equilibrium with little or no streaming pattern. We want to establish for triaxial initial configurations a constraint on the collapse factor $`𝒞`$ in (6) by solving the equations of motion for a pressure-less ellipsoidal fluid integrated over a timescale for complete relaxation suggested by $`N`$-body simulations. The motion of a triaxial uniform ellipsoid, of axes $`a_1>a_2>a_3`$, is governed by a set of harmonic equations (BT+87, Table 2.1)

$$\frac{\text{d}\stackrel{~}{v}_i}{\text{d}\stackrel{~}{t}}=_i\mathrm{\Phi }(\stackrel{\mathbf{~}}{𝒙})=\frac{\stackrel{~}{x}_i}{\stackrel{~}{x}_1^3}A_i(𝒂)$$
(15)
$$\frac{\text{d}\stackrel{~}{x}_i}{\text{d}\stackrel{~}{t}}=\stackrel{~}{v}_i$$
(16)
where $`\stackrel{~}{x}_i(\stackrel{~}{t})=a_i(\stackrel{~}{t})/a_i(0)`$, $`\stackrel{~}{v}(\stackrel{~}{t})=\text{d}\stackrel{~}{x}_i/\text{d}\stackrel{~}{t}`$ are dimensionless functions of the dimensionless time $`\stackrel{~}{t}t/t_{\mathrm{ff}}`$. The coefficients $`A_i(𝒂)`$ are known from potential theory. For instance we have

$$A_1(𝒂)2\frac{a_2(0)a_3(0)}{a_1^2(0)}\frac{F(\theta ,k)E(\theta ,k)}{k^2\mathrm{sin}^3\theta }$$

with similar definitions for $`A_2,A_3`$, and

$$k(\stackrel{~}{t})\left(\frac{a_1^2a_2^2}{a_1^2a_3^2}\right)^{1/2};\theta (\stackrel{~}{t})\mathrm{cos}^1(\frac{a_3}{a_1}).$$
(17)
Here $`F,E`$ are incomplete elliptical integrals (see BT+87 for details). We may identify the most relevant configurations by inspecting the gravitational energy. The self-gravitating potential energy $`W`$ is known for uniform ellipsoids from

$$W(𝒂,t)=\frac{3}{5}\frac{M^2}{a_1(t)}\frac{F(\theta ,k)}{\mathrm{sin}\theta }$$
(18)
with the definitions (17). The energy $`W`$ in (18) diverges when the ellipsoid collapses to a rod (or spindle) which is the case when $`𝒂(a_1,0,0)`$. It remains finite when two axes are non-zero, which includes collapse to a disc. Presumably the force and tidal fields are maximised when ellipsoids develop spindles, together with $`W`$. To constrain the tidal field, it is therefore sufficient to determine when an ellipsoid forms a spindle, or, generally, what maximum value $`W`$ may reach during evolution. We were not able to determine these analytically and have resorted to a numerical integration of the equations of motion. We found it useful to introduce the parameter $`\tau `$ defined as

$$\tau \frac{a_2a_3}{a_1a_3}0.$$
(19)
Thus axisymmetric prolate spheroids ($`a_2=a_3`$) all have $`\tau =0`$, whereas axisymmetric oblate spheroids ($`a_1=a_2`$) have $`\tau =1`$. Triaxial structures assume intermediate values.

## 3 Solutions for Pressure-less Ellipsoids.

### 3.1 Method & Tests

For given axes $`𝒂=(a_1,a_2,a_3)`$ we may integrate (15) and (16) subject to the initial conditions $`\stackrel{~}{x}_i(0)=1,\stackrel{~}{v}_i=0.`$ Note that the equations are singular when any of the axes vanishes. To integrate through such singularities we enforce time-symmetry by reversing the flow: $`\stackrel{~}{v}_i\stackrel{~}{v}_i`$ whenever $`\stackrel{~}{x}_iϵ`$, with $`ϵ`$ a free parameter. We set up a fourth-order Runge-Kutta integrator (Press et al. 1992) and varied $`ϵ`$ from $`10^3`$ down to $`1\times 10^6`$ without appreciable differences in the integrated global quantities such as maximum $`W`$ and time (see Table 1). However details of the fluid configurations converged to good accuracy only when $`ϵ5\times 10^5`$ or less. To quantify the quality of orbit-integration, we computed both axial lengths and velocity components for the representative case where $`(a_1,a_2,a_3)=(1,\frac{1}{2},\frac{1}{3})`$; then initially $`\tau =\frac{1}{4}`$ from (19). We computed a similar quantity $`\tau _v`$ from the velocity components $`(\dot{a}_1,\dot{a}_2,\dot{a}_3)`$ which we evaluated at $`t=t_{max}`$, when $`W/W(0)`$ is maximum.

The results are listed in Table 1 for various values of $`ϵ`$. By comparing the runs of individual components $`(a_1,\dot{a}_1)`$ and those of $`(\tau ,\tau _v)`$ with decreasing $`ϵ`$, we may conclude that both geometric and velocity ellipsoids vary little with $`ϵ`$. This does not hold for individual components, such as major axis, $`a_1`$. This last quantity must reach a minimum a few times larger than $`ϵ`$ to ensure that the dynamics is resolved properly. Drawing from the results in the Table, we set $`ϵ=2\times 10^5`$ or $`1.5\times 10^5`$ in all our integrations as a minimal condition to accurate integration.

In conjunction with the value of $`ϵ`$, the choice of timestep is crucial: we adopted a time-adaptative scheme tailored to the instantaneous free-fall time (1) at each step. This allowed us to resolve in time increases in potential energy by factors up to $`340`$, while keeping errors below the 1% level (though not for axisymmetric systems, see below).

When integrating equations (15) and (16), care must be taken that the indices $`(1,2,3)`$ are circulated between each axis to identify the current major and minor axes properly, and allow the correct evaluation of the force field. All integrations were done in three dimensions but we found it necessary to enforce symmetry in the potential when treating initially axisymmetric configurations in order to prevent large numerical errors. With enforced symmetry, we computed correctly the growth of singular axisymmetric potentials for a collapse factor reaching $`40`$. As a test, we integrated through the singularities formed through major-axis collapse of axisymmetric spheroids with initial axes $`(1,\frac{1}{2},\frac{1}{2})`$ and $`(1,1,\frac{1}{2})`$. The error in binding energy at the end of integration reached $`5\%`$ and $`6\%`$ in each case respectively, which we take as reference later when sampling the space of triaxial initial configurations. Note that these errors were accrued during a single integration through singularity: the slow divergence of $`W\mathrm{}`$ as axes vanish makes agreement with theory intractable, however the important quantity for these cases must be the time $`t`$ when $`W`$ diverges (since the spindle morphology is known), which could be evaluated to high accuracy.

In order to determine the sensitivity of our integration scheme to the symmetry of the system, we integrated a triaxial configuration with axes $`=(1,0.99,\frac{1}{2})`$, i.e. nearly oblate axisymmetric. In this case integration yielded a maximum collapse factor of $`26`$ and total error accrued for the system binding energy of $`0.5\%`$ with the same set-up used for the strictly axisymmetric case discussed before. We concluded from this that the integrator resolves relatively small departures from axisymmetry, of the order of one percent in axial ratios; and that these small departures are sufficient to avoid spindles and large energy errors, while still collapsing by appreciable factors.

We have monitored the three axes $`\stackrel{~}{x}_i`$ as a function of time to determine whether a spindle forms, which, in view of our approximations, is the case when

$$𝒂(a_1,ϵ,ϵ);a_1>ϵ:\mathrm{definition}\mathrm{of}\mathrm{a}\mathrm{spindle}.$$
(20)
If this occurs, then a self-gravitating body may assume an even larger collapse factor $`𝒞`$ defined in (6) since our scheme does not resolve the dynamics below that scale. The identification of configurations of high potential energy is made difficult due to the slow divergence of $`W`$ with vanishing minor axes $`a_2,a_3`$ (spindle), or all three axes simultaneously (spherical case). By contrast, the case when only one axis vanishes (by definition not the major axis) is integrable to high accuracy. We therefore computed the logarithmic averaged length $`l(a_1a_2a_3)^{1/3}`$, which we used together with (20) to ensure that the maximal potential energy computed occurred at the same time as the minimum of $`l`$. As a final precaution, we rejected any integrated solution that accumulated errors in binding energy exceeding 10%, in view of our tests with axisymmetric configurations. When (20) does not occur, the collapse factor reaches a finite maximum: this maximum is then the limit any $`N`$-body collapse calculation for this initial geometry may reach, since discreteness effects will only increase the growth of kinetic energy and slow down collapse.

### 3.2 Results

Anticipating the results from $`N`$-body calculations of Section 5, we integrated (15) and (16) up to $`t=t_{\mathrm{ff}}(\stackrel{~}{t}=1)`$, corresponding to one free-fall time (1) in spherical symmetry. We then explored the parameter space of $`(a_2(0),a_3(0))`$ by sampling the parameter $`\tau `$ uniformly in the interval , for a total of 900 pairs $`(a_2,a_3)`$, while fixing the major axis to an initial value $`a_1(0)=1`$. We evaluated (18) and kept the largest value, $`W_{max}`$, found during integration. The results are displayed on Fig. 1. The bottom panel graphs the surface of maximum energy for all pairs ($`a_2(0),a_3(0)`$). It is striking that large islands exist where Max($`W)O(10)`$, whereas all axisymmetric configurations with $`\tau =1`$ or 0 (i.e., $`a_2=1`$ and $`a_2=a_3`$, respectively) must develop spindles and infinite $`W`$. The vertical axis has been capped to $`W_{max}/W(0)=25`$ for clarity. Larger values were not prohibited in the course of integrating numerically. The presence of fragmented regions with large $`W`$ in the plane $`(a_2,a_3)`$ are indicative of the formation of spindles or very compact configurations, in the sense of our equation (20). We note the presence of islands of as few as a single point were $`W_{max}/W(0)25`$, suggesting a complex topology. Details of the topology of the energy surface $`W`$ are not important to the main argument and will not be pursued further.

Cases where spindles formed from initially triaxial configurations turn out to be exceptional. For these cases, a repeat of the integration with a smaller $`ϵ`$ lead to larger Max($`W`$) at the time the spindle formed, requiring careful step-wise integration over a small time interval. The remainder of the integrated solution, however, was left largely unchanged, giving confidence that singularities were correctly identified and cured. As stated earlier, we rejected all runs that accumulated errors in binding energy through full time-integration exceeding 10%: in total some 62 (i.e. 7.1%) of all initial configurations were rejected on this basis. A graph of $`\delta E/E_o`$ vs Max($`W/W(0)`$) showed no clear trend with increasing Max($`W/W(0)`$) for these 62 cases, as would have been the case if a single integration through singularity at high density was accountable for the bulk of the error budget: instead, large numerical errors develop owing to repeated integration through singularities, which may occur in rapid succession if the initial configuration is significantly non-axisymmetric. Indeed of the 62 cases with significant energy errors, 37 initially had axis ratios $`a_3/a_1=1/10`$ or lower; two more showed near sphericity, with $`a_2/a_1`$ and $`a_3/a_1>0.999`$. The remaining 23, however, showed no peculiarities in their initial values, or in the maximum $`W`$ computed. They were, nevertheless, left out of the analysis.

From Fig. 1 we may quantify the fractional area in ($`a_2,a_3`$) leading to large potential energy and collapse factor. The inset gives the projected isocontours in the parameter space. If we assume fair sampling, then the total area covered by the highest-level contours is estimated easily using a rectangular grid to cover the contoured area. In this way we compute a net fraction of $`30\%`$ of the total area exceeding an increase in $`W/W(0)`$ of 25. Note that practically all triaxial configurations with $`a_2\mathrm{}<0.4`$ and $`a_3\mathrm{}<0.2`$ likely develop large $`W`$. These considerations are quantified more accurately by sorting $`W_{max}/W(0)`$ in increasing order for all pairs $`(a_2,a_3)`$. This gives an integrated distribution of collapse factor $`𝒞`$. Figure 2 plots the integrated fraction of 836 solutions as function of $`W_{max}/W(0)`$. Two-thirds of these solutions reached a collapse factor $`𝒞<22.4`$, and 86% have $`𝒞<40`$. We sought a correlation between the morphological parameter $`\tau `$ and the maximum potential reached. On Fig. 2 we also plot the mean $`\tau `$ computed at $`t=0`$ for each ten-percentile interval, in increasing order of $`W_{max}/W(0)`$. We find a broad trend such that ellipsoidal initial conditions with larger $`\tau `$ tend to collapse to deeper potentials. We may identify for the first of these bins, which has $`<\tau >=0.24`$, the broad low-$`W_{max}`$ valley seen on Fig. 1. The non-monotonic relation of $`<\tau >`$ with $`W_{max}/W(0)`$ may be guessed at if we look at a scatter plot of this quantity versus $`\tau `$ directly, as shown on the top panel of Fig. 3. The points clustered to the bottom left corner of the graph clearly reflect the trend of small $`\tau `$ to yield small $`W_{max}/W(0)`$. The broad trend we compute for $`<\tau >`$ can be guessed from shifting a horizontal ruler vertically up the $`W`$ axis: there is a suggestion of a gap in the data which account for the dip in $`<\tau >`$ in the 50-60% interval (when $`W_{max}/W(0)20`$). Our conclusions concerning the significance of this gap must be moderated by the large deviations about the mean values. A more robust signature of dynamical evolution is a shift of the distribution of $`\tau `$ towards lower values during infall. This may be measured by computing the shape parameter $`\tau `$ at the time when $`W`$ is maximum, and compared to the initial value, $`\tau _0`$, for a given configuration. The result is shown on the bottom panel of Fig. 3. At maximum $`W/W(0)`$, the mean $`\tau =0.29`$ approximately, down from 0.50 for the initially uniform distribution. By the definition (19), this implies a roundening of the two minor axes, in the same fashion as occurs for axisymmetric spheroids, that is $`a_3/a_2`$ increases during infall (see Boily et al. 1999 for a discussion of this issue). If the systems were allowed to virialise, as will happen in an $`N`$-body calculation allowing orbit mixing and violent relaxation, in contrast to the fluid model presented here, the effect measured during infall of the increasing ratio $`a_3/a_2`$, so enhancing axisymmetry, would be offset by the onset of radial orbit instability, which is known to enhance triaxiality (see e.g. Aguilar & Merritt 1990). We have not addressed here the question of which effect would prevail in determining the equilibrium of haloes or galaxies. This issue would require more realistic density profiles than used here and is beyond the scope of the present paper.

### 3.3 Summary

We sum up the results for pressure-less self-gravitating ellipsoidal collapse as follows: 1) The bulk of initially triaxial figures does not reach a collapse factor $`𝒞=W/W(0)`$ exceeding 40. We find that 86% of the 836 configurations examined collapsed to smaller values; 2) Triaxial figures that collapse by small factors (small increase of $`W/W(0)`$) tend to have preferentially small values of the parameter $`\tau `$ (i.e., are more axisymmetric, see Fig. 2); 3) the axis ratio $`a_3`$ to $`a_2`$ increases during infall, leading to more axisymmetric configuration at maximum collapse (Fig. 3); 4) since no orbit mixing or fragmentation mode develop in the fluid solutions, which would contribute to boost velocity dispersion and stop infall , we deduce that the collapse factors $`𝒞`$ obtained are absolute maxima.

In order to apply correctly these results to galaxies and haloes, we must first quantify the impact of discreteness effects of finite-$`N`$ systems as discussed in Section 2. We proceed empirically with $`N`$-body numerical calculations. We turn first to the task of reproducing the theoretical expectation (7) for spherically symmetric systems. This is followed by a set of calculations of non-spherical initial configurations, from which we seek trends with particle number to compare with the results for perfect fluids of Fig. 1.

## 4 $`N`$-body calculations: Set-Up and Tests

### 4.1 Numerical code and units

Our intention is to cover as wide a range in particle number as possible to seek out correlations applicable to galaxies and haloes. The nested-grid code superbox was used (Fellhauer et al. 2000). This is an FFT Poisson solver with Cartesian grids, of uniform resolution for each cell of a given grid. There are three grid levels, each within one another, thus enabling enhanced resolution where it is required. This is a crucial feature for the problem at hand, since the structures collapse by large factors and the density increases accordingly. For computational purposes the units were chosen such that the total mass and initial radius of the system = 1, however the gravitational constant $`G=2`$. The free-fall time (1) for uniform spheres is therefore

$$t_{\mathrm{ff}}=\frac{\pi }{4}0.785398\mathrm{}.$$
(21)
A version of the code was set up where the integration timestep of a leap-frog scheme scales down with the inverse square root of the local density maximum during the simulation, in order that the timestep $`\delta t`$ remains in the same ratio to the instantaneous dynamical time (or, $`t_{\mathrm{ff}}`$ evaluated from  but with $`<\rho >`$ now the density at time $`t`$). The three levels of resolution allowed by the code were set such that when the system reaches a minimum size, the number of particles per cell is of the order of a few, hence interactions between particles are resolved approximately on the smallest scales. We found in practice that grid resolution smoothes out forces between particles and hence leads to a less deep potential minimum at the centre of gravity. As we show below, the net effect of gridding can be easily brought to low error levels.

### 4.2 Tests in spherical symmetry

We carried out checks of our set-up in spherical symmetry before conducting a survey of particle number and geometry. All simulations of violent relaxation proceed from zero velocity initially. Our results are summarised in Table 2. The analysis of Section 2.1 suggests a clear relation (Kepler’s) between system radius and the time interval $`\delta t_{\mathrm{ff}}`$ before collapse. From (5) and (6),

$$r(t_{\mathrm{ff}}t)r(0)\mathrm{sin}^2ϵ(t_{\mathrm{ff}}t)^{2/3}.$$

We recover this relation for large particle number simulations shown on Fig. 4. This graphs the evolution of four constant mass (Lagrangian) shells for a uniform density sphere of one million particles. The initial radius $`r(0)=1`$. Because the free-fall time is independent of position, the Lagrangian radii must remain in the same relative ratio to one-another; each must converge to the Keplerian regime near full collapse. The time of collapse $`(tt_{\mathrm{ff}}=0)`$ is off-scale on the right-hand side of the logarithmic abscissae on the figure. Two set-ups are illustrated, of low- (left-hand) and high-resolution (right-hand) grids. The linear high-resolution achieved, $`l=`$ grid size / number of cells $`=0.05/64=0.0008`$, is still large when compared with the mean inter-particle distance $`l_{int}`$ expected from (7),

$$\mathrm{Volume}\mathrm{at}\mathrm{bounce}=N(l_{int}/2)^3=$$

$$\mathrm{Initial}\mathrm{volume}÷\mathrm{N}=r(0)^3/N,$$

or $`l_{int}2\times N^{2/3}=2\times 10^4`$. If we count the average number of particles in a (cubic) cell at the bounce, we find $`(l/4l_{int})^3O(1)`$ particles. The mass distribution is therefore well sampled, and as a result both the constant ratios of Lagrangian radii and their match of the Keplerian relation are reasonably well recovered. By comparison, for runs with reduced number of mesh points, from 128 to 32, we find $`100`$ particles in each cell at the bounce. Poor resolution of the mass distribution leads to artificial deviations from the Keplerian tracks (see left-hand panels on Fig. 4). The Lagrangian radii spread out, which results in shell crossing at the centre while the outer shells continue to fall in: this causes artificial orbit mixing and a gentler collapse.

Our computational strategy must therefore ensure that the mass profile is well resolved at all times. Because the one-dimensional spherical collapse provides the strictest numerical test of our numerical set-up, we first recover the scaling (7) for large particle number to refine the code’s grid resolution. We do this for values of $`N`$ ranging from $`10^4`$ to $`16\times 10^6`$. The mass distribution is mapped accurately by Lagrangian radii sorted on concentric shells. We have measured the collapse factor (6) using both the ratio of gravitational radius, $`r_g`$, and two shells enclosing 30% and 60% of the total mass. As can be seen from Table 2, the collapse factor found for given grid resolution and particle number varies considerably depending on the choice of Lagrangian radii, as well as with the ratio of gravitational energies. However the free-fall time is recovered to 1.6% or better. These different values obtained for the collapse factor were used to define error bars on the averaged quantities. The trend with particle number is displayed on Fig. 5(a). Results in spherical symmetry are plotted as circles on the figure: the filled circles represent results from this paper, while open circles are taken from ALP+88. Each point is the average of data listed in Table 2 for given $`N`$, but excluding the low-resolution runs. The scaling relation (7) shown as the solid line on the figure provides a good fit to all data points.

## 5 Results for $`N`$-body spheroidal collapse $`(\tau =0`$ or 1)

Our results for pressure-less fluids show that initially axisymmetric distributions (with $`\tau =0`$ or 1) develop spindles and divergent gravitational energy as they collapse. These configurations may as a result show significant dependencies on particle number in an $`N`$-body realisation of the solution. We therefore explore the case of the collapse of spheroidal distributions first.

### 5.1 Oblate and prolate spheroids

To construct spheroidal distributions, we squeezed or stretched the axes of spheres to achieve the sought geometry. We consider two cases, a prolate $`\tau =0`$ spheroid of initial axes $`=(2,1,1)`$ and an oblate one with axes $`=(2,2,1)`$. We then performed $`N`$-body calculations with $`10^4,10^5`$ and $`10^6`$ particles and compared the outcome with the pressure-less fluid solutions for the same initial configurations.

The results are displayed on Fig. 6. This graphs the gravitational energy $`W`$ as a function of time. The dimensionless axes $`\stackrel{~}{x}_i`$ of the fluid solution are also displayed, where we have indicated the formation of discs or spindles, according to whether a single or two axes vanished at the time indicated. Spindle or disc each gives rise to sharp features in the profile of $`W`$. The numerical solutions for $`\tau =0`$ or 1 show clear dependencies on particle numbers, in the sense that the larger-$`N`$ calculations map the features of the fluid solution more closely. Better agreement with the fluid solutions are to be expected as we increase particle number, since the fluid solutions corresponds to $`N\mathrm{}`$. This is difficult to assess quantitatively for the solutions as a whole. However we may isolate features that support this view. For instance, as $`N`$ is increased from $`10^4`$ to $`10^6`$ particles, the maximum potential energy achieved in both cases displayed increases and corresponds to ‘spindles’ in the fluid solution. For the case of the $`\tau =1`$ (oblate) fluid spheroid, a spindle forms following collapse of the major-axis at $`t0.623`$. The time of maximum potential is $`t=0.646,0.642`$ and 0.638 respectively for the $`N=10^4,10^5`$ and $`10^6`$ runs, when these maxima shifts upward with $`N`$, from 9.57 to 13.2, and 15.7 for the largest-$`N`$ run. The value for the fluid solution $`\mathrm{}`$ formally. Note that the $`N=10^6`$ run is the only one with a rapid re-collapse to a disc singularity (see inset at $`t0.64`$) similar to the fluid solution. A similar comment can be made for the $`\tau =0`$ (prolate) fluid spheroid, where the sequence of singularities is reversed: a spindle forms first, followed by disc and spindle singularities. The two numerical calculations both follow the fluid solution relatively well, with one important difference: at the time of the first spindle, $`t0.370`$, the $`10^4`$ particle run shows an increase in $`W/W(0)`$ much reduced compared with the $`10^6`$ particle run (3.77 to 5.27, or 70% as much); at $`t0.525`$ a second spindle forms but now the two $`N`$-body runs develop very similar extrema in $`W/W(0)`$: 5.42 and 5.61, respectively. This means that following rebound through the firsts spindle, both calculations suffer a comparable degree of orbit mixing, which then smoothes out the second singularity seen in the fluid. As for the $`\tau =1`$ case, the phase of the $`N`$-body curves tunes up to the fluid solution as $`N`$ increases: thus the first singularity at $`t0.370`$ is found at $`t=0.385(N=10^4)`$ and $`t=0.376(N=10^6)`$, which differs with the fluid solution by only 1.6%.

The two cases $`\tau =0`$ and 1 displayed on Fig. 6 both develop spindles which will reach arbitrarily large potential energy as $`N\mathrm{}`$. Both behave in a qualitatively identical way in this respect. We decided therefore to investigate in more details only the relation of the solution for the oblate $`\tau =1`$ case to the number of particles, before exploring triaxial initial configurations.

### 5.2 Series of oblate spheroids

We constructed oblate spheroidal distributions as indicated above. The equator of the spheroids lies in the XY plane. The diameter was kept fixed and only one axis resized to achieve the desired aspect ratio, hence the gravitational energy $`W`$ is magnified with decreasing $`a_3(0)/a_1(0)`$. In the limit $`a_3(0)0`$ we compute a radial free-fall time

$$t_{\mathrm{ff}}(a_3[0]=0)=\sqrt{\frac{4}{3\pi }}\times t_{\mathrm{ff}}(\mathrm{Eq}.[\text{21}])0.5116.$$

The time of maximum contraction would therefore lie between this and the value (21) for spheres.

The results are listed in Table 3. The errors on the collapse factor $`𝒞=W/W(0)`$ are computed from variations about the mean value for fixed number of particles. The aspect ratios initially lie between 1/6 and 9/10, for particle numbers ranging from $`10^4`$ to $`10^7`$. All results are graphed as triangles on Fig. 5(b). (The effect of different initial aspect ratios on $`𝒞`$ is discussed in §5.3.) We have added points obtained from simulations with direct-summation codes by Boily et al. (1999) for small-$`N`$ systems to those of the present study (filled triangles on the figure). For large-$`N`$ calculations ($`N\mathrm{}>10^5`$ and beyond) the black triangles mark a gently increasing trend, well matched with a power-law dependence $`𝒞N^\alpha `$ with $`\alpha 1/6`$. For $`N10^4`$ or smaller, the fit remains good but note the large scatter for points obtained with a direct-summation code (open triangles on the figure). The range of collapse factors measured for the 10,000 particle runs listed in Table 3 allows for a multiplicative factor of 3/2 between maximum and minimum values of $`W/W(0)`$. The results for the direct-summation runs would allow a somewhat larger range, of perhaps 5/2. Whether this is cause for concern is debatable because of the small number of runs in this bin; it may be that particle-particle interactions, better resolved in the direct-summation scheme, increase the scatter somewhat, though not the mean values, which we recover well with the FFT scheme.

The trend with particle number $`N`$ for spheroidal distributions is never well fitted with the scaling $`N^{1/3}`$ of spherical distributions, though the data differ by only small factors for small-$`N`$ systems. The results of Section 3 suggest that any increasing trend of collapse factor $`𝒞`$ with particle number would be a sensitive function of the symmetry of the initial distribution, or of evolution towards axisymmetry in the course of evolution.

### 5.3 Another look at the LMS flow

We investigated the role played by the Lin-Mestel-Shu instability in numerical $`N`$-body calculations of violent relaxation. For homogenous systems, the LMS instability develops as an aspherical system collapses from rest first down its shortest axis. For oblate spheroids, the collapse down the minor axis $`z`$ occurs in a time

$$t(z=0)\sqrt{\frac{a_3(0)}{a_1(0)}}t_{\mathrm{ff}}(\mathrm{Eq}.\text{21}),$$
(22)
where $`t_{\mathrm{ff}}`$ is the free-fall time for spheres. For the initially flattest spheroids in our sample, of $`a_3(0):a_1(0)=1:6`$, we compute $`t=0.41t_{\mathrm{ff}}0.32`$. The data given in Table 3 and Fig. 6 & 8 show that maximum collapse (or, $`W`$) occurs at later times, as collapse down the major axis sets in. This applies to all simulations. Thus self-induced LMS-type of instabilities are not sufficient by themselves to halt collapse since the kinetic energy dispersion $`\sigma `$ grows anisotropically. This anisotropy persists up to the time of maximum $`W`$ and beyond, even for low-$`N`$ particle number calculations, and imprints the virialised equilibrium that follows (Boily et al. 1999).

A question remains which concerns the relative importance of the LMS instability compared to the fragmentation modes of instability that control collapse in spherical symmetry (cf. Section 2.1). Both types of instability will develop during the collapse of a spheroid, however the growth rate of the LMS instability deduced from (22) will be higher for spheroids with initially small aspect ratios. In Section 2.2, we have argued that the growth rate of the LMS instability is always more rapid than the fragmentation instability if the initial aspect ratio is less than about 0.75. We would, therefore, expect a signature of this instability in the form of a stream at a later stage of collapse, namely the bounce. We seek out evidence for this in our sample of runs of Table 3. Below we refer to the stream as an ‘LMS flow’, which should not be confused with the instability described by Lin and co-workers.

Since (22) shows a relation between minor-axis collapse time and initial aspect ratio, we expect a similar relation between maximum collapse factor $`𝒞=W_{max}/W(0)`$ (the bounce) and initial aspect ratio for the cases when an LMS flow drives the dynamics at that time. Crucial to the argument is the relative phase of the velocity components at the bounce. For the chosen spheroidal initial conditions, the motion can be divided in cylindrical $`z`$ and $`R`$ components, and we may align the minor axis with the $`z`$-component of the reference frame. The expectation for LMS flows is that when both $`<v_z>`$ and $`<v_R>`$ are negative inward (in phase) at the bounce, the value of $`W`$ achieved should be larger than when $`<v_z>`$ and $`<v_R>`$ are of opposite signs (out of phase), i.e. one inward, the other outward. Unless fragmentation or other types of instability manage to erase the signature LMS flow, the relative phase of the velocity components at the bounce will be set by the initial system aspect ratio.

We may compare the ensemble of calculations of Table 3, first by normalising individual values of $`W/W(0)`$ for given $`N`$ to the mean value for that series; in this way we remove the scaling $`N^{1/6}`$ between simulations with different particle numbers. We then sort the normalised values by increasing order of the initial aspect ratio, $`a_3(0)/a_1(0)`$. The results are displayed on Fig. 7. The sinusoidal pattern of the data is unmistakable. The data may be fitted with an sinusoid of amplitude $`1/2`$, which is comparable or larger than the intrinsic scatter of the points at given aspect ratio. Thus the LMS flow has equal or more impact on the potential and the system configuration than other factors which predict no dependence with initial aspect ratio, such as the growth of velocity dispersion by internal fragmentation modes (see Section 2). At the bottom of Fig. 7 we have sketched the relative phase of the velocity components obtained for the analytic fluid solution for oblate spheroids. The step-wise histogram indicates in-phase (high step) motion or out-of-phase (low step) motion. The arrows indicate the polarity of the motion. The correspondence with the numerical data is only suggestive: the pattern itself appears somewhat out of phase with the histograms. Since the pressure-less fluid solution from which the histogram was constructed does not suffer from any type of instability, the poor agreement with the data would suggest that instabilities other than the LMS instability are not completely negligible to set the system properties at the bounce.

We note that for $`N=10^4`$ particle runs and initial aspect ratio $`>1/2`$ or so, the data lie near the normalised $`𝒞=1`$ value (cf. Table 3 and Fig. 7, open squares), hence are not caught in the sinusoidal pattern driven by the streaming motion. However our data also indicates that $`10^4`$ particle runs with initial aspect ratios $`<1/2`$ or so match the range of values obtained with larger-$`N`$ runs. Thus fragmentation modes of instability may still be the dominant factor controlling the bounce for simulations with particle numbers $`\mathrm{}<10^4`$ and initially near-spherical distributions, but not so when the initial aspect ratio is sufficiently small.

## 6 Results for $`N`$-body triaxial collapse $`(\tau 0`$ or 1)

We now extend our study to triaxial configurations. We repeated the exercise of Section 5, but this time varying two initial aspect ratios as parameters. The set-ups used to obtain numerical and analytic solutions were as before. The results are illustrated on Fig. 8. We considered three triaxial configurations with $`\tau =1/3,1/2`$ and 2/3, so covering both prolate and oblate structures. For these cases the analytic fluid solution did not develop spindles and hence the potential energy remained finite for the duration of integration. As for the spheroidal calculations, the numerical $`N`$-body calculations come ever closer to the analytic solution with increasing particle number. For example, for $`N=10^6`$ particle runs, the potential energy already comes within 20% of the fluid solution at maxima. Higher particle numbers would only bring modest differences and convergence as $`N\mathrm{}`$ is therefore very slow. The $`10^4`$ particles runs remains approximately $`50\%`$ out of step with the fluid values, and hence the quantities involved with such low-$`N`$ calculations of relaxation processes are to be treated with caution in applications to galaxy or halo formation problems.

## 7 Discussion and conclusion

We have sought to constrain the tidal field developing around galaxies and haloes as they form. To do this we studied the growth of gravitational energy during the violent relaxation of ellipsoidal bodies. We used both an analytic pressure-less gas model and $`N`$-body numerical integration to set absolute limits on calculations of galaxy and halo formation involving $`N`$ point masses.

We found using the pressure-less gas model that close to 9/10 (86%) of all ellipsoidal triaxial configurations increase their gravitational energy by at most a factor 40 (cf. Figs. 1 and 2). We confirmed this with $`N`$-body calculations using up to one million particles.

We studied axisymmetric and spherical uniform distributions. We extended the scaling of collapse factor with particle number for homogeneous spheres, $`𝒞N^{1/3}`$, to $`N=16`$ millions with the code superbox. We noted that axisymmetric spheroidal distributions also show increasing collapse factor with particle number: The run of data points is well fitted by the power-law $`𝒞N^{1/6}`$, much gentler than for spherically symmetric distributions (see Fig. 5\[b\]). The scatter in the data for spheroids as such does not allow to fall back on the $`𝒞N^{1/3}`$ scaling for spheres, even when we try matching runs with different initial aspect ratio (cf. Table 3). We pointed out that the extrema of binding energy and hence tidal forces met by such systems depend sensitively on the formation of prolate structures (spindles) and are much gentler otherwise.

The growth of velocity dispersion during collapse can be attributed both to global fragmentation modes and Lin-Mestel-Shu-type of pancaking. We presented evidence to the effect that in calculations involving more than $`N=10^4`$ particles, the Poissonian seeds of fragmentation modes leads to growth in kinetic energy such that the sum remains smaller than, though not negligible before, the growth of kinetic energy attributable to the surface mode (Lin-Mestel-Shu) of instability. For $`N\mathrm{}<10^4`$ and initial aspect ratio $`\mathrm{}>1/2`$, the data suggests that fragmentation modes play an equally important rôle for up to $`N=10^4`$ particles (see Fig. 7). For $`N>10^5`$, the LMS flow at the bounce sets the system properties both in terms of potential depth and velocity components for the full range of the initial conditions studied here, where aspect ratios where taken in the range from 1/6 to 9/10. For this range of particle number and more, the collapse of systems with different initial morphologies can be distinguished without ambiguity on Fig. 5, which gives confidence on that the gravatitational infall has been properly resolved.

Taken together, these results imply that the formation of axially- or spherically-symmetric haloes and galaxies lead to deeper potentials during the violent relaxation phase and hence to more pronounced tidal fields than for non-symmetric ones. The stronger tidal fields would in turn reduce the rate of survival of sub-condensations or satellites orbiting within them. Consequently, we would expect galactic halo morphology and satellite populations to be correlated, in the sense that galaxies with rounder massive haloes would harbour fewer satellites, while triaxial haloes would harbour many more satellites, all other parameters being equal. This does not take into account the long-term fate of galactic satellites: Tidal forces do not subside after virialisation, and eventually will cause the disruption of all galactic satellites after a period of time (e.g. Ibata et al. 1994, Klessen & Kroupa 1998, Bullock et al. 2001). Thus the above statement refers to the time of formation only, once the system has virialised.

Another direct consequence of our results is that spherically symmetric haloes should be more centrally concentrated than non-spherical ones in virial equilibrium. We may expect this to bear on the kinematics of observed galaxies. Unfortunately the current observational constraints on the halo shapes are not sufficiently precise to allow us to test our prediction. Halo axial ratios have been measured for hardly over a dozen galaxies. What is more worrisome, however, is that the different techniques seem to give systematically different results. Merrifield (2002) summarises nicely the situation for disc galaxies. Their halos appear to be axisymmetric and oblate, with their axes of symmetry coalligned with the disk axes. The most reliable measurements for the minor to major axis ratios come from polar rings, but galaxies having such structures may not be a representative sample of disk galaxies, because they might be the results of recent mergers. Measurements from the flaring of the HI disk give systematically smaller values than those of polar rings measurements. An application to our own Galaxy (Olling and Merrifield 2000) shows that such values may be valid only if the value of the distance from the sun to the Galactic center and the local Galactic rotation speed are smaller than what is currently believed. Certainly some progress is necessary before the measurements attain the precision we need for testing our prediction.

Gravitational lensing offers some hope by constraining the distribution of total (dark + visible) gravitational mass inside a given volume from the symmetry of the lensed image, which will not respect the centre of mass of the system if it is not spherically symmetric. For instance, Maller et al. (2000) have applied such a lensing technique to the spiral B1600+434. The deconvolution procedure however does suggest that the shape of the halo deduced remains dependent on the choice of halo density profile (isothermal or otherwise) and symmetry. It may be that the systematic application of such techniques to sufficiently large samples would reveal a correlation in the sense that we indicated above.

The results obtained for uniform-density distributions should be contrasted with results obtained for non-uniform initial distributions. Theis & Spurzem (1999) investigated the morphological evolution of initially cold Plummer distributions. The collapse factors they obtained are given on Fig. 5(a) and 5(b). These confirm earlier results by ALP+88 of lower values of $`𝒞`$ for non-uniform systems and our own arguments of Section 2. The curve fitting the Theis & Spurzem data is shown on Fig. 5(b) shifted down with respect to the one obtained for uniform spheres. This new curve now intersects with the collapse factors obtained for aspherical distributions and $`N10^3`$ particles. A Plummer model shows an extended envelope of mass density $`\rho r^{5/2}`$. Thus for the same particle number, the collapse factor of a Plummer model is reduced in comparison with a uniform sphere, presumably due to shell crossing taking place near the centre. However, Plummer models with larger particle numbers also collapse by larger factors, and hence other systems with initially steep profiles will, too. In spherical symmetry, the collapse of mass distribution with radial dependence (e.g. $`\rho r^\alpha `$) has bearing on accretion problems, since the mass shells reach the centre a various rates in time. The currently favoured road to galaxy formation would have many clumps converging to the centre of mass. Since the scaling we have obtained for uniform-density profiles may be extended to non-uniform profiles, as shown with the Plummer model, we may hope that the relation of gravitational gradient to initial morphology will also find application to cosmological models of galaxy and galactic halo formation.

To recover the physics of collapsing systems adequately in $`N`$-body calculations requires a sufficiently large number of particles in order to disentangle effects of mass distribution and morphology. The results of Fig. 5(b) suggest a fiducial number $`N100,000`$ particles as a clean demarcation (where a large gap appears between the lower dotted and solid line on the figure) between initially spherical, axisymmetric and triaxial distributions. Furthermore, the effect of varying the initial density profile appears only to shift the zero-point of the curves, and hence does not affect the relation of maximum collapse factor $`𝒞=W_{max}/W(0)`$ to particle number.

## Acknowledgements

Thanks are due to Albert Bosma and Rachel Somerville for comments on a draft version of this paper; and to Albert Bosma for his hindsight on dark matter haloes. We are grateful to an anonymous referee for a detailed and constructive report. CMB was funded by the Sonder Forschungsbereich (SFB) 439 programme in Heidelberg. An E.G.I.D.E. grant awarded to CMB in 2001 by the French Ministère des Affaires Étrangères helped nbring this project to completion.