text:
# Universal persistence exponents in an extremally driven system

## Abstract

The local persistence $`R(t)`$, defined as the proportion of the system still in its initial state at time $`t`$, is measured for the Bak–Sneppen model. For 1 and 2 dimensions, it is found that the decay of $`R(t)`$ depends on one of two classes of initial configuration. For a subcritical initial state, $`R(t)t^\theta `$, where the persistence exponent $`\theta `$ can be expressed in terms of a known universal exponent. Hence $`\theta `$ is universal. Conversely, starting from a supercritical state, $`R(t)`$ decays by the anomalous form $`1R(t)t^{\tau _{\mathrm{ALL}}}`$ until a finite time $`t_0`$, where $`\tau _{\mathrm{ALL}}`$ is also a known exponent. Finally, for the high dimensional model $`R(t)`$ decays exponentially with a non–universal decay constant.

PACS numbers: 05.70.Ln, 64.60.Ht

Extremally driven systems represent an important class of non–equilibrium models, with applications including sheared granular media, biological macroevolution and interface depinning . Their defining characteristic is that they are driven in the vicinity of an extremal quantity, often the minimum or maximum of an inhomogeneous scalar field. Such deterministic dynamics can be naturally realised in the limit of vanishing (thermal) noise . Moreover, many extremal models are also ‘critical’ in an analogous sense to a continuous phase transition in equilibrium systems .

However, the temporal evolution of these systems is in many cases not fully understood. This is clearly demonstrated by recent work on one of the simplest, and therefore most studied, extremally driven systems, namely the Bak–Sneppen model . This has been shown to age in a similar manner to glassy systems , suggesting that it does not reach a steady state in the timeframe of current simulation methods. This is in contrast to the usual claim that stationarity is reached ‘after an extensive transient’ . It now appears that this contradiction arises because most studies only measure functions of a single time variable, such as critical exponents etc., that can appear to become constant after a finite time even for a non–stationary system; a true test of stationarity (such as the aging result mentioned above) require the measurement of two–time correlation functions .

The aim of this paper is to investigate one aspect of the temporal behaviour of the Bak–Sneppen that has not yet been studied, in the hope that this will further elucidate the time evolution of extremally driven systems in general. The quantity that we consider is the local persistence $`R(t)`$, defined as the proportion of the system at a time $`t`$ that has not yet changed from its initial state. This quantity has been the focus of considerable attention in recent years, at least in part because it is often found to decay algebraically with a non–trivial persistence exponent $`\theta `$, $`R(t)(\mathrm{const})t^\theta `$, that (for equilibrium systems) cannot be related to either the static or dynamic exponents . The universality of $`\theta `$ remains an open subject. In many cases it is not universal; that is, it may depend on the lattice coordination number or the precise choice of interaction term . However, it has been found to be universal in directed percolation . We find this to also be the case for the Bak–Sneppen model, and confirm our numerical finding by deriving an expression for $`\theta `$ in terms of the known universal exponent $`\gamma `$ . We further find that, when starting from a supercritical initial configuration, $`R(t)`$ decays in an anomalous, non–algebraic manner. Finally, $`R(t)`$ always decays exponentially the infinite dimensional case.

The model: The Bak–Sneppen model is defined as follows . A scalar quantity $`f_i[0,1]`$ is independently assigned to every site $`i`$ of a $`d`$–dimensional hypercubic lattice, which contains a total of $`N`$ sites. For every unit time $`t`$, the smallest $`f_i`$ in the system is found, and it and its nearest neighbours are assigned new values, again drawn uniformly from $`[0,1]`$. Other initial conditions and interaction terms will be considered below. To generate the random values of the $`f_i`$ , we used L’Ecuyer’s 64 bit combined multiple recursive pseudo–random number generator MRG63k3a, which has a quoted period of $`2^{377}`$ . The typical number of operations required to find the global minimum was reduced from $`𝒪(N)`$ to $`𝒪(\mathrm{ln}N)`$ by using a binary search tree, which also constrains the system size to be a power of 2 .

Typical results for a periodic 1D lattice are presented in Fig. 1. It can clearly be seen that the curves of $`R(t)`$ for different system sizes collapse after rescaling the time scale by $`N`$, for all $`N>2^{10}10^3`$. For the largest system $`N=2^{21}`$, the simulation time was long and only one run was performed; however, for smaller systems $`R(t)`$ was averaged over 10-100 different initial configurations. Also plotted in the figure is the theoretical prediction $`R(t)t^{1.59}`$, which will be derived below. The numerical $`R(t)`$ appear to be converging towards this prediction at late times, but there is still a slight curvature to the lines (on log–log axes) even for the largest $`t`$ we were able to simulate, making direct measurement of the exponent difficult.

A similar picture emerges in two dimensions, as shown in Fig. 2. The theoretical prediction in this case $`R(t)t^{2.43}`$ (again, see below for details), which is consistent with the simulation data, although as before there is still some slight curvature to $`R(t)`$ for the largest $`t`$ available to our current computing resources.

The universality of the exponent was tested by altering the microscopic details of the interactions. Three alternative rules were considered : updating the second nearest neighbours as well as the nearest neighbours (“2nd NN’s”); changing the value of the nearest neighbours from $`f_i`$ to $`\frac{1}{2}(f_i+r)`$, where the annealed random variable $`r`$ is uniformly distributed on  (“Half NN”); and setting the minimum value to 1 but changing the nearest neighbours as normal (“Min to 1”). In all cases $`\theta `$ was found to be approximately the same. This has been tested in both one and two dimensions; the 1D case is shown in Fig. 3.

Theory: We shall now explain how the persistence exponent can be theoretically derived. This argument utilises known results for the Bak–Sneppen model, which for reasons of space will not be repeated here; the reader is instead referred to  for their justification. It is known that the location of the minimum or ‘active’ site initially jumps around the lattice in an uncorrelated manner, but as the system evolves the lattice variables $`f_i`$ become spatially correlated, and the active site tends to remain in a localised region of the lattice for a finite amount of time before jumping to another, uncorrelated region of the lattice. Each such period of localised activity is known as an ‘avalanche.’

Let the mean duration of an avalanche at time $`t`$ be denoted by $`S(t)`$, and the mean number of sites changed at least once during the avalanche by $`n_{\mathrm{cov}}(t)`$, where the angled brackets ‘$`\mathrm{}`$’ refer to averages over different initial configurations. By taking a suitably large system size $`N`$, these quantities will increase by an arbitrarily small amount between avalanches, allowing the continuum limit to be taken. It is known that both quantities increase algebraically in time ,

$`S(t)`$ $``$ $`A\left({\displaystyle \frac{t}{N}}\right)^{\gamma /(\gamma 1)},`$ (1)
$`n_{\mathrm{cov}}(t)`$ $``$ $`B\left({\displaystyle \frac{t}{N}}\right)^{1/(\gamma 1)}.`$ (2)
Here, $`\gamma `$ is a universal exponent that depends on the lattice dimensionality and the symmetries of the interactions. We assume $`\gamma >1`$; The case of $`\gamma =1`$, realised in high dimensions, will be discussed later.

By defintion, $`R(t)=1`$ at $`t=0`$ and decreases by an amount $`\delta R=1/N`$ whenever a lattice site changes its value for the first time. The average rate at which sites are changed by the avalanches is $`n_{\mathrm{cov}}(t)/S(t)`$; however, only a fraction $`R(t)`$ of sites have not already changed from their initial state. Recalling that the active site jumps to a random part of the system inbetween avalanches, then $`R(t)`$ obeys the following equation

$$\frac{\mathrm{d}R(t)}{\mathrm{d}t}=\frac{n_{\mathrm{cov}}(t)R(t)}{S(t)N}.$$
(3)
Subsituting (1) and (2) into (3) gives

$$R(t)(\mathrm{const})t^\theta $$
(4)
with the persistence exponent $`\theta =B/A`$.

Since $`A`$ and $`B`$ are non–universal constant prefactors to the algebraic growth of $`S(t)`$ and $`n_{\mathrm{cov}}(t)`$ respectively, one would expect that their ratio, and therefore $`\theta `$, is also non–universal; however, a peculiar feature of this model is that $`B/A`$ is universal. This can be seen by combining together two exact equations, referred to as the gap and gamma equations in , which gives the following equation relating $`S(t)`$ and $`n_{\mathrm{cov}}(t)`$,

$$\frac{\mathrm{d}S(t)}{\mathrm{d}t}=\frac{n_{\mathrm{cov}}(t)}{N}.$$
(5)
Substituting (1) and (2) into (5) gives the persistence exponent to be

$$\theta =\frac{B}{A}=\frac{\gamma }{\gamma 1}.$$
(6)
Since $`\gamma `$ is universal, then so is $`\theta `$. In 1 and 2 dimensions, $`\gamma =2.70`$ and 1.70 respectively , giving predicted values of $`\theta 1.59`$ and 2.43, which were the values used in Figs. 13. As a further check on this analysis, note that (1) and (2) predict that the quantity $`tn_{\mathrm{cov}}(t)/(NS(t))`$ will approach the constant value $`B/A=\theta `$ at late times. This is confirmed by the numerical results presented in Fig. 4.

The preceding argument breaks down when $`\gamma =1`$, which arises for dimensions above the upper critical dimension $`d_\mathrm{c}`$ , for which the most recent estimate is $`d_\mathrm{c}=4`$ . In this case, $`S(t)`$ and $`n_{\mathrm{cov}}(t)`$ both increase exponentially in time rather than algebraically . By repeating the same steps as above, it is straightforward to show that $`R(t)`$ should now decay exponentially with a non–universal decay constant $`C`$,

$$R(t)(\mathrm{const})\mathrm{e}^{Ct}$$
(7)
This is confirmed by the numerical results given in Fig. 5 for a system in which the minimum and another $`K1`$ randomly chosen elements are assigned new values, which is the usual way of simulating an infinite dimensional system for this model. It is clear that the decay of $`R(t)`$ is exponential with a decay constant that depends strongly on $`K`$.

Finally, we consider the effects of initial conditions. All of the preceding analysis and numerical results has assumed an initial configuration in which every $`f_i`$ is uniformly distributed over the range . This can be generalised by choosing an initial state in which the $`f_i`$ are uniformly distributed over the range $`[f_0,1]`$ and are spatially uncorrelated. It is known that the Bak–Sneppen model has a ‘critical value’ $`f_\mathrm{c}(0,1)`$ such that the system is subcritical if $`f_0<f_\mathrm{c}`$, and supercritical if $`f_0f_\mathrm{c}`$ . (The value of $`f_\mathrm{c}`$ is non–universal and depends on the lattice structure and the microscopic interaction term). We have checked that any subcritical initial state gives rise to the same persistence exponent in both one and two dimensions, as predicted by our earlier argument. However, for a supercritical initial configuration, the system enters into a single, infinite avalanche, and all of the preceding analysis fails. In fact, $`R(t)`$ in this situation is simply related to the growth of a single avalanche in time, which is already known to increase algebraically in time with a universal exponent $`\tau _{\mathrm{ALL}}`$ . Thus $`R(t)`$ does not decay algebraically but rather according to the expression

$$1R(t)(\mathrm{const})t^{\tau _{\mathrm{ALL}}}$$
(8)
until a late time $`t_0`$ when the avalanche ‘touches’ itself through the periodic boundary conditions. Confirmation of this scaling is given in Fig. 6 for one and two dimensions. In high dimensions, $`R(t)`$ decays exponentially for both classes of initial state.

Summary: We have found that the persistence $`R(t)`$ in the Bak–Sneppen model can decay in one of three different ways. In low dimensions, it decays with a universal exponent $`\theta `$ when starting from a subcritical state, but by the anomalous form $`1R(t)t^{\tau _{\mathrm{ALL}}}`$ when starting from a supercritical state. For high dimensions the decay is always exponential. We also provided a simple theoretical prediction for the value of $`\theta =\gamma /(\gamma 1)`$, which, when combined with the numerical results, gives strong confirmation of the universality of $`\theta `$.

It is not yet clear how this relates to the question of the existence of a steady state of the model, although it is interesting to note that $`\theta >1`$ always, which is significantly higher than typical persistence exponents, and means that the number of untouched sites will rapidly decay to zero. Thus the observed non–stationarity must be due to temporal correlations within an avalanche, rather than the spread of avalanches throughout the system.

Acknowledgements: The author would like to thank the hospitality of the Vrije Universiteit, Amsterdam, where part of this work was carried out. This work was funded by EPSRC (UK) grant no. GR/M09674.


format: text

text:
# Evidences of Bolgiano scaling in 3D Rayleigh-Bénard convection

\[

## Abstract

We present new results from high-resolution high-statistics direct numerical simulations of a tri-dimensional convective cell. We test the fundamental physical picture of the presence of both a Bolgiano-like and a Kolmogorov-like regime. We find that the dimensional predictions for these two distinct regimes (characterized respectively by an active and passive role of the temperature field) are consistent with our measurements.

\]

The dimensional theory of homogeneous and isotropic turbulence has been definitely settled long ago by the work of A. N. Kolmogorov . On the other hand, it is still missing a clear theoretical picture for the strong fluctuations in the energy dissipation field that lead to intermittency effects (i.e. non gaussian behaviour of probability distribution functions). Phenomenological theories have been proposed  but no systematic theory for computing experimentally measured numbers has been successful so far. The situation of “non-ideal” turbulence is even more controversial, already at the level of dimensional expectations. A typical realization, the one we will address in this paper, is the tri-dimensional (3D) Rayleigh-Bénard cell, described, in the Boussinesq approximation , by the following set of equations:

$`_t𝒗`$ $`+`$ $`\left(𝒗\right)𝒗=p+\nu ^2𝒗+\alpha gT\widehat{𝒛}`$ (1)
$`_tT`$ $`+`$ $`\left(𝒗\right)T=\chi ^2T`$ (2)
with isothermal boundary conditions on the upper and lower planes of a cell of height H: $`T(z=0)=+\mathrm{\Delta }T/2`$ and $`T(z=H)=\mathrm{\Delta }T/2`$. As usual, $`𝒗(𝒙,t)`$ is the velocity field and $`T(𝒙,t)`$ the temperature field. Kinematic viscosity and thermal conductivity are respectively $`\nu `$ and $`\chi `$, while the thermal expansion coefficient is $`\alpha `$ and gravity acceleration is $`g`$. In the following we will mainly focus on the longitudinal structure functions of $`v`$ and $`T`$: $`S_p(r)=\left[\left(𝒗(𝒙+𝒓)𝒗(𝒙)\right)\widehat{𝒓}\right]^p`$ and $`T_p(r)=\left[T(𝒙+𝒓)T(𝒙)\right]^p`$.

In this work we present some tests of the predictions for the structure functions defined above that can be derived in the scenario proposed years ago by Bolgiano  to describe convective turbulence.

Despite much research on the subject, sound evidence of the validity of the Bolgiano scenario and the recovery of Kolmogorov scaling at small scales is still missing. Good quality confirmation of the dimensional Bolgiano scenario was recently shown in a two dimensional numerical simulations . This result cannot be directly related to the 3D case because of the strong differences in the properties of the velocity field in 2D. Statistical properties of the velocity field were recently used with even simpler models (shell models for turbulence) . This approach is even less probing since these models were built precisely in order to implement Bolgiano scaling. Finally, a recent experiment even questioned the behaviour the the Bolgiano length inside a convective cell . Of course from an experimental point of view it can be difficult, if not impossible at all, to measure all relevant quantities and it can be even harder to have access to information at several (not just a few) positions inside the volume. The results presented in this letter show good consistency with the idea of a Bolgiano regime (i.e. a range of length scales where temperature driven buoyancy effects are dominant) but, because of the still limited resolution of “state of the art” numerical simulations, we will have to resort to somehow indirect tests.

This letter is organized as follows: a brief review of phenomenological expectations, details of our numerical simulations, data analysis and then concluding remarks.

Starting from equations (1) and (2), if one uses dimensional analysis and assumes homogeneous scaling for velocity and temperature differences (inside the inertial range), one ends up with two distinct scaling regimes. At small scales (Kolmogorov-like scenario), $`rL_B`$,

$`\delta v(r)`$ $``$ $`\epsilon ^{1/3}r^{1/3}`$ (3)
$`\delta T(r)`$ $``$ $`N^{1/2}\epsilon ^{1/6}r^{1/3}`$ (4)
while at large scales (Bolgiano-like scenario), $`rL_B`$,

$`\delta v(r)`$ $``$ $`\left(\alpha g\right)^{2/5}N^{1/5}r^{3/5}`$ (5)
$`\delta T(r)`$ $``$ $`\left(\alpha g\right)^{1/5}N^{2/5}r^{1/5}`$ (6)
The Bolgiano length, $`L_B`$, is an estimate of the distance at which the dissipative and buoyancy terms on the right hand side of equation (1) balance, valid under the assumption of a scaling behaviour. In the following we will use the $`z`$ dependent version of $`L_B(z)`$ introduced in  that uses averages of $`\epsilon `$ and $`N`$ defined at a given height $`z`$, $`\epsilon (z)=\nu /2_{ij}\left(_iv_j\right)^2_z`$ and $`N(z)=\chi /2_i\left(_iT\right)^2_z`$:

$$L_B(z)=\frac{\epsilon (z)^{5/4}}{N(z)^{3/4}(\alpha g)^{3/2}}$$
(7)
Our analysis has been performed on data coming from Direct Numerical Simulations (DNS) employing a by now standard Lattice Boltzmann scheme , on a massively parallel computer. The resolution of the numerical simulation was $`240^3`$ and the Rayleigh number ($`Ra=\alpha g\mathrm{\Delta }TH^3/\left(\nu \chi \right)`$) was approximately $`3.510^7`$. The Prandtl number was equal to unity. We performed a stationary simulation extending over approximately $`500`$ recirculation times and stored nearly $`400`$ independent configurations with full information on all velocity components and the temperature field. Boundary conditions were periodic in the $`x`$ and $`y`$ directions (in order to maintain homogeneity on horizontal planes) and isothermal at the top and bottom planes of the cell ($`z=0`$ and $`z=H`$).

Free slip boundary conditions (i.e. satisfying only the incompenetrability constraint) were used on the top/bottom planes for the velocity field. This choice was made in an attempt to reduce the effects of a viscous boundary layer close to the horizontal walls. Indeed, as will be clear from the following, thermal effects are dominant near the isothermal walls, so using no-slip boundary conditions might produce effects on the velocity statistics interfering with the ones coming from pure buoyancy. More details on this point will be given later on.

The most direct way to test the dimensional validity of the Bolgiano picture would be to measure structure functions and to check whether they scale with exponents close (apart from intermittency corrections) to the ones predicted by the sets of equations (3-4) and (5-6).

Unfortunately this cannot be done directly because of the limited resolution of our DNS.

Fig. 1 substantiates this comment, by plotting the velocity structure function of order 3, $`S_3(z,r)=\left(v(x+r,z)v(x,z)\right)^3`$ at two different position: one very close to the wall ($`z=10`$), and the other at the center of the cell ($`z=120=H/2`$).

As it can be seen from Fig. 1, no evident scaling range can be detected, even if a steepening of $`S_3(r)`$ is clearly detected as one comes close to the walls.

Given this state of affair, in order to test the consistency with the Bolgiano picture we performed two distinct, even if somehow less direct, tests.
The first test consists in checking that the Bolgiano length actually keeps track of the scales at which the buoyancy term balances the dissipative term in equation (1).

First, we measure $`L_B(z)`$ in terms of $`\epsilon (z)`$ and $`N(z)`$ and plot our results in Figure 2. From a first look at the behaviour of $`L_B(z)`$ we learn that Bolgiano effects should be measurable, if at all, near to the isothermal walls (where $`L_B(z)`$ is of the order of $`10^1`$).

Close to the center of the cell, $`L_B`$ is of the order of $`10^2`$ and Kolmogorov like behaviour is expected to be measurable at almost all scales.

We then measure, directly and independently from the previous quantity, the scale at which dissipation and buoyancy effects balance, i.e. we look for the scale $`\stackrel{~}{L}_B(z)`$ for which:

$$\epsilon (z)\alpha g\delta v(\stackrel{~}{L}_B)\delta T(\stackrel{~}{L}_B)_z$$
(8)
We provisionally consider $`\stackrel{~}{L}_B`$ a modified definition of the Bolgiano length. In Fig. 2 we plot $`L_B(z)`$ and $`\alpha \stackrel{~}{L}_B(z)`$ with a constant which was tuned to be $`\alpha 3.1`$.

As it can be seen, the two definitions yield the same behaviour apart from the multiplicative factor, $`\alpha `$. The reason for the multiplicative factor (of order unity) is due to the fact that the two definitions are dimensional estimates so they can miss a numerical prefactor. Considering this, the fact that the two definition behave in the same way after rescaling has to be considered, in our opinion, an excellent agreement.

Here we want to underline two points which, we believe, add relevance to this finding. First of all the two definitions are of course linked but definitely different. The “traditional” definition of $`L_B`$ (see also ), as from definition (7), comes from supposing that the two scaling laws in equations (3) and (5) merge at $`L_B`$, hence comes from solving the equation $`\epsilon (z)^{1/3}L_B(z)^{1/3}=\left(\alpha g\right)^{2/5}N(z)^{1/5}L_B(z)^{3/5}`$. The second definition, $`\stackrel{~}{L}_B(z)`$ comes instead from a direct measurement of the strength of the dissipative and forcing term in equation (1).

The second important point consist on the fact that our cell is not homogeneous; this adds strength to the equality between the two different definitions of $`L_B(z)`$.

As a consequence of this test we can claim that the Bolgiano scenario and the expected power law behaviour are consistent with the value measured for the terms appearing on the right hand side of equation (1). Furthermore we fully confirm (with higher accuracy) our former results for $`L_B`$ . We like to underline that extracting the behaviour of the Bolgiano length using measured quantities that are not those appearing in its definition could introduce large errors. For example in  a Bolgiano length was extrapolated as the scale at which there is a change of slope in a particular structure function. In using such a procedure one is dominated by strong finite size effects present on the structure functions.

We now proceed to our second test. Under the hypothesis of validity of Bolgiano scaling, and if enough resolution were available, one would expect to see power law behaviour in the inertial range in Fig. 1 with two distinct slopes $`1`$ and $`9/5`$, corresponding respectively to Kolmogorov and Bolgiano scaling. With available resolution, we are not able to detect a clear power law behaviour from Fig. 1, although we see a clear steepening of the structure function as we come close to the wall. What we are going to do in the following is to try to quantify as well as possible this change of slope.

We adopt the following procedure. We focus on the plot of $`S_6(r)`$ versus $`S_3(r)`$ and apply Extended Self Similarity (ESS, see ) in order to detect a trustable plateau in the local slopes. We found this plateau to correspond roughly to distances in the interval $`^v=[25,40]`$ for the velocity structure function and $`^T=[15,30]`$ for the structure functions of the temperature field. We then define two other interval slightly shifted to the left $`_{}^v=[20,35]`$, $`_{}^T=[10,25]`$ and to the right $`_+^v=[30,45]`$, $`_+^T=[20,40]`$. These interval were of course shifted by a reasonable amount, i.e. were still possible highest or lowest estimate for the same plateau.

We finally perform a power law fit to extract a scaling exponent on the structure functions for the velocity and for the temperature in the three intervals defined before.

We define the scaling exponents for the structure functions of interest as follows:

$`\delta v(z,r)^p`$ $``$ $`r^{\zeta _p(z)}`$ (9)
$`\delta T(z,r)^p`$ $``$ $`r^{\chi _p}(z)`$ (10)
$`\left(\delta v(z,r)\delta T(z,r)^2\right)^{\frac{p}{3}}`$ $``$ $`r^{\rho _p(z)}`$ (11)
Our fits provide a central value for the exponents and two, respectively, higher and lower estimates (corresponding to the shifted intervals). The procedure adopted in estimating the errors reflects the fact that the largest source of systematic error is connected with the choice of the fitting ranges and not with statistical accuracy.

In the Kolmogorov regime we expect the exponents of equations (9-11) to take the following dimensional values: $`\zeta _p=p/3`$, $`\chi _p=p/3`$ and $`\rho _p=p/3`$. In the Bolgiano regime, on the other hand, we expect the following dimensional values: $`\zeta _p=3/5p`$, $`\chi _p=p/5`$ and $`\rho _p=p/3`$.

In Figure 3 we plot the behaviour of the measured central exponents and their errors as a function of $`z`$ for $`\zeta _3(z)`$ and $`\chi _3(z)`$. The behaviour is qualitatively and quantitatively consistent with the expected scaling exponents: we observe a smooth transition from a Bolgiano dominated regime (near to the wall, small $`z`$) to a Kolmogorov regime (approaching the center of the cell).

In Figure 3 we have not shown the behaviour of $`\rho _3(z)`$ as it is consistent (within error bars) with the constant value $`1`$ and plotting it on the same figure would have made it unreadable.

Here we, once again, like to underline the consistency of our findings. From the theoretical picture we expect to see the crossover to the Kolmogorov scaling when $`L_B(z)H`$. This is indeed the case as $`L_B(z)H`$ for $`z40`$ and indeed we see $`\zeta _3(z)1`$ in the same range of $`z`$.

Another interesting question concerns the “real” statistics of the Bolgiano regime. Very recently interest is growing along this line of research because of the desire to understand the differences of the statistical properties of an active with respect to a passive scalar . In this work we focused on the gross features, i.e. on the dimensional behaviour. If we try to look at intermittency, by means of ESS we find a strong increase of intermittency for the velocity field approaching the isothermal walls.

Unfortunately we believe that with our simulation we are not in a position to make any definite statements about intermittency in the Bolgiano dominated regime.

Indeed a study of intermittency in the Bolgiano regime would involves positions nearby the isothermal walls (as only there the Bolgiano length scale is small enough to have an inertial range largely dominated by buoyancy effects). Recently it was found that intermittency increases in the velocity structure functions inside a viscous boundary layer . In an attempt to reduce the viscous boundary layer we decided to apply free-slip velocity boundary conditions to the isothermal walls free slip for the velocity. Unfortunately our idea cannot be successful because of two reasons. First, the free-slip boundary conditions does not completely suppress the boundary layer. A boundary layer thickness can be defined through an extrapolation of $`v_z(x)^2`$, and in our case turned out to be of the order of roughly $`15`$ grid spacings. Second, it was recently realized that mechanical turbulence in a free-slip channel also presents an enhancement of intermittency near to the walls .

A procedure to disentangle buoyancy and planar effect is clearly needed to make any definite statement on intermittency in the Bolgiano regime. Otherwise one cannot decide whether a change of intermittency is related to Bolgiano dynamics or to the increase occurring nearby boundary layers.

In order to clarify this point it would be important to perform a simulation of a Rayleigh-Bénard-like system with periodic boundary condition in all directions (i.e. an homogeneous Rayleigh-Bénard cell). We suggest to extend to 3D the study made in 2 dimensions in .

In order to do that one could write temperature field as the sum of a linear profile plus a fluctuating part, $`T(x,y,z)=T_{\text{lin}}(z)+T^{}(x,y,z)`$, with $`T_{\text{lin}}=\mathrm{\Delta }T/2(12z/H)`$ and obtain:

$`_t𝒗`$ $`+`$ $`\left(𝒗\right)𝒗=p+\nu ^2𝒗+\alpha gT\widehat{𝒛}`$ (12)
$`_tT^{}`$ $`+`$ $`\left(𝒗\right)T^{}=\chi ^2T^{}{\displaystyle \frac{\mathrm{\Delta }T}{H}}v_z`$ (13)
Choosing the parameters in order to have a Bolgiano length as small as possible, one would benefit of a wide range of scales were to study the buoyancy dominated flow. Further advantage of homogeneity would be the natural increase of statistics and also the applicability of tools like $`SO(3)`$ decomposition, to disentangle anisotropic terms . A study of this kind is in progress.

Concluding we have performed a number of basic tests in order to validate the scenario of Bolgiano-Kolmogorov scaling in a convective cell, within the limitations but also advantage of $`z`$ non homogeneity of our cell. We were able to confirm the transition between the two expected scenarios.

We acknowledge useful discussions with R. Benzi, S. Succi and R. Verzicco. All numerical simulations were performed on the APEmille computer at INFN, Sezione di Pisa.


text:
# On the detectability of distant Compton-thick obscured quasars

## 1 Introduction

Most of the 0.5–2 keV soft X-ray Background (XRB) has been resolved by ROSAT (Hasinger et al 1999; Lehmann et al 2001) and recently most of the harder 2–7 keV XRB has been resolved by Chandra (Mushotzky et al 2000; Giacconi et al 2001; Barger et al 2001; Brandt et al 2001) and XMM (Hasinger et al 2001). The source identifications are incomplete at the present time but indicate that about one third of the sources are associated with normal quasars, another third with optically-bright galaxies and the last third with optically-faint galaxies (Mushotzky et al 2000; Barger et al 2000; Alexander et al 2001). X-ray hardness ratios indicate that the fainter sources are generally harder and consistent with absorbed sources. This confirms obscuration models for the XRB (Setti & Woltjer 1989; Madau, Ghisellini & Fabian 1994; Comastri et al 1995; Gilli, Salvati & Hasinger 2001).

Detailed studies, including X-ray spectroscopy of the brighter sources (e.g. Crawford et al 2001) have however shown that the column densities of absorbing material lie in the range of $`10^{22}`$ to a few times $`10^{23}\mathrm{cm}^2`$. Such columns are insufficient to produce the $`\nu I_\nu `$ peak in the XRB at about 30 keV (Marshall et al 1982). In order to match this peak, sources with columns exceeding $`1.5\times 10^{24}\mathrm{cm}^2`$ are required, i.e. sources which are Compton thick.

Here we address this issue and ask whether such sources should have been observed in Chandra and XMM images, and whether their apparent absence indicates a problem for the obscuration model of the XRB. Compton-thick sources are certainly found in the local Universe, with two of the nearest active galactic nuclei (AGN), the Circinus galaxy and NGC4945 (Matt et al 2000, and references therein) being in this class. Good examples of powerful Compton-thick AGN are NGC6240 (Vignati et al 1999) which is at low redshift and IRAS 09104+4109 which is at redshift $`z=0.442`$ (Franceschini et al 2000; Iwasawa et al 2001). The main issue is whether Compton-thick quasars are common.

Sources which are mildly Compton-thick ($`10^{24}<N_\mathrm{H}<10^{25}\mathrm{cm}^2`$) show direct emission peaking above 10 keV. When they are heavily Compton thick ($`N_\mathrm{H}>10^{25}\mathrm{cm}^2`$, this emission peak is down-scattered by Compton recoil and is therefore highly degraded. How much is seen in the 2–8 keV band from low redshift Compton-thick objects depends almost completely on how much emission is scattered into our line of sight. For example, the direct line of sight to the well-known obscured AGN NGC 1068 has a column density exceeding $`10^{25}\mathrm{cm}^2`$ (Matt et al 1997), the X-ray emission seen is all scattered flux. The $`z=0.92`$ hyperluminous IRAS galaxy F15307, known to host a quasar from observed scattered (polarised) optical broad lines, shows no X-ray flux detectable by ASCA and ROSAT (Ogasaka et al 1997; Fabian et al 1996). The difference between these sources is probably determined by the geometry of the obscuring medium, which may be torus-like in the case of NGC 1068 and more spherical for F15307. Since synthesis models for the XRB spectrum indicate that most (85 per cent) of the AGN power is absorbed (Fabian & Iwasawa 1999), a spherical geometry is more likely to be relevant to the sources dominating the hard XRB than a torus one. The scattered flux fraction is likely to be small. It is the absorbed direct flux which we are interested in here.

Redshift does help to shift so render mildly Compton-thick quasars detectable, leading to an inverse K-correction effect (Wilman & Fabian 1999). However the sources still need to be luminous to be detectable. For the spectrum of the XRB, it is the mildly Compton-thick objects which dominate at the $`\nu I_\nu `$ XRB peak. The lowest redshift members of the Compton-thick class contribute most to the highest energy part of that peak.

To assess whether Compton-thick quasars ought to be detected, we have taken one simple obscured AGN model for the XRB (from Wilman, Fabian & Nulsen 2000; WFN), and predicted the number of Compton-thick sources expected per square degree. We have then passed the predicted spectra through the response curves of Chandra and XMM to predict the number expected to be detected. Our results indicate that some Compton-thick sources should be detectable, particularly in the Chandra 1 Ms images (Alexander et al 2001; Rosati et al 2001). Chandra does not have the collecting area above 8 keV to detect such hard sources readily; XMM has too much internal background (relevant because of its much larger PSF) at those energies. There may yet be renewed claims to resolve the XRB when distant Compton-thick sources are finally found in reasonable numbers. We argue that our results are not strongly dependent on the model adopted.

## 2 The WFN model

For the sake of completeness, we provide here an outline of the WFN model. At its heart is the semi-analytic galaxy formation code developed by Nulsen & Fabian (1995), in which the Cole & Kaiser (1988) block model is used to simulate the hierachical growth of clustering, and the scheme of Nulsen & Fabian (1995) is used to treat the behaviour of the gas within a collapsed halo. In short, gas within a radius $`R=R_{\mathrm{CF}}`$, where the free-fall time is less than the cooling time, rapidly forms stars, producing supernovae which can then expel some (for normal galaxies) or all (in the case of dwarf galaxies) of the remaining gas from the system. Any gas at $`R>R_{\mathrm{CF}}`$ participates in a cooling flow (CF). Newly formed normal galaxies are ellipticals unless all the hot gas is able to cool before the present or the next hierarchical collapse, in which case they are spirals. Any collapse with at most one infalling normal galaxy forms a normal galaxy, with the stars of dwarf galaxies contributing to the spheroid of the new system. Collapses with more than one normal galaxy form a group or a cluster, as mergers between normal galaxies are ignored.

Each of the smallest block model units contains a seed black hole of mass $`1.6\times 10^6`$$`\mathrm{M}_{}`$, and when a block collapses all the black holes associated with its merging subblocks are assumed to merge into a single black hole at the centre of the new galaxy. Such nuclear black holes are then fed by Bondi accretion of hot gas from the cooling flow atmospheres which surround them, in accordance with the model of Nulsen & Fabian (2000). Whilst thus accreting, the black hole is assumed to radiate with an efficiency of 10 per cent, with a 2-10 keV luminosity equal to 3 per cent of the bolometric one, with a power law plus reflection spectrum. This intrinsic spectrum is then absorbed by the isothermal distribution of cold dusty clouds deposited by the cooling flow, as suggested by Fabian (1999). The accretion is terminated by wind-driven gas expulsion when the mass of the black hole reaches a critical fraction of that in the surrounding spheroid, thereby accounting for the observed correlation between the mass of the remnant black hole and its host spheroid. Thereafter the object shines as an optical quasar for $`9\times 10^7`$ yr. For full details of the model, see WFN.

### 2.1 Wind or radiation pressure?

The critical mass of the black hole was determined by invoking a wind from the central engine (Fabian 1999; see also Silk & Rees 1998 who use an energy argument rather than the force one used here). We note here that the radiation pressure of the absorbed radiation alone yields a similar limit (i.e. if the force due to the radiation $`L_{\mathrm{abs}}/c`$, where the absorbed power $`L_{\mathrm{abs}}L_{\mathrm{Edd}}`$).

To pursue this in detail, we adopt the scenario of Fabian (1999), where a fraction $`f`$ of the mass in the core of an isothermal galaxy $`M(<r)=2v^2r/G`$ is in cold, X-ray absorbing gas. The gas within radius $`r_1`$ has been accreted into a black hole of mass $`M_{\mathrm{BH}}=f2v^2r_1/G,`$ leaving a column density $`N=fv^2/2\pi Gm_\mathrm{p}r_1`$ beyond. Rearranging these formulae means that

$$M_{\mathrm{BH}}=f^2v^4/G^2\pi m_\mathrm{p}N.$$

Much of the radiation from the black hole is absorbed by the gas beyond $`r_1`$ giving rise to a force on the column of $`f_1L_{\mathrm{abs}}/4\pi r^2c,`$ where $`f_1`$ accounts for how much the matter traps and reradiates the energy. $`f_10.5`$ for an isolated cloud but may be several for the envisaged optically-thick cloud surrounding the nucleus. Then at the limit where the outward force due to radiation balances the inward gravitational force on the column,

$$L_{\mathrm{abs}}=4\pi GMNm_\mathrm{p}cf_1^1.$$

(The formula assumes that the column is all at one radius; if it extends to $`r_{\mathrm{max}}`$ then a factor of $`\mathrm{log}(r_{\mathrm{max}}/r)`$ enters into the inward force expression.)

Let $`L_{\mathrm{abs}}=f_2L_{\mathrm{Edd}}=f_24\pi GMm_\mathrm{p}c/\sigma _\mathrm{T},`$ where $`\sigma _\mathrm{T}`$ is the Thomson cross section. Then

$$N=f_1f_2/\sigma _\mathrm{T}.$$

Substituting this into the above black hole mass formula then gives

$$M_{\mathrm{BH}}=\frac{v^4\sigma _\mathrm{T}}{\pi G^2m_\mathrm{p}}\frac{f^2}{f_1f_2},$$

or

$$M_{\mathrm{BH}}=2.3\times 10^9\frac{f^2}{f_1f_2}\left(\frac{v}{200\mathrm{km}\mathrm{s}^1}\right)^4.$$

This is the mass of the black hole when the gas (which fuels it) is blown away. The normalization changes slightly, by losing a factor of $`f/\pi ,`$ if $`r_1`$ is redefined as the Bondi accretion radius (as in WFN) for gas at the galaxy virial temperature.

It compares well with the results of Ferrarese & Merritt (2000) and Gebhardt et al (2000) for reasonable values of $`f,`$ $`f_1`$ and $`f_2`$:

$$M_{\mathrm{BH}}=1.4\times 10^8\left(\frac{f}{0.25}\right)^2\left(f_1f_2\right)^1\left(\frac{v}{200\mathrm{km}\mathrm{s}^1}\right)^4\mathrm{M}_{}.$$

The above result demonstrates that a wind is not essential for expelling the gas surrounding a growing black hole.

In the present model it appears more as if the galaxy potential determines the black hole mass, rather than the other way round. Note however that the gravitational binding energy of a galactic bulge, where the velocity dispersion of the bulge is $`300v_{300}\mathrm{km}\mathrm{s}^1,`$ is $`E_{\mathrm{bulge}}2\times 10^6v_{300}^2M_{\mathrm{bulge}}c^2`$. The energy from the central black hole $`E_{\mathrm{AGN}}5\times 10^4M_{\mathrm{bulge}}c^2`$. So only one per cent of that energy can have a major effect on the formation of that bulge and that effect may have occurred when both the black hole and galaxy were young.

## 3 The yield of Compton-thick sources

For the above model with the free parameters set to the values given in section 3 of WFN, we present in Figs. 1–4 results on the yield of Compton-thick quasars in deep exposures with Chandra and XMM. Firstly, Fig. 1 shows as a function of redshift the total areal density of Compton-thick AGN in the model with intrinsic 2–10 keV luminosities in excess of $`10^{44}`$$`\mathrm{erg}\mathrm{s}^1`$. As discussed by WFN, however, such sources are obscured by $`N_\mathrm{H}15N_\mathrm{T}`$ when the obscured phase begins (where $`N_\mathrm{T}=1.5\times 10^{24}`$$`\mathrm{cm}^2`$, the column density above which a source becomes Compton-thick), and by $`N_\mathrm{H}3N_\mathrm{T}`$ at its end. Since there is no scattered flux included in the model spectra and the direct emission is almost completely suppressed for $`N_\mathrm{H}>10N_\mathrm{T}`$ (especially below 10 keV), the vast majority of the sources contributing to Fig. 1 will not be visible to Chandra or XMM. This is borne out by Figs. 2, 3 and 4, described below, which were produced by convolving the object spectra in the observed frame with ‘response functions’ giving the effective area of the telescope and detector combination under consideration.

In Fig. 2 we show the areal density as a function of redshift of both the Compton-thick and Compton-thin AGN (the latter also includes the unobscured quasars), which would produce more than 30 counts in the 5–10 keV band in a 220 ks exposure with XMM on the EPIC pn chip. This is repeated in Fig. 3 for a 100 ks exposure. Fig. 4 shows the analogous plot for the 2–8 keV band with the Chandra ACIS S3 chip, with the count threshold reduced to 10 counts owing to the smaller PSF (and thus lower internal background) compared to XMM. The majority of the detected objects contributing to the yields are quite close to the count thresholds. For example, if the XMM threshold is increased to 40 counts, the yield of Compton-thick (Compton-thin plus unobscured) AGN falls over all redshift by a factor 2–3 (1.5–2). Similarly, increasing the Chandra detection threshold to 15 counts, results in a fall by a factor of 2.5 in the yield of Compton-thick objects, whilst the number of Compton-thin and unobscured objects are relatively unchanged.

Most of the predicted Compton-thick sources should become detectable in exposures of 1Ms (Fig. 5), such as are now being done with Chandra (e.g. Alexander et al 2001; Rosati et al 2001). In an $`8.4\times 8.4`$ arcmin2 region and 1 Ms exposure we predict 182 sources (Compton thick and thin) in the 2–8 keV band whereas Alexander et al (2001) find 102. We do not consider that this discrepancy in numbers is serious, since most of the sources would be close to the detection threshold. That many of the sources are seen in both soft and hard energy bands (Alexander et al 2001) is also not a major problem since a few per cent of the primary soft radiation may be scattered into our line of sight.

The important issue of the redshift distribution of the sources is considered in the next Section. where we also discuss the implications of these results for the interpretation of deep Chandra and XMM surveys.

## 4 Discussion

We see that the surface density of Compton-thick sources expected in a 220 ks exposure (equal in length to the initial Chandra observation of the HDF-N and its vicinity by Hornschemeier et al. 2001) with Chandra ACIS chip is about 200. A chip is $`8\times 8`$ arcmin or 0.0178 sq deg which means that 3–4 faint hard sources are expected in such an exposure. If the ACIS-I array is use, with 4 times the area, then the number rises to about 10 (given the degradation of sensitivity off axis: most of the photons detected are above 6 keV). The sources have redshifts of 3–7 and are basically detectable because of the inverse K-correction. Most of the Compton-thick above $`z=3`$ become detectable in deeper exposures of 1 Ms.

In an XMM pn chip exposure of similar duration, concentrating on the 5–10 keV band, we have about 40 sources detected per square degree, mostly in the redshift range of 2–5. The field of view is about 0.2 sq deg (30 arcmin diameter) leading to about 8 sources detected per 220 ks exposure, similar to the yield from Chandra. Source confusion should not be a problem at such high photon energies.

The adopted WFN model thus predicts that the direct absorbed emission from some Compton-thick quasars is detectable by Chandra and XMM with similar yields. The absorbed primary emission from these sources will be very hard and most will be detected close to the threshold for detection (10 counts for Chandra, 30 for XMM), as discussed in section 2.1. Scattered emission may however allow the sources to be detected also in softer energy bands. The appearance of such objects in the optical and near infrared bands depends on the rate of star formation. If this rate is low then they should appear as redshifted early-type bulges, without strong emission lines, and if high they will be brighter and of much later spectral type, probably with emission lines.

So far only one source with a spectrum consistent with being Compton thick has been reported (Norman et al 2001); its redshift $`z=3.8`$. The observed X-ray emission in this case is not the absorbed direct radiation but a reflected component. Only a guess can be made as to the power of the direct component. It is this direct absorbed component of such sources which is important for the spectral peak of the XRB.

Our model predicts that more than half the sources detected in a Ms exposure with Chandra will have redshifts ranging from 2–8. Most of them would have $`I`$ band magnitudes exceeding 24, i.e. are optically faint according to the definition of Alexander et al (2001). The precise magnitudes would depend upon the star formation rate in the host galaxy. If we assume that the obscured black hole growth phase coincides with the major star formation phase in the galaxy, then the star formation rate must be tens $`\mathrm{M}_{}\mathrm{yr}^1`$. Then the $`R`$ (and approximately $`I`$) magnitude of the hosts will be about 24 for objects at $`z3`$ (e.g. Pettini et al 2001). Objects at $`z6`$ will be about 1.5 mag fainter. Alexander et al (2001) have 47 objects fainter than $`I=24`$ (15 with $`I>25.3`$), which is about half the number we predict.

At this stage we consider that our model remains viable, with a factor of two in number densities arrangeable by a small alteration of some of the parameters (e.g. the gas fraction and/or the gas metallicity: note that most of the sources are predicted to be close to the detection threshold, even in 1Ms). Further work on the optically-faint population in the Ms Chandra fields will resolve this issue by comparison with the redshift distribution in Fig. 5.

We now ask how robust our conclusions might be, since they depend on a specific model. First we note that they should be a reasonable (factor of two) approximation to any model which has the obscured phase before the ’normal’, unobscured, quasar phase, since the model has that emission peaking between $`z=12`$ as observed. It also produces a mean local black hole density in agreement with that inferred from observation (Merritt & Ferrarese 2001). The higher redshift ($`z>3`$) Compton-thick sources contribute 10–20 per cent of the intensity of the XRB spectrum in the 5-12 keV band, and so are not of great importance to the spectrum; the lower redshift sources with $`1<z<2`$ contribute far more and dominate the spectrum in the 20–30 keV band. These last objects are completely undetectable with Chandra or XMM, unless they have scattered soft emission. The lowest redshift Compton-thick sources contribute most at the highest energies.

We conclude that distant Compton-thick quasars, which may dominate the XRB intensity over the 5–20 keV band and represent the growth phase of massive black holes, are just detectable with present instruments. A handful may be found in 100 ks exposures at flux levels close to the detection threshold. 50–100 distant Compton-thick sources are predicted to be detectable in Chandra deep fields of 1 Ms or more. The sources should be optically faint and probably appear as high redshift starburst galaxies.

## 5 Acknowledgements

ACF and CSC thank the Royal Society for support.


text:
# Global spiral modes in multi-component disks

## 1 Introduction

Global spiral arms observed in disk galaxies are the result of density waves propagating in the galactic disks (e.g. Binney & Tremaine binney87 (1987)). Starting from pioneering works by B. Lindblad (lindblad41 (1941), lindblad42 (1942)), a number of authors extensively investigated this idea, and a variety of techniques have been developed to study the stability properties of gravitating disks. Most of these studies concentrated on one-component gaseous or collisionless stellar disks. Galactic disks, however, are multi-component systems whose properties depend on all the constituents. Lin & Shu (1966) first addressed the question showing the relative importance of gas in providing the spiral gravitational field. Jog & Solomon (jog84 (1984)) (hereafter JS) considered the influence of a cold component on the stability properties of gravitating disks. They derived a local dispersion relation for a two-fluid instability of a star-gas disk and discussed a relative contribution of two components in the stability properties of the disk. In further investigations the results of JS were generalized by Bertin & Romeo (1988) who studied the stability properties of the discrete global spiral modes by applying asymptotic WKB-methods. They also introduced the concept of the effective $`Q`$-parameter of the two-component disks. Elmegreen (1995) re-considered the original equations derived by JS and reduced them to a single effective $`Q`$-parameter. Romeo (1992, 1994) included the effects of the finite thickness of stellar and gaseous disks. Recently, Rafikov (rafikov01 (2001)) gave the analysis of the axisymmetric gravitational stability of the thin rotating disk consisting of several components.

The non-asymptotic global modal stability properties of multi-component disks were studied to a lesser extent. Noh et al. (noh91 (1991)) studied the global stability properties of a thin two-fluid protoplanetary disk composed of gas and dust coupled by gravity and friction and found that the stability properties of these disks are strongly affected by a small amount of dust. Kikuchi et al. (kikuchi97 (1997)) studied global stability properties of a three-phase disk with phases coupled by self-gravity and interchange processes between the phases. Korchagin & Theis (korchagin99 (1999)) numerically studied the behavior of global modes in multi-component disks at the non-linear stage of global instability.

In this paper, we return to the problem of the global stability properties of multi-component disks. We consider the simplest possible system of that kind – a two-fluid thin disk coupled by the self-gravity of the components. We compare the results of our non-linear two-dimensional numerical simulations with the local stability criteria derived by JS and study the applicability of the Jog-Solomon criterion for the prediction of the global stability properties of multi-phase gravitating disks. We also study the non-linear behaviour of the unstable global modes in multi-phase disks determining the dependence of the saturation levels of the unstable global modes on the mass fraction of a cold component.

Sect. 2 describes the basic equations and the numerical code used in simulations. In Sect. 3 we present the equilibrium model used in our simulations. In Sect. 4 we briefly describe the two-fluid local stability criterion. Sect. 5 presents the results of our simulations and gives a comparison of the results of the numerical simulations with the two-fluid analytical stability criterion.

## 2 Basic Equations

To model the behaviour of the global perturbations in a gravitationally coupled multi-component disk we use the hydrodynamic approximation, modelling both a stellar and a gaseous component as polytropic fluids:

$$P_{s,g}=K_{s,g}\sigma _{s,g}^{\gamma _{s,g}}$$
(1)
Here $`P_{s,g}`$ are the vertically integrated pressures of each component, $`K_{s,g}`$ are the polytropic constants and $`\gamma _{s,g}`$ are the polytropic indices for stellar and gaseous disks, respectively. In this approach, stellar and gaseous ‘fluids’ differ by the values of their polytropic indices $`\gamma _{s,g}`$. We choose the value $`\gamma _s=2`$ for the stellar disk which reproduces the empirical “square root” proportionality between the velocity dispersion and the surface density found in galactic disks (Bottema bottema92 (1992)). Additional arguments justifying the description of the collisionless stellar disks by a fluid dynamical approach are given by Kikuchi et al. (kikuchi97 (1997)). The polytropic index for the gaseous disk is assumed to be $`\gamma _g=1.66`$.

The behaviour of a two-fluid gravitationally interacting system is described by the continuity and momentum equations for each component:

$$\frac{\sigma _{s,g}}{t}+\frac{1}{r}\frac{}{r}\left(r\sigma _{s,g}u_{s,g}\right)+\frac{1}{r}\frac{}{\varphi }\left(\sigma _{s,g}v_{s,g}\right)=0$$
(2)
$`{\displaystyle \frac{u_{s,g}}{t}}+u_{s,g}{\displaystyle \frac{u_{s,g}}{r}}+{\displaystyle \frac{v_{s,g}}{r}}{\displaystyle \frac{u_{s,g}}{\varphi }}{\displaystyle \frac{v_{s,g}^2}{r}}=`$
$`{\displaystyle \frac{1}{\sigma _{s,g}}}{\displaystyle \frac{P_{s,g}}{r}}{\displaystyle \frac{}{r}}\left(\mathrm{\Psi }+\mathrm{\Psi }_H\right),`$ (3)
$`{\displaystyle \frac{v_{s,g}}{t}}+u_{s,g}{\displaystyle \frac{v_{s,g}}{r}}+{\displaystyle \frac{v_{s,g}}{r}}{\displaystyle \frac{v_{s,g}}{\varphi }}+{\displaystyle \frac{v_{s,g}u_{s,g}}{r}}=`$
$`{\displaystyle \frac{1}{\sigma _{s,g}r}}{\displaystyle \frac{P_{s,g}}{\varphi }}{\displaystyle \frac{1}{r}}{\displaystyle \frac{}{\varphi }}\left(\mathrm{\Psi }+\mathrm{\Psi }_H\right).`$ (4)
Additionally, all components are coupled by their gravity, which can be expressed by the Poisson equation

$$\mathrm{\Psi }(r,\varphi )=_{R_{\mathrm{in}}}^{R_{\mathrm{out}}}_0^{2\pi }\frac{(\sigma _s(r^{},\varphi ^{})+\sigma _g(r^{},\varphi ^{}))r^{}dr^{}d\varphi ^{}}{\sqrt{r^2+r^22rr^{}\mathrm{cos}(\varphi \varphi ^{})}}.$$
(5)
Here $`u_{s,g}`$ and $`v_{s,g}`$ are the radial and azimuthal velocities of stars and gas within the disk, and $`\sigma _{s,g}`$ is the surface density of each component. The potential $`\mathrm{\Psi }`$ derived from the self-gravity of a two-fluid disk, the explicit contribution $`\mathrm{\Psi }_H`$ from the rigid halo, and the pressure gradient determine the behaviour of the disk. All of the dependent variables are functions of the radial coordinate $`r`$, the azimuthal angle $`\varphi `$, and the time $`t`$. The variables and the parameters of the model are expressed in units in which the gravitational constant $`G`$ is equal to unity, the unit of mass is equal to $`10^{10}M_{}`$ and the unit of length is $`2`$ kpc. This choice determines the velocity unit $`V_{\mathrm{unit}}=149.1`$ km s<sup>-1</sup> and the time unit $`T_{\mathrm{unit}}=1.32\times 10^7`$ yr.

To solve the two-component hydrodynamical equations (2) – (2) we use the multi-component fluid-dynamical code developed by Korchagin & Theis (korchagin99 (1999)). The code is a multi-phase realization of the Eulerian ZEUS-2D code (Stone & Norman stone92 (1992)) with the Van-Leer advection scheme. The two-fluid hydrodynamical equations are solved using equally spaced azimuthal and logarithmically spaced radial zones. In the numerical simulations discussed here we employed a grid of 256 $`\times `$ 256 cells. The Poisson equation is solved by the 2D Fourier convolution theorem in polar coordinates. We look for the global modes starting from initial random perturbations at the level $`10^6`$.

## 3 Equilibrium model

The equilibrium disk taken in our simulations is based on the equilibrium properties of real galaxies. Specifically, we applied the model of the spiral galaxy NGC 1566 by Korchagin et al. (korchagin00 (2000)). The disk of NGC 1566 is observationally well studied. Its rotation curve taken from Bottema (bottema92 (1992)) can be approximated by

$$v_0(r)=\frac{V_1r}{(r^2+R_1^2)^{3/4}}+\frac{V_2r}{(r^2+R_2^2)^{3/4}}.$$
(6)
The radial dependence of the vertical velocity dispersion of the stars in the disk of NGC 1566 has an exponential distribution

$$c_z=c_z(r)=c_{z0}\mathrm{exp}(r/h_{c_z})\text{.}$$
(7)
This uniquely determines the radial dependence of the axisymmetric surface density distribution assuming that the disk is self-gravitating in the vertical direction:

$$\sigma (r)=c_z^2(r)/\pi Gz_0\text{.}$$
(8)
Here $`G`$ is the gravitational constant, and $`z_0`$ is the effective thickness of the disk. According to equations (7) and (8), the radial distribution of the stellar surface density can be expressed by

$$\sigma _s(r)=\sigma _{0_s}\mathrm{exp}(r/h_\sigma )[1(r/R_{\mathrm{out}})^2]^5\text{,}$$
(9)
where $`h_\sigma `$ is radial scale length for the stellar surface density distribution. For numerical purposes we use in equation (9) a multiplier $`[1(r/R_{\mathrm{out}})^2]^5`$ which results in a vanishing surface density at the outer boundary of the disk. We further assume that the surface density profile of the gaseous disk and its thickness are the same as those in the stellar disk. However, the total mass of the gaseous disk is smaller compared to the mass of the stellar disk and varies in our models.

The radial velocity dispersion of the stellar and gaseous components can be derived from the equation of state by

$$c_{r_{s,g}}=(\gamma _{s,g}K_{s,g}\sigma _{s,g}(r)^{\gamma _{s,g}1})^{1/2}\text{ .}$$
(10)
The free parameters $`K_{s,g}`$ can be fixed by the equations (7) and (8) and by the assumption of a constant ratio of the vertical and radial velocity dispersions. For the latter, to get the disk with weak global modes we have chosen $`c_{sz}/c_{sr}=0.9`$. Our choice can be justified as follows. According to Bottema (1992), the measured velocity dispersion in the disk of NGC 1566 gives the real dispersion perpendicular to the galactic disk. In our one-component model, we fix the $`z`$-component of the velocity dispersion and hence the surface density and the total mass of the disk, and allow for a variation of the radial velocity dispersion. Increasing the ratio $`c_{sz}/c_{sr}`$ we decrease the radial velocity dispersion, hence destabilizing the disk. Simulations show that for the measured $`z`$-dispersion and rotation curve in the disk of NGC 1566, the disk is stable for the usually assumed ratio $`c_{sz}/c_{sr}=0.6`$. Therefore, we choose in our studies the ratio $`c_{sz}/c_{sr}=0.9`$ which leads to a weakly unstable disk. This ratio is larger than the measured ratio of the velocity dispersions $`c_{sz}/c_{sr}=0.50.6`$ in solar neighborhood, and close to the upper limit measured in other galaxies. E.g., the determination of this ratio in the disk of galaxy NGC 488 undertaken by Gerssen et al. (1997) gives the value $`0.7\pm 0.19`$. With our choice, however, we do not aim to model any particular galaxy.

The equilibrium rotation of the disk given by the equation (6) is balanced by the self-gravity, the pressure gradient of the disk and the external potential of a rigid halo:

$$\frac{v_0^2}{r}=\frac{1}{\sigma }\frac{dP_s}{dr}+\frac{d}{dr}\left(\mathrm{\Psi }+\mathrm{\Psi }_H\right).$$
(11)
By this, we can fix the halo properties. Knowing the initial surface density of the disk we get the disk’s self-gravity by solving the Poisson equation numerically. Equation (11) can then be used to calculate the gradient of the halo potential and the halo mass distribution. The general parameters of our equilibrium model are listed in Table 1.

## 4 Stability Criteria

The local stability of a thin gravitating disk can be expressed in terms of Toomre’s $`Q`$ parameter which for the gaseous disk has the form (Safronov safronov60 (1960), Toomre toomre64 (1964)):

$$Q\frac{c(r)\kappa (r)}{\pi G\sigma (r)}.$$
(12)
Here, $`c(r)`$ is the sound speed of the gas, $`\kappa (r)`$ is the epicyclic frequency and $`\sigma (r)`$ is the surface density of the disk. A gravitating disk is locally unstable when $`Q<1`$.

The stability criterion (12) was derived for tightly wound perturbations. However, it can be successfully used as an indicator for the stability of global modes in gravitating disks. Numerical simulations show that a locally stable stellar or gaseous disk can be globally unstable if the minimum value of $`Q`$ does not exceed $`2`$. For example, a one-component disk with the rotation and surface density distribution given by equations (6) and (9) can be globally unstable, if the minimum value of $`Q`$ does not exceed $``$ 1.6 (Korchagin et al. korchagin00 (2000)).

Disks of real galaxies consist of stellar and gaseous components. The stellar component comprises the largest fraction of the total mass of the disk. However, colder gaseous components can play a significant role for the stability properties of multi-component disks. JS derived a local stability criterion for a two-component fluid disk which can be written in the form $`Q_{\mathrm{JS}}<1`$ where

$$Q_{\mathrm{JS}}=\left(\frac{2}{Q_s}\frac{q}{1+q^2}+\frac{2}{Q_g}\frac{c_g}{c_s}\frac{q}{1+q^2c_g^2/c_s^2}\right)^1.$$
(13)
Here $`Q_s`$ and $`Q_g`$ are Toomre’s $`Q`$-parameters determined for the ‘stellar’ and gaseous fluids, correspondingly. $`q`$ is the dimensionless radial wavenumber of perturbations $`q=kc_s/\kappa `$. The relative contribution $`\mathrm{\Gamma }`$ of gas to stars is given by the ratio of the second to the first term of the stability parameter (13):

$$\mathrm{\Gamma }=\frac{\sigma _g}{\sigma _s}\frac{1+q^2}{1+q^2c_g^2/c_s^2}.$$
(14)
It can be seen from equation (13) that the relative contribution of the components is always larger than the ratio of the surface densities of the components, provided the gas is dynamically colder than the stars. In the limiting case of large wavenumbers $`q\mathrm{}`$ the contribution of the components is given by the ratio $`\sigma _gc_s^2/\sigma _sc_g^2`$ \- the result, obtained by Marochnik & Suchkov (marochnik74 (1974)).

We will use the criterion given by equation (13) for the analysis of instability of the global modes in two-component gravitating disks.

## 5 Results

To study the influence of a cold component on the overall dynamics of a two-fluid gravitating disk, we start our analysis with a one-component ‘stellar’ disk which we choose to be close to marginal stability. Fig. 1 shows the equilibrium properties of this disk. The local stability parameter $`Q`$ is well above unity with a minimum value $`Q_{\mathrm{min}}1.6`$.

The radial behavior of the $`Q`$-parameter is one of the key assumptions in the modal analyses of spiral density waves. In previous studies it was usually assumed that the profile of the stability parameter $`Q`$ has a barrier at the central regions of the disk and is equal to unity in the outer regions of the disk (e.g. Bertin et al. 1989). Such a behavior of an ‘effective’ $`Q`$-parameter is supposed to model the absence of cold gas as well as the transition in geometry from the disk to the nuclear bulge in the central regions of a galaxy, and to describe the increasing role of a cold gaseous component in the outer regions of a galactic disk. In this paper we do not introduce an effective $`Q`$-parameter, but we try to base our studies on observational properties of disk galaxies. We choose the basic state in our fiducial model based on the properties of the nearby late type spiral galaxy NGC 1566. The ‘stellar’ $`Q`$-parameter shown in Fig. 1 has an outer barrier which is typical for disks with exponentially decreasing radial surface density distributions and exponentially decreasing velocity dispersions.

Linear global modal analysis of this model reveals a few slowly growing modes with the fastest growing spiral $`m=2`$. Nonlinear 2D simulations confirm this result giving an exponential growth rate of the $`m=2`$ mode of approximately 0.11 in dimensionless units. Linear theory predicts a value of 0.12 which is in good agreement with non-linear simulations. Fig. 2 illustrates the evolution of the $`m=14`$ global modes in the one-component ‘stellar’ disk given in terms of the global Fourier amplitudes

$$A_m\frac{1}{M_d}\left|_0^{2\pi }_{R_{\mathrm{in}}}^{R_{\mathrm{out}}}\sigma (r,\varphi )r𝑑re^{im\varphi }𝑑\varphi \right|\text{ .}$$
(15)
The most unstable $`m=2`$ mode grows linearly until $`t110`$, and then saturates at a level $`A_20.2`$ (Fig. 2). The linear regime corresponds to an exponential growth of perturbations.

In the next step we split the one-component disk into a ‘stellar’ and a ‘gaseous’ subsystem keeping the total mass of the disk constant. We consider a set of models with a cold component consisting of 5%, 10% and 20% of the total mass of the disk. In the first series of simulations we keep the ‘stellar’ and gaseous equations of state the same for all the models. As a result, the gas velocity dispersion, which depends on the gas fraction, is not constant. Thus, the velocity dispersion of gas is about 20% of the velocity dispersion of the stars in the model with 5% gas admixture. This ratio increases up to 36% for the model with 20% gas admixture.

The strong influence of the gas fraction is shown in Fig. 3 which displays the growth rate of the most unstable $`m=2`$ mode in a one-component ‘stellar’ disk compared to the modal growth rate of the two-component models with 5%, 10%, and 20% gas admixtures. Already a five percent gas fraction enhances the modal growth rate by a factor of 1.7 compared to the one-component disk. A 20% gas admixture results in an increase of the growth rate by a factor of 3.7. Table 2 gives the masses and the dispersion ratio of the components in the computed models together with the growth rates of the most unstable modes, saturation levels and the minimum values of the two-component stability criterion $`Q_{\mathrm{JS}}`$ calculated by equation (13). The growth rate of the dominant mode increases in the models with higher gas content, being the same for the gaseous and the stellar components. However, the saturation level of the most unstable mode is systematically higher for the gaseous disk. In general, the saturation level slightly decreases in the models with higher gas content. A possible explanation of this behavior is that spiral shocks play a more important role in the saturation of exponentially growing modes when more cold gas is present.

The ratio of the velocity dispersions of the subsystems can also affect the stability of the disk. The influence of the velocity dispersions on the growth rates of the global modes was studied in a series of numerical experiments in which the fraction of the cold component was fixed to 10% of the total mass of the disk, while the velocity dispersion ratio was allowed to change. The results of these simulations are shown in Fig. 4. The growth rate of the most unstable global mode is only slightly affected by the ratio of the velocity dispersions. The saturation amplitude remains approximately at the same level. Therefore, the mass fraction of the cold component is a more important parameter for the growth rate of perturbations than the ratio of the velocity dispersions. This is in agreement with results of Bertin & Romeo (1988). They demonstrated (see their Fig. 5) that the criterion of marginal stability of a two-component disk also depends slowly on the ratio of the velocity dispersions of the components and more strongly on the ratio of their densities.

Disks with gaseous admixture show qualitatively the same evolution as the pure stellar disk. However, they evolve on a shorter time scale. The appearance of the global modes in all models is quite similar despite the considerable difference in their growth rates. Fig. 5 shows the logarithmically spaced surface density contour plots taken at the moments when the global amplitude of the $`m=2`$ mode in all four models is $`A_20.06`$. All the models develop an open two-armed spiral of comparable winding. Note, however, that the spiral arms are more narrow, more tightly wound and less extended when more gas is present. Analysis of the global modes with help of the local dispersion relation (Bertin & Romeo 1988) also shows a shift towards tighter spirals with a smaller corotation radius in models with higher gas content.

To use the stability criterion (13), one needs to know the wavenumber of the perturbations. Fig. 6 shows the radial dependence of the radial wavenumber for $`m=2`$ perturbations computed for the models I - III. The wavenumbers were determined numerically by calculating the pitch angles of the spiral arms shown in Fig. 5 at different radii, i.e. we used the relation

$$\mathrm{cot}i=\left|\frac{kR}{m}\right|.$$
(16)
To find a pitch angle of the spiral at a particular radius, we determine the line of the maximum amplitudes of the spiral arm and calculate the angle between the tangents to the line of maxima and the corresponding circle.

Fig. 7 shows the radial profile of the dimensionless Jog-Solomon stability parameter $`Q_{\mathrm{JS}}`$ computed with the wavenumbers taken from Fig. 6. Fig. 7 further illustrates the destabilizing role of a cold gas component: the minimum value $`Q_{\mathrm{JS}}`$ of the Jog-Solomon parameter decreases with increasing content of cold gas. The minimum value of the Jog-Solomon parameter $`Q_{\mathrm{JS}}`$ is about 1.36 for the two-component system with 20% gas admixture versus a minimum value of 2.2 found for the two-component model with 5% gas fraction. In all three models, the Jog-Solomon stability criterion exceeds unity, and all the models are locally stable. Such a behaviour is quite similar to the behaviour of one-component disks which are known can be globally unstable even if the minimum value of Toomre’s $`Q`$ parameter is above unity.

Decreasing the total mass of the disk results in its stabilization which consequently leads to an overall growth of the Jog-Solomon parameter $`Q_{\mathrm{JS}}`$. Fig. 8 illustrates this by displaying the growth rate of the dominant mode $`m=2`$ as a function of the minimum value of the Jog-Solomon parameter $`Q_{\mathrm{JS}}`$. The latter was obtained from numerical simulations of two-component disks with 10% gas admixture. These disks are stabilized when the minimum value of the Jog-Solomon parameter exceeds $`2.5`$.

## 6 Summary

In this paper we studied the applicability of the two-fluid local stability criterion derived by Jog & Solomon (jog84 (1984)) for the onset and growth of non-axisymmetric perturbations in two-component gravitating disks with an open spiral structure. We performed three sets of numerical experiments: one set with fixed total mass of the disk and fixed polytropic equation of state, but variable gas fraction. Another set of simulations was performed with a fixed total mass of the disk and a fixed ratio of the stellar and gaseous velocity dispersions. A final set of calculations deals with a variable disk mass. Our numerical experiments have shown that the mass contribution of the cold component strongly influences the stability properties of the disk. Even a low gas fraction of $`5\%`$ destabilizes the disk substantially, e.g. the growth rate of $`m=2`$ perturbations increases by a factor 1.7. All our disk models are unstable to $`m=2`$ spiral modes. More unstable disks develop narrower and more tightly wound spiral arms. The saturation levels are almost constant for all simulations, with a slightly enhanced level for the gaseous component. A comparison of the results of the direct two-dimensional numerical simulations with the JS stability criterion allows us to conclude that similarly to Toomre’s stability parameter $`Q`$ derived for one-component disks, the Jog-Solomon stability parameter $`Q_{\mathrm{JS}}`$ can be used for the prediction of the global stability properties of two-component gravitating disks which greatly simplifies the analysis of disk stability.

###### Acknowledgements.

We thank the referee Dr. A.B. Romeo for his valuable comments on this work. NO and VK gratefully acknowledge support by the Deutsche Forschungsgemeinschaft under grant RUS 17/75/00(S).


text:
# On the “initial” Angular Momentum of Galaxies

## 1 Introduction

The high frequency of disk galaxies and their exponential profiles over many scale radii is perhaps one of the most striking global features of observed galaxies. Theoretical studies investigating the origin of this angular momentum in these galaxies have a long history (e.g. Weizsäcker 1951). In most modern investigations the Ansatz of Mestel (1963) is used where one assumes that any parcel of fluid to retains its initial specific angular momentum. In hierarchical structure formation scenarios the luminous parts of galaxies form from gas that is cooling within a dark matter (DM) halos that continuously are merging to build larger and larger objects. These assumption are made also in most recent models of disk galaxies (e.g. Mo, Mao and White 1998). The typical initial angular momentum of halos may be estimated from linear theory (Peebles 1969). However, its distribution is typically measured from N–body simulations that follow the formation and evolution of structure in the dark matter only. In the following we derive a rough estimate of the the magnitude of pressure forces that contribute significantly to the torquing forces.

## 2 Forces and Torques

The critical density and the total mass, $`M`$, in a spherical volume of comoving radius $`R`$ with an over–density $`\delta =\rho /\overline{\rho }_m`$ of the mean matter density ($`\overline{\rho }_m=\mathrm{\Omega }_m\rho _c`$) are given by,

$`\rho _c`$ $`=`$ $`{\displaystyle \frac{3H_0^2}{8\pi G}},M=\delta {\displaystyle \frac{4\pi }{3}}R^3\mathrm{\Omega }_m\rho _c,`$ (1)
where the symbols have their usual meaning (see e.g. Peebles 1993). Let $`f`$ be the fraction of the torquing force from surrounding structures on this spherical over-density to the radial force of the sphere. I.e. the torquing force is $`F_t=fGM\rho /r^2`$ where $`r`$ is a proper radius of the sphere. Here $`f`$ is necessarily somewhat smaller than the radial force of its own gravity111Otherwise the fluid parcel would not collapse with the sphere.. To get a rough idea of what pressure gradient may change the torquing we simply compare

$`{\displaystyle \frac{1}{\rho _b}}{\displaystyle \frac{P}{r}}`$ $`=`$ $`f{\displaystyle \frac{GM}{r^2}}`$ (2)
Using difference instead of differential, substituting $`r=R/(1+z)`$ and using $`\mathrm{}R`$222This $`\delta R`$ is taken smaller than the total radius since for a realistic density distribution most of the mass is at large radii.to denote the comoving length scale over which the pressure gradient is taken we find,

$`\mathrm{}P`$ $``$ $`f{\displaystyle \frac{3H_0^4}{16G\pi }}\delta ^2\mathrm{\Omega }_b\mathrm{\Omega }_m(1+z)^4R\mathrm{}R.`$ (3)
Assuming constant density this translates to the very small temperature gradient of

$`\mathrm{}T`$ $``$ $`f{\displaystyle \frac{\mathrm{}P\mu }{k_B\rho _b}}=\delta {\displaystyle \frac{H_0^2m_P}{2k_B}}\mathrm{\Omega }_m(1+z)R\mathrm{}R,`$ (4)
$`\mathrm{}T`$ $``$ $`2.6\times 10^3\text{K}\delta {\displaystyle \frac{f}{0.1}}(1+z)R_{Mpc}{\displaystyle \frac{\mathrm{}R_{Mpc}}{0.3}}.`$ (5)
Such small pressure gradients are expected to originate from well known processes such as ionization fronts from reionization, galactic outflows driven by black holes, and galactic feedback as required in semi–analytical models of galaxy formation (e.g. Kauffmann et al. 2000). The larger the mass of the object the larger the pressure gradient has to be to alter the angular momentum acquisition of the gas as compared to the dark matter.

## 3 Simulations

To investigate the relation between angular momentum of the baryons and the DM we have performed two different cosmological smoothed particle hydrodynamics simulations using GADGET (Springel, Yoshida, and White 2000). Both simulations use 1283 for each the DM and the gas in a periodic volume of 10 comoving Mpc side-length. Both are evolved from redshift 60 to $`z=3`$ with initial conditions drawn from a realization appropriate for a spatially flat universe with $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_bh^2=0.02`$ and $`h=0.7`$. Both simulations follow only gravity and adiabatic gas physics. In one run we have artificially increased abruptly the temperature of the gas to $`4\times 10^6`$ Kelvin at redshift 7 to mimic a strong feedback case.

We use HOP (Eisenstein and Hut 1998) to identify centers of halos. The we go out radially and define a halo as all the mass within the radius where the over–density in the dark matter corresponds to 200 times the mean density. Within these halos we then measure the total and specific angular momentum vectors. In Figure 1 we show the distributions of the alignment angles between the dark matter and the baryonic component for all halos that have no other overlapping halo. Clearly for the purely adiabatic case their is remarkable agreement between the directions of the angular momenta with a median alignment angle of 20 degrees. However, as soon as pressure forces are introduces as in the ad–hoc feedback the respective angles become uncorrelated. The imperfect alignment between the DM and the gas in the purely adiabatic run comes most likely from non–linear angular momentum exchange at virialization. Such non-linear effects may well be aggravated in simulations which also include the cooling of the gas.

## 4 Discussion and Conclusions

In this short contribution we have pointed out that large scale pressure forces have significant influence on the formation of galaxies. Until the physical processes acting on the gas in the intergalactic medium are fully understood it is difficult to asses up to which mass scale such large scale non-gravitational forces may play a role. Because of space constraints we could only remark here that the collisionless and collisional components will exchange angular momentum during merging events. This will further diminish and possible correlation of the spins of DM and gas. Our results have important implications for models which associate the dark matter angular momentum vector and their correlations with observed ellipticities of galaxies (see contributions by Natarayan, Crittenden and Heavens in this volume). Since the directions of the angular momentum of radial shells within halos varies as function of radius (Moore this volume) and the luminous part of galaxies may be re–oriented from the torquing their non-spherical parent halos (Moore and Sellwood, this volume) we are led to conclude that current models of the formation and evolution of disk galaxies are oversimplified. However, one interesting exception may be low redshift massive elliptical galaxies. If most of their stars have formed at early times in smaller higher redshift objects the shapes and angular momentum distributions of gas and DM may be strongly correlated. This would be a natural explanation of the strong lensing observations of Kochaneck (this volume) which indicate that “light traces mass” for these galaxies.

## References


text:
# Thermalization and elliptic flow at RHIC

## 1 Introduction

A highly debated issue in ultrarelativistic heavy ion collisions is whether or not the colliding system reaches local thermal equilibrium at some stage of its evolution. Equilibrium with respect to inelastic collisions strongly constrains the ratios of particle abundances  (“chemical” equilibrium), as well as phase-space densities, which are obtained by combining informations from momentum spectra and two-particle HBT correlations; on the other hand, equilibrium with respect to elastic collisions constrains momentum distributions, and implies in particular that they are isotropic in the local rest frame. This is the “kinetic” equilibrium, on which we concentrate here.

Kinetic equilibrium itself has (at least) two facets. One is the equilibration between longitudinal and transverse degrees of freedom, i.e., the implication that in the local rest frame, longitudinal and transverse momenta are of the same order of magnitude. This aspect of thermalisation can be discussed from first principles at the partonic level, and there is now a vast literature on this subject. But experimental signatures deal in fact rather with equilibration among the two transverse degrees of freedom. Due to the high Lorentz contraction at ultrarelativistic energies, the typical transverse scale is much larger than the longitudinal scale, so that this “transverse equilibrium” is probably easier to achieve than “longitudinal-transverse equilibrium”.

In Sec. 2, we discuss some experimental signatures of transverse equilibrium, especially elliptic flow, and the relevance of hydrodynamical models in this context. In Sec. 3, we present a new method recently developed to obtain more reliable measurements of elliptic flow, which are required in order to draw definite conclusions on the issue of thermalisation.

## 2 Elliptic flow: a signature of transverse thermalisation

### 2.1 When are hydrodynamical models useful?

If the colliding system thermalises, its evolution follows the laws of hydrodynamics; it is therefore natural to use hydrodynamical models in order to define signatures of thermalisation. There are two ingredients in these models: i) an equation of state; ii) initial conditions, i.e., energy density, baryon density and fluid velocity on a space-like hypersurface (typically, at some initial time).

Which parameters are under control? The initial time $`t_0`$, at which the system thermalises, is hopefully short but to a large extent unknown, so that hydrodynamics has little predictive power concerning observables which depend strongly on $`t_0`$, such as thermal photon production. Similarly, there is no serious motivation for studying the very first stages of the collision within the framework of hydrodynamics, either within a one-fluid model  or a multifluid model, and the relevance of signatures based on such parametrisations is questionable. Finally, there is no well-defined prescription concerning the initial longitudinal fluid velocity and density, for which a variety of parametrisations exist, either “Landau”  or “Bjorken”  type.

The situation is much clearer concerning transverse degrees of freedom: the transverse collective velocity must be initially zero since each nucleon-nucleon collision populates the transverse momentum space randomly; the initial density profile in the transverse plane is strongly constrained by observed multiplicity distributions. Therefore, reliable signatures of thermalisation should rather be sought for in observables associated with transverse degrees of freedom: $`p_T`$ spectra, transverse radii  and azimuthal anisotropies, in particular elliptic flow, on which we concentrate here.

### 2.2 Why elliptic flow?

Elliptic flow is defined as a correlation between the azimuth $`\varphi `$ of an outgoing particle and the azimuth $`\mathrm{\Phi }_R`$ of impact parameter:

$$v_2=e^{2i(\varphi \mathrm{\Phi }_R)},$$
(1)
where brackets denote a statistical average. At ultrarelativistic energies, $`v_2`$ is positive for noncentral collisions. In a hydrodynamical picture, it results from anisotropic pressure gradients in the transverse plane, due to the almond-shaped region of the overlap region between the two nuclei. Microscopically, $`v_2`$ is created by rescattering among the produced particles, which makes it a sensitive probe of final state interactions: if there are none, it vanishes, while other observables such as $`p_T`$ spectra may still look “thermal”.

Furthermore, predictions of hydrodynamical models for $`v_2`$ are very stable: since $`v_2`$ is created by transverse pressure gradients, it strongly depends on the initial density transverse profile, which is well controlled as discussed above; it also depends significantly on the equation of state, on which it may thus provide valuable information; on the other hand, $`v_2`$ depends only weakly on arbitrary parameters, such as initial time and longitudinal velocity. Quite remarkably, simple hydrodynamical parametrisations are able to reproduce simultaneously the measured $`p_T`$ spectra, HBT radii and elliptic flow at RHIC.

### 2.3 Predictions for $`v_2`$

The momentum anisotropy $`v_2`$ calculated in hydro models is roughly equal to the anisotropy of the almond-shaped overlap region. This purely geometrical effect dominates the centrality dependence of $`v_2`$, which decreases linearly with the number of participants. Deviations from this behaviour can be used to signal a phase transition  or a departure from thermal equilibrium. The latter is expected to occur for the most peripheral collisions, where $`v_2`$ should be smaller than the hydro prediction. This is indeed observed in several specific transport models like UrQMD  (which however predicts a much too small value of $`v_2`$), QGSM, and AMPT. With a more systematic study, one could relate the observed centrality dependence of $`v_2`$ to the degree of thermalisation of the system. Experimental results vary significantly depending on the method used to analyse elliptic flow. We come back to this issue in Sec. 3.

Hydrodynamical calculations  were also able to predict the $`p_T`$ dependence of $`v_2`$ for identified hadrons, in remarkable agreement with experimental results: $`v_2`$ is almost linear in $`p_T`$ for pions and significantly smaller for protons. However, these non-trivial features are also reproduced by transport models. In addition, the latter predict a saturation  at high $`p_T`$ which is not seen in hydro, suggesting that many elastic collisions are necessary to build the flow at high $`p_T`$. This saturation, which is seen in the data, could also be related to hard physics.

## 3 Analysing elliptic flow with multiparticle correlations

### 3.1 Flow from azimuthal correlations

Since the reaction plane $`\mathrm{\Phi }_R`$ is unknown, $`v_2`$ cannot be derived directly from (1). It must be inferred from azimuthal correlations. The standard flow analysis  relies on the key assumption that all azimuthal correlations are due to flow, i.e., that angles relative to the reaction plane $`\varphi \mathrm{\Phi }_R`$ are statistically independent. This allows one to write the two-particle correlation as

$$e^{2i(\varphi _1\varphi _2)}=e^{2i(\varphi _1\varphi _R)}e^{2i(\varphi _R\varphi _2)}=(v_2)^2.$$
(2)
where brackets denote an average over pairs of particles belonging to the same event. One could also use higher order correlations, such as the four-particle correlation, which would give

$$e^{2i(\varphi _1+\varphi _2\varphi _3\varphi _4)}=(v_2)^4.$$
(3)
However, such equations are not quite correct, since they neglect various other sources of azimuthal correlations (“nonflow correlations”), which is no longer justified at ultrarelativistic energies.

### 3.2 Simple illustration of nonflow correlations

In order to illustrate nonflow correlations, we consider the following example: assume that in each event, $`M/2`$ pairs of particles are emitted, where both particles in a pair have collinear momenta, but pairs are emitted with random orientations. Since pairs are emitted randomly, there is no flow ($`v_2=0`$), but there are azimuthal correlations. In each event, there is a total of $`M(M1)/2`$ particle pairs, among which $`M/2`$ are correlated, hence the two-particle correlation

$$e^{2i(\varphi _1\varphi _2)}=\frac{1}{M1}.$$
(4)
A similar reasoning yields the four-particle correlation:

$$e^{2i(\varphi _1+\varphi _2\varphi _3\varphi _4)}=\frac{2M(M2)}{M(M1)(M2)(M3)}=\frac{2}{(M1)(M3)}.$$
(5)
Applying Eqs. (2) and (4), or (3) and (5), one would obtain $`v_21/\sqrt{M}`$, although there is no flow: this is the typical order at which nonflow correlations spoil the standard flow analysis.

### 3.3 Subtraction of nonflow correlations

The contribution of nonflow correlations can be greatly reduced by combining the informations from two- and four-particle correlations. In Eq. (5), a 4-uplet of particles gives a nonvanishing contribution if it consists of two correlated pairs, either (1,3) and (2,4) or (1,4) and (2,3). Subtracting the corresponding contributions, one obtains

$$e^{2i(\varphi _1+\varphi _2\varphi _3\varphi _4)}2e^{2i(\varphi _1\varphi _2)}^2=\frac{4}{(M1)^2(M3)}.$$
(6)
The l.-h. s. of this equation defines the cumulant of four-particle correlations, $`c_2\{4\}`$. The r.-h. s. is the contribution of nonflow correlations, of order $`1/M^3`$, i.e., much smaller than the corresponding contribution to the four-particle correlation (5), of order $`1/M^2`$. On the other hand, the contribution of flow remains of the same magnitude: from Eqs. (2), (3) and (6), one obtains $`c_2\{4\}=(v_2)^4`$. The subtraction therefore reduces the relative contribution of nonflow effects, and yields a more accurate estimate of the flow.

This method was recently applied to STAR data. The resulting value of $`v_2`$ are smaller than those obtained with the standard analysis, in particular for the most peripheral collisions: this is precisely where nonflow effects are expected to give the largest contribution since the multiplicity $`M`$ is smaller. The centrality dependence obtained with this method suggests that departures from thermalisation at RHIC may be larger than was previously thought.

This cumulant expansion can be worked out to arbitrary orders, and allows one to extract the genuine 4-, 6-particle correlations and beyond. The practical implementation of the method is described in detail elsewhere. Flow, which is essentially a collective phenomenon, contributes to all orders, while the relative contribution of nonflow correlations decreases as the order increases. Higher order cumulants therefore provide a unique possibility to check quantitatively that azimuthal correlations are indeed of collective origin.

## Acknowledgments

J.-Y. O. thanks Raimond Snellings and Raju Venugopalan for discussions.