text:
# Time-Decreasing Hazard and Increasing Time until the Next Earthquake

\[

## Abstract

The existence of a slowly always-decreasing probability density for the recurrence times of earthquakes implies that the occurrence of an event at a given instant becomes more unlikely as time since the previous event increases. Consequently, the expected waiting time to the next earthquake increases with the elapsed time, that is, the event moves away to the future fast. We have found direct empirical evidence of this counterintuitive behavior in two worldwide catalogs as well as in diverse local catalogs. Furthermore, the phenomenon can be well described by universal scaling functions.

\]

Many probability distributions have been proposed to account for the recurrence time of earthquakes , which is the time interval between successive earthquakes in a certain region. When aftershocks and mainshocks are considered together as a part of essentially one unique process , we have determined that a universal scaling law describes the probability density $`D(\tau )`$ of the recurrence time. In this way, for events of magnitude $`M`$ above a certain threshold $`M_c`$ in a given spatial area (whose limits do not need to depend on the tectonic background), $`D(\tau )`$ scales with the rate of seismic activity $`R`$ in the area as

$$D(\tau )=Rf(R\tau ),$$

where $`R`$ is defined as the mean number of earthquakes (with $`MM_c`$) per unit time and $`f`$ is a universal scaling function. This scaling is in fact the hallmark of the self-similarity of seismicity in the space-time-magnitude domain. Among several possible scaling functions, the best fit is obtained from a (truncated) gamma distribution,

$$f(\theta )=\frac{C}{a\mathrm{\Gamma }(\gamma )}\left(\frac{a}{\theta }\right)^{1\gamma }e^{\theta /a},$$

with $`\gamma `$ the shape parameter, $`a`$ the scale parameter, $`C`$ a correction to normalization (due to the truncation of the distribution close to zero), and $`\mathrm{\Gamma }(\gamma )`$ the gamma function. When $`\gamma <1`$, $`f`$ turns out to be a decreasing power law, accelerated by an exponential factor in the long-time limit.

We analyze many regions and $`M_c`$values from two worldwide catalogs (the NEIC-PDE and the Significant Earthquake Database from NOAA at NEIC ) and from several regional catalogs (Southern California, Japan, New Zealand, New Madrid (USA), the Iberian Peninsula, and the British Islands ). Each analyzed region is delimited by two meridians and two parallels , with linear size spanning from $`0.16^{}`$ (about 18 km) to the whole globe ($`2010^3`$ km), covering a large variety of tectonic environments, whereas the considered magnitudes range from larger than 1.5 to larger than 7.5 (this is about a factor $`10^9`$ in the minimum released energy).

Except for a $`360^{}\times 180^{}`$region, which covers the entire globe, the rest of the regions are defined by a window of $`L`$ degrees in longitude and $`L`$ degrees in latitude. The coordinates $`(x,y)`$ of the west-south corner of these regions can be obtained from the vector $`(k_x,k_y)`$ at Fig. 1’s labels as $`x=x_{min}+k_xL`$, $`y=y_{min}+k_yL`$, with $`(x_{min},y_{min})=(180^{},90^{})`$, $`(123^{},30^{})`$, $`(127^{},27^{})`$, $`(160^{},60^{})`$, $`(91^{},35^{})`$, $`(20^{},30^{})`$, and $`(10^{},45^{})`$ for the worldwide, Southern California, Japan, New Zealand, New Madrid, Iberian Peninsula, and British Island catalogs, respectively. The periods of study are (in years A.D. including the last one) 1973-2002 for the NEIC catalog, 1897-1970 for the NOAA one, 1988-1991, 1995-1998, and 1984-2001 for Southern California (denoted as SC88, SC95, and SC84), and 1995-1998, 1996-2001, 1975-2002, 1993-1997, and 1991-2001, for the rest of catalogs (in the same order as before). The regions and times of observation are selected in order that a period of stationary seismic activity is included, this means that aftershock sequences should not have too much weight in the seismicity of the region. When this is not the case (i.e., for very large aftershock activity) our analysis is still valid, but the scaling with the mean rate has to be replaced by a scaling with the instantaneous rate.

A maximum-likelihood estimation of the parameters using the rescaled recurrence time $`R\tau `$ for all the regions and $`M_c`$’s studied gives $`\gamma =0.74\pm 0.05`$, and $`a=1.23\pm 0.15`$, which yields a coefficient of variation $`CV1.2`$. The constant $`C`$ is determined from the normalization condition given the minimum value for which the gamma distribution holds; for $`\theta >0.05`$, $`C=1.10\pm 0.10`$ (see Fig. 1 (a)).

The results of the fit are shown in Fig. 1 (a) using the survivor function, which is defined as $`S(\tau )\text{Prob}[\tau ^{}>\tau ]=_\tau ^{\mathrm{}}D(\tau ^{})𝑑\tau ^{}`$ (where $`\tau ^{}`$ is a generic label for the recurrence time, while $`\tau `$ refers to a particular value of the same quantity). It is straightforward to obtain that, in our case, $`S(\tau )`$ should also verify a scaling relation, $`S(\tau )=G(R\tau )`$, with $`G(\theta )=CQ_\gamma (\theta /a)`$, and $`Q_\gamma (\theta /a)`$ the complement of the incomplete gamma function . The total agreement between these equations and the measured distributions is clear from the data collapse and the fitting curve in the plot, for intermediate and long recurrence times. The accuracy of the scaling law and the gamma fit is guaranteed as the seismic activity is stationary in this range of recurrence times. On the contrary, short times are usually not free of disturbances of the stationariness, due to the triggering of aftershock sequences, which destroy the universal scaling behavior. In order to treat all the distributions in the same way, we calculate the rate $`R`$ only for events in the scaling region, i.e., short recurrence times are not considered in the rate.

The knowledge of the probability distribution of the recurrence times allows one to answer two important questions about the temporal occurrence of earthquakes. First, for a certain region and for $`MM_c`$, we can study the probability per unit time of an immediate earthquake given that there has been a period $`\tau `$ without activity, using the hazard rate ,

$$\lambda (\tau )\frac{\text{Prob}[\tau <\tau ^{}\tau +d\tau |\tau ^{}>\tau ]}{d\tau }=\frac{D(\tau )}{S(\tau )},$$

where the symbol “$`|`$” accounts for conditional probability. From the previous formulas we get that $`\lambda (\tau )`$ scales as $`\lambda (\tau )=Rh(R\tau )`$, with

$$h(\theta )=\frac{1}{a\mathrm{\Gamma }(\gamma )}\left(\frac{a}{\theta }\right)^{1\gamma }\frac{e^{\theta /a}}{Q_\gamma (\theta /a)}.$$

This function turns out to be monotonically decreasing, tending as a power law to the value $`1/a`$ as $`\theta \mathrm{}`$. So, contrary to common belief, the hazard does not increase with the elapsed time since the last earthquake, but just the opposite; this is precisely the more direct characterization of long-term clustering .

Also, one can wonder about the expected time till the next earthquake, given that a period $`\tau _0`$ without earthquakes (in the spatial area and range of magnitudes considered) has elapsed,

$$ϵ(\tau _0)E[\tau \tau _0|\tau >\tau _0]=\frac{1}{S(\tau _0)}_{\tau _0}^{\mathrm{}}(\tau \tau _0)D(\tau )𝑑\tau .$$

This function can be referred to as the expected residual recurrence time  and again we find a scaling form for it, which is $`ϵ(\tau _0)=e(R\tau _0)/R`$, with the scaling function

$$e(\theta )=a\gamma \frac{Q_{\gamma +1}(\theta /a)}{Q_\gamma (\theta /a)}\theta .$$

This is an increasing function of $`\theta `$, which reaches an asymptotic value $`e(\theta )a`$. Therefore, the residual time until the next earthquake should grow with the elapsed time since the last one. Notice the counterintuitive behavior that this represents: if we decompose the recurrence time $`\tau `$ as $`\tau =\tau _0+\tau _f`$, with $`\tau _f`$ the residual time to the next earthquake, the increase of $`\tau _0`$ implies the increase of the mean value of $`\tau _f`$, but the mean value of $`\tau `$ is kept fixed. In fact, this is just a more dramatic version of the classical waiting-time paradox .

For the particular case of earthquakes this is even more paradoxical, since one would say that the longer the time one has been waiting for an earthquake, the closer it will be, due to the fact that as time passes stress increases on the corresponding fault and the next earthquake becomes more likely. (Nevertheless, one needs to distinguish between earthquakes on a given fault and earthquakes over a certain area.) The question was originally put forward by Davis et al. , who pointed out that if a lognormal distribution is a priori assumed for the recurrence times, the expected residual time increases with the waiting time. (However, the increase here was associated to the update of the distribution parameters as the time since the last earthquake, which was taken into account in the estimation, increased; we deal with a different concept of increasing residual time.) Sornette and Knopoff  showed that the increase (or decrease) depends completely on the election of the distribution, and studied the properties of a number of them. We are going to see that the observational data provide direct evidence against the simple picture of the next earthquake approaching in time.

Indeed, our mathematical predictions can be contrasted with the catalogs; both the hazard rate and the expected residual recurrence time can be directly measured with no assumption about their functional form. Following their definitions, these quantities are estimated as

$$\lambda (\tau )=\frac{n(\tau ,\tau +\mathrm{\Delta }\tau )}{n(\tau ,\mathrm{})\mathrm{\Delta }\tau },ϵ(\tau _0)=\frac{_{i|\tau _i>\tau _0}(\tau _i\tau _0)}{n(\tau _0,\mathrm{})},$$

where $`n(\tau _1,\tau _2)`$ denotes the number of quakes with recurrence time in the range $`(\tau _1,\tau _2)`$ and the sum in $`ϵ(\tau _0)`$ is computed only for earthquakes $`i`$ such that $`\tau _i>\tau _0`$ (and of course $`MM_c`$). From the results displayed in Figs. 1 (b) and 1 (c) it is apparent that the hazard rate decreases with time whereas the expected residual recurrence time increases, as we predicted. Moreover, both quantities are well approximated by the proposed universal scaling functions.

On the other hand, the part of the recurrence-time distribution that accounts for short times displays a typical behavior $`f(\theta )=K_1/\theta ^{1\alpha }`$, and the corresponding functions for the survivor function, hazard rate, and expected residual return time turn out to be:

$$G(\theta )=K_1(K_2\theta ^\alpha /\alpha ),h(\theta )=\frac{1}{K_2\theta ^{1\alpha }\theta /\alpha },$$

$$e(\theta )=\frac{K_3\theta ^{1+\alpha }/(1+\alpha )}{K_2\theta ^\alpha /\alpha }\theta ,$$

where the constants $`K_2`$ and $`K_3`$ depend on the rest of the distribution. An example for these functions with $`\alpha 0.2`$ is also represented in Fig. 1, showing the appropriate decreasing or increasing tendency in each case.

For the sake of concreteness, let us consider worldwide earthquakes with $`M7.5`$, which occur at a rate $`R=6`$ per year, roughly. In the days immediately after one event of this type, the expected time to the next one (anywhere in the world) is about 2 months (for $`\tau =6`$ days, we have $`R\tau =0.1`$, and $`e(0.1)1`$, see Fig. 1 (c)). If after 2 months the quake has not come, the expected residual time not only does not decrease but increases to 2.2 months ($`e(1)1.1`$, this would lead to $`4.2`$ months between both events), and if the elapsed time reached 1 year (which is unlikely but not impossible), the expected waiting time would further increase to 2.4 months ($`e(6)1.2`$). In the same way, the hazard rate would drop from 0.7 to 0.5 and to 0.4 month<sup>-1</sup> ($`h(\theta )1.4,1`$, and $`0.85`$), respectively. The same process is reproduced at all magnitude and spatial scales in a self-similar manner. An intuitive explanation of this phenomenon is that when the elapsed time since the last earthquake is large, the system enters into a long “drought period” in which the recurrence time is likely to be very large. Notice however that there is no fundamental difference between these drought periods and the rest of recurrence times, since all of them are governed by the same smooth distribution.

The universality of this behavior demands further explanation; nevertheless, it suggests the existence of a simple mechanism in which, as time passes, the variable that triggers rupture runs away from the rupture threshold (on average). The “excursions” of this variable would keep the memory of the last event stored in the system up to very long times to generate the negative aging observed.

The considerations reported here should be at the core of any research regarding earthquake-occurrence modeling  and predictability . Finally, in order to account for the self-similarity of these processes, the concept of self-organized criticality provides the most appealing framework up to now .

This author has benefited a lot from the original perspectives and deep insights of the late Per Bak. The author also thanks M. Boguñá, D. Sornette’s criticisms, the Ramón y Cajal program of the Spanish MCyT, and all the people at the Statistical Physics Group of the Universitat Autònoma de Barcelona, as well as those institutions that have made their catalogs available on the Internet.


format: text

text:
# Effects of Small-Scale Fluctuations of Neutrino Flux in Supernova Explosions

## 1 Introduction

It has been turned out that most simulations of core-collapse supernova explosions with spherical symmetry fail to produce a successful explosionLie01 . In addition, observations suggest that the ejecta of core-collapse supernova explosions are aspherical(e.g., Wa02 ). These facts lead us to multidimensional simulations. In 2-D and 3-D simulationsMiWiMa93 ; HeBeHiFrCo94 ; BuHaFr95 ; JaMue96 ; Me98 ; FrHe00 ; ShEbSaYa01 ; FrWa02 ; KiPlJaMu03 , it has been shown that multidimensional effects, such as convection inside the proto-neutron star and convective overturn around the neutrino-heated region, increase the explosion energy and can trigger a successful explosionHeBeHiFrCo94 ; JaMue96 ; KeJaMue96 .

When a proto-neutron star rotates, the neutrino flux is expected to be enhanced along the rotational (polar) axis. Janka and MönchmeyerJaMoe89a ; JaMoe89b  first discussed the possibility of aspherical neutrino emission from a rapidly rotating inner core. They argued that a neutrino flux along the polar axis might become three times greater than that on the equatorial plane.

Shimizu et al.ShYaSa94 ; ShEbSaYa01  proposed that the anisotropic neutrino radiation should play a crucial role in the explosion mechanism itself. They carefully investigated the effects of anisotropic neutrino radiation on the explosion energy. They found that only a few percent enhancement in the neutrino emission along the pole is sufficient to increase the explosion energy by a large factor, and leads to a successful explosion. They also found that this effect saturates around a certain degree of anisotropy. It should be noted here that the assumed rotational velocity of the inner core is very different between Janka et al.JaMoe89a ; JaMoe89b  and Shimizu et al.ShYaSa94 ; ShEbSaYa01 .

In the work of Shimizu et al.ShEbSaYa01 , they considered only a form of global anisotropy; the maximum peak in the neutrino flux distribution was located at the pole and the minimum at the equatorial plane. On the other hand, Burrows et al.BuHaFr95  have suggested that the neutrino flux can fluctuate with angle and time due to gravitational oscillation on the surface of the proto-neutron star. In this work, we introduce such small-scale fluctuations in the neutrino flux in our numerical code by modifying the angular distribution of the neutrino flux. We aim to study the effects of these small-scale fluctuations on the shock position, the explosion energy, and the asymmetric explosion. Details are found in ref.MaShMo03 .

## 2 Numerical Code

We perform 2-D simulations by solving hydrodynamic equations in spherical coordinates. A generalized Roe’s method is employed to solve the hydrodynamic equations with general equations of state (EOSs). The details of our numerical technique, together with the EOS and the initial condition used, are described in the previous articleShEbSaYa01 . In our study, we have improved the numerical code of Shimizu et al.ShEbSaYa01 ; the cells in the $`\theta `$-direction were shifted by half of the cell sizeSh95  in order to avoid a numerical error near the pole, although the error was not serious for the investigation of the explosion energy.

In the present work, the local neutrino flux is assumed to beMaShMo03

$$l_\nu (r,\theta )=\frac{7}{16}\sigma T_\nu ^4c_1\left(1+c_2\mathrm{cos}^2(n_\theta \theta )\right)\frac{1}{r^2},$$
(1)
where $`\sigma `$ is the Boltzmann constant, and $`T_\nu `$ is the temperature on the neutrinosphere. In equation (1), $`n_\theta `$ represents the number of waves in the $`\theta `$-direction. The case of $`n_\theta =1`$ corresponds to the global anisotropy, namely, no fluctuation. $`c_2`$ is a parameter which is related to the degree of anisotropy in the neutrino radiation. We see in equation (1) that the neutrino fluxes in the $`x`$ (equatorial) and $`z`$ (polar) directions become $`l_xl_\nu (r,\theta =90^{})c_1`$ and $`l_zl_\nu (r,\theta =0^{})c_1(1+c_2),`$ respectively. The degree of anisotropy $`l_z/l_x`$ is then represented as

$$\frac{l_z}{l_x}=1+c_2.$$
(2)
The value of $`c_1`$ is calculated from $`c_2`$ and $`n_\theta `$ so as to adjust the total neutrino flux to that in the spherical model at the same $`T_\nu `$.

It should be noted here that the amplitude of fluctuations in the neutrino flux distribution for an observer far from the neutrinosphere and that on the neutrino-emitting surface are different. When we observe the neutrino flux far from the neutrinosphere, the local neutrino flux is seen as equation (1). On the other hand, the neutrino flux on the neutrino-emitting surface has a profile similar to equation (1) but the amplitude is different. In the latter, $`c_2`$ is replaced by $`a`$ where $`a`$ is a parameter which represents the degree of anisotropy of the neutrino flux on the neutrinosphere. It is preferable that we compare the results for the same value of $`a`$, since $`a`$ is more directly related to explosion dynamics. The value of $`c_2`$, therefore, is calculated from a given $`a`$, depending on $`n_\theta `$. Although it is difficult to calculate the exact relationship between $`c_2`$ and $`a`$, we can estimate it by assuming that the strength of the neutrino flux on the neutrinosphere is approximated by a profile of step function. This makes it possible to relate the value of $`c_2`$ to the value of $`a`$ for each $`n_\theta `$. For detail, see Madokoro et al.MaShMo03 .

We set $`a=0.31`$ in this work. This value of $`a`$ is chosen in such a way that the value of $`l_z/l_x`$ for the global model ($`n_\theta =1`$) becomes 1.10. The values of $`c_2`$ for each fluctuation model ($`n_\theta =3,5`$) are accordingly calculated. These are summarized in Table 1. The neutrino temperature on the neutrino-emitting surface $`T_\nu `$ is assumed to be 4.70 MeV. In our simulation, we have 500 nonuniform radial zones which cover from 50 to 10000 km in radius. For $`\theta `$, we have 62 uniform angular zones from $`\theta =0^{}`$ to $`\theta =90^{}`$ with equatorial symmetry for $`n_\theta =1`$ and $`3`$. Note that we use a double value of the angular zones for $`n_\theta =5`$, which is different from that used in ref.MaShMo03 .

## 3 Results

Figure 1 depicts the contour maps of the dimensionless entropyShEbSaYa01  distribution with the velocity fields for the three models at $`t250`$ ms after the shock stall. The shock front is represented by the crowded contour lines at $`r20003000`$ km for the model A1-T470, $`r16002200`$km for the model A3-T470, and $`r1400`$km for the model A5-T470. We see in Figure 1 that the shock front is largely distorted in a prolate form for the model A1-T470. This is because the neutrino heating along the pole is more intensive than that on the equatorial plane. Due to increased pressure in the locally heated matter near the polar axis, the shock front along the pole is pushed up, resulting in a prolate deformation. We find that the shock position becomes less extended than that of the global anisotropy, and the shock front approaches the spherical shape when the mode number of fluctuation increases. This trend is especially remarkable for the model A5-T470 in which the shock front is almost spherical.

Figure 2 shows the evolution of the explosion energy, as well as the thermal, kinetic, and gravitational energies for the three models. The energy gain for the case of $`n_\theta =1`$ is the highest among others at all stages of the explosion. It is also seen in Figure 2 that the explosion energy decreases as the mode number of fluctuations in the neutrino flux increases and finally approaches that of the spherical explosion.

Thus, we found that there are remarkable differences in the explosion energy depending on the mode number of the fluctuations. We also found that larger number of modes in the fluctuations makes the result closer to that of the spherical explosion. This is because any small-scale fluctuations on the neutrinosphere are greatly averaged out when the neutrino emission is observed far enough from the neutrino-emitting surface. Moreover, we found that a certain broad space is needed to be heated by neutrinos to revive the stalled shock wave rigorously and hence the global anisotropy ($`n_\theta =1`$) is the most effective to increase the explosion energy. These results support the claim by Shimizu et al.ShEbSaYa01 .

## 4 Conclusion

We have investigated the effects of small-scale fluctuations in the neutrino flux on the core-collapse supernova explosion. The profile of the neutrino radiation field was specified taking its geometric effects into account. Since the small-scale fluctuations are averaged out for radiative and hydrodynamic reasons, the results of the fluctuation models become closer to that of the spherical explosion. Consequently, the global anisotropy is the most effective mechanism in increasing the explosion energy when $`L_\nu `$ is given. This supports the claim made by Shimizu et al.ShEbSaYa01 .

## Index


text:
# Return to return point memory

## Abstract

We describe a new class of systems exhibiting return point memory (RPM) that are different from those discussed before in the context of ferromagnets. We show numerically that one dimensional random Ising antiferromagnets have exact RPM, when configurations evolve from a large field. However, RPM is violated when started from some stable configurations at finite field, unlike in the ferromagnetic case. This implies that the standard approach to understanding ferromagnetic RPM systems will fail for this case. We also demonstrate RPM with a set of variables that keep track of spin flips at each site. Conventional RPM for the spin configuration is a projection of this result, suggesting that spin flip variables might be a more fundamental representation of the dynamics. We also present a mapping that embeds the antiferromagnetic chain in a two dimensional ferromagnetic model, and prove RPM for spin-exchange dynamics in the interior of the chain with this mapping.

Many systems exhibit the remarkable phenomenon of “Return Point Memory” (RPM) RPMexpt . This is most easily demonstrated by looking at hysteresis loops for ferromagnets, Figure 1, where an external field $`H`$ is lowered from saturation to a minimum value $`H_{min}`$ and then raised by some intermediate amount to $`H_{max}`$ before being lowered again to $`H_{min}`$. In systems with RPM, the state of the system the second time the field reaches $`H_{min}`$ is identical to the first. The result generalizes to more complicated variations in $`H,`$ under the constraints shown in Figure 1.

One of the most well known demonstrations of this phenomenon is Barkhausen noise Barkhaus  where the noise observed in changing $`H`$ of a ferromagnet is highly reproducible under repeated cycling of the field Barkexpt . These experiments show that not only does the system return to a state having the the same total magnetization but also gives evidence that the domain wall configurations are identical.

Many systems do not exhibit RPM. For example spin glass Hamiltonians with both ferromagnetic and antiferromagnetic bonds, violate return point memory and instead exhibit subharmonic limit cycles under the application of a periodic field Subharmonic . Therefore from a theoretical point of view, it is of interest to try to find the conditions for a system to exhibit RPM. A major advance in this direction came from a proof of return point memory for systems with purely ferromagnetic interactionsSethna . For a broad class of ferromagnetic models with zero temperature dynamics, RPM was elegantly proven Sethna  by building on the earlier no-passing theorem for charge density waves Middleton .

The essential ingredient in the proof Sethna  of RPM is that a partial ordering of spin configurations is preserved under the application of $`H`$. By partial ordering, we mean that if one has two spin configurations $`\alpha ,\beta `$ corresponding to spins $`s_1^\alpha ,\mathrm{},s_1^\alpha `$ and $`s_1^\beta ,\mathrm{},s_1^\beta `$ then one says that $`\alpha \beta `$ if $`s_i^\alpha s_i^\beta `$ for all $`i`$. Preservation of this ordering means that if a time dependent field $`H_\alpha (t)`$ is applied to $`\alpha `$, and a smaller field $`H_\beta (t)H_\alpha (t)`$ (for all time $`t`$) is applied to $`\beta `$, then $`\alpha \beta `$ for all times. That is, $`\alpha `$ cannot “pass” $`\beta `$. This no-passing constraint was proved to be valid for ferromagnetic systems with single spin flip dynamicsSethna .

An interesting question to ask is whether or not all systems satisfying RPM require similar conditions as were needed to prove the ferromagnetic caseSethna . In more detail, is it necessary for every pair of two states, that satisfy $`\alpha \beta `$, to satisfy the no-passing requirement? Here we show that, rather surprisingly, the answer to this question is no. This is shown by examining fully antiferromagnetic Ising chains in one dimension, with zero temperature (deterministic) dynamics identical to those used for ferromagnetic systems. We find that when started from a large $`H`$ and fully saturated magnetization, the system always satisfies RPM. However, unlike the ferromagnetic case, if one starts in a random state that is stable at some large field $`H_s`$, and lowers the field to $`H_{min}`$, then raising the field to $`H_{max}<H_s`$ and returning it to $`H_{min}`$ changes the state. From a practical perspective, this is not a severe restriction, since if $`H`$ is saturated in the distant past, RPM is valid for any subsequent evolution of $`H(t)`$ foot:temperature . However, for the purpose of proving RPM, this implies that there is no definition of the bounding operator “$``$” (at least none for which $`\alpha \alpha `$ is true of all $`\alpha `$) with which a proof along the lines of the ferromagnetic case can be constructed.

Even for states descended from the saturated state, for which RPM is satisfied, we find that no-passing is violated for the spin configuration: in Figure 1, the spin configuration on branch (3) is not bounded above and below by (1) and (2) respectively. However we have been able to construct a new “spin flip” variable that does satisfy no-passing when starting from a high field. (This is not true starting from a random configuration; the remarks of the previous paragraph apply to any variable and bounding operator.) No-passing for the spin flip variable implies that it also satisfies RPM. (However, as with the spin configuration, RPM for a variable does not imply no-passing.) Since the spin configuration can be obtained as a projection of the flip state, this version of RPM is stronger, suggesting that this new variable may be a more fundamental way of understanding these systems.

We also investigate whether if one moves beyond single spin flip dynamics Glauber  it is possible for the anti-ferromagnetic chain to show RPM in the same sense as the ferromagnetic case. We were able to show that with spin exchange dynamicsKawasaki  that, aside from the ends, conserve magnetization, antiferromagnetic chains do show RPM in the same configuration-independent way that ferromagnets do. We did this by embedding this one dimensional antiferromagnetic problem in a two dimensional ferromagnetic system which has single spin flip dynamics and therefore shows return point memory. However this mapping must fail for single spin dynamics in a rather interesting way; because under this mapping, single spin flips become nonlocal and the standard proofSethna  does not applyfoot:itshouldn't  Therefore it is not just the Hamiltonian, but the dynamics as well, that determine whether or not a system satisfies RPM.

We consider the random antiferromagnetic Ising model:

$$=\underset{i}{}[J_is_is_{i+1}+(h_i+H)s_i]$$
(1)
where the bonds $`J_i`$ and local fields $`h_i`$ are independent random variables. All the $`J_i`$’s are negative, and the $`h_i`$’s are equally likely to be positive and negative. $`H`$ is the externally applied field. Initially, $`H`$ is large and positive, and all the spins point up. Thereafter, the field is changed adiabatically. At any field, a spin is flipped if doing so reduces the energy $``$ of the system. This spin flip can render other spins unstable, in which case the process is repeated till there are no more spins to flip. If several spins are unstable, the one whose flipping reduces the energy the most is flipped. However, because an avalanche propagates outwards from the original site, and the left and right propagating directions are disjoint, the same results would be obtained if all unstable spins were flipped simultaneously. For all the numerical results reported in this paper, $`10^7`$ random choices of $`\{J_i,h_i\}`$ were tested.

Figure 2 shows a typical hysteresis loop, with random bond disorder but no random fields ($`h_i=0`$). The bonds are drawn from a distribution uniform over $`[1,0].`$ Return point memory is seen at $`H=1.4`$ the hysteresis loop. RPM is also found when the $`h_i`$’s are drawn from a distribution uniform over $`[1,1],`$ if $`J_i=1`$ for all $`i.`$ In both cases, although it cannot be shown in the figure, RPM exists for the full spin configuration rather than just the overall magnetization. However, if the $`J_i`$’s are not equal and the random fields are non-zero, we find that RPM fails if $`\delta h0.01`$ and $`\delta J0.01`$ in the random bond and random field cases respectively. We therefore conclude that either $`\delta h`$ or $`\delta J`$ must be zero for RPM. The results are the same for open and periodic boundary conditions caveat . (Even when $`\delta h`$ and $`\delta J`$ are both non-zero, the deviation from RPM is quite small, and hard to detect if one averages the hysteresis loop over realizations of randomness. A similar phenomenon was observed earlier for Sherrington-Kirkpatrick spin glasses Pazmandi .)

An important difference between the ferromagnetic and antiferromagnetic cases is that a spin at a single site can flip several times while the magnetic field is varied monotonically. Thus if the field is lowered, a spin pointing up can be triggered and flip down; if its neighbors have already flipped down, they can then be pushed back up by the new spin flip. As a result of this, the magnetization does not vary monotonically with $`H`$. This can be seen in the plot for a single realization of randomness in Figure 2. In more detail, it is possible observe that i) an avalanche that starts from a site and destabilizes both its neighbors is only possible for a configuration $``$ going to $``$ (or its mirror image), where the initial site is in the middle. The next nearest neighbors are stabilized, and the avalanche only covers three sites. ii) an avalanche that starts from a site and destabilizes only one neighbor is only possible for a configuration $``$ going to $``$ (or its mirror image), where the initial site is at the end. The avalanche only covers two sites. Thus as $`H`$ is varied, the chain evolves through single spin flips, two-site avalanches with $`\mathrm{\Delta }M=0`$ and three-site avalanches which have $`\mathrm{\Delta }M=1`$ for decreasing $`H`$ and $`\mathrm{\Delta }M=1`$ for increasing $`H.`$ These results and more have been proved earlier with random field disorder (without bond disorder), for the major hysteresis loop Shukla ; the full shape of the hysteresis loop is found analytically Shukla . We have extended the results of Shukla  to prove ii) for the major loop and i) for the entire hysteresis curve unpub .

Motivated by the observation of retrograde variation of the magnetization with $`H,`$ we construct an alternative representation of the dynamics in terms of spin flips instead of the spin configuration. Initially, when all the spins point up, the flip variable is zero at each site. Thereafter, each time a spin at site $`i`$ is reversed, the flip variable $`l_i`$ is increased by 1 if this happens when the field $`H`$ is increasing, and decreased by 1 if this happens when $`H`$ is decreasing. Clearly, along any branch of the hysteresis loop, while $`H`$ varies monotonically, so must each $`l_i.`$ Also, $`s_i=12[l_imod2],`$ and if two configurations $`\alpha `$ and $`\beta `$ satisfy the condition that $`l_i^\alpha l_i^\beta `$ is even for all $`i,`$ they correspond to the same spin state. In our numerical simulations, we find that for the cases when RPM is valid, it also holds for the flip configuration. Since the configuration $`\{s_i\}`$ is a projection of $`\{l_i\},`$ this is a stronger result than RPM, and suggests that the underlying dynamics in terms of $`\{l_i\}`$ is fundamental to random antiferromagnetic chains.

With $`m_i=_{ji}l_j,`$ we find that no-passing is satisfied: if $`H`$ is decreased from $`H_{max}`$ to $`H_{min},`$ increased to $`H_{max}`$ and then returned to $`H_{min},`$ for any $`H`$ and any site $`i`$ the value of $`m_i`$ on the third segment of this path is bounded below and above by the corresponding $`m_i`$’s on the first and second segments (see Figure 1). This is not true for the spin variables $`s_i`$ foot:li . However, as emphasized earlier, RPM is not satisfied if one starts from an arbitrary initial state at some $`H`$ instead of the saturated state, so that unlike the ferromagnetic case Sethna  the proof of RPM must take into account the ancestry of a state.

In order to see to what extent RPM is influenced by the dynamics used for the model, we now consider spin-exchange dynamics Kawasaki  instead of single spin flip Glauber . A pair of neighboring spins that are oriented opposite to each other are exchanged if it is energetically favorable to do so. Since such a move does not change the overall magnetization, in order for there to be a response to a magnetic field, we allow single spin flips at the two ends of the chain (only open chains are considered).

This problem can be solved by embedding the antiferromagnetic chain in a two dimensional ferromagnetic model. We first consider the case when there is only random bond disorder. Figure 3 shows a two dimensional square lattice of spins, rotated by an angle $`\pi /4.`$ Ferromagnetic bonds connect next nearest neighbors, but (without random fields) not nearest neighbors. As shown in the figure, the vertical bonds are all zero, and the horizontal bonds are identical within each vertical strip. At the top and bottom boundaries, the boundary conditions force all the spins to be up and down respectively. Free boundary conditions are used on the side walls. Thus in its ground state, there is one horizontal domain wall across the system. As shown in the figure, we adopt a convention in which the domain wall consists of line segments oriented at $`\pm \pi /4,`$ i.e. along the principal directions of the square lattice. The mapping from the two dimensional system to the one-dimensional chain is as follows: if any line segment of the two dimensional domain wall is oriented at $`\pi /4`$ or $`\pi /4,`$ the corresponding spin in the antiferromagnetic chain is 1 or -1 respectively. The two dimensional ground state corresponds to alternating spins in the chain, as is appropriate when $`H`$ is zero.

For a general shape of the domain wall, whenever two successive segments point in the same direction, a (horizontal) bond is broken, whereas this does not happen when they point in opposite directions. By choosing the horizontal bond strengths to be $`2J_1,2J_2,2J_3\mathrm{},`$ correlated vertically, the energy of the antiferromagnetic chain is increased by $`2J_i`$ when spins $`i`$ and $`i+1`$ point in the same direction compared to when they are opposite, as desired for an antiferromagnetic chain. The magnetic field $`H`$ couples to $`_is_i`$ for the chain, which is equivalent to the difference in height between the ends of the two dimensional domain wall. This is equivalent to a magnetic field $`H`$ on the rightmost column of the two dimensional system, with the left end of the domain wall tethered. It is also possible to generalize the model to include random bond disorder for the chain: nearest neighbor bonds of strength $`2|h_i|`$ are introduced in the $`i`$’th column, oriented at $`\pi /4`$ if $`h_i`$ is positive and $`\pi /4`$ if $`h_i`$ is negative.

With this construction, all bonds are ferromagnetic for the two dimensional system. Further, the fields at the side boundaries vary monotonically with $`H.`$ Further, spin-exchange for the chain is equivalent to single spin flips in the two-dimensional lattice. The results of Ref. Sethna  can therefore be invoked. We conclude that, with these dynamics, RPM is valid for all configurations, and is valid for simultaneous random field and random bond disorder. As we have seen, neither of these statements is valid for single spin flip dynamics for the chain; the two dimensional analog of spin flip at a site on the chain is to move the entire domain wall to the right of the site up or down by one unit if the spin flips up or down foot:untether .

In this paper we have shown that the hysteresis loop for a random Ising antiferromagnetic chains at zero temperature exhibits return point memory (RPM). For spin flip dynamics, the result is history dependent, being valid only for configurations that start from saturated magnetization and a large magnetic field. This is unlike the result for ferromagnets, where the result is valid for all configurations, indicating that the mechanism for RPM is different from the ferromagnetic case. (Also, RPM is only valid if either random field or random bond disorder is present, but not both, a restriction that does not apply to ferromagnets.) For spin exchange dynamics, we have proved RPM by mapping to a two-dimensional ferromagnetic model, and have therefore shown that it is as general: valid for all configurations, and with simultaneous random field and bond disorder. This implies that RPM depends on the Hamiltonian and the dynamics used.

We thank Peter Young for useful comments. AD acknowledges support from NSF grant DMR 0086287.


text:
# The CDF Calorimeter Upgrade for Run IIb

## 1 THE CENTRAL PRESHOWER AND CRACK DETECTORS

The CDF Central Preshower (CPR) and Central Crack (CCR) detectors will be replaced at the time the silicon detector is replaced for Run IIb.

In 1992 the CDF Collaboration installed gas detectors on the front face of the central calorimeter in order to sample early showers and cover the $`\varphi `$-cracks between calorimeter wedges, as shown in Fig. 1. The CPR has been extensively used in electron identification (ID), providing about a factor 2-3 more rejection of charged pions that pass all other cuts. This extra rejection has been crucial in soft electron ID for b-jet tagging, as was shown in the first top evidence paper . The CPR has been used in several publications involving photon ID. By using conversion rates, which are energy independent, it extended the QCD measurement of direct photons by more than 100 GeV in photon transverse momentum $`P_T`$ .

The CCR, located after a 10 radiation length thick tungsten bar, has been checked for large pulse heights in all the rare events CDF has observed in Run I. The $`\varphi `$-cracks cover about 8% of the central detector, and in events with multiple electromagnetic objects, the possibility of one object hitting the crack is quite high.

The present slow CPR and CCR gas

detectors will suffer the luminosity increase foreseen for Run IIb. In order to maintain the same Run I capabilities, they will be replaced by scintillator counters read out by Wave-Length Shifting (WLS) fibers. The new CPR will also have a better segmentation and will be used to improve the jet energy resolution by both correcting for energy loss in the dead material in front of it and adding its information in jet algorithms incorporating charged tracking.

### 1.1 The Detector Design fron Run IIb

The new CPR will be based on 2cm thick scintillator tiles segmented in $`\eta `$ and $`\varphi `$ and read out by a 1mm diameter WLS fiber running into a groove on the surface of each tile. Six tiles (12.5x12.5 cm2 each) will cover the front face of each calorimeter tower, and the tiles will be assembled in 48 modules like the one shown in Fig. 2 covering the 48 central calorimeter wedges. After leaving the tiles, the WLS fibers will be spliced to clear fibers which will terminate into plastic connectors at the higher $`\eta `$ edge of each module. There $``$5m long optical cables will transmit the light to 16-channel PhotoMultiplier Tubes (PMTs) at the back of the wedge. A current prototype, consisting in scintillator tiles provided by the CDF JINR Dubna group and Pol.Hi.Tech fibers, provided a light yield of $``$20 (12) photoelectrons at the exit of the tile (after all the optical chain), exceeding the design requirement.

The new CCR will use the same technique but the available space will limit the scintillator thickness to 5mm. Ten tiles, $``$5cm wide, will cover each $`\varphi `$-crack with the same calorimeter segmentation of 10 towers/wedge.

## 2 THE EM TIMING PROJECT

The CDF Collaboration is adding timing information into the readout of the Central and Plug electromagnetic calorimeters (CEM and PEM) using a technique similar to the hadron TDC system. This upgrade would significantly improve the potential of the CDF detector to do high-$`P_T`$ searches for new physics in data samples with photons in the final state by: 1) reducing the the cosmic ray and beam halo sources of background; 2) checking that all photons in unsual events are from the primary interaction. With sufficient calibration data, there is even the possibility of searching for very long-lived particles which decay (1-10ns) into photons.

### 2.1 The Hardware Project

The signal from the PMT goes to a Transition Board on the back of a calorimeter readout VME crate. All the lines are passed through the VME backplane into an Amplifier Shaper Discriminator (ASD) which effectively turns the signal into an LVDS digital pulse suitable for use by a TDC.

While the 960 PEM PMTs already have a dynode output designed into the base, a custom splitter is used for the output of the 960 CEM PMTs. The splitter is a fully passive element, completely connectorized, which works by inductively coupling the primary line (for the energy measurement) to the secondary output. The primary output loses a negligible amount of the charge and the secondary line only takes $``$15% of the output voltage for use to fire the ASD/TDC system.

Several tests of the splitter, both on the test bench and on the detector itself, show no difference between the input shape and height before and after inserting the splitter into the system. For energies $`>`$ 4 (2) GeV for the CEM (PEM), the splitter fires the ASD/TDC system with 100% efficiency. The intrinsic timing resolution, measured using the CEM LED and a splitter into a TDC channel of the hadron calorimeter, is $``$1.1ns, dominated by the 1ns TDC resolution.


text:
# Dissipationless collapse of a set of 𝑁 massive particles

## 1 Introduction

It is intuitive that the gravitational collapse of a set of $`N`$ masses is directly related to the formation of astrophysical structures like globular clusters or elliptical galaxies (the presence of gas may complicate the pure gravitational $`N`$-body problem for spiral galaxies). From an analytical point of view, this problem is very difficult. When $`N`$ is much larger than 2, direct approach is intractable, and since Poincaré results of non analyticity, exact solutions may be unobtainable. In the context of statistical physics, the situation is more favorable and, in a dissipationless approximation111 The dissipationless hypothesis is widely accepted in the context of gravitational $`N`$-body problem because the ratio of the two-body relaxation time over the dynamical time is of the order of $`N`$. For a system composed of more than $`10^4`$ massive particles a study during a few hundreds dynamical times can really be considered as dissipationless, the unique source of dissipation being two-body encounters., leads to the Collisionless Boltzmann Equation (hereafter denoted by CBE)

$$\frac{f}{t}+𝐩.\frac{f}{𝐫}+m\frac{\psi }{𝐫}.\frac{f}{𝐩}=0$$
(1)
where $`f=f(𝐫,𝐩,t)`$ and $`\psi =\psi (𝐫,t)`$ are respectively the distribution function of the system with respect to the canonically conjugated $`(𝐫,𝐩)`$ phase space variables and the mean field gravitational potential. As noted initially by Hénon (1960), this formalism holds for such systems if and only if we consider $`N`$ identical point masses equal to $`m`$. This problem splits naturally into two related parts: the time dependent regime and the stationary state. We can reasonably think that these two problems are not completely understood. The transient time dependent regime was investigated mainly considering self-similar solutions (Lynden-Bell & Eggleton (1980), Henriksen & Widrow (1995), Blottiau et al. (1988) and Lancellotti & Kiessling (2001)). These studies conclude that power law solutions can exist for the spatial dependence of the gravitational potential (with various powers). Nevertheless, there is no study which indicates clearly that the time dependence of the solutions disappears in a few dynamical times, giving a well defined equilibrium-like state. On the other hand, applying Jeans theorem (e.g. Binney & Tremaine (1987) hereafter BT87, p. 220), it is quite easy to find a stationary solution. For example, every positive and integrable function of the mean field energy per mass unit $`E`$ is a potential equilibrium distribution function for a spherical isotropic system. Several approaches are possible to choose the equilibrium distribution function. Thermodynamics (Violent Relaxation paradigm: Lynden-Bell (1967), Chavanis (2002), Nakamura (2000)) indicate that isothermal spheres or polytropic systems are good candidates. Stability analyses can be split into two categories. In the CBE context (see Perez & Aly (1998) for a review), it is well known that spherical systems (with decreasing spatial density) are generally stable except in the case where a large radial anisotropy is present in the velocity space. This is the Radial Orbit Instability, hereafter denoted by ROI (see Perez & Aly (1998), and Perez et al. (1998) for a detailed analytic and numeric study of these phenomena) which leads to a bar-like equilibrium state in a few dynamical times. In the context of thermodynamics of self-gravitating systems, in a pioneering work by Antonov (1962), it was shown that an important density contrast leads to the collapse of the core of system (see Chavanis (2003) for details).
In all these studies there is no definitive conclusion, and the choice of the equilibrium distribution remains unclear. Introducing observations and taking into account analytical constraints, several models are possible: chronologically, we can cite (see for example BT87, p. 223-239) the Plummer model (or other polytropic models), de Vaucouleurs $`r^{1/4}`$ law, King and isochrone Hénon model or more recently the very simple but interesting Hernquist model (Hernquist (1990)) for spherical isotropic systems. In the anisotropic case, Ossipkov-Merritt or generalized polytropes can be considered. Finally for non spherical systems, there also exists some models reviewed in BT87 (p. 245-266). Considering this wide variety of possibilities, one can try to make accurate numerical simulations to clarify the situation. Surprisingly, such a program has not been completely carried on. In a pioneering work, van Albada (1982) remarked that the dissipationless collapse of a clumpy cloud of $`N`$ equal masses could lead to a final stationary state that is quite similar to elliptical galaxies. This kind of study was reconsidered in an important work by Aguilar & Merritt (1990), with more details and a crucial remark concerning the correlation between the final shape (spherical or oblate) and the virial ratio of the initial state. These authors explain this feature invoking ROI. Some more recent studies (Cannizzo & Hollister (1992), Boily et al. (1999) and Theis & Spurzem (1999)) concentrate on some particularities of the preceding works. Finally, two works (Dantas et al. (2002) and Carpintero & Muzzio (1995)) develop new ideas considering the influence of the Hubble flow on the collapse. However, the problem is only partially depicted.
The aim of this paper is to analyse the dissipationless collapse of a large set of N Body systems with a very wide variety of ‘realistic’ initial conditions. As we will see, the small number of particles involved, the numerical technique or the specificity of the previous works did not allow their authors to reach a sufficiently precise conclusion. The layout of this paper is as follows. In section 2 we describe in detail the numerical procedures used in our experiments. Section 3 describes the results we have obtained. These results are then interpreted in section 4, where some conclusions and perspectives are also proposed.

## 2 Numerical procedures

### 2.1 Dynamics

The Treecode used to perform our simulations is a modified version of the Barnes & Hut (1986) Treecode, parallelised by D. Pfenniger using the MPI library. We implemented some computations of observables and adapted the code to suit our specific problems. The main features of this code are a hierarchical $`O(N\mathrm{log}(N))`$ algorithm for the gravitational force calculation and a leap-frog algorithm for the time integration. We introduced an adaptative time step, based on a very simple physical consideration. The time step is equal to a fraction $`n_{ts}`$ of the instantaneous dynamical time $`T_d`$ (222The fraction $`n_{ts}`$ is adapted to the virial parameter $`\eta `$ and ranges roughly from $`n_{ts}=300`$ when $`\eta =90`$ to $`n_{ts}=5000`$ when $`\eta =08`$. The dynamical time we used is given by
$$T_d=\frac{{\displaystyle \underset{i=1}{\overset{N}{}}}\sqrt{x_i^2+y_i^2+z_i^2}}{{\displaystyle \underset{i=1}{\overset{N}{}}}\sqrt{vx_i^2+vy_i^2+vz_i^2}}$$
) , i.e. $`\mathrm{\Delta }t=T_d/n_{ts}`$ . The simulations were run on a Beowulf cluster (25 dual CPU processors whose speed ranges from 400MHz to 1GHz).

### 2.2 Initial Conditions

The initial virial ratio is an important parameter in our simulations. The following method was adopted to set the virial ratio to the value $`V_{initial}`$. Positions $`𝐫_i`$ and velocities $`𝐯_i`$ are generated. We can then compute

$$V_p=\frac{2K}{U}$$
(2)
where

$$K=\underset{i=1}{\overset{N}{}}\frac{1}{2}m_i𝐯_i^2$$
(3)
and

$$U=\frac{G}{2}\underset{ij}{\overset{N}{}}\frac{m_im_j}{\left(max((𝐫_i𝐫_j)^2,ϵ^2)\right)^{1/2}}.$$
(4)
In this relation $`ϵ`$ is a softening parameter whose value is discussed in section 2.3.2. As the potential energy depends only on the positions, we obtain a system with a virial ratio equal to $`V_{initial}`$ just by multiplying all the particle velocities by the factor $`\left(V_{initial}/V_p\right)^{1/2}`$. For convenience we define

$$\eta =\left|V_{initial}\right|\times 10^2$$
(5)
#### 2.2.1 Homogeneous density distribution (Hη)

As we study large $`N`$-body systems, we can produce a homogeneous density by generating positions randomly. These systems are also isotropic. We produce the isotropic velocity distribution by generating velocities randomly.

#### 2.2.2 Clumpy density distribution (C$`{}_{}{}^{n}{}_{\eta }{}^{}`$)

A type of inhomogeneous systems is made of systems with a clumpy density distribution. We first generate $`n`$ small homogeneous spherical systems with radius $`R_g`$. Centers of these subsystems are uniformly distributed in the system. The empty space is then filled using a homogeneous density distribution. In the initial state, each clump contains about 1% of the total mass of the system and have a radius which represents $`5\%`$ of the initial radius of the whole system. These systems are isotropic.

#### 2.2.3 Power law $`r^\alpha `$ density distribution (P$`{}_{}{}^{\alpha }{}_{\eta }{}^{}`$)

We first generate the $`\phi `$ and $`z`$ cylindrical coordinates using two uniform random numbers $`u_1`$ and $`u_2`$:

$$(z,\phi )=(2u_11,2\pi u_2).$$
(6)
Using the inverse transformation method, if

$$r=RF^1\left(u\right)\text{with }F\left(r\right)=\frac{1}{S}_\iota ^rx^{2\alpha }𝑑x$$
(7)
where $`R`$ is the radius of the system, $`u`$ is a uniform random number, $`\iota 1`$ and

$$S=_0^1x^{2\alpha }𝑑x,$$
(8)
then the probability density of $`r`$ is proportional to $`r^{2\alpha }`$, and the mass density $`\rho `$ is proportional to $`r^\alpha `$. Finally, one gets

$$𝐫=\left[\begin{array}{c}r\sqrt{1z^2}\mathrm{cos}\phi \\ r\sqrt{1z^2}\mathrm{sin}\phi \\ rz\end{array}\right]$$
(9)
These systems are isotropic.

#### 2.2.4 Gaussian velocity distribution(G$`{}_{}{}^{\sigma }{}_{\eta }{}^{}`$)

Most of the systems we use have a uniform velocity distribution. But we have also performed simulations with systems presenting a gaussian initial velocity distribution. These systems are isotropic, but the x-, y- and z-components of the velocity are generated following a gaussian distribution. Using a standard method we generate two uniform random numbers $`u_1`$ and $`u_2`$, and set

$$v_i=\sqrt{2\sigma ^2\mathrm{ln}u_1}\mathrm{cos}(2\pi \sigma ^2u_2)i=x,y,z$$
(10)
where $`\sigma `$ is the gaussian standard deviation.

#### 2.2.5 Global rotation (R$`{}_{}{}^{f}{}_{\eta }{}^{}`$)

Some of our initial systems are homogeneous systems with a global rotation around the $`Z`$-axis. The method we choose to generate such initial conditions is the following. We create a homogeneous and isotropic system (an H-type system). We then compute the average velocity of the particles.

$$\overline{v}=\frac{1}{N}\underset{i=1}{\overset{N}{}}𝐯_i$$
(11)
We project the velocities on a spherical referential, and add a fraction of $`\overline{v}`$ to $`v_\varphi `$ with regard to the position of the particle. We set

$$v_{i,\varphi }=v_{i,\varphi }+f\frac{\rho _i\overline{v}}{R}$$
(12)
where $`f`$ is a parameter of the initial condition, $`\rho _i`$ is the distance from the particle to the $`Z`$-axis and $`R`$ the radius of the system. The amount of rotation induced by this method can be evaluated through the ratio:

$`\mu =K_{rot}/K,`$ (13)
where $`K`$ is the total kinetic energy defined above, whereas $`K_{rot}`$ is the rotation kinetic energy defined by Navarro & White (1993):

$$K_{rot}=\frac{1}{2}\underset{i=1}{\overset{N}{}}m_i\frac{(𝐉_i\widehat{𝐉}_{tot})^2}{[𝐫_i^2(𝐫_i\widehat{𝐉}_{tot})^2]}$$
(14)
Above, $`𝐉_i`$ is the specific angular momentum of particle i, $`\widehat{𝐉}_{tot}`$ is a unit vector in the direction of the total angular momentum of the system. In order to exclude counter-rotating particles, the sum in equation (14) is actually carried out only over those particles which verify the condition $`(𝐉_i\widehat{𝐉}_{tot})>0`$.

#### 2.2.6 Power-law initial mass function(M$`{}_{}{}^{k}{}_{\eta }{}^{}`$)

Almost all the simulations we made assume particles with equal masses. However, we have created some initial systems with a power-law mass function, like

$$n(M)=\alpha M^\beta $$
(15)
The number of particles whose mass is $`MmM+dM`$ is $`n(M)dM`$. In some models, the value of $`\alpha `$ and $`\beta `$ depends on the range of mass that is considered. We have used several types of mass functions, among them the initial mass function given by Kroupa (2001) $`(k=I)`$, the one given by Salpeter (1955) $`(k=II)`$ and an $`M^1`$ mass function $`\left(k=III\right)`$. In order to generate masses following these functions, we first calculate $`\alpha _k`$ to produce a continuous function. We then can calculate the number of particles whose mass is between $`M`$ and $`M+dM`$. We generate $`n(M)`$ masses

$$m_i=M+udM1in(M)$$
(16)
where $`0u1`$ is a uniform random number. In the initial state, these systems have a homogeneous number density, a quasi homogeneous mass density and they are isotropic.

#### 2.2.7 Nomenclature

We indicate below the whole set of our non rotating initial conditions.

* Homogeneous Hη models: H88, H79, H60, H50, H40, H30, H20, H15 and H10
* Clumpy C$`{}_{}{}^{n}{}_{\eta }{}^{}`$ models: C$`{}_{}{}^{20}{}_{67}{}^{}`$, C$`{}_{}{}^{20}{}_{65}{}^{}`$, C$`{}_{}{}^{20}{}_{61}{}^{}`$, C$`{}_{}{}^{20}{}_{48}{}^{}`$, C$`{}_{}{}^{20}{}_{39}{}^{}`$, C$`{}_{}{}^{20}{}_{29}{}^{}`$, C$`{}_{}{}^{20}{}_{14}{}^{}`$, C$`{}_{}{}^{20}{}_{10}{}^{}`$, C$`{}_{}{}^{20}{}_{07}{}^{}`$ and C$`{}_{}{}^{03}{}_{10}{}^{}`$
* Power Law P$`{}_{}{}^{\alpha }{}_{\eta }{}^{}`$ models: P$`{}_{}{}^{2.0}{}_{50}{}^{}`$, P$`{}_{}{}^{2.0}{}_{09}{}^{}`$, P$`{}_{}{}^{1.0}{}_{50}{}^{}`$, P$`{}_{}{}^{0.5}{}_{50}{}^{}`$, P$`{}_{}{}^{1.0}{}_{10}{}^{}`$, P$`{}_{}{}^{1.5}{}_{08}{}^{}`$ and P$`{}_{}{}^{1.5}{}_{40}{}^{}`$
* Gaussian velocity profiles G$`{}_{}{}^{\sigma }{}_{\eta }{}^{}`$ models: G$`{}_{}{}^{1}{}_{50}{}^{}`$, G$`{}_{}{}^{2}{}_{50}{}^{}`$, G$`{}_{}{}^{3}{}_{50}{}^{}`$, G$`{}_{}{}^{4}{}_{12}{}^{}`$ and G$`{}_{}{}^{5}{}_{50}{}^{}`$
* Mass spectra M$`{}_{}{}^{k}{}_{\eta }{}^{}`$ models: M$`{}_{}{}^{I}{}_{50}{}^{}`$, M$`{}_{}{}^{II}{}_{50}{}^{}`$, M$`{}_{}{}^{III}{}_{51}{}^{}`$, M$`{}_{}{}^{I}{}_{35}{}^{}`$, M$`{}_{}{}^{II}{}_{25}{}^{}`$, M$`{}_{}{}^{III}{}_{15}{}^{}`$ and M$`{}_{}{}^{I}{}_{07}{}^{}`$

For all these models we ran the numerical simulations with 30 000 particles (see § 3.1)

### 2.3 Observables

#### 2.3.1 Units

Our units are not the commonly used ones (see Heggie & Mathieu (1986)). We did not set the total energy $`E`$ of the system to $`0.25`$ because we wanted to prescribe instead the initial virial ratio $`V_{initial}`$, the size $`R`$ of the system and its mass $`M`$. We thus have $`M=1`$, $`R=10`$ and $`G=1`$, and values of $`V_{initial}`$ and $`E`$ depending on the simulation. We can link the units we have used with more standard ones. We have chosen the following relationships between our units of length and mass and common astrophysical ones:

$$M=10^6M_{\mathrm{}}\text{ and }R=10pc$$
(17)
Our unit of time $`u_t`$ is given by:

$$1u_t=\sqrt{\frac{R_c^3G_sM_s}{R_s^3G_cM_c}}\mathrm{4.72\hspace{0.17em}\hspace{0.17em}10}^{11}s=\mathrm{1.50\hspace{0.17em}\hspace{0.17em}10}^4yr$$
(18)
where variables $`X_s`$ are expressed in our simulation units and variables $`X_c`$ in standard units.

#### 2.3.2 Potential softening and energy conservation

The non conservation of the energy during the numerical evolution has three main sources.
The softening parameter $`\epsilon `$ introduced in the potential calculus (cf. equation 4) is an obvious one. This parameter introduces a lower cutoff $`\mathrm{\Lambda }`$ in the resolution of length in the simulations. Following Barnes & Hut (1989), structural details up to scale $`\mathrm{\Lambda }10\epsilon `$ are sensitive to the value of $`\epsilon `$. Moreover, in order to be compatible with the collisionless hypothesis, the softening parameter must be greater than the scale where important collisions can occur. Still following Barnes & Hut (1989), this causes

$$\frac{\epsilon }{10}\frac{Gm}{v^2}$$
(19)
In our collapse simulation with $`310^4`$ particles, this results in $`\epsilon \eta 2/3`$. The discretization of time integration introduces inevitably another source of energy non conservation, particularly during the collapse. The force computation also generates errors. The choice of the opening angle $`\mathrm{\Theta }`$, which governs the accuracy of the force calculation of the Treecode, is a compromise between speed and accuracy. For all these reasons, we have adopted $`\epsilon =0.1`$. This choice imposes $`\eta 6`$ (for $`3010^4`$ particles). This trade-off allowed to perform simulations with less than $`1\%`$ energy variation without requiring too much computing time. For each of our experiments, the total CPU time ranges between 3 to 24 hours for $`3000u_t`$ and $`310^4`$ particles. The total agregated CPU time of all our collapse experiments is approximately 6 months.
We have tested two other values of the softening parameter ($`\epsilon =0.03`$ and $`\epsilon =0.3`$) for several typical simulations. These tests did not reveal significant variations of the computed observables.

#### 2.3.3 Spatial indicators

As indicators of the geometry of the system, we computed axial ratios, radii containing $`10\%`$ ($`R_{10}`$), $`50\%`$ ($`R_{50}`$) and $`90\%`$ ($`R_{90}`$) of the mass, density profile $`\rho \left(r\right)`$ and equilibrium core radius. The axial ratios are computed with the eigenvalues $`\lambda _1`$, $`\lambda _2`$ and $`\lambda _3`$ of the (3x3) inertia matrix $`I`$, where $`\lambda _3\lambda _2\lambda _1`$ and, if the position of the particle $`i`$ is $`𝐫_i=(x_{1,i};x_{2,i};x_{3,i})`$

$$\{\begin{array}{ccc}I_{\mu \nu }=\underset{i=1}{\overset{N}{}}m_ix_{\mu ,i}x_{\nu ,i}\hfill & & \text{for }\mu \nu =1,2,3\hfill \\ & & \\ I_{\mu \mu }=\underset{i=1}{\overset{N}{}}m_i(r_i^2x_{\mu ,i}^2)\hfill & & \text{for }\mu =1,2,3\hfill \end{array}$$
(20)
The axial ratios $`a_1`$ and $`a_2`$ are given by $`a_2=\lambda _1/\lambda _21.0`$ and $`a_1=\lambda _3/\lambda _21.0`$.
The density profile $`\rho `$, which depends only on the radius $`r`$, together with the $`R_\delta `$ ($`\delta =10,50,90`$) have a physical meaning only for spherical or nearly spherical systems. For all the spatial indicators computations we have only considered particles whose distance to the center of mass of the system is less than $`6\times R_{50}`$ of the system. This assumption excludes particles which are inevitably ‘ejected’ during the collapse333The number of excluded particles ranges from 0% to 30% of the total number of particles, depending mostly on $`\eta `$. For example, the number of excluded particles is 0% for H80, 3% for C$`{}_{67}{}^{}{}_{}{}^{20}`$, 5% for H50, 22% for C$`{}_{10}{}^{}{}_{}{}^{20}`$ and 31% for H10..
After the collapse a core-halo structure forms in the system. In order to measure the radius of the core, we have computed the density radius as defined by Casertano and Hut (see Casertano & Hut (1985)). The density radius is a good estimator of the theoretical and observational core radius.
We have also computed the radial density of the system. The density is computed by dividing the system into spherical bins and by calculating the total mass in each bin.

#### 2.3.4 Statistical indicators

When the system has reached an equilibrium state, we compute the temperature of the system

$$T=\frac{2K}{3Nk_B}$$
(21)
where $`K`$ is the kinetic energy of the system, $`k_B`$ is the Boltzmann constant (which is set to 1) and the notation $`A`$ denotes the mean value of the observable $`A`$, defined by

$$A=\frac{1}{N}\underset{i=1}{\overset{N}{}}A_i$$
(22)
In order to characterise the system in the velocity space we have computed the function

$$\kappa \left(r\right)=\frac{2v_{i,rad}^2_{rr_i<r+dr}}{v_{i,tan}^2_{rr_i<r+dr}}$$
(23)
where $`v_{i,rad}`$ is the radial velocity of the ith particle, and $`v_{i,tan}`$ its tangential velocity. For spherical and isotropic systems ($`a_1a_21`$ and $`\kappa \left(r\right)1`$), we have fitted the density by
1- a polytropic law

$$\rho =\rho _0\psi ^\gamma $$
(24)
which corresponds to a distribution function

$$f\left(E\right)E^{\gamma 3/2}$$
(25)
2- an isothermal sphere law

$$\rho =\rho _1e^{\psi /s^2}$$
(26)
which corresponds to a distribution function

$$f\left(E\right)=\frac{\rho _1}{\left(2\pi s^2\right)^{3/2}}e^{\frac{E}{s^2}}$$
(27)
Using the least square method in the $`\mathrm{ln}(\rho )\mathrm{ln}(\psi )`$ plane we get $`(\gamma ,\mathrm{ln}\left(\rho _0\right))`$ and $`(s^2,\mathrm{ln}\left(\rho _1\right))`$.

## 3 Description of the results

We have only studied systems with an initial virial ratio corresponding to $`\eta [7,88]`$. In such systems, the initial velocity dispersion cannot balance the gravitational field. This produces a collapse. After this collapse, in all our simulations the system reaches an equilibrium state characterised by a temporal mean value of the virial ratio equal to $`1`$, i.e $`\eta =100`$, and stationary physical observables. These quantities (defined in section 2) are presented in a table of results in the appendix of this paper. The following results will be discussed and interpreted in section 4.

### 3.1 Relevant number of particles

In all previous works on this subject (van Albada (1982), Aguilar & Merritt (1990), Cannizzo & Hollister (1992) and Boily et al. (1999)), the authors did not really consider the influence of the number of particles on their results. In the first two and more general works, this number is rather small (not more than a few thousands in the largest simulations). The two other studies are more specific and use typically $`10^4`$ and, in a few reference cases, $`2.10^4`$ particles.

In order to test the influence of the number of particles on the final results, we have computed several physical observables of some collapsing systems with various numbers of particles. The results are presented in Figure 1. In order to check the influence of $`N`$ in the whole phase space, we have studied positions and velocities related observables: $`a_1`$, $`a_2`$, $`R_{10}`$, $`R_{50}`$ and $`R_{90}`$ and parameters of isothermal and polytropic fit models namely $`\gamma `$ and $`s^2`$. Moreover, in order to be model-independent, we have studied three representative initial conditions: H80, H50 and H10, i.e. respectively initially hot, warm and cold systems. The number of particles used in each case ranges from $`10^2`$ to $`10^5.`$ We can see in Figure 1 that some observables are N-dependent when $`N<3.10^4`$. In particular, $`R_{50}`$ and $`R_d`$ present a monotonic variation larger than $`50\%`$ when $`N`$ varies from $`10^2`$ to $`10^5`$ and the ellipticity of the final state is overestimated for small values of $`N`$. As a conclusion of this preliminary study, we claim that the relevant number of particles for collapse simulations is $`N3.10^4`$. As all simulations have been completed with a total energy loss smaller than 1%, we state moreover that this result is independent of the numerical scheme used (Treecode or Direct N-Body). As a consequence, the simulations presented hereafter have been performed using $`N=3.10^4`$ particles.

### 3.2 Morphological segregation

An important study by Aguilar & Merritt (1990) shows that, in the case of an initial density profile $`\rho r^1`$, the shape of the virialized state depends on $`\eta `$: a very small $`\eta `$ leads to a flattened equilibrium state, when a more quiet collapse produces a spherical one. Our investigations concern a wide range of different initial conditions (homogeneous, clumpy,…) and show that the influence of $`\eta `$ depends on the properties of the initial system444 The particular case of rotating initial conditions is discussed in a special section.. Figure 2 shows the axial ratios of the equilibrium state reached by our simulations. In fact, only a few simulations produced a final state with an ellipticity greater than E1. Every homogeneous initial condition (i.e. Hη, G$`{}_{\eta }{}^{}{}_{}{}^{\sigma }`$ and M$`{}_{\eta }{}^{}{}_{}{}^{k}`$) resulted in a spherical equilibrium state independently of the values of $`\eta `$ we tested. Cold clumpy systems have a weakly flattened equilibrium state. The only final systems with an ellipticity significantly greater than E1 are those produced by the collapse of cold P$`{}_{\eta }{}^{}{}_{}{}^{\alpha }`$.
Previous studies invoked ROI to explain this morphological segregation. However, it seems that ROI requires inhomogeneities near the center to be triggered.

### 3.3 Characteristic size segregation

In addition to the morphological segregation, presented in the previous section, we discovered a finer phenomenon.
In the Figures 3, 4, 5, 6 and 7, we have plotted the mass density of all equilibrium states produced by the collapse of our initial conditions as a function of the ratio $`r/R_{50}`$. These plots represent the density at the end of our simulations (after about 100 crossing times). These functions do not significantly evolve after the collapse except for M$`{}_{}{}^{I}{}_{07}{}^{}`$. For this special case, a comparative plot is the subject of Figure 8.

All equilibrium states we obtain clearly fall into two categories:

* Flat Core Systems
All these systems present a core halo structure, i.e. a large central region with a constant density and a steep envelope. These systems are typically such that $`\mathrm{ln}\left(R_{50}/R_d\right)<0.5`$ and $`\mathrm{ln}\left(R_{10}/R_d\right)<0.05`$.
* Small Core Systems
For such systems, the central density is two order of magnitude larger than for Flat Core systems. There is no central plateau and the density falls down regularly outward. These systems are typically such that $`\mathrm{ln}\left(R_{50}/R_d\right)>0.7`$ and $`\mathrm{ln}\left(R_{10}/R_d\right)>0.1`$.

The diagram $`\mathrm{ln}\left(R_{10}/R_d\right)`$ vs $`\mathrm{ln}\left(R_{50}/R_d\right)`$ is the subject of Figure 9. One can see in this figure that each equilibrium state belongs to one or the other family except in a few particular cases. In the Flat Core family we found all Hη, G$`{}_{}{}^{\sigma }{}_{\eta }{}^{}`$ and M$`{}_{}{}^{k}{}_{\eta }{}^{}`$ systems except M$`{}_{}{}^{I}{}_{7}{}^{}`$, and two P$`{}_{}{}^{\alpha }{}_{\eta }{}^{}`$ systems namely P$`{}_{}{}^{0.5}{}_{50}{}^{}`$ and P$`{}_{}{}^{1}{}_{50}{}^{}`$. These systems are all initially homogeneous or slightly inhomogeneous (e.g. P$`{}_{}{}^{0.5}{}_{50}{}^{}`$ and P$`{}_{}{}^{1}{}_{50}{}^{}`$ systems). In the Small Core family, we found all C$`{}_{}{}^{20}{}_{\eta }{}^{}`$ and the P$`{}_{}{}^{2.0}{}_{50}{}^{}`$ and P$`{}_{}{}^{2.0}{}_{09}{}^{}`$ systems. These systems are all initially rather very inhomogeneous. Finally, there are 5 systems in-between the two categories: C$`{}_{}{}^{03}{}_{10}{}^{}`$, P$`{}_{}{}^{1.5}{}_{40}{}^{}`$, P$`{}_{}{}^{1.5}{}_{08}{}^{}`$, P$`{}_{}{}^{1}{}_{10}{}^{}`$ and M$`{}_{}{}^{I}{}_{07}{}^{}`$. This last model is the only one which migrates from Flat core set (when $`t10Td`$) to the edge of the small core region (when $`t100Td`$).

### 3.4 Equilibrium Distribution Function

In order to compare systems in the whole phase space, we fitted the equilibrium state reached by each system with two distinct isotropic models, e.g. polytropic and isothermal (see equations 24-27 or BT87, p.223-232). Figure 10 shows these two fits for the P$`{}_{}{}^{0.5}{}_{50}{}^{}`$ simulation. The technique used for the fit is described in section 2 of this paper. The result obtained for this special study is the following: the equilibrium states reached by our initial conditions can be fitted by the two models with a good level of accuracy. As long as $`\eta <70`$, the polytropic fit gives a mean value $`\gamma =4.77`$ with a standard deviation of $`\mathrm{\hspace{0.17em}2.48\hspace{0.33em}10}^1`$. This deviation represents $`5.1\%`$ of the mean value. This value corresponds typically to the well known Plummer model for which $`\gamma =5`$ (see BT87 P.224 for details). When the collapse is very quiet ( typically $`\eta >70`$ ) polytropic fit is always very good but the value of the index is much larger than Plummer model, e.g. $`\gamma =6.86`$ for H79 and $`\gamma =7.37`$ for H88. The corresponding plot is the subject of the Figure 12. All the data can be found in appendix. As we can see on the example plotted in Figure 10, the isothermal fit is generally not as good as the polytropic one. On the whole set of equilibria, isothermal fits give a mean value $`s^2=\mathrm{2.5\hspace{0.17em}\hspace{0.17em}10}^2`$ with a standard deviation of $`\mathrm{1.6\hspace{0.17em}\hspace{0.17em}10}^2`$ ($`60\%`$). The corresponding plot is the subject of Figure 11. All the data can be found in appendix.

In fact both isothermal and polytropic fits are reasonable: as long as the model is able to reproduce a core halo structure the fit is correct. The success of the Plummer model, which density is given by

$$\rho \left(r\right)=\frac{3}{4\pi b^3}\left[1+\left(\frac{r}{b}\right)^2\right]^{5/2}$$

can be explained by its ability to fit a wide range of models with various ratio of the core size over the half-mass radius. The adjustment of this ratio is made possible by varying the free parameter $`b`$. We expect that other core halo models like King or Hernquist models work as well as the Plummer model. As a conclusion of this section, let us say that as predicted by theory there is not a single universal model to describe the equilibrium state of isotropic spherical self-gravitating system.

### 3.5 Influence of rotation

We saw in section 3.2 a source of flattening for self-gravitating equilibrium. Let us now show the influence of initial rotation, which is a natural candidate to produce flattening. The way we have added a global rotation and the significance of our rotation parameters $`f`$ and $`\mu `$ are explained in section 2.2.5.

The set of simulations performed for this study contains 31 different elements. The initial virial ratio ranges from $`\eta =10`$ to $`\eta =50`$, and the rotation parameter from $`f=0`$ (i.e. $`\mu =0`$) to $`f=20`$ (i.e. $`\mu =0.16`$ when $`\eta =50`$). As a matter of fact, equilibrium states always preserve a rather important part of the initial rotation555We observed that $`\mu `$ is always smaller in the equilibrium state than in the initial one, typically each rotating systems conserves 65% of the initial $`\mu `$ and, observed elliptical gravitating systems generally possess very small amount of rotation (see e.g. Combes et al. (1995)). We thus exclude large values of $`f`$ .
Our experiments exhibit two main features (see Figure 13): on the one hand, rotation produces a flattened equilibrium state only when $`f`$ exceeds a triggering value (typically $`f=f_o4`$). On the other hand, we have found that for a given value of $`\eta `$, the flatness of the equilibrium is roughly $`f`$independent, provided that $`f>f_o`$.

### 3.6 Thermodynamical segregation

As we study isolated systems, the total energy $`E`$ contained in the system is constant during the considered dynamical evolution . This property remains true as long as we consider collisionless evolutions. For gravitational systems, this means that we cannot carry out any simulation of duration larger than a few hundred dynamical times. We have obviously taken this constraint into account in our experiments. All systems which experience a violent relaxation reach an equilibrium state which is stationary in the whole phase space. Spatial behaviour like morphological segregation produced by ROI was confirmed and further detailled thanks to our study. A new size segregation was found in section 3.3. Now let us consider another new segregation appearing in the velocity space. Each equilibrium state is associated with a constant temperature $`T`$, calculated using equation $`\left(\text{21}\right)`$. More precisely, we have calculated the temporal mean value 666 The temporal mean value is computed from the time when the equilibrium is reached until the end of the simulation of the temperatures, evaluated every one hundred time units. As we can see in Figure 14, after the collapse and whatever the nature of the initial system is, the temperature is a very stable parameter.
Figure 15 shows the $`ET`$ diagram of the set of all non rotating simulations. It reveals a very interesting feature of post-collapse self-gravitating systems.

On the one hand, the set of systems with a total energy $`E>0.054`$ can be linearly fitted in the $`ET`$ plane. We call this set Low Branch 1 (hereafter denoted by LB1, see Figure 15). On the other hand, the set of the systems with a total energy $`E<0.054`$ splits into two families. The first is an exact continuation of LB1. Hence we named it Low Branch 2 (hereafter denoted by LB2). The second can also be linearly fitted, but with a much greater slope (one order of magnitude). We label this family High Branch (hereafter denoted by HB).
In LB1 or LB2, we find every H, G, and M systems with $`\eta >25`$, every P and every C with $`n>10`$. In HB, we find C$`{}_{}{}^{03}{}_{10}{}^{}`$ and every H, G and M systems with $`\eta <25`$. This segregation thus affects violent collapses (cold initial data): on the one hand, when $`\eta >25`$ all systems are on LB1, on the other hand for $`\eta <25`$, initially homogeneous or quasi-homogeneous (e.g. C$`{}_{10}{}^{}{}_{}{}^{03}`$) systems reach HB when initially inhomogeneous systems stay on LB2 instead.

## 4 Interpretations, Conclusions and Perspectives

Let us now recapitulate the results we have obtained and propose an interpretation:

1. The equilibrium state produced by the collapse of a set of $`N`$ gravitating particles is $`N`$independent provided that $`N>\mathrm{3.0\hspace{0.17em}10}^4`$.
2. Without any rotation, the dissipationless collapse of a set of gravitating particles can produce two relatively distinct equilibrium states:

* If the initial set is homogeneous, the equilibrium has a large core and a steep envelope.
    * If the initial set contains significant inhomogeneities ($`n>10`$ for clumpy systems or $`\alpha >1`$ for power law systems), the equilibrium state has only a small core around of which the density falls down regularly.

The explanation of this core size segregation is clear: it is associated to the Antonov core-collapse instability occurring when the density contrast between central and outward region of a gravitating system is very big. As a matter of fact, if the initial set contains inhomogeneities, these collapse much more quickly than the whole system777Because their Jeans length is much more smaller than the one of the whole system. and fall quickly into the central regions. The density contrast becomes then very large and the Antonov instability triggers producing a core collapse phenomenon. The rest of the system then smoothly collapses around this collapsed core. If there are no inhomogeneities in the initial set, the system collapses as a whole, central density grows slowly without reaching the triggering value of the Antonov instability. A large core then forms. Later evolution can also produce core collapse: this is what occurs for our M$`{}_{}{}^{I}{}_{07}{}^{}`$ system (see Figure 8). This is an initially homogeneous system with Kroupa mass spectrum which suffers a very strong collapse. As the mass spectrum is not sufficient to bring quickly enough a lot of mass in the center of the system, Antonov instability does not trigger and a large core forms. As the collapse is very violent, an increasing significant part of particles are progressively ”ejected” and the core collapse takes progressively place. This is the same phenomenon which is generally invoked to explain the collapsed core of some old globular clusters (e.g. Djorgowski et al. (1986)): during its dynamical evolution in the galaxy, some stars are tidally extracted from a globular cluster, to compensate this loss the cluster concentrate its core, increasing then the density contrast, triggering sooner or later the Antonov instability.
3. Without any rotation, the collapse (violent or quiet) of an homogeneous set of gravitating particles produces an E0 (i.e. spherical) isotropic equilibrium state. There are two possible ways to obtain a flattened equilibrium:

* Introduce a large amount of inhomogeneity near the center in the initial state, and make a violent collapse ($`\eta <25`$).
    * Introduce a sufficient amount ($`f>4`$) of rotation in the initial state.

These two ways have not the same origin and do not produce the same equilibrium state.
In the first case, one can reasonably invoke the Radial Orbit Instability: as a matter of fact, as it is explained in a lot of works (see Perez et al. (1998) for example) two features are associated to this phenomenon. First of all, it is an instability which needs an equilibrium state from which it grows. Secondly, it triggers only when a sufficient amount of radial orbits are present. The only non rotating flattened systems we observed just combine this two conditions: sufficient amount of radial orbits because the collapse is violent and something from which ROI can grow because we have seen in the previous point that inhomogeneities collapse first and quickly join the center. The fact that cold P$`{}_{}{}^{\alpha }{}_{\eta }{}^{}`$ systems are more flattened than C$`{}_{}{}^{\alpha }{}_{\eta }{}^{}`$ ones is in complete accordance with our interpretation: as a matter of fact, by construction, power law systems have an initial central overdensity, whereas clumpy systems create (quickly but not instantaneously) this overdensity bringing the collapsed clumps near the center. The ROI flattening is oblate ($`a_21`$ and $`a_1<1`$).
The rotational flattening is more natural and occurs when the centrifugal force overcomes the gravitational pressure. The rotational flattening is prolate ( $`a_2>1`$ and $`a_11`$). We notice that initial rotation must be invoked with parsimony to explain the ellipticity of some globular clusters or elliptical galaxies. As a matter of fact, these objects are very weakly rotating systems and our study has shown that the amount of rotation is almost constant during the collapse.
4. Spherical equilibria can be suitably fitted by both isothermal and polytropic laws with various indexes. It suggests that any distribution function of the energy exhibiting an adaptable core halo structure (Polytrope, Isothermal, King, Hernquist,…) can suitably fit the equilibrium produced by the collapse of our initial conditions.
5. There exists a temperature segregation between equilibrium states. It concerns only initially cold systems (i.e. systems which will suffer a violent collapse): for such systems when $`\eta `$ decreases, the equilibrium temperature $`T`$ increases much more for initially homogeneous systems than for initially inhomogeneous systems. On the other hand, whatever their initial homogeneity, quiet collapses are rather all equivalent from the point of view of the equilibrium temperature: $`T`$ increases in the same way for all systems as $`\eta `$ decreases. This feature may be the result of the larger influence of the dynamical friction induced by the primordial core on the rapid particles in a violent collapse.

All these properties may be directly confronted to physical data from globular clusters (see Harris catalogue Harris (1996)) or galaxies observations.
As a matter of fact, in the standard ”bottom-up” scenario of the hierarchical growth of structures, galaxies naturally form from very inhomogeneous medium. Our study then suggests for the equilibrium state of such objects a potential flattening and a collapsed core. This is in very good accordance with the E0 to E7 observed flatness of elliptical galaxies and may be a good explanation for the presence of massive black hole in the center of galaxies (see Schodel et al. (2002)).
On the other hand Globular Clusters observations show that these are spherical objects (the few low flattened clusters all possess a low amount of rotation), and that their core is generally not collapsed (the collapsed core of almost 10% of the galactic Globular Clusters can be explained by their dynamical evolution through the galaxy). Our study then expect that Globular Clusters form from homogeneous media.
These conclusions can be tested using the $`ET`$ plane. As a matter of fact, we expect that an $`ET`$ plane build from galactic data would not present any High Branch whereas the same plane build from Globular Clusters data would.

## Acknowledgments

We thank the referee for the relevance of his remarks and suggestions. We thank Joshua E. Barnes, who wrote the original Treecode. We also particularly thank Daniel Pfenniger for the use of the parallel Treecode. The simulations were done on the Beowulf cluster at the Laboratoire de Mathématiques Appliquées from the École Nationale Supérieure de Techniques Avancées.
This is a preprint of an article accepted for publication in Monthly Notices of The Royal Astronomical Society ©2003 The Royal Astronomical Society.

## Appendix A Tables of results