text:
# Density profile of a strictly two-dimensional Bose gas at finite temperature

## 1 Introduction

Recent experiments have realized quasi two-dimensional (2D) cold Bose systems by tuning the anisotropy of the trapping potential  and have stimulated interest in studying Bose-Einstein condensation in a low-dimensional regime. An important difference between the macroscopic 2D fluid and the 3D one is that at finite temperature phase fluctuations in 2D destroy the long-range order, in agreement with the Mermin-Wagner theorem . A Kosterlitz-Thouless transition to a superfluid state still occurs on cooling through the binding of vortices of opposite vorticity  and leads to an algebraic decay of the one-body density matrix. However, in a trapped 2D fluid the modifications of the density of states caused by the confining potential allow a true condensate to exist even at finite temperature.

Another important consequence of lowered dimensionality is that the T-matrix for two-body collisions in vacuo at low momenta and energy, which should be used to obtain the collisional coupling parameters to lowest order in the particle density, vanishes in the strictly 2D limit as the s-wave scattering length becomes larger than the width of the axial trapping . It is then necessary to evaluate the scattering processes between pairs of Bose particles by taking into account the presence of a condensate and a thermal cloud through a many-body T-matrix formalism .

This formalism has already been used in a number of studies of low dimensional Bose gases, dealing in particular with phase fluctuations and the Kosterlitz-Thouless transition in a variational approach , with a mean-field evaluation of the breathing-mode frequency in a trapped 1D gas , and with the equilibrium density profile of a 2D condensate . In the limit of zero temperature the condensate-condensate coupling parameter has been related to the two-body T matrix by considering that, when two particles in the condensate collide at zero momentum, they both require an energy equal to the chemical potential $`\mu `$ to be excited out of the condensate . With increasing temperature the population of the excited states becomes non-negligible and a microscopic theory which also takes into account the depletion of the condensate is required . An additional coupling parameter to describe the scattering processes between an atom in the condensate and an atom in the thermal cloud has been introduced by Stoof and co-workers  through the two-body T-matrix at energy $`\mu `$.

The present paper is the finite-temperature extension of the work by Tanatar et al. , who studied the cross-over from a 3D to a 2D regime in the equilibrium density profiles at zero temperature. Here we limit ourselves to the strictly 2D case, where the scattering length has become larger than the axial thickness of the cloud. We give a semi-analytical expression for the many-body T-matrix elements corresponding to the condensate-condensate and to the condensate-thermal cloud coupling parameters ($`g_2`$ and $`g_1`$, respectively) and we use them in a two-fluid model  to evaluate the condensate fraction and the density profiles of the bosonic gas at increasing temperature. Congruently with the work of Tanatar et al.  we focus on the experimental parameter of Görlitz et al.  for the case of largest anisotropy.

The paper is organized as follows. Section. 2 summarizes the general treatment of the many-body T-matrix for two-body scattering processes in the presence of a condensate and a thermal cloud. In Sec. 3 we evaluate a simplified form of the two-fluid model for a 2D Bose gas at finite temperature. Finally, we present and discuss our results in Sec. 4.

## 2 Coupling parameters in a mixture of condensed and thermal bosons

The scattering processes between pairs of atoms in a gas consisting of a Bose-Einstein condensate and a thermal cloud are described by the many-body T-matrix $`T^{MB}(E)`$ as a function of the energy $`E`$. The coupling parameters are given by the matrix elements $`𝐤^{}|T^{MB}(E)|𝐤T^{MB}(𝐤,𝐤^{},𝐊;E)`$, taken in the limit of zero energy and momenta. Here $`𝐤`$ and $`𝐤^{}`$ are the incoming and outgoing relative momenta of the pair of center-of-mass momentum $`𝐊`$. We shall consider only the condensate-condensate and condensate-thermal cloud couplings and neglect the scattering between thermally excited atoms in the present case of a dilute gas.

Before discussing the many-body T-matrix, however, we shall first recall the behavior of the two-body T-matrix that describes collisions between pairs of particles in vacuo. We shall then discuss how the two T-matrices are related in the appropriate limit for a two-fluid system.

### 2.1 The two-body T matrix

The two-body T-matrix is the solution of the Lippmann-Schwinger equation,

$`𝐤^{}|T^{2B}(\overline{E})|𝐤`$ $`=`$ $`𝐤^{}|V(|𝐫_1𝐫_2|)|𝐤`$ (1)
$`+`$ $`{\displaystyle \underset{𝐪}{}}𝐤^{}|V(|𝐫_1𝐫_2|)|𝐪{\displaystyle \frac{1}{\overline{E}2ϵ_𝐪}}𝐪|T^{2B}(\overline{E})|𝐤,`$
with $`V(|𝐫_\mathrm{𝟏}𝐫_\mathrm{𝟐}|)`$ being the interparticle potential. The center-of-mass energy for the pair of atoms is $`\overline{E}`$ and each atom of mass $`m`$ has single-particle excitation energy $`ϵ_𝐪=\mathrm{}^2q^2/2m`$, as collision takes place in free space. In Eq. (1) the collision is described by a single-loop interaction between the two atoms plus contributions involving all possible transition routes from state $`|𝐤`$ to state $`|𝐤^{}`$ via intermediate states $`|𝐪`$.

In the case of hard-disk potential of strength $`V_0=4\pi \mathrm{}^2/m`$ and in the dilute limit $`𝐤a,𝐤^{}a1`$ where $`a`$ is the 2D scattering length, the solution of Eq. (1) is

$$𝐤^{}|T^{2B}(\overline{E})|𝐤\frac{4\pi \mathrm{}^2/m}{\mathrm{ln}|4\mathrm{}^2/\overline{E}ma^2|}.$$
(2)
In this case the T-matrix is independent of the momenta, but depends on the logarithm of the energy and vanishes as $`\overline{E}`$ approaches zero. Therefore, the presence of the surrounding gas must be taken into account in the collisional processes. Following the proposal of Morgan et al. , this is done by setting $`\overline{E}=2\mu `$ for the mutual scattering of two condensate particles, this being the energy required for them to be excited out of the condensate. By a similar argument Al Khawaja et al.  set $`\overline{E}=\mu `$ for the scattering between a boson in the condensate and a boson in the thermal cloud.

### 2.2 The many-body T-matrix

According to the above argument, a first approximation for the many-body T-matrix at zero momenta and energy is

$$T_n^{MB}(0,0,0;0)=T^{2B}(0,0,0;n\mu )$$
(3)
with $`n`$ = 1 for condensate-thermal cloud scattering and $`n`$ = 2 for condensate-condensate scattering. Further many-body effects can enter the scattering processes attended by the presence of a condensate from attributing a Bogoliubov spectrum to the intermediate states . At low momenta the many-body T-matrix is the solution of the integral equation,

$`𝐤^{}|T^{MB}(E)|𝐤`$ $`=`$ $`𝐤^{}|V(|𝐫_\mathrm{𝟏}𝐫_\mathrm{𝟐}|)|𝐤+{\displaystyle \underset{𝐪}{}}𝐤^{}|𝐕(|𝐫_\mathrm{𝟏}𝐫_\mathrm{𝟐}|)|𝐪`$ (4)
$`\times {\displaystyle \frac{1+2N(\omega _𝐪)}{E2\mathrm{}\omega _𝐪}}𝐪|T^{MB}(E)|𝐤`$
where $`\mathrm{}\omega _𝐪ϵ_𝐪+\mu `$ are the Bogoliubov excitation energies in the Hatree approximation $`\mu ϵ_𝐪`$ and $`N(\omega _𝐪)\{\mathrm{exp}[\beta (ϵ_𝐪+\mu )]1\}^1`$ are the corresponding population factors. For a contact interaction potential Eq. (4) yields

$$T^{MB}(0,0,0;E)=\left[\frac{1}{V_0}\underset{𝐪}{}\frac{1+2N(\omega _𝐪)}{E2\mathrm{}\omega _𝐪}\right]^1.$$
(5)
As discussed by Stoof and Bijlsma , Eq. (5) is still affected by an ultraviolet divergence, which can be remedied by replacing $`V_0`$ in favor of the two-body T-matrix. The final result is

$$T_n^{MB}(0,0,0;0)=T^{2B}(0,0,0;n\mu )\left[1+T^{2B}(0,0,0;n\mu )\underset{𝐪}{}\frac{N(\omega _𝐪)}{\mathrm{}\omega _𝐪}\right]^1,$$
(6)
where $`T^{2B}(0,0,0;n\mu )`$ is given by Eq. (2) with $`\overline{E}=n\mu `$.

### 2.3 Calculation of coupling parameters

The sum in the RHS of Eq. (6) can be evaluated analytically by replacing the sum over intermediate states by an integral over momentum. Setting $`N(\omega _𝐪)=_{s=1}^{\mathrm{}}\mathrm{exp}(s\beta \mathrm{}\omega _𝐪)`$, we get

$$\underset{q}{}\frac{N(\omega _q)}{\mathrm{}\omega _q}=\frac{m}{2\pi \mathrm{}^2}\underset{s=1}{\overset{\mathrm{}}{}}Ei(s\beta \mu )$$
(7)
where $`Ei`$($`x`$) is the exponent-integral function. In the asymptotic low-temperature regime ($`\beta \mu 1`$) this thermal-population term can be approximated by

$$\underset{s=1}{\overset{\mathrm{}}{}}Ei(s\beta \mu )(\beta \mu )^1\mathrm{ln}[1\mathrm{exp}(\beta \mu )].$$
(8)
A numerical illustration is given in Fig. 1 for values of the system parameters appropriate to the experiments on 23Na by Görlitz et al.  (see Sec. 4).

In summary, the coupling parameters in the 2D Bose gas are given by

$$g_n^{MB}=\frac{4\pi \mathrm{}^2/m}{\mathrm{ln}|4\mathrm{}^2/(nm\mu a^2)|2_{s=1}^{\mathrm{}}Ei(s\beta \mu )},$$
(9)
with $`n`$=2 for collisions between pairs of condensate atoms and $`n`$=1 for collisions between an atom in the condensate and a thermally excited atom. These parameters depend on temperature both through the chemical potential and through the excited-state population factor given in Eq. (7) and asymptotically approximated by Eq. (8). If this population factor is dropped, one obtains from Eq. (3) the corresponding “two-body” coupling parameters as

$$g_n^{2B}=\frac{4\pi \mathrm{}^2/m}{\mathrm{ln}|4\mathrm{}^2/(nm\mu a^2)|}.$$
(10)
In the next Section we evaluate the chemical potential and hence the coupling parameters in Eqs. (9) and (10) through a self-consistent evaluation of the density profiles in a two fluid-model.

## 3 Equilibrium properties in the two-fluid model

Let $`n_c(r)`$ and $`n_{nc}(r)`$ be the particle density profiles for the condensate and for the thermal cloud in a 2D Bose gas which is radially confined inside a isotropic planar trap described by the external potential $`V_{ext}(r)=m\omega ^2r^2/2`$. The two-fluid model  combines a solution of the Gross-Pitaevskii equation for the condensate with a Hatree-Fock model of the thermal cloud, which is treated as an ideal gas subject to an effective potential $`V_{eff}(r)`$. In the present case,

$$V_{eff}(r)=V_{ext}(r)+2g_1n_c.$$
(11)
As already noted, we are neglecting the collisions between pairs of bosons belonging to the thermal cloud. It is well known that a full numerical solution of the Gross Pitaevskii equation can be avoided when the kinetic energy term in it can be neglected . This yields the condensate density in the Thomas-Fermi approximation as

$$n_c(r)=(1/g_2)[\mu V_{ext}(r)2g_1n_{nc}(r)]\theta (\mu V_{ext}(r)2g_1n_{nc}(r)).$$
(12)
Since we are interested in examining qualitative behaviors rather than in attaining high numerical accuracy, we have adopted Eq. (12) for the condensate density and found that discontinuities occurring in the density profiles at the Thomas-Fermi radius can be eliminated by the simple expedient of introducing momentum cut-off in the expression for $`n_{nc}(r)`$. That is, we calculate $`n_{nc}(r)`$ from

$$n_{nc}(r)=\frac{m}{2\pi \mathrm{}^2\beta }\mathrm{ln}\left\{1\mathrm{exp}\left[\beta \left(\mu V_{eff}(r)\frac{p_0^2}{m}\right)\right]\right\},$$
(13)
where we take $`p_0=\sqrt{2mg_1n_{nc}}`$ . The model is then evaluated by solving self-consistently Eqs. (11)-(13) together with the condition that the areal integral of $`n_c(r)`$+$`n_{nc}(r)`$ should be equal to the total number $`N`$ of particles.

Before turning to a presentation of our numerical results, let us point out the wide range of validity of the present model. A mean-field treatment is valid when the diluteness condition $`n_ca^21`$ holds and if the temperature of gas is outside the critical region. With regard to the thermal cloud, no significant differences have been found between the predictions of the Hatree-Fock and Popov approximations in the regime $`n_ca^21`$, except at very low temperature where the thermal cloud is becoming negligible .

## 4 Results and discussion

For a numerical illustration, we have taken the values of particle number, the radial trap frequency, and the scattering length as appropriate for 23Na atoms in the experiment of Görlitz et al.  ($`N=5\times 10^5,\omega =188.4`$Hz, $`a=2.8`$ nm). Whereas in their experiment collisions are in a 3D regime, we focus on a strictly 2D regime that could be reached experimentally by increasing either the trap anisotropy parameter or the scattering length.

In Fig. 2 we compare the temperature dependence of the coupling parameters $`g_n^{2B}`$ (long dashes for $`n`$=2 and short dashes for $`n`$=1) with that of the coupling parameters $`g_n^{MB}`$ (full and dotted lines, respectively). It is evident that the many-body screening of the interactions due to the occupancy of excited states is quite large and rapidly increasing with temperature.

Such many-body screening has, however, very little effect on equillibrium properties of the gas for our choice of system parameters. Figure 3 reports the condensate fraction $`N_0/N`$ for the $`g_n^{MB}`$ model as a function of temperature, in comparison with that of an ideal Bose gas at the same values of the system parameters. As is well known, the transition temperature and the condensate fraction are lowered by the interactions. However, the $`g_n^{2B}`$ model gives results that are practically the same as the $`g_n^{MB}`$ one.

Finally, the panels in Fig. 4 show the evolution of the density profiles for the condensate and for the thermal cloud with increasing temperature from near absolute zero to the critical temperature $`T_c0.95T_0`$ with $`T_0=(\sqrt{6}/\pi )\mathrm{}\omega \sqrt{N}`$ being the critical temperature of the ideal gas . The results in Fig. 4(a) are in good agreement with those of Tanatar et al. , except that the tails of the profile are missed in the Thomas-Fermi approximation. Again the screening of collision from the occupancy of the excited states is very small and becomes barely visible at $`T`$ 0.75 $`T_0`$.

In conclusion, the use of Eq. (3) to describe the many-body effects in two-body scattering processes in 2D Bose-condensed gas appears to be very good in regard to equillibrium properties at large values of particle number. A decrease in the number of particles lowers the chemical potential and may lead to observable effects for $`N`$ 103.

This work was partially supported by INFM under the project PRA-Photonmatter. Two of us (K.K.R. and M.P.T.) wish to thank the Condensed Matter Section of the Abdus Salam International Center for Theoretical Physics for their hospitality during the preparation of this work.

Figure captions

Figure 1. The correction term $`A(\beta \mu )`$=$`_{n=1}^{\mathrm{}}Ei(n\beta \mu )`$ from excited-state occupancy (dashed line) and its approximate form from Eq. (8) (solid line) as functions of $`\beta \mu `$. In the inset a zoom of the region $`0.2<\beta \mu <1`$ is shown.

Figure 2. Interaction strengths (in units of $`\mathrm{}^2/m`$) as functions of temperature $`T`$ (in units of ideal-gas critical temperature $`T_0`$).

Figure 3. Condensate fraction $`N_0/N`$ as a function of temperature $`T`$ (in units of $`T_0`$) from the $`g_n^{MB}`$ model and (dashed line) compared with the non-interacting gas (solid line).

Figure 4. Density profile of the condensate (solid line) and the thermal cloud (dashed line) in the $`g_n^{MB}`$ model (in units of $`N/a_{ho}^2`$ with $`a_{ho}=(\mathrm{}/m\omega )^{1/2}`$) as functions of radial distance $`r`$ (in units of $`a_{ho}`$) at various values of the temperature ($`T/T_0`$=0, 0.3, 0.5, 0.8 and 0.95 from (a) to (e)).


format: text

text:
# Magnetorotational Instability in a Couette Flow of Plasma

## 1 Introduction

A central problem in the theory of accretion disks in astrophysics is understanding the fundamental mechanism of angular momentum transfer. A robust anomalous outward angular momentum transport must operate in order for accretion to occur Frank, King & Raine (2002). A phenomenological theory of turbulent angular momentum transport ($`\alpha `$-disks) proposed by Shakura Shakura (1972) was in the basis of our understanding of how accretion occurs and still remains viable to date. The puzzle of the origin of turbulence in hydrodynamically stable disks was resolved when magnetorotational instability (MRI), originally discovered in Refs. Velikhov (1959) and Chandrasekhar (1960), was applied to accretion disks by Balbus and Hawley Balbus & Hawley (1991). MRI causes MHD turbulence to develop in an initially weakly magnetized fluid.

Despite its importance, MRI has never been observed in the laboratory. Recently, two experiments have been proposed to test MRI in a differentially rotating flow of liquid metal (Couette flow) between two rotating cylinders Ji, Goodman & Kageyama (2001); Noguchi et al. (2002). A great deal of theoretical work on investigating MRI in Couette flow of liquid metals has been done confirming that MRI can be excited for magnetic Reynolds number, Rm, exceeding a few Goodman & Ji (2002); Rüdiger & Zhang (2001); Rüdiger et al. (2003); Willis & Barenghi (2002). A particular attention was given to why MRI was not observed in previous experiments with hydromagnetic Couette flows of liquid metals Goodman & Ji (2002); Rüdiger et al. (2003). Because of the very small magnetic Prantdl number of metals, $`\text{Pm}10^6`$, the rotation needs to be very fast to achieve $`\text{Rm}>1`$. Indeed, kinetic Reynolds number needs to exceed $`\text{Re}=\text{Rm}/\text{Pm}10^6`$. In previous experiments the rotation speed was not high enough to achieve $`\text{Re}10^6`$. For such high Reynolds numbers the flow in the experiment is likely to become turbulent even without a magnetic field for Rayleigh stable rotation profiles. The instability may have a nonlinear hydrodynamic nature Richard & Zahn (1999) or may be due to the Ekman circulation induced by the end plates Noguchi et al. (2002). The presence of such turbulence can affect the conditions for the excitation of MRI Noguchi et al. (2002).

Here we consider plasma as alternative to liquid metals to use in MRI experiment. By changing temperature and density of the plasma by a few times around values $`T5\text{eV}`$ and $`n10^{14}\text{cm}^3`$ one can vary Pm in the range of about $`10^2`$ to $`10^2`$. For velocities of plasma of the order of the thermal speed of ions and the typical size of the apparatus of about $`50\text{cm}`$ Rm is in the range from $`10^2`$ to $`10^3`$, while Re is in the range from $`5`$ to $`10^4`$. These plasma parameters are readily achievable in the laboratory Wang et al. (2002). Therefore, in plasma experiment it can be relatively easy to have high enough Rm allowing for MRI to grow, while keeping Re modest and the flow laminar. Laminar character of the flow can allow a detailed study of the structure of the unstable mode and the secondary flow without intervening noise from turbulence. For higher $`\text{Re}10^4`$, the transition from laminar to turbulent flows in hydromagnetics can be investigated. Moreover, the effects of wide variations of Pm can be studied with a plasma MRI experiment.

## 2 Description of the Experiment

The basic setup of possible plasma MRI experiment is illustrated in Fig. 1. This experiment is now under construction at Los Alamos National Laboratory. It can be also used for observing laminar plasma dynamos Wang et al. (2002). Tentative set of specifications and typical parameters of plasma for that experiment is shown in Table 1. The plasma is produced by the electric discharge in hydrogen and is injected into the space between two coaxial conducting cylinders. Mean free path for Coulomb collisions is much smaller than the radii of the cylinders, and MHD description of plasma is appropriate. Cylindrical coordinate system $`r`$, $`\varphi `$, $`z`$ is used with the axis of symmetry coinciding with the central axis of the cylinders. Axial uniform magnetic field $`B_{0z}`$ is produced by the coils around the cylinders and is preexistent of the discharge. Therefore, plasma is coupled to the axial magnetic field lines already at the formation and can slide along the axial magnetic field lines filling the space between two cylinders. The outer cylinder is grounded and the voltage $`\mathrm{\Phi }_0`$ is applied to the inner cylinder during and after the plasma injection, in order to maintain the rotation. Neglecting the boundary and sheath modification of any applied electric field, plasma rotates with $`𝐄\times 𝐁`$ drift,

$$𝐄_0=\frac{\mathrm{\Phi }_0}{\mathrm{ln}(R_2/R_1)}\frac{1}{r}𝐞_r=\frac{𝐕_0\times 𝐁_0}{c},$$
(1)
where $`𝐄_0`$ is the equilibrium electric field, $`𝐁_0=B_{0z}𝐞_z`$ is the equilibrium axial magnetic field, $`R_1`$ is the radius of the inner cylinder, $`R_2`$ is the radius of the outer cylinder, and $`c`$ is the speed of light.

Since plasma is highly collisional and is thermalized before it reaches to the experimental region of the chamber, we further make an assumption that the stationary state is isothermal:

$$\mathrm{\Gamma }P_0/\rho _0=C_s^2=\text{constant},$$
(2)
where $`P_0(r)`$ is the equilibrium pressure, $`\rho _0(r)`$ is the equilibrium density, $`C_s`$ is the sound speed, and $`\mathrm{\Gamma }`$ is the ratio of specific heats ($`\mathrm{\Gamma }=5/3`$ for hydrogen).

The Larmor radius of electrons in the experiment is also much smaller than the radii of the cylinders, but is comparable to the Coulomb mean free path. Despite this, here we assume for simplicity that the conductivity of the plasma is isotropic and is given by non-magnetized expression. The coefficient of magnetic diffusivity, $`\eta `$, and dynamic viscosity coefficient $`\rho \nu `$ are independent of density and depends only on temperature. Then, $`\eta `$ and $`\rho \nu `$ remain constant throughout the volume of the plasma in the isothermal approximation considered here.

In the equilibrium state of the steady rotation the magnetic field is $`B_{0z}=\text{constant}`$ and the velocity is $`𝐕_0=(0,r\mathrm{\Omega }(r),0)`$, with the angular velocity profile

$$\mathrm{\Omega }(r)=\frac{\mathrm{\Phi }_0}{B_{0z}\mathrm{ln}(R_2/R_1)}\frac{c}{r^2}.$$
(3)
Balance between the centrifugal force and the gradient of pressure in the equilibrium state is

$$r\mathrm{\Omega }^2\rho _0=\frac{dP_0}{dr}.$$
(4)
Using isothermal condition, Eq. (4) can be solved to give the radial dependence of $`\rho _0`$

$$\rho _0=\rho _0(R_1)\mathrm{exp}\left[_{R_1}^r\frac{r\mathrm{\Gamma }\mathrm{\Omega }^2}{C_s^2}𝑑r\right]\text{.}$$
(5)
When the rotation speed $`\mathrm{\Omega }(R_1)R_1`$ is comparable to the sound speed, the centrifugal force is significant enough to cause the compression of the plasma toward the outer cylinder.

## 3 The Instability in a Radially Stratified Plasma

We consider perturbations of the equilibrium state described above: $`𝐕=𝐕_0+𝐯`$, $`𝐁=𝐁_0+𝐛`$, $`\rho =\rho _0+\rho _1`$, $`P=P_0+P_1`$, and linearize MHD equations in small perturbations. We idealize the problem by considering cylinders of infinite length and take the dependence of perturbations on $`t`$, $`\theta `$, and $`z`$ as $`\mathrm{exp}[i(\omega t+m\theta +k_zz)]`$. If the equilibrium density $`\rho _0(r)`$ varies significantly between $`R_1`$ and $`R_2`$, general perturbations of the plasma are compressible. However, the phenomenon of MRI is due to the stretching of the magnetic field lines by the differential rotation coupled with the action of centrifugal force. Boussinesq approximation neglects the changes of the volume of the displaced parcel of plasma as it gets quickly adjusted to the new pressure equilibrium at a new radial location in the stratified plasma. It allows to capture centrifugal force acting on a displaced parcel of plasma but excludes compressible modes, which are not essential for the development of MRI. This simplifies calculations substantially. Here we adopt Boussinesq approximation. Linearized continuity equation becomes equivalent to two equations:

$`{\displaystyle \frac{\rho _1}{t}}+(𝐕_0)\rho _1+𝐯\rho _0=0\text{,}`$ (6)
$`𝐯=0\text{.}`$ (7)
Other linearized MHD equations are:

$`𝐛=0\text{,}`$ (8)
$`{\displaystyle \frac{𝐛}{t}}=\times (𝐯\times 𝐁_0)+\times (𝐕_0\times 𝐛)+\eta ^2𝐛\text{,}`$ (9)
$`\rho _0{\displaystyle \frac{𝐯}{t}}+\rho _0(𝐯)𝐕_0+\rho _0(𝐕_0)𝐯+\rho _1(𝐕_0)𝐕_0=P_1+\rho _0\nu _0^2𝐯`$
$`{\displaystyle \frac{1}{4\pi }}(B_{0z}b_z)+{\displaystyle \frac{1}{4\pi }}(𝐛)𝐁_0+{\displaystyle \frac{1}{4\pi }}(𝐁_0)𝐛\text{,}`$ (10)
where $`\nu _0`$ is the unperturbed value of kinematic viscosity. Note, that $`\rho \nu =\rho _0\nu _0`$, so that the viscosity coefficient $`\rho \nu `$ remains unperturbed. The following reductions are done. First, we note that any solution of Eq. (9) satisfies Eq. (8) automatically. This means that out of four equations provided by (8) and (9) together, one should be omitted. We choose to omit $`z`$-component of Eq. (9) and use Eq. (8) to express $`b_z`$ via $`b_r`$ and $`b_\theta `$ in the $`r`$ and $`\theta `$ components of Eq. (9) as well as in all components of Eq. (10). We also use Eq. (7) to express $`v_z`$ via $`v_r`$ and $`v_\theta `$. Eq. (6) is used to express $`\rho _1`$ via $`v_r`$ and substitute into momentum Eq. (10). As for Eq. (10) we replace its $`z`$-component by the equation obtained by taking the divergence of Eq. (10). Using Eqs. (8) and (7) the divergence of Eq. (10) can be reduced to the equation with the only second order radial derivative being the $`^2\mathrm{\Pi }/r^2`$, where $`\mathrm{\Pi }=P_1+B_{0z}b_z/(4\pi )`$ is the perturbation of the total pressure. In summary, we obtain five second order linear differential equations in $`r`$ for five variables, $`\mathrm{\Pi }`$, $`b_r`$, $`b_\theta `$, $`v_r`$, and $`v_\theta `$: the divergence of Eq. (10), $`r`$ and $`\theta `$ components of induction Eq. (9), and $`r`$ and $`\theta `$ components of momentum Eq. (10)

An ideal conductor boundary is a good approximation, since $`\eta `$ for metallic walls is $`<800\text{cm}^2\text{s}^1`$, which is always much smaller than that of plasma, so even a thin metallic wall can be considered as a good conductor. Boundary conditions for conducting walls are given by $`b_r=0`$, $`db_\theta /dr+b_\theta /r=0`$, $`v_r=0`$, $`v_\theta =0`$, $`dv_r/dr=0`$ at both $`r=R_1`$ and $`R_2`$.

We discretize the equations and boundary conditions on a uniform grid in rescaled variable $`x=\mathrm{ln}(r/R_1)`$ with $`N=200`$ grid points. The problem is reduced to finding values of $`\omega `$ such that the determinant of the matrix $`5N\times 5N`$ is zero. We calculate the determinant by using LU-decomposition and search for zeros by using Newton iteration method. After an eigenvalue of $`\omega `$ is found, the corresponding eigenfunction is calculated by using back substitution.

We find $`m=0`$ axisymmetric unstable mode in a wide range of parameters achievable in the experiment. Let us focus on one case with fixed temperature and density given in Table 1. Then, $`C_s(T)`$ and $`\text{Pm}(\rho ,T)`$ are also fixed. The remaining parameters, which can be adjusted in the experiment, are $`\mathrm{\Phi }_0`$ and $`B_{0z}`$. For angular velocity profile (3) the value of $`\mathrm{\Phi }_0`$ is uniquely related to the value of $`\mathrm{\Omega }_1=\mathrm{\Omega }(R_1)`$. The velocity of the plasma should be subsonic everywhere for Boussinesq approximation to be a reasonable approximation. At supersonic speeds, plasma is pushed close to the outer wall leaving very rarefied region in the bulk between cylinders. In this case, it is more appropriate to discuss buoyant or Parker instabilities than MRI. The velocity of rotation is maximum near the inner cylinder, therefore, our calculations are valid for $`\mathrm{\Omega }_1R_1<C_s`$. The maximum value $`\mathrm{\Omega }_{1\mathrm{m}\mathrm{a}\mathrm{x}}=C_s/R_1`$ is given in Table 1.

We plotted the MRI growth rate $`\text{Im}(\omega )/\mathrm{\Omega }_1`$ for $`m=0`$, $`k_z=\pi /L`$ mode in Fig. 2. We also searched for the unstable non-axisymmetric modes, but did not find any. In a small region to the right in Fig. 2 with small contour separation the mode is purely growing ($`\text{Re}(\omega )=0`$), in the rest of the unstable region the mode is growing and oscillating ($`\text{Re}(\omega )>0`$). The MRI is absent for weak magnetic fields ($`v_{\mathrm{A}z}/C_s10^2`$) and small wave numbers ($`k_zL<2\pi `$) because driving force from magnetic field is too weak to overcome the viscosity and diffusivity. The MRI is also suppressed in high magnetic field region due to magnetic tension. For a given magnetic field, only finite number of modes in the $`k_z`$ direction will be excited.

The drawback of using plasmas is limited time available before plasma recombines enough to significantly reduce its coupling with the magnetic field. During the confinement time of the order of $`10^3\text{s}`$ initial perturbations grow by $`10`$ times, and can be registered experimentally. This confinment time is within the range of planned experiment.

The authors are grateful to Zhehui Wang of Los Alamos National Laboratory for encouraging the present work and support. Discussions with Hantao Ji and Jeremy Goodman were very useful. KN acknowledges support from DOE through the LDRD-ER program at Los Alamos National Laboratory. VP acknowledges support from DOE grant DE-FG02-00ER54600.


text:
# 1 Introduction

## 1 Introduction

One of very few places where one has realistic hopes of detecting fundamental physics at or near the Planck scale is in the Cosmological Microwave Background Radiation (CMBR). The reason for this is that inflation expands microscopic fluctuations into macroscopic seeds for structure formation, and thereby possibly imprinting trans-Planckian physics on large scales. For some recent reviews with references see . Unfortunately our understanding of Planckian or stringy physics, and the way inflation is embedded into a fundamental theory, is not sufficiently well developed to allow for reliable predictions of what to be expected.

Various models of possible trans-Planckian physics have been proposed indicating the rough nature of possible effects. In  it was argued that quite independently of the exact mechanism, one can expect a modulated primordial spectrum with an amplitude linear in $`H/\mathrm{\Lambda }`$. The main idea was to mimic the unknown physics in terms of a non-standard vacuum choice at the fundamental scale. This is quite natural since the vacuum is not unique in an expanding universe. In  this was phrased in terms of a (complex) one parameter family of de Sitter invariant vacua, called the $`\alpha `$-vacua, .

The introduction of the $`\alpha `$-vacua has made it possible to analyze the trans-Planckian proposal in a systematic way. There has been a lively debate  as to whether a departure from the Bunch-Davies (or thermal) vacuum is at all possible and whether field theoretical arguments put severe constraints on physics that is expected to be inherently Planckian. In this paper we will not address this issue, but instead focus on the more general question of to what extent non-thermal departures from the Bunch-Davies vacuum are possible.

The question we want to discuss is whether, and how long, non-thermal traces of physics before a possible beginning of inflation can survive and affect the CMBR. A similar question was asked in  where non-thermal initial conditions were introduced at a specific time and their effects on the CMBR analyzed. Technically, as will be explained in the next section, the analysis is a direct parallel to what was done in  even though the physics is slightly different.

It is reasonable to expect that any non-thermal deviation from the Bunch-Davies vacuum will experience thermalization due to the Hawking radiation present during inflation. In this paper we are in particular interested in the actual thermalization time $`\tau `$, that is, the maximal time non-thermal deviations can survive in de Sitter space. We will argue, in two distinct but possibly related ways, that this time is of order $`\tau R^3/l_p^2`$, where $`R`$ is the horizon radius. This is substantially longer than the naive expectation of a time scale of the order $`1/T`$.

This, one could argue, is generic for quantum systems inside isolated boxes and does not have anything particularly to do with de Sitter space or Hawking radiation. With this in mind, we discuss in more general terms the meaning of the relaxation time mentioned above. In particular we try to relate it to another typical time scale relevant for finite isolated systems, namely the Poincaré recurrence time.

The plan of the paper is as follows. In section 2 we briefly review the trans-Planckian problem. Particular attention is given to generic signatures resulting from the assumption that something special happens at a definite scale, above which new physics is supposed to take over. The “transverse” case, where instead something special happens at a definite time during the rolling of the inflaton field, is also discussed. In section 3 the issue of thermalization is addressed. As already mentioned, the typical thermalization time for non-thermal excitations in de Sitter space is argued to be of order $`\tau R^3/l_p^2`$. The section ends by attempts to put this into a more general context. We conclude in section 4.

## 2 Signatures of new physics in the CMBR

### 2.1 Modified initial conditions on a fixed scale

It has been argued in several works that physics at or beyond the Planck scale (or string scale) might leave a non negligible imprint in the CMBR. The basic argument behind the claim is that the rapid expansion of the universe during inflation magnifies small scale physics and makes it accessible on large scales. Going back in time the quantum fluctuations responsible for the structure in the CMBR eventually reaches a linear scale of the order the Planck scale and one would expect that Planckian physics determine their form. In fact, in an expanding universe the notion of a vacuum is not unique. The standard procedure is to ignore the presence of the Planck scale and trace the quantum state back to the infinite past. In an accelerated universe (where the scale factor grows faster than the Hubble radius) the linear scale of the fluctuations can be made arbitrarily smaller than the Hubble scale and the vacuum becomes unique. In exact de Sitter space the vacuum is called the Bunch-Davies vacuum and happens to be thermal.

There is clearly something unsatisfactory with this procedure. A more conservative approach was advocated in  where it was argued that one should stop when the modes reached the Planck scale (or string scale) and impose initial conditions there. It was shown that a particular choice of initial conditions (the instantaneous Minkowski vacuum) leads to a corrected primordial spectrum of the form

$$P(k)=\left(\frac{H}{\stackrel{}{\varphi }}\right)^2\left(\frac{H}{2\pi }\right)^2\left(1\frac{H}{\mathrm{\Lambda }}\mathrm{sin}\left(\frac{2\mathrm{\Lambda }}{H}\right)\right).$$
(2.1)
It was also argued that the linear dependence on $`H/\mathrm{\Lambda }`$ in the correction term is rather generic. If, furthermore, we are not in exact de Sitter space but instead experience a slow roll where $`H`$ is changing with time (which is translated into scale on the sky) the primordial spectrum will have a characteristic modulation. See also . The precise amplitude and phase of the modulation depends in a detailed way on the initial conditions chosen, but the modulation is a generic effect and a typical sign of trans-Planckian physics. It is amusing to compare the effect with the acoustic peaks which are due to oscillations after the fluctuations have reentered the horizon. The modifications we are studying have their origin in oscillations before the fluctuations exit.

### 2.2 Modified initial conditions at a fixed time

In  an analysis very similar to the one in  was performed. The idea was to impose initial conditions different from the Bunch-Davies vacuum not at a fixed scale at all times, but at a fixed time on all scales. The main difference is that the fixed scale $`\mathrm{\Lambda }`$ in eq. (2.1) should be replaced by the physical scale of the mode at the initial time. This way of imposing initial conditions breaks de Sitter invariance and gives rise to a different kind of modulation. There is also a numerical factor, as explained in , depending on the change of the slowroll parameters at the initial time. It is useful to draw a diagram which illustrates the two possibilities (see Fig. 1).

In  it was argued that one should choose the Bunch-Davies vacuum at the fundamental scale. As a consequence, as argued in , any departure from the Bunch-Davies vacuum will be washed away on small scales. It is only on the very largest scales that traces of the initial conditions remain, and if inflation has proceeded longer than the minimal 60 or 70 e-foldings, the scales where the imprints are located has not yet reentered (and never will if the universe now starts to accelerate again).

The aim of the next section is to argue that traces of non-thermal physics much further back can remain hidden in the trans-Planckian regime and manifest themselves as modified initial conditions at the fixed scale.

## 3 Time scales of thermalization

### 3.1 Thermalization in de Sitter space

Having a non-thermal state in a thermal environment suggests that the construction has a limited lifetime. The question we will try to answer in this section is how long this lifetime really is. That is, for how long can inflation, with a non-thermal vacuum, continue before the vacuum relaxes and becomes thermal? A naive first estimate would suggest that the thermalization time in de Sitter space would be related to the temperature and be given by $`1/T`$, where $`T`$ is the Gibbons-Hawking temperature . It is easy to see, however, that this crucially depends on the strength of the interactions that are causing the thermalization. In an expanding universe the interactions must be faster than the expansion of the universe for the thermalization to happen in this way. Since the time scale of the expansion is $`H^1T^1`$ this depends in a detailed way on the actual interactions. In a slow-roll inflation, with its flat potential, the interactions are not fast enough and thermalization will not take place.111Note that this assumption is implicit also in  when it is argued that physics can remember initial conditions at a fixed time all the way back through many e-foldings. To summarize, after the fluctuations have grown larger than the horizon, there is no causal process which could affect them. As a consequence the fluctuations freeze and any non-thermal signatures are kept intact.222One way to see that a thermalization time of order $`1/TR`$ is far too short, is to consider our present stage of cosmic acceleration. If we consider this state as a perturbed de Sitter space, one finds that with a thermalization time of order $`R`$ the lifetime for the perturbations would be roughly $`R10`$ billion years. Now our sun for example (being a small perturbation of the present de Sitter phase) will even $`1000`$ billion years into the future still exist as a small perturbation, by then in shape of a cooling white dwarf.

This global perspective could, however, be misleading from the point of view of a local observer . It is important to note that taking a global perspective of the process ignores possible subtleties having to do with horizons. Since the physics of horizons, especially in de Sitter space, is largely mysterious one needs to be careful when they come into play. After all, the counterpart of black hole complementarity, , is not well understood in de Sitter space (see, e.g. ). Let us therefore take a local perspective and try to figure out what would be observed.

We start with a mode that is created at a definite scale and then expands to larger size. A local observer would see how the mode approaches the horizon but never crosses it. The closer the horizon it gets the more redshifted it becomes. From the local observer’s point of view anything approaching the horizon effectively enters into a region with ever higher temperature and will at some point in time be completely “thermalized” and turn into radiation. It therefore seems like there is indeed a finite lifetime for any excitations above the thermal state.

It is important to note that our first argument, where thermalization was excluded, was made using quantum field theory augmented with Hawking radiation. It is only when we take a further step, invoking physics related to the horizon – eventually described by quantum gravity, holography and complementarity – that we are able to find signs of thermalization.

Before presenting the conditions we believe need to be fulfilled in order to give a correct estimate of the thermalization time, let us comment on a proposal made in , where the same issue was discussed. There thermalization also was invoked based on effects appearing only once the physics of horizons was taken into account. The condition the authors of  claim is determining the thermalization time is the time it takes an object (initially starting close to the center of the static patch) to be of order a Planck length away from the horizon. Indeed, from the local observers point of view the temperature at that point is of order the Planck scale, naively implying that gravitational interactions quickly thermalize the object in question. This estimate of the thermalization time is easily shown to give $`\tau R\mathrm{ln}R`$.

We believe this way of arguing suffers from a number of different problems. First of all, it only concerns objects approaching the observers horizon, it has nothing to say about the observer herself and things bound to her. In particular, the time it takes for her to be thermalized, using this criterion, would be infinite. This is clearly wrong. In order to overcome this problem one could argue that, somehow, the time it takes for local objects to thermalize is also of order $`R\mathrm{ln}R`$. This, however, must also be wrong, as can be concluded from footnote 2 (and noting that $`\mathrm{ln}R`$ only gives a factor of order 100 for our present universe).

A second, perhaps more serious, objection concerns the actual validity of the argument even for objects approaching the horizon. If we focus on what an observer actually would see as an object recedes towards the horizon, one needs to be careful. Since the rate of the photons (emerging from the horizon) received by our observer is of order $`1/R`$ (see ), the time it would take for her to see the object being thermalized must be a much larger time than $`R\mathrm{ln}R`$, during which time only of order $`\mathrm{ln}R`$ photons could have been detected. In fact, the time it will take is of order $`R^3`$, as will be shown below.

The reason we find it important to focus on what actually would be seen is that, just like in inflation, one can imagine that the de Sitter phase is abruptly turned off, being followed by a more standard cosmological evolution. When this happens, the object will, of course, return to the observers causal patch at some time in the future . Since the situation between the object and the observer is symmetric, it is clear that the object will be in as good shape as the observer. If the observer has not been thermalized by then, then neither should the object be, implying that the estimated thermalization time in  is far too short.

Indeed, considering the symmetric situation we have between the observer and the object (being for example another observer) and the fact that they can meet again some time after the de Sitter phase has turned off, seems to imply that the estimated thermalization time should be the same for local objects as for those who approach the horizon, even from the perspective of one single observer.

With this last observation in mind let us make two independent estimates of the thermalization time, one considering objects approaching the horizon, and one considering local object, bound to the observer. We will find that they are indeed of the same order.

As mentioned above, in order to decide how fast an object (falling towards the horizon) is being thermalized, one should focus on what is actually seen by a local observer. From this point of view the thermalization process will be very slow considering the fact that the rate of photons received by the observer is of the order $`1/R`$, i.e. one photon every $`R`$ . Let us now try to find out what actually happens to the object (according to the observer). To do that we think of the horizon as an area consisting of $`R^2/l_{pl}^2`$ Planck cells, and remember that the photon has a wavelength of order Planck scale when emitted and can indeed resolve specific Planck cells.

Now first assume the object in question is something really simple, corresponding to an information content much smaller than the $`R^2`$ number of degrees of freedom of the horizon. This would mean that only a few of the Planck cells are involved in encoding the object. In the extreme case of an object with entropy of $`𝒪(1)`$, one would need to wait until of the order $`R^2`$ photons has been emitted to be reasonably sure to see a photon coming from the burning of the object. In the other extreme one can think of an object consisting of the order $`R^2`$ degrees of freedom. In this case it is clear that one has to wait until of the order $`R^2`$ photons has been emitted, in order for all parts of the object to have been burnt. And so regardless of the size of the object, one has to wait a time,

$$\tau \frac{R^2}{l_{pl}^2}R=\frac{R^3}{l_{pl}^2},$$
(3.1)
in order to actually see the destruction.333We note that precisely this time scale was found in , where in that paper $`R^3`$ could, in the context of black holes, be interpreted as the time, after some object had fallen into the black hole, an observer would need to wait in order to recover the corresponding information. At this point one could object and say that since the blueshift is so strong near the horizon, anything close to it would thermalize almost immediately, leading to a much shorter thermalization time. However, taking the blueshift into account also requires taking the time dilation into account. Indeed, this makes processes happening close to the horizon appear extremely slow. This is one way to see that the large blueshift is perfectly consistent with the relatively long time scale in eq. (3.1).

We therefore suggest that the maximal lifetime of non-thermal excitations in de Sitter space is given by $`\tau R^3/l_{pl}^2`$. It is important to emphasize that this relatively long time scale follows only if other non-gravitational interactions are frozen. If this is not the case, and those interactions occur faster than the expansion of the universe, the characteristic thermalization time will scale like the naive $`1/T`$. In that case we do not need any reference to holography or complementarity. Thermalization occurs regardless of perspective.

Now let us try to estimate the time it takes to thermalize a local object, bound to the observer. Let us reconsider the possibility that local interactions do give rise to a thermalization, but only if we take physics near the Planck scale into account. Again it is the Hawking radiation that we expect can do the job. The probability that the Hawking radiation, with wavelength of order $`R`$, will interact (and thereby thermalize) with physics at a fundamental length scale, say at the Planck length, is small, but obviously nonzero. The question is, therefore, how long do we need to wait in order for the probability of scattering of a given Planck cell with the de Sitter radiation to be of $`𝒪(1)`$? The claim is that this will give a rough estimate of the thermalization time of the underlying microphysics. The interaction rate is $`\mathrm{\Gamma }=\sigma nv`$, where for this particular process we have for the cross section $`\sigma l_{pl}^2`$, the number density of the radiation $`nT^31/R^3`$ and the relative velocity $`v=c=1`$. Then the time $`\tau `$ it takes for this process to be likely is given by the condition $`\mathrm{\Gamma }\tau 1`$ implying,

$$1\sigma nv\tau l_{pl}^21/R^31\tau \text{ }\text{ }\tau \frac{R^3}{l_{pl}^2},$$
(3.2)
which is seen to coincide, up to orders of one, with the previous result. Therefore, regardless of whether local objects or objects falling towards the horizon are concerned, the thermalization time will be the same. We argued above that this must be the case based on the symmetry between the observer and the object and by noting that, if the de Sitter phase is only temporary, they will meet again. We find it encouraging that the above results are in agreement with this assessment.

To summarize, we have made three attempts to estimate the thermalization time. The one where the effects of horizons was ignored led to an apparent contradiction with an estimate based on horizons and holography. A reconciliation of these two estimates need an understanding of complementarity. A reasonable conclusion is, though, that $`\tau `$ represents the longest time that a non-thermal deviation can survive from a local perspective.

We also found that the longer time scale can be obtained from Hawking radiation acting through Planck scale physics. It would be interesting to understand this connection better. The Planck scale argument suggests that the longer time scale could have global importance by affecting the initial conditions. That is, the vacuum, at Planck scale, would be thermalized after a time $`\tau `$.

Now let us see what the implications of these results are for inflation. In inflation the Hubble constant is constrained from observations to be no larger than $`H10^4\text{ }m_{pl}`$. With this input the thermalization time for non-thermal excitations ($`\alpha `$-vacua included) is found to be of order $`\tau R^3=1/H^310^{12}\text{ }t_{pl}`$. Comparing this with the time needed for the required number of e-foldings, which for 70 e-foldings is $`t_{infl}70/H710^5\text{ }t_{pl}`$, one can conclude that the thermalization time allows for visible effects of non-thermal behavior in the CMBR, with room to spare. On the other hand, if we imagine that inflation went on uninterrupted for a very long time, thermalization effects become important.444More generally, the thermalization time, in units of e-foldings, relevant for initial conditions on a scale $`\mathrm{\Lambda }`$ can be expected to be $`\left(\mathrm{\Lambda }/H\right)^2`$. If the desired amplitude for the corrections to the primordial spectrum is $`10^2`$ the number of allowed e-foldings become $`10^4`$.

### 3.2 General considerations of thermalization in cosmological like boxes

Let us investigate the thermalization a bit further. The time scale is in general given by

$$\tau \mathrm{\Gamma }^1E^2/T^3.$$
(3.3)
In the expanding de Sitter space, with a constant Hawking temperature, the expanding modes redshift like $`Ee^{Ht}\epsilon `$. Cutting of at the horizon scale, where everything freezes, we find $`ET`$ and as a consequence the naive $`\tau 1/T`$. A slightly more detailed calculation focuses on the integral $`𝑑t\text{ }\mathrm{\Gamma }`$. It is only if this integral diverges before horizon crossing – i.e. there is time for an infinite number of interactions – that we get exact thermalization. In our case we find

$$^{t_{cross}}𝑑t\text{ }\mathrm{\Gamma }\frac{T^3}{H\epsilon ^2}(e^{2Ht_{cross}}1)T/H1$$
(3.4)
where we have used $`e^{2Ht_{cross}}\epsilon H\epsilon /T1`$. Note that we have taken a global perspective for our argument. The only way for the integral to become large is if the interactions are truly strong, unlike what is the case in inflation, giving rise to a large dimensionless prefactor in the above expression.

A main point of our paper has been to argue that there could be another way to get thermalization. We first argued for this from a holographic point of view. We then found a local derivation which gave the same result. Let us now proceed along the second line of approach to see whether we can learn something more. We assume that Planckian physics can be modelled by a fixed $`E`$ of order $`\mathrm{\Lambda }`$ and that the Planckian physics not only affects what is going on at lower energies, but also suffers a possibly thermalizing back reaction from low energy physics. Under these assumptions the integral trivially diverges and the system, as we have argued in the previous subsection, eventually thermalizes. To be slightly more general, we can consider a temperature which is not constant (unlike the Gibbons-Hawking temperature in de Sitter space), but instead redshifts with the expansion of the universe. We then find

$$𝑑t\text{ }\mathrm{\Gamma }𝑑t\text{ }\frac{T^3}{\mathrm{\Lambda }^2}𝑑t\frac{1}{a^3\mathrm{\Lambda }^2},$$
(3.5)
which diverges only if the scale factor grows slower than $`t^{1/3}`$. For this to be the case, one needs an equation of state $`p=\sigma \rho `$ where $`\sigma >1`$ which is unphysical. In conclusion we find that we must put the system in a box that prevents the temperature from redshifting. This is precisely what happens in de Sitter space.555It is amusing to speculate whether there are astrophysical systems where the time scale eq. (3.1) could be important. In this context one may note that the hottest stars, with internal temperatures of the order 1 GK (or $`10^5`$ eV) give rise to a time scale of $`10^{17}`$ years, which clearly is much too long to be of any interest. However, there are several systems emitting radiation in the TeV range which translated into temperature becomes $`10^{16}`$ K. If this is inserted into eq. (3.1), one finds a relaxation time of about an hour. If our reasoning is correct, one could speculate that in a region where TeV temperatures are sustained for an hour or so, effects of quantum gravity could give rise to a complete thermalization of any physical system (including baryon number violation etc.).

### 3.3 Remarks on time scales in finite isolated systems

One should note that the time scale we have discussed, although long, is very much shorter than the Poincaré recurrence time (see Fig. 2).

The Poincaré recurrence time is the time scale over which macroscopic thermal fluctuations are expected to happen (for recent qualitative discussions on this in the context of de Sitter space, see for example ). The typical time scale is given by $`t_{rec}e^S`$, where $`S`$ is the entropy of the system. The second law is only meaningful near a macroscopic fluctuation and states that time always points towards increasing entropy. Our time scale eq. (3.1) is the maximum decay time (or, by time reversal symmetry, rise time) of a macroscopic fluctuation.666In  it was noted that there is an interesting geometrical connection between these time scales in a space time where inflation suddenly stops.

We suggest that this is a general behavior for finite isolated quantum systems, not limited to de Sitter space. Another such system of particular interest is large black holes in Anti de Sitter space (for discussions on long, and very long, time scales in this context, see for example ). Since these black holes are eternal they are particularly intriguing laboratories to test some of the ideas mentioned above. Interestingly, both time scales (maximal relaxation and recurrence) seems to tell us something about retrieval of information that has fallen into the black hole . The mechanism for this and the meaning of these time scales is, however, still obscure . Even though tempting, we do not to attempt to address these issues here.777It is important to point out though that in the discussion of thermalization in de Sitter space, we mean thermalization strictly in the thermodynamic sense. We do not mean actual loss of coherence at the microscopic level. However, these residual correlations are not likely to be detectable and so effectively the spectrum would be thermal.

## 4 Conclusions

In this paper we have investigated the problem of thermalization in de Sitter space. We have argued, based on holography, that the maximum thermalization time – relevant if all non-gravitational interactions are weak and slower than the expansion of the universe – is given by eq. (3.1).888This is provided that the trans-Planckian physics can be viewed as a closed system and there are no sources of negentropy. We have also argued that this allows trans-Planckian physics to retain a memory of physics preceding inflation for up to $`N=H\tau 1/RR^3=R^2/l_{pl}^2`$ e-foldings, which for $`R10^4\text{ }l_{pl}`$ gives $`N10^8`$.

We have also speculated that the time scale could have a more universal importance and not only be relevant for de Sitter space. In particular there seems to be two typical time scales involved in box-like configurations. One is the maximal relaxation time (which we interpreted in terms of thermalization in de Sitter space), which is the typical decay (or rise) time for macroscopic fluctuations. The other is the recurrence time, which is the typical time scale for which these macroscopic fluctuations are expected to happen.

## 5 Acknowledgments

UD is a Royal Swedish Academy of Sciences Research Fellow supported by a grant from the Knut and Alice Wallenberg Foundation. The work was also supported by the Swedish Research Council (VR).


text:
# STELLAR ORBITS AND THE INTERSTELLAR GAS TEMPERATURE IN ELLIPTICAL GALAXIES

## 1 Introduction

The anisotropy of stellar orbits in elliptical galaxies can be estimated from the temperature of the hot interstellar gas through which they move. These two quite different galactic attributes are intimately related by the Jeans equation for the stars and the condition for hydrostatic equilibrium in the gas. This simple relationship is valuable since both the anisotropy and gas temperature are difficult to extract from optical and X-ray observations, respectively.

Many massive elliptical galaxies are nearly spherical (Merritt & Trembly 1996) and slowly rotating, but their stellar velocity ellipsoids are not in general isotropic. The non-spherical nature of the stellar velocity dispersion is represented by the parameter

$$\beta =1\sigma _t^2/\sigma _r^2$$

where $`\sigma _r`$ is the radial stellar velocity dispersion and $`\sigma _t`$ is the dispersion in a transverse direction, i.e. $`\sigma _t^2=\sigma _\theta ^2=\sigma _\varphi ^2`$. If the orbits are predominantly radial, $`0<\beta <1`$, the line of sight velocity profile becomes more strongly peaked than a Gaussian profile with increasing projected radius $`R`$; if the orbits are mostly tangential, $`\mathrm{}\beta <0`$, the profile is more flat-topped and becomes broader with increasing radius. Both $`\beta (r)`$ and the galactic potential $`\mathrm{\Phi }(r)`$ can be determined from optical observations of the line of sight velocity dispersion $`\sigma (R)`$, the optical surface brightness distribution and the deviation of the stellar line profiles from a Gaussian as expressed by the line-symmetric coefficient $`h_4(R)`$ in a Gauss-Hermite expansion (e,g, van der Marel & Franx 1993).

Stellar line profiles observed in most massive E galaxies indicate a preference for radial orbits with $`\beta 0.3`$ (Bender, Saglia & Gerhard 1994; Gerhard et al. 1998; Saglia et al. 2000). Accurate $`\beta (r)`$ profiles from optical data require high quality data and considerable care in reduction and analysis. The radially anisotropic nature of stellar orbits in luminous elliptical galaxies provides important and otherwise unavailable information about the merger history of these galaxies (Naab, Burkert & Hernquist 1999), so an improved or independent determination of $`\beta (r)`$ would be desirable.

Similarly, the radial variation of the hot gas temperature $`T(r)T_{vir}10^7`$ K in elliptical galaxies depends on the spatial and spectral resolutions of X-ray detectors and the accuracy of the three-dimensional decomposition that converts $`T`$ as a function of projected radius $`R`$ to physical radius $`r`$. In addition, it is often unclear whether the hot gas at any radius has a single temperature or a multitude of temperatures as might be expected if the gas were cooling. It is unclear if the gas cools at all. For example, high resolution X-ray spectra with XMM-Newton of the large E galaxy NGC 4636 fail to show emission lines expected from gas at intermediate temperatures, such as the 0.574 keV OVII line (Xu et al. 2002), suggesting that the gas is not cooling below $`T_{vir}/3`$ as in classical cooling flows. However, Bregman, Miller & Irwin (2001) detected the OVI 1032,1038Å doublet in NGC 4636 emitted from gas at $`T3\times 10^5`$ K, implying that cooling near the expected rate may occur after all. As an additional source of confusion, the observed X-ray spectra can often be significantly improved by assuming two quite different discrete temperatures at each galactic radius (e.g. Buote 2002 for NGC 1399; Buote et al. 2002a and Tamura et al. 2003 for NGC 5044). While the origin and physical nature of gas with only two (or a limited continuum of) temperature phases may be conceptually problematical, it may be related to certain heating mechanisms. Nevertheless, there are good reasons for suspecting the accuracy or interpretation of gas temperature profiles $`T(r)`$ in elliptical galaxies.

In the following discussion we illustrate the close relationship between the radial variations of $`\beta (r)`$ and $`T(r)`$ using data for the well-observed giant elliptical NGC 4472. While currently available data for this galaxy are adequate to illustrate this relationship, more accurate optical and X-ray data will be necessary to verify $`\beta (r)`$ or $`T(r)`$ with confidence.

We recognize that NGC 4472 is not perfectly spherical; its optical image has been variously classified as E1 or E2, However, optical determinations of $`\beta (r)`$ for NGC 4472 and other nearly-spherical elliptical galaxies have been made under the assumption of spherical symmetry (e.g. Kronawitter et al. 2000). Indeed the very concept of the $`\beta `$ anisotropy parameter is only valid in spherical geometry. Similarly, the gas temperature $`T(r)`$ determined from X-ray observations of NGC 4472 and other nearly spherical elliptical galaxies are based on the assumption of spherical symmetry. One justification for this is that the stellar gravitational potential that confines the stars and hot gas within the half-light radius is always more spherical than the underlying stellar luminosity distribution (e.g. Brighenti & Mathews 1996).

Our discussion of the close relationship between stellar anisotropy and gas temperature is a thinly-veiled request for observers to concentrate on those luminous E galaxies for which the most accurate X-ray and optical data can be acquired.

## 2 Relating $`\beta (r)`$ and $`T(r)`$

The Jeans equation for the radial stellar velocity dispersion $`\sigma _r`$ is

$$\frac{1}{\mathrm{}_{}}\frac{d(\mathrm{}_{}\sigma _r^2)}{dr}+\frac{2\beta \sigma _r^2}{r}=g$$
(1)
where $`\mathrm{}_{}`$ is the stellar luminosity density corresponding to the stellar velocity dispersions $`\sigma _r`$ and $`\sigma _t`$. The equation for hydrostatic equilibrium in the hot interstellar gas is

$$\frac{1}{\rho }\frac{d(\rho c^2)}{dr}=g.$$
(2)
Here $`c^2=kT/\mu m_p`$ is the isothermal sound speed, $`m_p`$ is the proton mass and $`\mu =0.61`$ is the molecular weight. The uniformity of the stellar mass to light ratio within the effective radius of NGC 4472 is excellent evidence that the hot gas is in hydrostatic equilibrium (Brighenti & Mathews 1997). Both equations contain the same gravitational term $`g=GM(r)/r^2`$. Eliminating this term results in an equation explicitly independent of the confining mass:

$$\frac{d\sigma _r^2}{dr}+\frac{\sigma _r^2}{r}\left[\frac{d\mathrm{log}\mathrm{}_{}}{d\mathrm{log}r}+2\beta \right]=\frac{dc^2}{dr}+\frac{c^2}{r}\left[\frac{d\mathrm{log}\rho }{d\mathrm{log}r}\right].$$
(3)
This equation has been discussed previously, particularly in the context of the anisotropic orbits of galaxies in rich clusters that contain hot gas (e.g. Fabricant et al. 1989). In addition to hydrostatic equilibrium, the validity of this equation depends on the assumption that the gas pressure exceeds that of any other interstellar pressure such as turbulence, cosmic rays, etc.

If $`\beta (r)`$ and $`T(r)`$ (as well as $`\rho (r)`$ and $`\mathrm{}_{}(r)`$) are known from observations, then Equation (3) can be solved for $`\sigma _r(r)`$ and the line of sight velocity dispersion is then found from

$$\sigma ^2(R)=\frac{_R^{r_t}\sigma _r^2(r)\left[1(R/r)^2\beta \right]\mathrm{}_{}(r)(r^2R^2)^{1/2}r𝑑r}{_R^{r_t}\mathrm{}_{}(r)(r^2R^2)^{1/2}r𝑑r}$$
(4)
(e.g. Binney & Mamon 1982). In this equation $`R`$ is the projected radius and $`r_t`$ is some radius significantly beyond the limit of the available data (we assume $`r_t=3R_e`$). If $`\beta (r)`$ is known, Equation (3) can be solved for $`\sigma _r(r)`$ for various parameterized temperature profiles $`T(r)`$ until $`\sigma (R)`$ from Equation (4) agrees with observed velocity dispersions. By this means it is possible to derive the gas temperature profile $`T(r)`$ from a knowledge of $`\beta (r)`$. Conversely, if the gas temperature $`T(r)`$ is securely known from X-ray observations, then $`\beta (r)`$ can be determined by requiring that solutions $`\sigma ^2(R)`$ of Equations (3) and (4) pass through the observations. This close relationship between $`\beta (r)`$ and $`T(r)`$ is possible and valuable because $`\mathrm{}_{}(r)`$, $`\rho (r)`$ and $`\sigma ^2(R)`$ are much easier to determine from observations than either $`\beta (r)`$ or $`T(r)`$.

Except for the very central region, the stellar brightness in NGC 4472 is well-fit with a de Vaucouleurs profile. The stellar B-band luminosity density $`\mathrm{}_{}(r)\mathrm{}_{,deV}(r)`$ is determined by the total luminosity $`L_B=7.89\times 10^{10}`$ $`L_{B,}`$ and the effective radius: $`R_e=104^{\prime \prime }=1.733^{}=8.57`$ kpc, assuming a distance $`d=17`$ Mpc. However, within a small “break” radius, $`r_b=2.41^{\prime \prime }=200`$ pc, the stellar luminosity flattens to $`\mathrm{}_{,core}(r)=\mathrm{}_{,deV}(r_b)(r/r_b)^{0.90}`$ (Gebhardt et al. 1996; Faber et al. 1997). In Figure 1 we compare the (un-normalized) stellar luminosity density $`\mathrm{}_{}(r)`$ in NGC 4472 with the hot gas density profile $`\rho (r)`$ determined from X-ray observations. The remarkable proportionally $`\mathrm{}_{}(r)n_e^2(r)`$ was first discovered by Trinchieri, Fabbiano & Canizares (1986).

X-ray observations of the hot gas temperature in the central region of NGC 4472 are shown in Figure 2. Only three observations of the gas temperature are available in the region within 10 kpc where $`\beta `$ has been determined from optical observations. The filled circles show ROSAT deprojected gas temperature observations from Buote (2000). The linear temperature profile defined by these two temperature observations

$$T(r)=0.7480\times 10^7+7.251\times 10^5r_{kpc}$$
(5)
is shown in Figure 2 with a solid line. From Chandra data Soldatenkov, Vikhlinin & Pavlinsky (2003) have recently provided an additional temperature measurement of $`T/10^7\mathrm{K}=0.766\pm 0.023`$ at radius 1” = 0.0824 kpc that is shown with a filled square in Figure 2. The linear fit to the Buote data passes within the error bars of the Soldatenkov et al. observation, so all three observations are consistent with a linear temperature profile within the effective radius $`R_e=8.57`$ kpc. Although the X-ray data is sparse in the region of interest, it is in perfect qualitative agreement with linear gas temperature profiles observed in other bright elliptical galaxies in $`r\begin{array}{c}<\hfill \\ \hfill \end{array}R_e`$ (Fig. 1 of Brighenti & Mathews 1997).

### 2.1 Determining $`\beta (r)`$ from $`T(r)`$

The open circles in Figure 3 show observations of the line of sight stellar velocity dispersion in NGC 4472 (from Bender et al. 1994) that were used by Kronawitter et al. (2000) to determine $`\beta (r)`$ from optical observations alone. The filled circles in Figure 3 are the data of Fried & Illingworth (1994). Both sets of observations lie entirely within the effective (half-light) radius, $`R_e=104^{\prime \prime }=8.57`$ kpc. For some reason the normalization is systematically different in the two sets of optical observations. The Fried-Illingworth data have smaller error bars, less overall scatter and resemble in overall form the velocity dispersion profiles $`\sigma (R)`$ of other well-observed ellipticals \[e.g. Gerhard et al. (1998) for NGC 6703\]. The observational uncertainties apparent in the data sets shown in Figure 3 make it difficult to determine with confidence any radial dependence of the stellar anisotropy. Consequently, we assume that the orbital anisotropy $`\beta `$ in Equation (3) is constant over the range $`r>r_o=0.34`$ kpc.

In Figure 3 we compare solutions of Equations (3) and (4) with the observed line of sight stellar velocity dispersion in NGC 4472. Each solution curve $`\sigma (R)`$ is based on the linear gas temperature profile in Figure 2 and a particular value of $`\beta `$. More information about the solutions shown in Figure 3 can be found in Table 1. Of the three solid lines designed to fit the Fried-Illingworth data, solution 2 with initial $`\sigma _r(r_o)=430`$ km s-1 and $`\beta =0.7125`$ fits the data best. Solutions 1 and 3 show additional marginally acceptable fits although they fall respectively below and above the Fried-Illingworth data at smaller projected radii $`R`$. Evidently $`\beta =0.71\pm 0.15`$ is the most likely anisotropy based on the gas temperature observations. In some of our solutions of Equation (3) $`\sigma _r^2(r)`$ becomes unphysically negative at a radius $`r_z`$ less than the outer limit of our integration in Equation (4), $`r_t=3R_e`$; for these solutions we replace $`r_t`$ with $`r_z`$. In solutions 4 - 6, shown with dashed lines in Figure 3, we find that the anisotropy parameter must have values $`\beta =0.63\pm 0.15`$ to be in approximate agreement with the data of Bender et al. (1994). Once $`\sigma _r(r_o)`$ is selected for each of these dashed line solutions in Figure 3, the vertical position of the solution is extremely sensitive to $`\beta `$, i.e., $`\beta `$ must be specified to one part in 1000 to pass through the center of the Bender et al. data. One may doubt if Nature can be so discriminating. Nevertheless, for either set of optical data, the corresponding values of $`\beta `$ are about twice as large as the $`\beta (r)\begin{array}{c}<\hfill \\ \hfill \end{array}0.34`$ found by Kronawitter et al. (2000) from the Bender et al. data.

### 2.2 Determining $`T(r)`$ from $`\beta (r)`$

The preferred orbital anisotropy in NGC 4472 determined by Kronawitter et al. (2000) from the Bender et al. observations can be accurately represented with

$$\beta _K(r)=0.449r_{kpc}^{0.7333}\mathrm{exp}(r_{kpc}/2.472),$$

which peaks at $`r1.65`$ kpc where $`\beta _K0.34`$. In this section we assume that $`\beta _K(r)`$ is correct and solve Equation (3) for $`c(r)^2`$ and $`T(r)`$. As before we consider only linear temperature variations defined by two temperatures $`T_1`$ and $`T_2`$ evaluated respectively at $`r_1=2.18`$ kpc and $`r_2=8.55`$ kpc. We seek values of $`T_1`$ and $`T_2`$ for which solutions $`\sigma (R)`$ of Equation (4) pass through or near the Bender et al. optical data.

We begin with a set of solutions for which $`\beta =\beta _K(r)`$ and the temperature profile is constrained to agree with the observed linear $`T(r)`$ in Figure 2. When $`\mathrm{}_{}(r)`$ and $`\rho (r)`$ are taken from Figure 1 as before, the only remaining undetermined parameter required to solve Equations (3) and (4) is $`\sigma _r(r_o)`$. The four short dashed lines in Figure 4 show the line of sight velocity dispersion profile $`\sigma (R)`$ for $`269.10\sigma _r(r_o)269.25`$ km s-1, corresponding to solutions 7 - 10 in Table 1. Since we are using the Kronawitter et al. $`\beta _K(r)`$, it is appropriate to compare these solutions only with the data of Bender et al. from which $`\beta _K(r)`$ was derived. Of these, the best-fitting solution is 9, but the overall flat or positive slope \[$`d\sigma (R)/dR>0`$\] at larger $`R`$ does not resemble the $`\sigma (R)`$ profiles of most well-observed E galaxies. Since none of these solutions is a particularly good fit to the Bender et al. data, it appears that $`\beta _K(r)`$ found from optical data alone is somewhat inconsistent with the observed hot gas temperature profile.

To explore alternate X-ray temperature profiles that may agree better with the $`\sigma (R)`$ of Bender et al., we consider solution 11 in which the hot gas is isothermal at the mean value of the two Buote temperatures in Figure 2, $`T_1=T_2=1.137\times 10^7`$ K. This solution, shown with a dotted line in Figure 4, decreases slightly with $`R`$, in better agreement with data in $`R\begin{array}{c}<\hfill \\ \hfill \end{array}2`$ kpc. Solutions with negative $`dT/dr`$ may fit even better, but all known X-ray observations for E galaxies have positive $`dT/dr`$ for $`r<R_e`$. Finally, we consider solution 12 based on a temperature profile having the same slope as the observations in Figure 2, but with uniformly higher temperatures, $`T_1=1.306\times 10^7`$ K and $`T_2=1.768\times 10^7`$ K. In solution 12 (solid line in Figure 4) the radially decreasing line of sight dispersion $`\sigma (R)`$ follows the Bender et al. data better than any other curve in Figure 4, possibly suggesting that the gas temperature should be higher than the observations in Figure 2. Therefore, either the gas temperatures from X-ray observations in Figure 2 are too low or $`\beta _K(r)`$ and the $`\sigma (R)`$ data from which it was derived are inaccurate.

## 3 Summary and Final Remarks

We have shown that Equations (3) and (4) establish a useful correspondence between the anisotropy of stellar orbits in elliptical galaxies $`\beta (r)`$ and the temperature of the hot interstellar gas $`T(r)`$. In practice, the $`\beta T`$ relation is based on the assumption of hydrostatic equilibrium for the hot gas and requires reliable observations of the stellar and gas densities and the stellar line of sight velocity dispersions. This connection between $`\beta `$ and $`T`$ is possible only if the galaxies are nearly spherical, as emphasized by Maggorian & Ballantyne (2002).

If the gas temperature profile $`T(r)`$ in an E galaxy is securely known from X-ray observations, $`\beta (r)`$ can be determined from the observed stellar velocity dispersion using Equations (3) and (4). In principle, this method of determining $`\beta (r)`$ could be done at galactic radii much larger than $`R_e`$ where stellar observations of $`h_4(R)`$ may be difficult but where observations of $`\sigma (R)`$ and X-ray gas temperatures are generally reliable. Conversely, if $`\beta (r)`$ can be observed with high precision from optical observations, the hot gas temperature $`T(r)`$ can be found from Equations (3) and (4) by using the gas density from X-ray observations. This means of determining $`T(r)`$ may be particularly valuable in regions of the hot gas atmospheres where X-ray spectra indicate a two-temperature or limited multi-temperature thermal structure (e.g. Buote 2002b) that may complicate the determination of the average local gas temperature.

As a representative example, we have used the $`\beta `$$`T`$ relation to determine $`\beta (r)`$ in the bright elliptical galaxy NGC 4472 from the hot gas temperature profile $`T(r)`$ found from X-ray observations. The resulting anisotropy, $`\beta 0.71\pm 0.15`$, indicates that the stellar orbits are considerably more radially biased than implied by the smaller $`\beta _K(r)\begin{array}{c}<\hfill \\ \hfill \end{array}0.3`$ determined from optical observations alone (Kronawitter et al. 2000). These large $`\beta `$ for NGC 4472 are based on the Fried-Illingworth stellar velocity dispersions $`\sigma (R)`$. Somewhat lower anisotropies, $`\beta =0.63\pm 0.15`$ are consistent with the velocity dispersion observations of Bender et al., but these $`\beta `$ still exceed the $`\beta _K(r)\begin{array}{c}<\hfill \\ \hfill \end{array}0.3`$ found optically by Kronawitter et al. (2000). Alternatively, if the Kronawitter et al. $`\beta =\beta _K(r)`$ and the Bender et al. (1994) data are correct, then the gas temperature in NGC 4472 may be higher than previously thought from X-ray observations and possibly have a different slope $`dT/dr`$. These conclusions apply only to the limited region $`r<R_e`$ for which $`\sigma (R)`$ data currently exist for NGC 4472.

It is clear, at least for NGC 4472, that currently available observations are not sufficiently accurate to establish either $`\beta (r)`$ or $`T(r)`$ with complete confidence. Accurate optical observations of $`\sigma (R)`$ and stellar line profiles can be observed with 8 - 10 meter telescopes to several $`R_e`$ and $`T(r)`$ is also often known at these large radii. In addition, it may be possible to determine gas temperature profiles with more accurate or prolonged X-ray observations than those in Figure 2. The $`\beta (r)`$ determined from X-ray emission can serve to calibrate the accuracy of purely optical determinations or to provide a test for hydrostatic equilibrium. To accomplish this, it would be desirable to obtain high quality stellar dispersion and gas temperature observations of those galaxies – such as NGC 4472, NGC 4649, NGC 1399, M87, etc. – that are luminous at both optical and X-ray frequencies.

Studies of the evolution of hot gas in elliptical galaxies at UC Santa Cruz are supported by NASA grants NAG 5-8409 NAG 5-13275 and NSF grants AST-9802994 and AST-0098351 for which we are very grateful.


text:
# Keys to Cosmology – Clusters of Galaxies

## 1 Introduction

Clusters of galaxies are very versatile tools for various types of analyses. They can be used to determine cosmological parameters as well as physical processes on large scales and extreme environments. Out of this large field we select here mass and dark matter measurements, determination of the dynamical state and the investigation of the interaction between the cluster galaxies and the intra-cluster gas.

## 2 Mass Determination and Dark Matter

As clusters of galaxies are large structures which can be regarded as being representative for the universe as a whole their ratio of baryonic matter to dark matter is a crucial number. Dark matter makes up about 80% of the total cluster mass, with the total cluster mass being determined by different methods.

One of the methods uses X-ray observations. The X-ray emitting gas traces the total cluster potential and hence traces the total cluster mass. With this mass determination one finds the following mass fractions: mass in galaxies 3-5%, mass in the intra-cluster gas 15-20%.

For the X-ray mass determination two assumptions are required: hydrostatic equilibrium and spherical symmetry. Tests with a sample of X-ray clusters together with analytical models showed that even if clusters are not spherically symmetric, but elongated the determination of total masses and gas mass fractions is quite reliable (Piffaretti et al. 2003). Only if projected clusters masses are calculated for the comparisons with masses from gravitational lensing differences are visible (see Fig. 1). Therefore part of the observed discrepancies between X-ray masses and lensing masses can come from these projection effects.

The second assumption, hydrostatic equilibrium, was tested by hydrodynamic simulations and found to be a valid assumption for relaxed i.e. non-merger clusters (Evrard et al. 1996; Schindler 1996). Clusters in the process of a major merger are obviously not in hydrostatic equilibrium and therefore the mass determination can easily be wrong by a factor of two.

Another possible source for the discrepancy between X-ray and lensing masses found in some clusters can be non-thermal pressure. While in the centre of clusters the negligence of non-thermal pressure from cosmic rays and magnetic fields can lead to a significant mass underestimate (Colafrancesco et al. 2003), we found that in the outer part, where the mass is usually determined, the non-thermal pressure from magnetic fields is negligible (Dolag & Schindler 2000). Only if the clusters are in the process of merging, their magnetic field can temporarily be so strongly enhanced that a mass underestimate down to 50% can be found.

A second way to determine the mass is gravitational lensing. Light from background galaxies is deflected by the huge mass of a cluster. Therefore the images of these background galaxies are distorted and show up as arcs (= strong lensing). The arcs are very thin and typically very faint structures. Under non-ideal observing conditions (e.g. bad seeing) they are easily dispersed and disappear into the background. Even under ideal conditions they are not easy to detect, because they are often just above the background level. To remove the noise and make faint structures better visible usually smoothing is applied. Unfortunately in the case of such thin structures as arcs the smoothing procedure often leads to a dispersion of the few photons so that the arcs are not detectable. To prevent this dispersing we have developed an algorithm that automatically smooths only along the arcs and not perpendicular to them, by so-called “anisotropic diffusion” (see Fig. 2, Lenzen et al. 2003). The subsequently applied source finding procedure extracts all the information from the sources necessary to distinguish arcs from other sources. This new algorithm is much more efficient in finding gravitational arcs than existing source detection algorithms because it is optimised just for this purpose.

The matter density of the universe $`\mathrm{\Omega }_m`$ can be determined with mass determinations of galaxy clusters in two ways. One possibility is to measure total mass and baryonic mass and assume that the ratio between the two is representative for the universe as a whole. This methods yields values around $`\mathrm{\Omega }_m0.3`$ (Ettori et al. 2003, Castillo-Morales & Schindler 2003). Another way is to determine the mass function which yields similar results (Rosati et al. 2002, Reiprich & Böhriger 2002). In the latter method one must be careful with selection effects and source confusion (see Gil-Merino & Schindler 2003).

## 3 Dynamical State

The dynamical state of clusters depends sensitively on the matter density and slightly also on the dark energy. Therefore the determination of the dynamical states of many clusters is another independent way to constrain cosmological parameters.

As from observations one can obtain only snapshots of the cluster evolution, it is ideal to use numerical simulations to learn how different dynamical states appear in X-ray images or temperature maps. Therefore we perform hydrodynamic simulations of the intra-cluster medium in a collaboration between University of Innsbruck (Domainko, Kapferer, Kimeswenger, Mair, Schindler, van Kampen), University of Edinburgh (Mangete, Ruffert) and Max-Planck-Institut für Gravitationsphysik (Benger). We use a grid code (Piecewise Parabolic Method, Colella & Woodward 1984) with multiply nested grids (Ruffert 1992).

Preliminary results of the simulations are shown in Figs. 3 and 4. Subcluster mergers lead to numerous shocks, which are moving outwards. The shocks appear in the temperature as well as in the density distribution as steep gradients. It is possible that particles, which are responsible for the radio haloes in several clusters, are accelerated to relativistic energies in these shocks (Giovannini & Feretti 2002).

Even before the collision of subclusters the gas is heated. The gas between the two subclusters is compresses and therefore it is heated. Hence a hot region between the subclusters is visible in the temperature map of the intra-cluster gas. For more information, images and films see http://astro.uibk.ac.at/astroneu/hydroskiteam/index.htm.

Such a feature of a heated region is e.g. observed in the cluster CL0939+4713 (De Filippis et al. 2003). The X-ray image shows (apart from several point-like sources) two extended subclusters. The temperature map clearly shows a hot region between the two subclusters as it is expected for the early stage of an (almost) central collision (see Fig. 5).

## 4 Interaction between Galaxies and the Intra-Cluster Gas

Until recently galaxies and the intra-cluster gas have been treated mostly independently. But in the last years a lot of examples of interaction between these two components have been found. The evidence does not only come observations of galaxies in several wavelengths, but also from observations of the intra-cluster gas, in particular of the metal abundances. As the metal content of the intra-cluster gas is on average roughly one third of the solar value, i.e. the total amount of metals in the intra-cluster gas is about the same as in all galaxies together, a lot of gas must have been transported from the galaxies into the intra-cluster medium. Several transport processes have been suggested: ram pressure stripping (Gunn & Gott 1972), galactic winds (De Young 1978), galaxy-galaxy interaction, jets from AGN and others. A closer look at the different metal enrichment processes is therefore very important for the understanding of cluster formation and galaxy evolution. So far numerical simulations including only some of the processes gave quite discordant results on the efficiency of the different processes (e.g. Cen & Ostriker 1999; Aguirre et al. 2001; Metzler & Evrard 1994, 1997). To improve on this we have started a comprehensive project. We use the hydrodynamic simulations mentioned above and include all possible enrichment processes. We calculate simulated metallicity maps for direct comparison with observations.

Finally, the first observed metallicity maps are available because the X-ray satellites XMM and CHANDRA have sufficiently high sensitivity and provide the possibility for spatially resolved spectroscopy. Therefore the evolution of metals in the intra-cluster gas as well as the spatial distribution can be measured now and be compared with simulations. Recently, in a number of clusters variations in the overall metallicity distribution have been found – not only in the cluster centre (Perseus cluster: Schmidt et al. 2002, A2199: Johnstone et al. 2002, A3558 and 3C129: Furuzawa et al. 2002, A1060: Yamasaki et al. 2002, A3671: Hudaverdi et al. 2002, 2A0335+096: Tanaka et al. 2002, AWM7: Furusho 2002, Cl0939+4713: De Filippis et al. 2003).

First results of the simulations are shown in Figs. 6 and 7. The simulations shown here include so far only ram-pressure stripping. While many other simulations have concentrated on the effect of the stripping on the galaxies (Abadi et al. 1999; Mori & Burkert 2000; Quilis et al. 2000; Vollmer et al. 2001; Schulz & Struck 2001; Toniazzo & Schindler 2001; Otmianowska-Mazur & Vollmer 2003), we concentrate here on the effects of the stripping process on the intra-cluster gas.

The gas that is stripped off the galaxies is visible not only in the metallicity maps of the intra-cluster gas, but also in the temperature maps, because this gas is cooler than the surrounding intra-cluster medium. It takes a while for the gas from the galaxies to mix with the intra-cluster gas. Sometimes one can see even shock waves passing over these cool regions and the cool gas is not heated up immediately. Also the metallicity distribution is not expected to be homogeneous shortly after the stripping, because the gas is not mixed immediately. This expectation is indeed supported be the metallicity variations observed in the metallicity maps mentioned above. In the simulations we see the highest stripping rate and hence the highest metallicities in regions, where subclusters fall onto the main cluster.

For more information on the simulation method, first results in form of images and films see http://astro.uibk.ac.at/astroneu/hydroskiteam/index.htm.

### Acknowledgements

This work is supported partly by the Austrian Science Foundation FWF, Project P15868, and by the EU Exchange Programme “Training and Research on Advanced Computing Systems” (TRACS) between Innsbruck and Edinburgh.


text:
# Global Monte Carlo for fermions using ordered statistics

## 1 INTRODUCTION

The standard method for stochastic representation of the fermion determinant is the pseudofermion method . It allows to trade the computation of the determinant in favour of extra noise in the path integral of lattice QCD. The other alternative is to use the stochastic Taylor expansion of . This suffers from the so called probability sign problem which is created by the noise . The problem is less apparent if the method is applied to a fractional power of the determinant, a fact which is used by the Kentucky Monte Carlo algorithm .

In this talk I show that the sign problem created by the noise can be eliminated altogether if one uses a certain statistic from a sample of unbiased noisy estimators of the effective fermion action. As expected, a statistics which is an unbiased estimation of the effective fermion action cannot yield an unbiased estimation of the fermion determinant. The opposite logic suggests that a biased estimation of the former could lead in principle to an unbiased estimation of the latter if this is done ”correctly”. This is indeed the case if one chooses a certain order statistic from the sample. I show below how this can be realized using the central limit theorem and the asymptotic distribution of a central order statistic.

## 2 USING THE CENTRAL LIMIT THEOREM

Let $`X_1,X_2,\mathrm{},X_n`$ be unbiased noisy estimators of the fermion effective action $`S_{\text{eff}}=\text{Tr log}(D^{}D)`$, where $`D`$ is a generic Dirac. A particular way (and popular) to do it is using real (there is nothing special, one can use complex numbers as well) random vectors $`Z_1,Z_2,\mathrm{},Z_n`$ such that:

$$X_i=Z_i^{}[\text{log}(D^{}D)]Z_i$$
(1)
Computation of the right hand side can be done for example using Lanczos methods as described in detail in . It is easy to show that $`X_i`$ is an unbiased estimator of $`S_{\text{eff}}`$. Another useful result (that so far went unappreciated and wasn’t used) is the normality of the $`X_i`$ distribution. This can be easily shown for $`Z_2`$ noise vectors using the central limit theorem. Note that the whole idea relies on this result.

## 3 BASICS OF ORDER STATISTICS

Everyone has encountered at least once the notion of median as an unbiased estimator of the population mean. It is identified as the middle point after the sample of real numbers $`X_1,X_2,\mathrm{},X_n`$ drawn from the same continuous (to fix the idea) distribution $`f(x),xR`$ has been ordered in the ascending order:

$$X_{(1)},X_{(2)},\mathrm{},X_{(n)}$$
(2)
This is possible if the sample size $`n`$ is odd. In case $`n`$ is even the median is defined to be any number in the interval $`[X_{(n/2)},X_{(n/2+1)}]`$.

The median is a statistic as it is the arithmetic average. It is nothing particular to the median as an instance of a statistic calculated from the ordered data. In general, $`X_{(k)},k=1,\mathrm{},n`$ is called the $`k`$’th ordered statistic.

Distribution of $`X_{(k)}`$ can be found using a simple heuristic. First, let $`f_k(x)dx`$ be the probability that $`X_{(k)}`$ is in an infinitesimal interval $`dx`$ about $`x`$. This means that one of the sample variables is in this interval, $`k1`$ are less than $`x`$ and $`nk`$ are greater than $`x`$. The number of ways of choosing these variables is the multinomial coefficient $`C(n;k1,1,nk)`$. If $`F(x)`$ is the cumulative distribution function of the sample variables then one gets:

$$f_k(x)dx=[F(x)]^{k1}[1F(x)]^{nk}f(x)dx$$
(3)
This result can be formally proven. As a starting classic is the book of H.A. David . A recent introduction to the subject is provided in .

## 4 ASYMPTOTIC THEORY

The distribution of an order statistic can be of little use for direct analytical calculations of moments for most of known parent distributions. Of course one can use numerical methods to compute numerical values to the desired accuracy. However, as the sample size grows the order statistic distribution approaches some limiting distributions. The situation is analogous to the central limit theorem. In general the distribution is not necessarily asymptotically normal as one could have expected although this will be the case for the central order statistics that will be dealt with below.

For $`0<\alpha <1`$, let $`i=[n\alpha ]+1`$, where $`[n\alpha ]`$ represents the integer part of $`n\alpha `$. Let also $`F`$ be absolutely continuous with $`f`$ positive at $`F^1(\alpha )`$ and continuous at that point. Then the following result applies:

Theorem (asymptotic distribution of a central order statistic). As $`n\mathrm{}`$,

$$\sqrt{n}f[F^1(\alpha )]\frac{X_{(i)}F^1(\alpha )}{\sqrt{\alpha (1\alpha )}}N(0,1)$$
(4)
The proof can be found eg, in .

Example. Suppose that the distribution is symmetric around the population mean $`\mu `$. Let also assume that the variance $`\sigma ^2`$ is finite and $`f(\mu )`$ is finite and positive. For simplicity $`n`$ is taken to be odd. Then the median $`X_{(n+1)/2}`$ is an unbiased estimator of $`\mu `$ and asymptotically normal. Further, Var$`[X_{(n+1)/2}][4nf(\mu )^2]^1`$.

## 5 UNBIASED ESTIMATION OF THE FERMION DETERMINANT

Let the $`X_i,i=1,\mathrm{},n`$ be a normally distributed sample of unbiased estimators of the fermion effective action. I am interested in the expected value of $`e^{X_{(i)}}`$ with $`i=[n\alpha ]+1`$ and $`\alpha `$ as above. Let $`XN(\mu ,\sigma )`$. It is easy to show that:

$$<e^X>=e^{\mu +\sigma ^2/2}$$
(5)
Hence, an unbiased estimation of the effective action gives a biased estimation of the fermion determinant. Now let $`x_\alpha :=(F^1(\alpha )\mu )/\sigma `$. Then $`X_{(i)}`$ is asymptotically normal with mean $`\mu +x_\alpha \sigma `$ and variance $`2\pi \sigma ^2\alpha (1\alpha )e^{x_\alpha ^2}/n`$. For the ordered statistic one gets:

$$<e^{X_{(i)}}>=e^{\mu x_\alpha \sigma +\pi \sigma ^2\alpha (1\alpha )e^{x_\alpha ^2}/n}$$
(6)
To have an unbiased estimator of the fermion determinant in terms of an ordered statistic one should have $`<e^{X_{(i)}}>=e^\mu `$. This can be achieved if:

$$x_\alpha =\pi \sigma \alpha (1\alpha )e^{x_\alpha ^2}/n$$
(7)
The crucial question here is: given $`n`$ and $`\sigma `$ is there any solution of the above equality in terms of $`x_\alpha `$? To answer this question I express first $`\alpha `$ in terms of $`x_\alpha `$. Since $`(X\mu )/\sigma N(0,1)`$ then I have:

$$\alpha =P(X<x_\alpha )=\frac{1}{2}+\frac{1}{2}\text{erf}(\frac{x_\alpha }{\sqrt{2}})$$
(8)
Hence, the equation to be solved is:

$$x_\alpha =\frac{\pi \sigma }{4n}[1\text{erf}(\frac{x_\alpha }{\sqrt{2}})^2]e^{x_\alpha ^2}$$
(9)
Using Taylor expansion around $`x_\alpha =0`$ one gets:

$$x_\alpha =\frac{\pi \sigma }{4n}[1+O(x_\alpha ^2)]$$
(10)
For small $`x_\alpha `$ the $`O(x_\alpha ^2)`$ can be neglected one gets a unique solution:

$$x_\alpha =\frac{\pi \sigma }{4n},x_\alpha <<1$$
(11)
In this case the sample size should be much greater then $`\sigma `$ and $`i`$ can be calculated to give:

$$i=[n\alpha ]+1=[\frac{n}{2}+\sqrt{\frac{\pi }{2}}\frac{\sigma }{4}]+1$$
(12)
Thus, using order statistics it is possible to have an unbiased estimator of the fermion determinant from a biased high estimator of the fermion effective action.

However, the estimator is only asymptotically unbiased meaning that the sample size has to be large enough. In practice, at least from the examples below the asymptotic regime is reached for $`n=40`$. In many application to lattice fermions $`\sigma `$ is large anyway. Thus, requiring $`n>>\sigma `$ will satisfy the conditions of the asymptotic normality as well.

Note also that $`\sigma `$ is unknown and has to be estimated from the sample. This is a source of systematic error that may bias the stochastic estimation of the determinant. To compute the error one has to substitute $`\sigma `$ with its estimator $`S`$ in the formula for $`x_\alpha `$. Then the average error is given by:

$$<e^{x_\alpha \sigma +\frac{\pi \sigma ^2}{4n}}>=<e^{\frac{\pi \sigma }{4n}(\sigma S)}>$$
(13)
The error has to be computed in order to control the systematic errors. This can be done using for $`\sigma `$ large samples to guarantee small errors and then substitute in the above formula.

## 6 APPLICATION TO GLOBAL MONTE CARLO ALGORITHMS

The application of choice of the stochastic determinant is in Monte Carlo algorithms. I have tested here in case of lattice Schwinger model. To simulate the model I made new gauge field proposals using the gauge action which consists of 10000 sweeps of local Metropolis steps followed by a global accept/reject step to include two degenerate flavors of fermions. Figure 1 shows the square Wilson loops as a function of the linear size of the loop using three algorithms: 1) computing exactly the determinant ratios, 2) using the stochastic method with $`n=30`$ noisy estimators and 3) the same as in 2) but with $`n=40`$ noisy estimators. This example shows that already $`n=40`$ estimators are enough to reach the exact result.