text:
# A Three Dimensional Magnetohydrodynamic Pulse in a Transversely Inhomogeneous Medium

## 1 Introduction

Problems of heating of the open corona of the Sun and acceleration of the solar wind are closely related with the interaction of MHD waves with plasma inhomogeneities and, in particular, with coupling of compressible and incompressible components of the waves. Investigation of coupling mechanisms between the waves is of prime importance because, on one hand, the Alfvén waves are usual candidates of energy transport from the lower layers of the solar atmosphere to the corona, and, on the other, compressible perturbations are subject to much more efficient dissipation than the incompressible Alfvén waves (see, e.g., Narain & Ulmschneider (1996)). Also, the compressive waves, as opposed to the Alfvén waves, can transport energy across the magnetic field, due to the fact that their propagation in space is not constrained by the field lines.

The original idea of the phase mixing of incompressible Alfvén waves (Heyvaerts & Priest (1983)) was based on the following argument: when plasma has a density gradient perpendicular to the magnetic field, local Alfvén speed is a function of the transverse coordinate. Thus, when an Alfvén wave propagates along the field its perturbations on the adjacent field lines become out of phase. This stretching of Alfvén wave front creates progressively smaller spatial scales across the field. In turn, because the dissipation is proportional to the wave number squared, phase mixing leads to enhanced dissipation of the Alfvén wave directly. In the compressible plasma, as demonstrated by Malara, Petkaki & Veltri (1996); Nakariakov et al. (1997); Nakariakov et al. (1998); Botha et al. (2000); Tsiklauri, Arber & Nakariakov (2001); Tsiklauri, Nakariakov & Arber (2002), phase mixing of linearly polarized plane Alfvén waves leads to the enhanced nonlinear generation of fast magnetoacoustic waves. However, it was established, in 2.5D geometry, for the harmonic Alfvén wave Botha et al. (2000), and for a wide spectrum Alfvén pulse by Tsiklauri, Arber & Nakariakov (2001) that compressive perturbations, which are initially absent from the system, do not grow to a substantial fraction of the initial Alfvén wave due the destructive wave interference effect.

In this work we consider fully three dimensional geometry and study interaction of linear (with the non-linear effects totally ignored) MHD waves with a one-dimensional inhomogeneity of the plasma, taking into account compressibility of the plasma and the localization of the MHD pulse in the direction perpendicular to both the magnetic field and the inhomogeneity gradient. In particular, we shall study how the phenomenon of phase mixing is affected by these factors. A similar problem was studied by De Groof & Goossens (2000) in the context of resonant absorption of MHD waves in coronal loops as a heating mechanism. However, our treatment is more general as we perform a direct numerical 3D simulation without resorting to Fourier transform in time and $`y`$ direction (and, consequently, our study can easily be generalized to the case of 2D and 3D structuring) allowing us to consider an initial value problem. Besides, De Groof & Goossens (2000) used a harmonic form of the initial perturbation as opposed to our spatially localized, Gaussian, one. Recently, Hood, Brooks & Wright (2002) showed that the wide spectrum regime of phase mixing can be quite different from the harmonic one, as in fact, phase mixing of localized Alfvén pulses results in a slower, algebraic, damping as opposed to the standard exponential damping of harmonic Alfvén waves, which suggest that localized Alfvénic perturbations will transport energy higher into the corona than harmonic ones. Such Alfvénic perturbations could be generated e.g. by transient events such as solar flares, coronal mass ejections, etc. While most of the phase-mixing studies concentrated on harmonic perturbations (e.g. Heyvaerts & Priest (1983); Hood, Ireland & Priest (1997); Hood et al. (1997); Nakariakov et al. (1998); Ruderman, Nakariakov & Roberts (1998); DeMoortel et al. (1999); DeMoortel, Hood & Arber (2000); Grappin, Léorat & Buttighoffer (2000); Botha et al. (2000)), only few considered spatially localized ones (Nakariakov et al. (1997); Tsiklauri, Arber & Nakariakov (2001); Tsiklauri, Nakariakov & Arber (2002); Hood, Brooks & Wright (2002)).

An additional motivation to this study is connected with the growing interest to the problem of interaction of MHD waves with 2D and 3D plasma structures and irregularities, as the coronal and wind plasmas are observed to be structured in all three dimensions. 2D and 3D structuring can dramatically affect properties of MHD waves. In particular, as it was shown in the incompressible regime by Similon & Sudan (1989) and confirmed in numerical experiments performed by Petkaki, Malara & Veltri (1998) and Malara, Petkaki & Veltri (2000), the structuring dramatically increases the efficiency of wave dissipation. However, in those studies, the compressible effects were not taken into account. From this point of view, our investigation of the interaction of a 3D MHD pulse with a 1D plasma inhomogeneity in the compressible regime provides a necessary building element of the general theory of MHD wave interaction with plasma inhomogeneities.

The present model is based upon the MHD description of plasma and we focus on the problems that are of relevance to the heating of solar corona, acceleration and the dynamics of solar wind. While, similar studies exist that deal with, for instance, terrestrial auroral applications. Namely, an alternative to MHD, more relevant, particle-in-cell simulations (using a guiding center implicit code for the electrons) have been performed (Génot, Louarn & Le Quéau (1999); Génot, Mottez & Louarn (2000); Génot, Louarn & Mottez (2001)), which study the propagation and collisionless dissipation (via effective electron beam generation) of the Alfvén waves in the auroral density inhomogeneities (cavities).

In this work we study the propagation part of the problem, i.e. we consider an ideal plasma limit. The work is in progress to include finite plasma resistivity in order to investigate quantitatively the dependence of the decay of Alfvénic part of the MHD pulse upon the coupling to the existing compressive waves.

The paper is organized as follows: in section 2 we formulate our model. In section 3 we present the results of numerical simulation, while we close in section 4 with the discussion of main results.

## 2 The model

In our model we use equations ideal MHD

$$\rho \frac{V}{t}+\rho (V)V=p\frac{1}{4\pi }B\times \mathrm{curl}B,$$
(1)
$$\frac{B}{t}=\mathrm{curl}(V\times B),$$
(2)
$$\frac{p}{t}+Vp+\gamma pV=0,$$
(3)
where $`B`$ is the magnetic field, $`V`$ is plasma velocity, $`\rho `$ is plasma mass density, and $`p`$ is plasma thermal pressure. In what follows we use $`5/3`$ for the value of $`\gamma `$.

We solve equations (1)-(3) in Cartesian coordinates ($`x,y,z`$). Note that as we solve a fully 3D problem we retain variation in the $`y`$-direction, i.e. ($`/y0`$). The equilibrium state is taken to be an inhomogeneous plasma of density $`\rho _0(x)`$ and a uniform magnetic field $`B_0`$ in the $`z`$-direction. We consider a plasma configuration similar to the one investigated in Malara, Petkaki & Veltri (1996); Nakariakov et al. (1997); Nakariakov et al. (1998); Botha et al. (2000); Tsiklauri, Arber & Nakariakov (2001), i.e. the plasma has a one-dimensional inhomogeneity in the equilibrium density $`\rho _0(x)`$ and temperature $`T_0(x)`$. The unperturbed thermal pressure, $`p_0`$, is taken to be constant everywhere.

Next, we do usual linearization of the Eqs. (1)-(3) and write them in component form as following

$$\rho _0(x)\frac{V_x}{t}+\frac{p}{x}\frac{B_0}{4\pi }\left(\frac{B_x}{z}\frac{B_z}{x}\right)=0,$$
(4)
$$\rho _0(x)\frac{V_y}{t}+\frac{p}{y}+\frac{B_0}{4\pi }\left(\frac{B_z}{y}\frac{B_y}{z}\right)=0,$$
(5)
$$\rho _0(x)\frac{V_z}{t}+\frac{p}{z}=0,$$
(6)
$$\frac{B_x}{t}B_0\frac{V_x}{z}=0,$$
(7)
$$\frac{B_y}{t}B_0\frac{V_y}{z}=0,$$
(8)
$$\frac{B_z}{t}+B_0\left(\frac{V_x}{x}+\frac{V_y}{y}\right)=0,$$
(9)
$$\frac{p}{t}+\gamma p_0\left(\frac{V_x}{x}+\frac{V_y}{y}+\frac{V_z}{z}\right)=0.$$
(10)
It is useful to re-write Eqs. (4)-(10) in a form of three coupled wave equations as following

$$\left[_{tt}^2\left(c_s^2(x)+c_A^2(x)\right)_{xx}^2c_A^2(x)_{zz}^2\right]V_x$$
(11)
$$\left[\left(c_s^2(x)+c_A^2(x)\right)_{xy}^2\right]V_y\left[c_s^2(x)_{xz}^2\right]V_z=0,$$

$$\left[_{tt}^2\left(c_s^2(x)+c_A^2(x)\right)_{yy}^2c_A^2(x)_{zz}^2\right]V_y$$
(12)
$$\left[\left(c_s^2(x)+c_A^2(x)\right)_{xy}^2\right]V_x\left[c_s^2(x)_{yz}^2\right]V_z=0,$$

$$\left[_{tt}^2c_s^2(x)_{zz}^2\right]V_z\left[c_s^2(x)_{xz}^2\right]V_x$$
(13)
$$\left[c_s^2(x)_{yz}^2\right]V_y=0,$$

where $`c_A(x)=B_0/\sqrt{4\pi \rho _0(x)}`$ and $`c_s(x)=\sqrt{\gamma p_0/\rho _0(x)}`$ denote local Alfvén and sound speeds respectively.

We solve Eqs.(4)-(10) numerically after re-writing them in a dimensionless form using following normalization: $`B_{x,y,z}=B_0\overline{B}_{x,y,z}`$, $`(x,y,z)=a_{}(\overline{x},\overline{y},\overline{z})`$, $`c_A(x)=B_0/\sqrt{4\pi \rho _0(x)}=B_0/\sqrt{4\pi \rho _{}}/\sqrt{32\mathrm{tanh}(\lambda x)}=c_A^{}/\sqrt{32\mathrm{tanh}(\lambda x)}`$, $`t=(a_{}/c_A^{})\overline{t}`$, $`V_{x,y,z}=c_A^{}\overline{V}_{x,y,z}`$. Note, that $`c_s(x)=\sqrt{\gamma \beta /2}c_A(x)`$, where $`\beta `$ stands for the ratio of thermal and magnetic pressures $`\beta =p_0/(B_0^2/8\pi )`$. Here, $`\lambda `$ is a free parameter which controls the steepness of the density profile gradient. In our simulations we use $`\lambda =0.5`$. In what follows we omit bars on top of the physical quantities.

## 3 Numerical Results

In order to solve Eqs.(4)-(10) numerically we have written a new numerical code dt4dx10, which uses a high order finite difference scheme. Namely, it evaluates 10-th order centered spatial derivatives, and advances solution in time using 4-th order Runge-Kutta algorithm. Therefore, dt4dx10 is $`O[(\mathrm{\Delta }t/T)^4]`$ accurate in time and $`O[(max(\mathrm{\Delta }x,\mathrm{\Delta }y,\mathrm{\Delta }z)/L)^{10}]`$ accurate in space, with $`T`$ and $`L`$ denoting run time and linear size of the simulation domain. The use of such a high-order numerical approach was motivated by the very nature of the problem considered, as development of phase mixing leads to the generation of very steep profiles in the wave front.

The simulation cube size is set by the limits $`25.0x25.0`$, $`25.0y25.0`$ and $`25.0z25.0`$. Boundary conditions used in all our simulations are zero-gradient in all three spatial dimensions.

We have performed calculation on various resolutions in attempt to achieve convergence of the results. The graphical results presented here are for the spatial resolution $`128^3`$, which refers to number of grid points in $`x,y`$ and $`z`$ directions respectively. We have also performed calculation on the spatial resolution $`256^3`$ and we found that the results converge perfectly. This is understandable due to the high order of the scheme even for the $`128^3`$ resolution time-error, $`O[(\mathrm{\Delta }t/T)^4]7\times 10^8`$, and spatial-error, $`O[(max(\mathrm{\Delta }x,\mathrm{\Delta }y,\mathrm{\Delta }z)/L)^{10}]8\times 10^{22}`$. Since a $`256^3`$ resolution run takes about 8 hours on Compaq ES40 with eight EV6 500-MHz processors, while a $`128^3`$ resolution takes only 1 h on Compaq ES40 with four EV6 500-MHz processors producing the same results, we opted for the latter, less CPU-consuming, alternative.

In the numerical simulations the MHD perturbation is initially a plane (with respect to $`x`$-coordinate) pulse, which has a Gaussian structure in $`y`$ and $`z`$-coordinates

$$V_y(x,y,z,t=0)=c_A(x)\mathrm{exp}\left((\alpha _yy)^2(\alpha _zz)^2\right).$$
(14)
Here, $`\alpha _y`$ and $`\alpha _z`$ are free parameters which control the strength of gradients in $`y`$ and $`z`$ direction of the initial perturbation. As the problem considered is linear, the wave amplitude can be taken to be normalized to unity.

In the geometry considered, when the initial perturbation depends on $`y`$ it cannot be regarded as pure Alfvénic one. In fact, in such a pulse, all three waves – Alfvén, fast and slow magnetosonic waves – are inter-coupled so that there is no use of their separation per se. The pulse is set to be initially plane in the $`x`$-direction. This allows us to emphasize the effect of the inhomogeneity on the pulse evolution. From the point of view of applications, this simply means that the initial characteristic size of the pulse in that direction is greater than the scale of the inhomogeneity.

The particular choice of the initial condition Eq. (14) is motivated by the argument that when we set $`B_y`$ initially to zero we automatically guarantee fulfillment of $`\mathrm{div}B=0`$ when $`\alpha _y0`$, and as the system evolves it adjusts itself which mode to excite (depending whether $`\alpha _y`$ is zero or not). For instance, if we choose our initial conditions as $`V_y`$ given by Eq. (14), $`B_y(x,y,z,t=0)=\mathrm{exp}\left((\alpha _yy)^2(\alpha _zz)^2\right)`$, with $`\alpha _y=0`$ and the rest of physical quantities set to zero, we would excite a pure Alfvénic pulse traveling in the negative direction along $`z`$-axis. However, as long as $`\alpha _y0`$ we have to set $`B_y(x,y,z,t=0)=0`$ in order to fulfill $`\mathrm{div}B=0`$. Besides, our choice of the localized “kinematic” (the velocity is initially perturbed only, while the magnetic field is constant) perturbation Eq. (14) is well motivated by the fact that such perturbations can arise, e.g., during coronal mass ejections, solar flares, or other violent events (see, e.g., Roussev et al. 2001), etc.

In Figure 1 we present time evolution of the transverse component of the plasma velocity, $`V_y`$ and the density perturbation $`\rho ^{}/\rho _0(x)`$. The choice of these physical quantities is motivated by the following reasoning: the transverse velocity component $`V_y`$ is subject to Alfvén wave phase mixing and the density perturbation demonstrates the effect of the plasma compressibility. The left column of Fig. 1 shows three snapshots of $`V_y`$ and two snapshots ($`\rho ^{}/\rho _0(x)`$) through the $`y=0`$ cross-section, the middle column shows three snapshots of $`V_y`$ and two snapshots ($`\rho ^{}/\rho _0(x)`$) through the $`x=0`$ cross-section, and the right column shows three snapshots of $`V_y`$ and two snapshots ($`\rho ^{}/\rho _0(x)`$) through the $`z=0`$ cross-section. The snapshots of $`V_y`$ and $`\rho ^{}/\rho _0(x)`$ at times $`t=7.5`$ and 15 are placed such that it is easier to spot the differences. The snapshot of $`\rho ^{}/\rho _0(x)`$ at $`t=0`$ is not present in the figure as it is identically zero and it is generated later by the $`y`$-gradients in the function $`V_y`$. This compressible perturbation is a significant component of the pulse. What we gather from this graph is as following:

* $`y=0`$ cross-section: the initial MHD pulse, which is plane with respect to $`x`$-coordinate, is split in two D’Alambert’s solutions, with half of the amplitudes, traveling in two opposite directions along the magnetic field. Because the unperturbed density is inhomogeneous across the magnetic field (in the $`x`$-direction) the local Alfvén speed depends upon $`x`$. Thus, the Alfvénic pulse is phase-mixed and its initially plane front is continuously distorted creating transverse gradients. Note, that the right wing of the pulse is asymmetric to the left one, this is due to the fact that its amplitude drops from 1.0 in the $`x>0`$ domain to $`1/\sqrt{5}`$ in the $`x<0`$. The generated density perturbation evolves in the similar manner, however there is a notable difference – as plasma $`\beta =0.5`$ in this case, we observe that the density perturbation travels at a slower speed (slow magnetosonic wave speed) than the Alfvénic one. If we look at the right wing of $`V_y(x,0,z)`$ at $`t=15`$, we see that it traveled to z=15, which is consistent with the Alfvén speed equals 1 (recall that for $`x>0`$ $`c_A(x)=1`$). The density perturbation as it travels along the field should have a velocity of a slow magnetosonic wave, which for the field aligned propagation coincides with the local speed of sound. For the parameters considered $`c_s(x)=\sqrt{\gamma \beta /2}c_A(x)=0.65c_A(x)`$. Therefore, the right wing of the density perturbation at $`t=15`$ should have traveled to a position $`z9.7`$. That is what we actually see from the bottom panel of the left column in Fig.1.
* $`x=0`$ cross-section: the initial Alfvénic pulse which has an asymmetry as $`\alpha _y\alpha _z`$ Gaussian bell shape is split in two smaller amplitude bells traveling in two opposite directions along the magnetic field plus some semi-elliptic shaped wake expanding outwards in all directions. Note, that because the unperturbed density is homogeneous across the magnetic field in $`y`$-direction we do not see any distortions on the wave front with respect to the $`y=0`$ line as we saw in the previous case. The generated density perturbation evolves in the similar manner, however a notable difference (can be better seen from Fig. 2) is that the solution consists of positive and negative spikes plus semi-elliptic shaped wakes that are symmetric with respect to $`\rho ^{}/\rho _0(x)=0`$ plane. This can be understood readily as due to the conservation of total mass, integral of density perturbation over the volume should, indeed, be zero.
* $`z=0`$ cross-section: the MHD pulse evolution, viewed through this cross-section, is similar to the case of $`y=0`$ cross-section (see above). Again we see two phase-mixed D’Alambert’s solutions propagating in the opposite directions. The notable difference is that through $`z=0`$ cross-section the pulse looks wider as $`\alpha _y=0.6`$ while $`\alpha _z=1.0`$. The generated density perturbation evolves in the similar manner, as the Alfvénic one. However, in differ from $`y=0`$ cross-section, we see no difference in traveling speeds. This can be explained as following: the pulse was not excited as a “pure Alfvén wave”. In fact, should the perturbation $`V_y`$ be purely Alfvénic there would have been no propagation along $`y`$ axis at all (recall that the magnetic field is directed along the $`z`$-axis). Thus, what we see through the $`z=0`$ cross-section both for $`V_y`$ and $`\rho ^{}/\rho _0(x)`$ is the compressive part of the perturbation traveling at the fast magnetosonic wave speed. For the parameters considered, the fast magnetosonic wave speed is $`\sqrt{1+\gamma \beta /2}c_A(x)=1.2c_A(x)`$. Therefore, the right wing of both $`V_y`$ and density perturbation at $`t=15`$ should have traveled to a position $`y17.9`$. That is what we actually observe from the two bottom panels of the right column in Fig.1.

Also, the figure shows the obvious anisotropy in the pulse propagation, connected with the presence of two chosen directions: the direction of the straight magnetic field and the direction of the density gradient.

Perturbation of $`B_y`$ is initially absent from the system, but as it is a potential component of the kinetic counterpart, $`V_y`$, of the wave, it is soon generated and further evolves similarly to $`V_y`$ in a form of two negative and positive phase-mixed pulses traveling into two opposite directions. Also, the development of the pulse is accompanied by generation of compressible components. The structure of the density perturbations in the pulse at the time $`t=15`$ is shown in Figures. 2, 3 and 4. Note, that in all the figures, pre-factor $`ϵ`$ stands to underscore the fact that the linear problem, which we study, does not depend on the initial amplitude. Thus, we use $`A=1`$, while we have to bear in mind that linear approximation itself is valid for small amplitudes that is why we use small, arbitrary pre-factor $`ϵ`$ in our notations.

A fairly good quantity describing the relation between different components of the pulse is the maximum of absolute value over the whole simulation domain. We plot this quantity for all physical variables in Figure 5. In this figure, thick solid lines correspond to the resolution $`256^3`$, while thin solid lines to that of $`128^3`$, both for $`\alpha _y=0.6`$. It is remarkable that on all panels these two curves practically do overlap, which serves as a proof of convergence of our simulation results. The last figure in the bottom row presents the maximum of absolute value over the whole simulation domain of $`\mathrm{div}B`$, and we indeed observe almost perfect fulfillment of the fundamental law, $`\mathrm{div}B=0`$, which comes as a bonus of high-order, centered, finite difference numerical scheme. As expected, Fig. 5 shows that when $`\alpha _y=0`$, the pulse is perfectly Alfvénic (perturbing $`V_y`$ and $`B_y`$ only). This conclusion can be deducted either from analyzing Eqs. (11)-(13) or resorting to a classic mechanical analogy of coupled pendulums. In effect, Eqs. (11)-(13) also describe three inter-coupled mathematical pendulums which are located in $`xOz`$ plane. Thus, as long as we do not perturb these pendulums such that $`\alpha _y0`$, they will always oscillate in the $`xOz`$ plane. Thus, what we see in Fig. 5 when $`\alpha _y=0`$ (dash-dotted lines) is that initial kinematic Alfvén perturbation ($`V_y`$) is split in half-amplitude D’Alambert’s solutions (both for $`V_y`$ and $`B_y`$) and no other physical quantity is generated. However, when we switch on the coupling, $`\alpha _y0`$, compressible perturbations are generated.

In order to investigate the effect of phase mixing on the created quasi-steady MHD state, we produced time series of all physical quantities in several points of the simulation 3D box, Figures 7 and 8. Namely, at the points (x=0,y=-8,z=0, 4, 8), which are located in the phase mixing (spatial inhomogeneity) domain, and at (x=16,y=-8,z=0, 4, 8), which are far away from it. The three panels (from the left to the right), in effect, trace the dynamics in the $`z`$-direction, i.e. along the regular magnetic field lines. These types of data would be obtained by, for example, three satellites which are located in three different regions of solar wind.

There are two noteworthy features that can be gathered from these plots: First, we see that in the phase mixing region, $`V_x`$, which is associated with the fast magnetosonic component of the pulse, attains about 40 times larger values than in the region that is far from the phase mixing region. That is a sensible result, since it is known that phase mixing efficiently generates oblique fast magnetosonic waves. Second, the total magnetic field perturbation $`|BB_0|`$ is positively  correlated with the pressure (as well as density, which is proportional to the gas pressure) perturbation, which indeed is what was expected from the fast magnetosonic wave.

Also, Figs 5-8 demonstrate that, interacting with the plasma inhomogeneity, the longitudinally (parallel to the magnetic field) propagating part of the initial perturbation develops to an almost non-evolving (without change in amplitude) propagating state. In this quasi-steady MHD state all components of the pulse tend to propagate without change in amplitude, even though the effect of phase mixing is in action all the time ($`x`$-component of the inverse characteristic spatial scale tends to infinity).

In Figure 6 we investigate the energetics of our numerical simulation. In fact, we observe nearly perfect ($`\pm 0.2`$% error) conservation of total energy by dt4dx10 numerical code. The major conclusion which can be drawn from this graph is that when $`\alpha _y=0`$ there is no internal (compressive) energy generation, while with the increase of $`\alpha _y`$ its final, asymptotic levels increase progressively.

The parametric space of the problem is studied in Figure 9 by plotting the final levels of relative density perturbation, approached by the MHD pulse, as a function of the initial localization of the pulse in the $`y`$-direction $`\alpha _y`$, and plasma $`\beta `$. We gather from Fig. 9 that the achieved levels of density fluctuations depend weakly on $`\alpha _y`$, while there is somewhat stronger dependence on $`\beta `$.

## 4 Conclusions

In this paper we present results of numerical modeling of interaction of an initially localized MHD pulse with a transverse inhomogeneity of the plasma. The pulse is non-plane in all three spatial directions and is subject to the effect of phase mixing. An impulsively generated perturbation of the plasma velocity develops to an anisotropically propagating pulse, which has a significant compressible component. More specifically:

* The non-uniformness of the pulse in the homogeneous transverse (or, in other words, in the third, perpendicular to both the magnetic field and the inhomogeneity gradient) direction (non-zero $`\alpha _y`$) leads to appearance of compressive perturbations in the pulse.
* In the presence of plasma inhomogeneity, a non-uniform MHD pulse is essentially compressible. More specifically, when the pulse is initially plane in the homogeneous transverse direction, $`\alpha _y=0`$, there is no internal (compressive) energy generation, while with the increase of $`\alpha _y`$, the levels of the compressible energy increase progressively.
* A propagating MHD pulse asymptotically reaches a quasi-steady state. The final levels of density perturbation, which can be considered as a measure of the compressibility in the pulse, depend weakly on $`\alpha _y`$, while there is somewhat stronger dependence on $`\beta `$.
* A smooth 1D transverse inhomogeneity of the plasma supports propagation of compressible MHD pulses. This mechanism of wave guiding is different from the well-understood phenomenon of refraction of fast magnetosonic waves, as the profile of the Alfvén speed in the inhomogeneity does not have a minimum.
* Weakly non-plane MHD waves are subject to phase mixing and can qualitatively be considered in terms of the 2.5D theory of Alfvén wave phase mixing. Consequently, the 2.5D theory of Alfvén wave phase mixing remains relevant to the 3D case too. However, quantitative theories of the interaction of MHD waves with plasma inhomogeneities should include the compressibility of the plasma as a necessary ingredient.

Our main conclusion is that the effects of three-dimensionality, compressibility and inhomogeneity should be all together taken into account in the wave-based theories of coronal heating and solar wind acceleration, and as well as in the theories of MHD turbulence.

###### Acknowledgements.

The authors are grateful to Tony Arber for valuable advise when dt4dx10 code was written and for a number of valuable discussions. DT acknowledges financial support from PPARC. Numerical calculations of this work were performed using the PPARC funded Compaq MHD Cluster at St Andrews and Astro-Sun cluster at Warwick.


format: text

text:
# The Unique History of the Globular Cluster 𝜔 Centauri

## 1 Omega Centauri, the most massive globular cluster in the Galaxy

Omega Centauri is unique among globular clusters in the Milky Way in that its member stars exhibit a wide spread in metallicity (Freeman & Rodgers, 1975; Freeman & Norris, 1981; Mallia & Pagel, 1981). This has led to the general belief that $`\omega `$ Cen was chemically self-enriched, while the other clusters were not (Suntzeff & Kraft 1996; see Meylan & Heggie 1997 for a review). In this context, it is often noted that $`\omega `$ Cen is the most massive globular cluster in the Milky Way and thus possibly the most capable of retaining its own stellar ejecta (e.g., Ikuta & Arimoto 2000). The globular cluster G1 in the halo of the Andromeda galaxy has similar properties; it appears to be self-enriched and is also one of the most massive clusters in its host galaxy (Meylan et al., 2001).

Omega Centauri has several other intriguing and potentially related properties. The metal-rich stellar population is more centrally concentrated and has a lower velocity dispersion than the metal-poor population (Norris et al., 1997). Moreover, the metal-rich population has no detectable rotation, while the metal-poor population rotates relatively rapidly ($`7`$ km s-1), consistent with the apparent flattening of the cluster, which is among the highest of all globular clusters in the Galaxy (Frenk & Fall, 1982). The metal-rich stellar population also appears to be younger than the metal-poor population, consistent with the idea of self-enrichment over a period $`\mathrm{\Delta }\tau =35`$ Gyr (e.g., Hughes & Wallerstein 2000; Hilker & Richtler 2000; Lee et al. 1999). This is supported by the large enhancement in s-process elements in the metal-rich population relative to that in globular clusters with similar metallicities (Norris & Da Costa, 1995; Vanture, Wallerstein & Brown, 1994). These properties have led to the suggestion that $`\omega `$ Cen managed to retain the ejecta from the first generation of stars in the asymptotic giant branch (AGB) phase, which then formed a centrally concentrated, metal-enriched, younger generation of stars (e.g. Norris et al., 1996; Smith et al., 2000).

In this paper we demonstrate that $`\omega `$ Cen is not special in the ability to retain its own stellar ejecta. The escape speed in $`\omega `$ Cen is not the highest among the Galactic globular clusters. We find that more than half of the clusters could retain at least 20% of their AGB winds. It is therefore surprising that $`\omega `$ Cen is the only globular cluster in the Milky Way to show signs of self-enrichment.111 Sarajedini & Layden (1995) have inferred a dispersion of metallicity of 0.16 dex in M54, a globular cluster projected in the center of the Sagittarius dwarf galaxy. Also, Ree et al. (2001) have suggested that two other massive globular clusters (NGC 6388 and NGC 6441) may have an age and/or metallicity spread of 1.2 Gyr and 0.15 dex, respectively. These estimates, however, are not as large as those for $`\omega `$ Cen and may result from differential reddening or other systematic errors. Until these are confirmed, we assume that $`\omega `$ Cen is the only Galactic globular cluster with multiple stellar populations.  We conclude that this cluster could not have evolved in isolation on its present orbit in the Galaxy and that it must have been substantially different in the past. We discuss other models suggested previously, including the idea that $`\omega `$ Cen is the nucleus of a dwarf galaxy that lost its outer envelop by tidal stripping (Zinnecker et al., 1988; Freeman, 1993).

## 2 The escape velocity from the cluster

To compute the escape velocities of Galactic globular clusters, we use the photometric data from the catalog by Harris (1996). Out of 147 clusters in the catalog, we choose 141 for which the data are available on the distance, integrated magnitude, core radius, and concentration parameter. All clusters are fit by single-mass isotropic King models with a constant mass-to-light ratio, $`M/L_V=3`$ in solar units.

First, we check that the photometric models are consistent with the directly measured central velocity dispersions for a sample of 56 clusters from Pryor & Meylan (1993), their Table 2. The most reliable estimates of the dispersion are based on individual stellar radial velocities, but these estimates depend on the position of stars relative to the cluster center. The core velocity dispersion of $`\omega `$ Cen is 17 km s-1 (Merritt, Meylan & Mayor, 1997), with the indication that it rises in the center up to 25 km s-1 (van Leeuwen et al., 2000).

We thus avoid a direct comparison and show instead in Figure 1 two separate histograms of the central velocity dispersions as inferred from the photometric models and as measured directly. The distributions are qualitatively similar. Also, the photometric model of $`\omega `$ Cen predicts a central velocity dispersion of about 17 km s-1, the same as the observed value.

Fig. 1.— Distribution of the central line-of-sight velocity dispersion of the Galactic globular clusters, as inferred from the photometric models (left panel) and observed spectroscopically (right panel). The arrows point to the velocity dispersion of $`\omega `$ Cen.

We note that in the past 20 years there has been only one study (Dubath et al., 1997) providing homogeneous velocity dispersion measurements for a significant number of clusters. Now, with the availability of large-field survey spectrographs such data should be easier to obtain and would provide a very useful constraint on the mass models of globular clusters.

Having calibrated the photometric models, we can use them to determine the escape velocity. This is calculated using the reduced gravitational potential, $`v_{\mathrm{esc}}(2\mathrm{\Phi }_t2\mathrm{\Phi })^{1/2}`$, in the center (giving the maximum escape velocity) and at the half-mass radius (providing a threshold for most stars). Here, the tidal potential $`\mathrm{\Phi }_t`$ accounts for truncation at the tidal radius.

Figures 2 and 3 demonstrate that $`\omega `$ Cen is not special in its escape velocity. There are 11 clusters with higher $`v_{\mathrm{esc},0}`$ and 8 clusters with higher $`v_{\mathrm{esc},\mathrm{h}}`$.

We also find that the ratio of the central escape velocity to the central velocity dispersion can be fit by a simple relation

$$v_{\mathrm{esc},0}/\sigma _0=3.7+0.9(c1.4),$$
(1)
where $`c`$ is the concentration of the King models. The fit is accurate to 3% in the range $`c=0.42.8`$.

We now calculate the fraction of AGB winds retained by a cluster as a function of the wind velocity $`u_w`$. We integrate the distribution function of stars, $`f(E)`$, over the phase space volume within which the total velocity of the wind particles is below the escape speed at the star’s position:

$$f_w=\frac{1}{M}f(E)𝑑𝐫𝑑𝐯\frac{1}{4\pi }𝑑\mathrm{\Omega }_uh(v_{\mathrm{esc}}(r)\left|𝐯+𝐮_w\right|),$$
(2)
where $`\mathrm{\Omega }_u`$ is a solid angle of the wind velocity $`𝐮_w`$, and $`h(x)`$ is a unit step function.

Fig. 2.— Central escape velocity for 141 Galactic globular clusters vs their mass $`M`$. Omega Centauri is denoted by an asterisk. The escape velocity is calculated from the photometric single-mass King models with $`M/L_V=3`$.

Fig. 3.— Escape velocity at the half-mass radius for the Galactic globular clusters vs their mass $`M`$. Omega Centauri is denoted by an asterisk.

Figure 4 shows the result for clusters of various concentration parameters, from $`c=0.6`$ to 2.4. The fraction $`f_w`$ declines gradually and even the clusters with $`u_w3\sigma _0`$ can retain 10 to 20% of their winds. Since according to Loup et al. (1993) the characteristic terminal velocity of AGB winds is $`u_w15`$ km s-1 (although a velocity as high as 90 km s-1 has been measured for a field red giant by Dupree, Sasselov & Lester 1992), more than half of the Galactic clusters should be able to retain a sizable fraction of their winds.

Fig. 4.— Fraction of AGB winds retained by clusters of concentration $`c`$ (marked by numbers) as a function of the wind velocity in units of the central velocity dispersion. The top solid line shows the model with the highest binding energy ($`c=1.4`$), the dotted ($`c=0.6`$) and dashed ($`c=2.4`$) curves give the lower bounds. All models with intermediate concentrations fall between the plotted lines. For all concentrations, the fraction can be fit reasonably accurately by $`f_w(xu_w/\sigma _0)=1/[1+(x/2)^2+(x/4)^4+(x/3)^8]`$.

## 3 The orbit and the possibility of self-enrichment

The question now is: could $`\omega `$ Cen (and other clusters) retain the wind material throughout their evolution in the Galaxy? The oldest stars in the Galactic disk are dated to be at least $`9.5\pm 0.8`$ Gyr old (Oswalt et al., 1996; Leggett et al., 1998). The age of the Universe is determined fortuitously well from the coincidental CMB correlation (Knox et al., 2001) to be $`14.0\pm 0.5`$ Gyr. Thus, it is quite possible that the dense gaseous disk was already assembled when the Universe was about 4 Gyr old, equal to the age spread of stars in $`\omega `$ Cen. If this was the case, ram pressure of the disk would strip any gas left in $`\omega `$ Cen, as we show below.

The current Galactocentric velocity of $`\omega `$ Cen from the recently measured proper motion (Dinescu et al., 1999) is relatively low: $`V_R=31`$ km s-1, $`V_\theta =65`$ km s-1, $`V_z=4`$ km s-1, with the uncertainty $`\pm 10`$ km s-1. The speed of only 70 km s-1 makes the cluster strongly bound and thus, its present position at 6.3 kpc from the Galactic center is likely to be close to the orbital apocenter. Indeed, calculations of Dinescu et al. (1999) give the pericenter distance $`R_p=1.2`$ kpc and the apocenter distance $`R_a=6.4`$ kpc. Therefore $`\omega `$ Cen remains close enough to the Galactic center to pass through the disk twice in each orbital period, which is relatively short, $`P_{\mathrm{orb}}=1.2\times 10^8`$ yr.

We calculate the amount of gas that might have been present in $`\omega `$ Cen as follows. The total mass-loss rate of AGB stars depends on the slope of the IMF between 1 and 8 $`M_{\mathrm{}}`$, $`dN/dmm^\gamma `$, on the mass of the white dwarf remnant, $`m_{wd}`$, and on the main sequence lifetimes of stars, $`t_{ms}`$:

$$\frac{dM_w}{dt}=\frac{dN}{dm}(mm_{wd})\left(\frac{dt_{ms}}{dm}\right)^1,$$
(3)
where $`m=m_{ms}(t)`$ is the inversion of the function $`t_{ms}(m)`$. We consider two cases, $`\gamma =2.2`$ (following Kroupa, Tout & Gilmore 1993) and $`\gamma =3.3`$ (Scalo, 1986). Using the code of Hurley, Pols & Tout (2000), we calculate $`t_{ms}(m)`$ for the oldest stars in $`\omega `$ Cen, with the metallicity $`[Fe/H]=1.7`$, and find that in the mass range of interest it is well fit by $`t_{ms}5\times 10^9(m/M_{\mathrm{}})^{2.4}`$ yr. The masses of stellar remnants have been determined observationally in the nearby star clusters (Weidemann, 2000), and the result can be accurately fit by $`m_{wd}=0.444M_{\mathrm{}}+0.084m_{ms}`$. Combining all these expressions, we find the mass-loss rate at the age of 4 Gyr: $`dM_w/dt10^2M_{cl}`$ Gyr-1. Here $`M_{cl}3\times 10^6M_{\mathrm{}}`$ is the present mass of the cluster. (Interestingly, the mass-loss rate is almost the same for both values $`\gamma =2.2`$ and $`\gamma =3.3`$.) Thus, in half the orbital period the cluster would accumulate $`M_g=2\times 10^3M_{\mathrm{}}`$ of gas.

We integrate the orbit of $`\omega `$ Cen backward in time using an accurate Galactic potential that fits most observational constraints of the Milky Way (model $`A_1`$ of Klypin, Zhao & Somerville 2002). We model the evolution of the Galaxy starting with a spherical halo 14 Gyr ago, compress the baryons into the disk and the bulge over an exponential timescale of 3 Gyr and then turn them into stars over the next 3 Gyr. The bulge is gradually converted into the present flattened shape.

We find the locations where the cluster crosses the plane of the disk and then compare the ram pressure of the disk gas with the restoring gravitational force per unit area at the half-mass radius of $`\omega `$ Cen. In many passages $`\rho _{\mathrm{disk}}v_{\mathrm{rel}}^2>GM_{cl}M_g/(4\pi R_h^4)`$ and therefore the AGB ejecta should be stripped from the cluster. Also, the current mass of $`\omega `$ Cen is too small for dynamical friction to affect the orbit.

To summarize, we have shown that if (1) $`\omega `$ Cen had formed in isolation with its present mass, and (2) the Galactic disk was in place prior to the formation of the metal-rich stars, then the cluster would encounter the disk gas in its orbital motion and lose the gas accumulated from the enriched AGB winds. The same argument applies to all globular clusters within 10 kpc from the Galactic center. They could not have formed stars from the material enriched by s-process elements.

Alternatively, if all the stars in the oldest globular clusters form prior to the assembly of the Galactic disk, then $`\omega `$ Cen must be unrealistically special. Why did it retain its enriched gas and formed new populations of stars over 4 Gyr whereas a few dozen other clusters, capable of doing so, did not. As we have shown, even though $`\omega `$ Cen is the most massive cluster its potential well is similar to many others. So we are left again with the puzzle of why $`\omega `$ Cen is different from all other Galactic globular clusters.

We conclude that it is highly unlikely that $`\omega `$ Cen formed and evolved in isolation and enriched itself with heavy elements of the AGB winds.

## 4 In Closing: Towards a correct story of $`\omega `$ Cen

In light of the unusual properties of $`\omega `$ Cen, it is becoming a popular scenario that the cluster may be the core of a disrupted dwarf galaxy (e.g., Freeman 1993). The attraction is that a continuous infall of gas to the center of the dwarf galaxy may lead to a variable star formation history. Kormendy (1985) has shown that dwarf galaxies split on the plot of the central surface brightness vs core radius into dense dE galaxies (like M32) with the central velocity dispersion $`\sigma _0`$ in excess of 50 km s-1, and diffuse dSph galaxies (like Draco) with $`\sigma _0<20`$ km s-1. Tidal stripping of the outer envelope would decrease the total luminosity of the dwarf galaxy but would not significantly change its core properties, $`\mu _0`$ and $`R_c`$. Thus, the cores of dE would have too high a velocity dispersion and the cores of dSph would have too low density to count as globular clusters.

The counter-examples to this rule are nucleated dSph galaxies, sometimes classified as dE,N (Peterson & Caldwell, 1993). NGC 205 (Meylan, 2001) has a core velocity dispersion of 30 km s-1, dropping to 15 km s-1 in the bright compact nucleus (Carter & Sadler, 1990). The recently discovered Sagittarius galaxy apparently has a globular cluster M54 positioned at its center (Ibata, Gilmore & Irwin, 1994). Also, Drinkwater et al. (2000) detected five Ultra-Compact Dwarf galaxies in the Fornax cluster of galaxies, which are brighter than globular clusters but almost as compact. They could be the result of mergers of young star clusters, as suggested by Fellhauer & Kroupa (2002). It remains to be seen if the nuclei of these galaxies, when stripped of their envelopes, would have similar kinematic and chemical properties to $`\omega `$ Cen. Also, detailed orbital calculations are necessary to determine if the cluster can be brought in to its present location near the Galactic center (see Zhao 2001).

The uniqueness of $`\omega `$ Cen among the Galactic globular clusters still poses unanswered questions. In the hierarchical formation scenario we expect a few dozen progenitor galaxies to make up the Milky Way. If $`\omega `$ Cen did in fact form at the center of one of them, why did not other almost as massive clusters form in the other progenitors and have similar chemical enrichment? Were they entirely disrupted or was the history of $`\omega `$ Cen completely different? We intend to address these issues in forthcoming work.

We thank Chris Tout for discussions. OG is supported by the STScI Institute Fellowship, and JEP acknowledges the continuing support by the STScI visitors program. The electronic tables of the escape velocities, velocity dispersions, and main-sequence lifetimes are available at http://www.stsci.edu/$``$ognedin/gc/.


text:
# Probing the stellar populations of early-type galaxies: the SAURON survey

## 1. Introduction

In recent years very significant progress has been made in both the theoretical and observational study of the stellar populations of early-type galaxies. Stellar population synthesis models now rely on detailed stellar evolution tracks and much more complete stellar libraries at higher resolution (e.g. Vazdekis et al. 1996, Kodama & Arimoto 1997). Recent work on nearby samples of galaxies have often focused on the breaking of the age-metallicity degeneracy, following the widely used Lick/IDS indices to probe dynamically hot galaxies (Jørgensen 1999). Kuntschner (2000) emphasized the effect of non-solar abundance particularly when interpreting the observed line indices of early-type, metal rich, galaxies: Fornax ellipticals thus exhibit a relatively narrow spread in age, the distribution of galaxies in index/index diagrams being mostly attributed to differences in metallicity.

Although the conclusions reached by these studies do indeed converge towards a scenario consistent with the so-called hierarchical universe, the question originally posed by the Lick group two decades ago is still pertinent: when and how do early-type galaxies form? Most observational campaigns along these lines have considered line-strength gradients only, and where thus only able to probe a one-dimensional (long-slit) section of the galaxies. Two-dimensional spectroscopy has been mostly used for probing the kinematics of galaxies. Except in a handful of cases (e.g. Emsellem et al. 1996, Peletier et al. 1999, del Burgo et al. 2001), the available stellar absorption lines have only served a surprisingly restricted purpose.

One of the aims of the SAURON project is to provide the 2D spatial coverage which is missing. In this paper, we describe the approach the SAURON group is following, and present some preliminary results regarding the stellar populations for a few prototypical examples of galaxies in the SAURON sample.

## 2. From where, to where?

The essential ingredients of a good recipe for the formation and evolution of early-type galaxies should include a treatment of the star formation history, together with a scheme that links this to their dynamics and morphology. We first sum up briefly what we know about their chemical and dynamical properties.

### 2.1. What do we know?

Early-type galaxies are found preferentially in high local density regions. Their properties vary, from giant luminous ellipticals which are assumed to be triaxial, pressure supported objects with shallow nuclear brightness profiles, to low luminosity ellipticals, which exhibit disky isophotes, are consistent with isotropic rotators and include steep power law central cusps. In this context, the bulges of early-type spirals resemble low-luminosity ellipticals. Most if not all galaxies are thought to harbour central supermassive black holes, the masses of which correlate with the mean stellar dispersion of the spheroid (Ferrarese & Merritt 2000, Gebhardt et al. 2000).

Colours of ellipticals become redder with higher luminosity, an effect usually attributed to a parallel increase in metal enrichment and alpha element abundance. The star formation timescale in ellipticals is thought to be short, and ellipticals in general are viewed as old, coeval systems in which, as already emphasized, metallicity increases with luminosity. S0s have a larger spread in (luminosity weighted) age than ellipticals (Kuntschner & Davies 1998, Kuntschner 2000; see Fig. 1).

As we look back in time, star formation rates increase and interactions and mergers become more common. Already by $`z=0.5`$, the fraction of blue galaxies in rich clusters has increased dramatically compared to $`z=0`$, a trend known as the Butcher Oemler effect (Butcher & Oemler 1984). Furthermore HST studies (Dressler et al. 1997) have shown that these blue galaxies are spirals that are 2 to 3 more common in intermediate redshift clusters with lenticular galaxies being correspondingly less common. The signatures of mergers are observed in todays lenticulars and ellipticals as e.g. shells, kiloparsec scale decoupled cores, or counter-rotating gas disks.

### 2.2. What would we like to know?

Any scenario for the formation and evolution of early-type galaxies should explain the basic changes in their properties in terms of luminosity, Hubble type and environment. What information can we hope to dig out? Galaxies (besides our own) are observed projected on the sky, leaving us two spatial dimensions, and only one (line-of-sight) velocity.

The first obstacle is therefore to solve the inverse problem of determining the intrinsic shapes of these galaxies. In dynamical terms, we can push things further by determining the dominant orbit families which constitute their skeleton. These certainly hold traces of the involved formation and evolution processes, including passive secular evolution. Then questions, more directly related to the hierarchical merging scenario, follow: how common are decoupled cores, how common and massive are black holes in galactic nuclei?

Other integrated quantities such as colours and line indices can then help us to examine the stellar population content of these galaxies. How common are the disks in ellipticals, and is their stellar population different (younger?) from the rest of the galaxy. Can we constrain the stellar populations of the decoupled cores, and more generally of individual orbit families? To determine the origin of the different morphological/dynamical components we therefore need to reconcile the kinematics and stellar populations.

### 2.3. SAURON: the goals

The SAURON project was motivated by the questions summarized above. The minimal requirements included the two-dimensional mapping of the kinematics and line strengths (including at least one age and one metallicity indicator) of a reasonable sample of early-type galaxies up to at least one half $`R_e`$. We first designed a dedicated integral field spectrograph, to cover a field of view of 33$`\times `$41 arcsec2, and to work in the spectral domain around the Mg triplet at 5170 Å that includes Mg and Fe and H$`\beta `$ features as well as the \[OIII\] emission lines (see Bacon et al. 2001 for a description of the instrument). We then defined a representative sample of 72 galaxies, 24 for each E, S0 and Sa classes, evenly distributed in cluster and field type environments (de Zeeuw et al. 2002).

The SAURON observational campaign, representing more than 50 nights at the William Herschel Telescope (La Palma), is producing data with unprecedented uniformity and quality, and represents a unique database to unravel the stellar population - kinematics connection.

## 3. Age and composition of a galaxy

In practice, how do we determine the stellar population content of a galaxy, when it cannot be resolved into individual stars? The first step, adopted for stellar synthesis modeling, is to simulate the spectrum of a cluster whose stars share a single age and metallicity. Such a single burst spectrum can be built assuming an initial mass function (IMF), and following the stars along their evolution in the HR diagram. Each stellar spectrum is then matched using corresponding observed spectra. The final composite spectrum is the sum of the individual spectra weighted according to their respective contributions. This can be repeated for a grid of different ages and metallicities, and even different IMFs, to provide a set of template synthesized spectra. The observed galaxy spectrum can then be fitted, providing values for the luminosity-weighted mean age and metallicity.

Apart from the fact that a galaxy spectrum includes the contributions of stars with different ages and metallicities (integrated along the line of sight), the fitting of the full spectrum is sometimes not optimal as it may amplify errors in the flux calibration, does not necessarily focus on the relevant information, and requires very high signal-to-noise data. Since the age and metal enrichment are directly reflected in the depths of the absorption lines, these can be used as tracers of the stellar populations.

### 3.1. The Lick system

The Lick/IDS system of indices was thus designed with the purpose of focusing on important absorption features. These indices are easily measured, and can be compared with the values obtained from stellar synthesis models. Some indices are well known to be more age sensitive (e.g. hydrogen Balmer lines), or metalliticity sensitive (e.g. Fe lines).

However the Lick/IDS system was defined about 2 decades ago, and has a few limitations:

* The spectral resolution used for the Lick system is $`8.4`$ Å, diluting the available spectral signatures.
* The standard stellar spectra used to calibrate the Lick system were not flux calibrated.
* As defined, the Lick indices are rather sensitive to the required correction for velocity broadening.
* The measurements of the Lick stars are based on a non-linear detector.

These factors arise largely from the increase in efficiency and precision of modern instruments. Although we now see their limitations the Lick indices have played a pivotal role in the study of early-type galaxy populations for almost 30 years.

### 3.2. New models, indices and libraries

However, recent stellar synthesis models, such as those developed by Vazdekis and collaborators (Vazdekis 1999), do not just provide line indices, but also predict the spectrum over the SAURON range at a spectral resolution of 1.8 Å! This gives us the ability to define new, better adapted, indices, and to operate thorough spectral fitting of the galaxy data, even allowing for mixtures of different stellar populations.

In the context of the SAURON survey, we have therefore designed a new index, based on the original Fe5270 Lick one: it is less dependent on the broadening correction, and more sensitive to metallicity. New age-metallicity grids have been derived using the Vazdekis models which include up-to-date isochrones and stellar libraries. We have also observed a significant number of stars with SAURON, including high-metallicity ones, in order to compile a consistent library of stellar templates (Bacon et al. 2001).

## 4. Line strength maps with SAURON

### 4.1. Validation of the instrument

The first step before interpreting the SAURON line-strength maps is to validate the data itself by checking whether our line-strength measurements agree with previously published values. We have done this for a number of galaxies for which line-strengths are available in the literature, and Fig. 2 shows such a comparison in the case of NGC 5813. Overall, the agreement is excellent up to the edge of the SAURON field of view. There are some slight offsets between the different data sets, but these are within the quoted error bars.

### 4.2. Emission lines and how to get rid of them

One important issue when deriving line-strengths is to make sure that emission lines are not contaminating the spectra, as this can sometimes be a major source of error particularly for the H$`\beta `$ and Mgb indices (Goudfrooij & Emsellem 1996). In order to correct for the contribution of emission lines, we have designed a simple three-step procedure. We first fit the spectra using a complete spectral library which includes both stars and galaxies without nebular emission. This fit is achieved by masking the spectral regions where emission lines may be present. The fit obviously requires the knowledge of the stellar kinematics (velocity and dispersion), also obtained by avoiding the contaminated regions. We then subtract the fits to obtain pure emission line spectra, which are themselves fitted using different components for each emission line (H$`\beta `$, \[OIII\], \[NI\]). These emission-line fits are finally subtracted from the original spectra, (ideally) leading to emission-line free spectra. Line-strengths are measured from the emission-corrected spectra.

This procedure is illustrated in the case of NGC 5813, the central region of which exhibits the presence of emission line gas. Fig. 3 shows the H$`\beta `$ map before and after correction for the emission line contributions. In the former, there is a clear structure of lower H$`\beta `$ values which corresponds to a gaseous filament. As shown in the top right panel of Fig. 2, the H$`\beta `$ and \[OIII\]$`\lambda `$4959 emission lines are barely visible in the original spectra, but appear clearly in the residual pure emission-line spectra. After correction, the H$`\beta `$ map is flat over the full SAURON field of view, and shows no evidence for the presence of the decoupled core which is apparent only in the stellar velocity field.

## 5. SAURON results: two decoupled cores

In this Section, we present preliminary results of the SAURON observations of two decoupled cores, to illustrate the importance of two-dimensional spectroscopy in probing the stellar populations of early-type galaxies.

### 5.1. NGC 4365: the classical case

NGC 4365 is an elliptical known to harbour a kinematically decoupled core (Surma & Bender 1995). The kinematic major-axis of the core was assumed to lie at 90 degrees from that of the outer prolate-like body. However the SAURON kinematical maps revealed a more complex morphology with the zero velocity curve having a slowly varying position angle in the outer part (Davies et al. 2001). This strongly hints at a triaxial geometry the details of which can only be constrained by the full two-dimensional coverage provided by SAURON.

The line-strength maps lead to more surprising results. The Mgb increases towards the centre, with isocontours which closely resemble the galaxy isophotes. The H$`\beta `$ map is flat over the whole field, except perhaps for a slight unresolved peak at the very centre of the galaxy (Fig. 4). The stellar population of the core and the main body are found to be indistinguishable in age, metallicity and abundance ratios (Fig. 5). This strongly suggests that the decoupled core of NGC 4365 is the result of an early accretion, the core and the outer body of the galaxy being roughly coeval. It also tells us that the structure observed in NGC 4365 corresponds to a stable dynamical configuration.

### 5.2. NGC 4150: the post-starburst

The case of NGC 4150, an S0 galaxy, is at first sight even more puzzling. The stellar velocity field shows the presence of a small counter-rotating core with a low velocity amplitude of about 20 km/s peak to peak. The Mgb map of NGC 4150 shows a central dip, and an enhanced H$`\beta `$ in the central few arcseconds! When projected on a H$`\beta `$ versus \[MgFe5270\] model grid, this corresponds to a slightly decreasing metallicity and strongly decreasing age towards the centre (Fig. 6). It thus seems that the central region of NGC 4150 was the scene of a recent gaseous accretion event followed by a burst of star formation.

## 6. Conclusions

In this paper, we presented the cases of NGC 4365 and NGC 4150 to illustrate the importance of two-dimensional spectroscopy for studying the stellar populations of galaxies. However, the complexity and richness in the structures of early-type galaxies will only be fully revealed when the complete set of SAURON datacubes has been analysed.

The SAURON campaign will end in April ’02. We will then have a unique data set to probe the stellar content of early-type galaxies. We are also gathering new stellar libraries and state-of-the-art stellar synthesis models. As emphasized in this paper, this is a critical step if we wish to fully exploit the SAURON datacubes. The efforts of the SAURON team are therefore not only directed to acquiring the observations for the representative sample of 72 early-type galaxies, but are also focused on the development of new analysis tools and models. These tools will allow us to link the chemical and dynamical history of ellipticals and early-type spirals, and to ultimately constrain the formation and evolution scenarios.

## References

Bacon, R., Copin, Y., Monnet, G., et al. 2001, MNRAS, 326, 23

del Burgo, C., Peletier, R. F., Vazdekis, A., Arribas, S., Mediavilla, E. 2001, MNRAS, 321, 227

Butcher, H., Oemler, A. Jr., ApJ, 285, 426

Davies, R., Kuntschner, H., Emsellem, E., et al. 2001, ApJL, 548, 33

Dressler, A., Oemler, A. Jr., Couchm W. J., et al., 1997, ApJ, 490, 577

Emsellem, E., Bacon, R., Monnet, G., Poulain, P. 1996, A&A, 312, 777

Ferrarese, L., Merritt, D. 2000, ApJ 539, L9

Gebhardt, K., Bender, R., Bower, G., et al. 2000, ApJ, 539, L13

Goudfrooij, P., Emsellem, E. 1996, A&A, 306, 45

Kodama, T., & Arimoto, N. 1997, A&A, 320, 41

Kuntschner, H., Davies, R. 1998, MNRAS, 295, 29

Kuntschner, H., 2000, MNRAS, 315, 184

Jørgensen, I. 1999, MNRAS, 306, 607

Peletier, R. F., Vazdekis, A., Arribas, S., et al., 1999, MNRAS, 310, 863

Surma, P., Bender, R. 1995, A&A, 298, 405

Vazdekis, A., Casuso, E., Peletier, R. F., Beckman, J. E. 1996, ApJS, 106, 307

de Zeeuw, T., P., Bureau, M., Emsellem, E., 2002, MNRAS 329, 513


text:
# The continuum limit of quark number susceptibilities

## I Introduction

With the RHIC now in its second year of running, it is important to pin down lattice predictions for the high temperature phase of QCD. In this context, fully non-perturbative measurements of quark number susceptibilities gott ; gavai ; qnch ; 2fl  are important for four reasons. Firstly, they are directly related to experimental measurements of event-to-event fluctuations in particle production koch . Secondly, earlier results qnch ; 2fl  showed a strong jump in the susceptibility across the phase transition, but indicated a significant departure from weak-coupling behaviour, and it is important to check whether this persists into the continuum. Third, resummed perturbative computations of this quantity have now become available toni , thus making it possible to accurately test the importance of non-perturbative contributions to this quantity. Finally, with continuum extrapolated results in hand, one can address the question of whether the strangeness production seen in heavy-ion collisions can be quantitatively explained as a signal of the quark-gluon plasma.

We have recently presented systematic results for quark number susceptibilities in quenched QCD qnch  as well as for QCD with two flavours of light dynamical quarks 2fl . In these studies a large range of temperatures, $`T`$, was covered at a series of different quark masses, $`m`$. These computations were done at a fixed cutoff with the lattice spacing, $`a=T/4`$ while keeping the finite volume effects under control so that the thermodynamic limit could be taken reliably. There was a 3–5% difference between the quenched and dynamical computations. In this paper we examine the cutoff dependence of the quenched results and report their continuum (zero lattice spacing) limit.

The partition function of QCD with three quark flavours is

$$Z(T,\mu _u,\mu _d,\mu _s)=𝒟U\mathrm{exp}[S(T)]\underset{f=u,d,s}{}detM(T,m_f,\mu _f),$$
(1)
where the temperature $`T`$ determines the size of the Euclidean time direction, $`S(T)`$ is the gluonic part of the action and the determinants of the Dirac matrices, $`M`$, contain as parameters the quark masses, $`m_f`$, and the chemical potentials $`\mu _f`$ for each of the flavours $`f=u,d,s`$. We also define the chemical potentials $`\mu _0=\mu _u+\mu _d+\mu _s`$, $`\mu _3=\mu _u\mu _d`$ and $`\mu _8=\mu _u+\mu _d2\mu _s`$, which correspond to the diagonal flavour $`SU(3)`$ generators. Note that $`\mu _0`$ is the usual baryon chemical potential and $`\mu _3`$ is an isovector chemical potential.

Quark number densities are defined as

$$n_i(T,\mu _u,\mu _d,\mu _s)=\frac{T}{V}\frac{\mathrm{log}Z}{\mu _i},$$
(2)
and the susceptibilities as

$$\chi _{ij}(T,\mu _u,\mu _d,\mu _s)=\frac{T}{V}\frac{^2\mathrm{log}Z}{\mu _i\mu _j},$$
(3)
where $`V`$ denotes the spatial volume. The subscripts $`i`$ and $`j`$ are either of the index sets $`f`$ or $`\alpha =0`$, 3 and 8. We lighten the notation by writing the diagonal susceptibilities $`\chi _{ii}`$ as $`\chi _i`$. We determine the susceptibilities at zero chemical potential, $`\mu _f=0`$. In this limit, all $`n_i(T)=0`$. Since we work with $`m_u=m_d=m<m_s`$, we also have $`\chi _{03}=0`$.

Flavour off-diagonal susceptibilities such as

$$\chi _{ud}=\left(\frac{T}{V_3}\right)\mathrm{tr}M_u^1M_u^{}\mathrm{tr}M_d^1M_d^{}$$
(4)
are given entirely in terms of the expectation values of disconnected loops. Since $`m_u=m_d`$, we obtain $`\chi _{us}=\chi _{ds}`$ with each defined by an obvious generalization of the formula above. Of the flavour diagonal susceptibilities we shall use

$$\chi _s=\left(\frac{T}{V_3}\right)\left[\left(\mathrm{tr}M_s^1M_s^{}\right)^2+\mathrm{tr}\left(M_s^1M_s^{\prime \prime }M_s^1M_s^{}M_s^1M_s^{}\right)\right].$$
(5)
$`\chi _u=\chi _d`$ are given by a generalization of this formula. Numerically, the simplest quantity to evaluate is the diagonal iso-vector susceptibility

$$\chi _3=\frac{1}{2}\left(\frac{T}{V_3}\right)\mathrm{tr}\left(M_u^1M_u^{\prime \prime }M_u^1M_u^{}M_u^1M_u^{}\right).$$
(6)
Two more susceptibilities are of interest. These are the baryon number and electric charge susceptibilities,

$$\chi _0=\frac{1}{9}\left(4\chi _3+\chi _s+4\chi _{ud}+4\chi _{us}\right)\mathrm{and}\chi _q=\frac{1}{9}\left(10\chi _3+\chi _s+\chi _{ud}2\chi _{us}\right).$$
(7)
Note that $`\chi _0`$ is the baryon number susceptibility for three flavours of quarks. As a result, this expression differs from the iso-singlet quark number susceptibility for two flavours, defined in gott , both in overall normalization and by terms containing strangeness. In our numerical work we have chosen to use staggered quarks. Hence, to normalise to one flavour of continuum quarks and compensate for fermion doubling on the lattice, we have to multiply each of the traces in eqs. (46) by a factor of $`1/4`$.

The quark mass appears in two logically distinct places— first in the operators which define the susceptibilities in eqs. (46), and secondly in the determinant in the partition function of eq. (1) which defines the weight for the averaging of these operators. The first defines a valence quark mass and the second the sea quark mass. In this paper we adopt the quenched approximation, whereby the determinants in eq. (1) are set equal to unity. However, we shall vary the valence quark masses. As shown earlier, unquenching the sea quarks changes results by 3–5% qnch ; 2fl .

Note that all the flavour off-diagonal susceptibilities are exactly zero in an ideal gas. We shall denote by $`\chi _{FFT}^3`$ the ideal gas value for $`\chi _3`$. On an $`N_t\times N_s^3`$ lattice with $`N_c`$ colours and lattice spacing $`a`$, we find

$$a^2\chi _{FFT}^3=\frac{N_c}{2N_tN_s^3}\underset{p}{}\left\{D^2\mathrm{sin}^2p_0\mathrm{cos}^2p_0+D^1\left(\mathrm{sin}^2p_0\frac{1}{2}\right)\right\},$$
(8)
where the spectrum of momenta is $`p_0=(2\pi /N_t)(n+1/2)`$ with $`0n<N_t`$ and $`p_i=2\pi n/N_s`$ with $`0n<N_s`$, and $`D=(ma)^2+_\mu \mathrm{sin}^2p_\mu `$. For a given $`m/T_c`$ and $`T/T_c`$, the value of the quark mass in lattice units is

$$ma=\left(\frac{1}{N_t}\right)\left(\frac{m}{T_c}\right)\left(\frac{T_c}{T}\right).$$
(9)
The formula in eq. (8) is normalised such that there is exactly one copy of each flavour of quarks. The values of $`\chi _s`$, $`\chi _0`$ and $`\chi _q`$ for an ideal gas can then be simply obtained as linear combinations of $`\chi _{FFT}^3`$ for different valence quark masses.

## II Lattice results

We have previously reported measurements of $`\chi _3`$ in quenched QCD with 4 time slices ($`N_t=4`$) at temperatures of $`1.1T_c`$, $`1.25T_c`$, $`1.5T_c`$, $`2T_c`$ and $`3T_c`$ gavai . We have extended these measurements to higher temperatures through two simulations on $`4\times 20^3`$ lattices at $`\beta =6.6`$ and $`6.7`$. These correspond to temperatures of $`4.13T_c`$ and $`4.7T_c`$. All these computations are performed at fixed bare quark masses of $`m/T_c=0.03`$, $`0.3`$, $`0.5`$, $`0.75`$ and 1 (the data for $`m/T_c=0.5`$ are new).

The continuum limit has been taken by going to several smaller lattice spacings. One set of computations is performed with a lattice spacing $`a=1/6T`$ which is 33% smaller than the lattice spacing used earlier. We have taken data on $`6\times 20^3`$ lattices at $`\beta =6.0625`$, $`6.3384`$ and $`6.7`$ corresponding to $`T/T_c=1.333`$, 2 and $`3.133`$, and on a $`6\times 24^3`$ lattice at $`\beta =6.8`$ corresponding to $`T/T_c=3.532`$. One further computation was made on an $`8\times 18^3`$ lattice at $`\beta =6.55`$. This corresponds to $`T/T_c=2`$ with lattice spacing $`a=1/8T`$. The quark mass $`m/T_c`$ is kept independent of $`a`$ and $`T`$, and hence $`ma`$ decreases with increasing $`T`$ at fixed $`N_t`$ or with increasing $`N_t`$ at fixed $`T`$ according to eq. (9).

The lattice scale has been set using the plaquette measurements of bielef  and the analysis performed in myold . The conversion of $`\beta `$ to $`T`$ using data for $`N_t=4`$ and 8 give the two curves in Figure 1. The difference between the two curves is within the statistical uncertainty of 5% on the scale assignment. There is the same degree of statistical uncertainty in the bare quark masses and in all other scales. We shall not display this inherent scale uncertainty in the measurements, but it should be kept in mind.

The simulations have been performed with a Cabbibo-Marinari pseudo-heatbath technique with 3 $`SU(2)`$ subgroups updated on each hit. An initial 1000 sweeps have been discarded for thermalisation. On the $`N_t6`$ lattices we have used 80 configurations separated by 1000 sweeps for the measurements. On the $`N_t=8`$ lattice we have used 55 configurations separated by 500 sweeps.

We have previously checked that the spatial lattice size, $`L=N_sa`$, is quite irrelevant to the values of $`\chi `$ measured as long as it is sufficiently large compared to the inverse pion mass gavai ; 2fl . In fact, at all the couplings and quark masses we have used, $`m_\pi L5`$. The most important other constraint on $`L`$ comes from the requirement that it should be large enough to prevent spatial deconfinement saumen . This is ensured by taking $`N_s/N_tT/T_c`$ in all our simulations.

The traces are evaluated by the usual stochastic technique,

$$\mathrm{Tr}A=\frac{1}{2N}\underset{i=i}{\overset{N}{}}R_i^{}AR_i,$$
(10)
where $`R_i`$ are a set of $`N`$ uncorrelated vectors with components drawn independently from a Gaussian ensemble with unit variance. Each vector has three colour components at each site of the lattice. Since we use a half lattice version of the Dirac operator for staggered Fermions, the number of components of each vector $`R_i`$ is one and a half times the number of lattice points. $`(\mathrm{Tr}A)^2`$ is obtained by dividing the set of $`R_i`$ into disjoint blocks, constructing $`\mathrm{Tr}A`$ in each block, taking all possible distinct pairs of such estimates, multiplying them and then averaging over the pairs. For $`N_t=4`$ it is possible to get accurate results with $`N10`$, although we have chosen to use $`N=80`$ for each gauge configuration 2fl .

It is easy to check that the variance of the estimator above is proportional to $`\mathrm{Tr}(A^2)`$. Since the diagonal element of the Fermion matrix is proportional to the quark mass $`ma`$, whereas the off-diagonal elements are bounded by unity, with increasing $`T`$ or decreasing $`a`$, $`\mathrm{Tr}A`$ decreases linearly, while its variance remains constant. Hence, for sufficiently small $`ma`$ the number of vectors has to increase quadratically with $`ma`$. We found that $`N_v=100`$ was sufficient for $`N_t=8`$ at $`T=2T_c`$, although at larger $`N_t`$ or $`T`$ significantly larger values of $`N_v`$ are needed. A further numerical problem arises from the cancellation of the two matrix elements which give $`\chi _3`$ (see eq. 6). While increasing $`N_t`$ and keeping the physical size of the lattice fixed, each of these terms increases quadratically but mutually cancel to give a number which decreases quadratically with $`N_t`$. At fixed word length, this implies a reduction in accuracy by a factor of $`N_t^4`$.

In Figure 2 we show how the continuum limit of $`\chi _3`$ is reached. There is an 8% decrease in $`\chi _3/\chi _{FFT}^3`$ in going from $`N_t=4`$ to $`N_t=6`$ lattices. However, within measurement errors, the continuum limit is reached already at $`N_t=6`$. This implies a cancellation of the $`N_t`$-dependence of $`\chi _3`$ and $`\chi _{FFT}^3`$. The reason for the atypical behaviour for $`N_t=4`$ is not in the numerator, which is non-perturbative, but in the denominator, which is the ideal gas.

For $`N_t=4`$, the spectrum of $`p_0`$ consists of $`\{\pm \pi /4,\pm 3\pi /4\}`$. As a result, $`\mathrm{sin}^2p_0=1/2`$ and the second term in the sum on the right of eq. (8) vanishes identically. This accident occurs only for $`N_t=4`$, and for all other values of $`N_t`$ this term contributes a non-vanishing value. It is easy to see that in the limit $`N_t\mathrm{}`$ this term contributes to the leading part of $`\chi _{FFT}^3`$. This is the reason the values for $`\chi _3/\chi _{FFT}^3`$ for $`N_t=4`$ lie well away from the continuum limit.

Our complete results for the dependence of $`\chi _3/\chi _{FFT}^3`$ on $`T/T_c`$ at several different quark masses for $`N_t=4`$ and 6 are shown in Figure 3. In view of the data shown in Figure 2, our results for $`N_t=6`$ are the continuum limit for the quenched theory: these are summarised in Table 1. Note that even at temperature as high as $`3.5T_c`$ the result is far below the ideal gas limit. Closer to $`T_c`$ there is an even more dramatic fall in the value of $`\chi _3/\chi _{FFT}^3`$. We have also measured the off-diagonal susceptibility $`\chi _{ud}`$. As already seen for $`N_t=4`$ this is consistent with zero. In Table 1 we show that the continuum limit of $`\chi _{ud}`$ also vanishes within errors of $`10^7T^2`$. It is therefore much smaller than the 3-loop perturbative result given in toni .

Note that the continuum limit of our measurements lie 20–30% below the ideal gas result. They also differ significantly from HTL as well as the skeleton graph resummed results of toni . At $`T=3T_c`$, HTL predicts $`\chi _3/\chi _{FFT}^3=0.90`$$`0.94`$ on varying the scale of $`\alpha _s`$ between $`\pi T`$ and $`4\pi T`$, and the resummed computation gives $`\chi _3/\chi _{FFT}^3=0.95`$$`0.97`$. Our measurement shows $`\chi _3/\chi _{FFT}^3=0.80`$ at $`3T_c`$ (see Table 1). A second computation of the HTL result is available for $`N_f=2`$ mustafa . Since this agrees with the HTL result of toni  for $`N_f=2`$, we believe the HTL result for $`\chi _3/\chi _{FFT}^3`$ is under control. There is thus a genuine discrepancy between these lattice results and existing perturbative computations.

The quark number susceptibilities are closely related to the screening correlator of a one-link separated quark bilinear operator which corresponds to a $`\rho `$ meson at zero temperature. The $`T>0`$ transfer matrix of the problem mixes this with a quark bilinear that corresponds to the $`\pi `$ old . We have earlier shown evidence for $`N_t=4`$ that $`\chi _3`$ is closely related to the only known non-perturbative quantity among the screening correlators, that is the screening mass $`M_S`$ coming from the two degenerate correlators that descend from the $`\pi `$ and $`\sigma `$ operators of the zero temperature theory 2fl . In view of the strongly non-perturbative character of $`\chi _3`$, as revealed by the comparison with HTL and resummation explained above, examining this correlation in the continuum becomes more significant. In Figure 4 we show data obtained on $`8\times 18^3`$ lattices which indicate that this correlation survives into the continuum.

## III Strangeness

Measurements of hadron yields in relativistic heavy-ion collisions have been analysed extensively heidel ; strange  and the observed enhancement of strangeness has been claimed to be evidence for quark-gluon plasma formation. The chemical composition is observed at the hadron freezeout temperature, which has been found to be very close to the value of $`T_c`$ in full QCD. It has been argued that this chemical composition cannot arise due to hadronic rescattering heinz . It would be interesting if it could be directly determined whether or not the composition is characteristic of the plasma for $`TT_c`$. We argue below that this question can be answered by a lattice QCD measurement such as we have presented above.

The quantity which seems to be directly connected to lattice measurements is the Wroblewski parameter wrob , $`\lambda _s`$, which measures the ratio of newly created primary strange to light quarks—

$$\lambda _s=\frac{2s\overline{s}}{u\overline{u}+d\overline{d}}.$$
(11)
A fluctuation-dissipation theorem relates the rate of production of quark pairs by the plasma in equilibrium to the imaginary part of the full complex quark number susceptibilities fw . If the inverse of the characteristic time scales of the QCD plasma are not close to the typical energy scales for the production of strange and lighter quarks, then the ratio of their production rates is just the ratio of the static susceptibilities that we have measured. Then, if the observed chemical composition is created in equilibrium, we should have

$$\lambda _s=\frac{2\chi _s}{\chi _u+\chi _d}\frac{\chi _s}{\chi _u},$$
(12)
where the last equality holds in the limit of equal $`u`$ and $`d`$ quark masses. These susceptibilities have to be evaluated on the lattice at the temperature and chemical potential, $`\mu _0`$, relevant to the collision. However, it has been shown recently cleymans  that $`\lambda _s`$ is insensitive to $`\mu _0`$ for SpS and RHIC energies. Furthermore a lattice measurement shows that the screening mass, $`M_S`$, does not change rapidly with increasing chemical potential miyamura . Then, in view of the correlation between $`M_S`$ and $`\chi `$ shown in Figure 4, we do not expect rapid changes with chemical potential in the ratio above. Both these arguments show that as a first approximation one can take the ratio in eq. (12) at zero chemical potential.

From our measurements of $`\chi `$ reported here, we can form the lattice “prediction” of $`\lambda _s`$ created in the plasma at $`T_c`$. This has been done by taking our measurements at $`T/T_c=4/3`$ and 2, and extrapolating them linearly to $`T_c`$. For $`\chi _u`$ we use data taken for $`m/T_c=0.03`$, which is light enough for this purpose. Since scaled quantities feel the smallest correction for unquenching 2fl , we take $`\chi _s`$ to correspond to $`m/T_c=1`$, which is appropriate to full QCD. In Figure 5 the resultant prediction is compared with the data collected in strange ; cleymans . It can be seen that there is a fair agreement with the data. However, this is subject to several assumptions—

1. The foremost assumption is that the characteristic time scales of the plasma are not close to the inverse energy scale of the production processes. It is not possible to test this assumption in an Euclidean computation. However, it has been suggested that these characteristic time scales could be observed in real or virtual photon production as spikes in the spectrum spikes . If these are seen and they lie in the energy range for the production, the assumption would be falsified.
2. Another important assumption is that chemical equilibration takes place in the plasma. If chemical equilibrium is not achieved but the light degrees of freedom achieve energy equipartition, then the effective temperature could be higher.
3. We have assumed slow variation of the ratio $`\chi _s/\chi _u`$ with chemical potential. While this seems to be a reasonable estimate in view of the results of cleymans ; miyamura  and our own observations about the relation between the lowest screening mass and the susceptibilities, one should keep in mind possible changes in the ratio due to violation of this assumption.
4. Our continuum results are obtained in quenched QCD, whereas the data corresponds to full QCD. As mentioned already, we have taken $`m/T_c=1`$, as appropriate to full QCD, in order to correct for this. However, an additional 5–10% shift in our results could arise when unquenching 2fl . This is not shown in Figure 5.
5. Our measurements are made away from $`T_c`$ and extrapolated down to this temperature. In principle this can be corrected by a computation made directly at $`T_c`$. However, due to the enormous cost of extracting physical quantities at $`T_c`$, we postpone this to the time when we attempt a full QCD simulation.

It might seem surprising that $`\lambda _s`$ corresponds to the chemical composition at $`T_c`$. However, a moment’s thought shows that this is entirely in accord with the assumption of chemical equilibrium. Even if the system is thermalized at a much higher temperature, it remains in equilibrium until it reaches $`T_c`$. Departure from equilibrium as the system cools through the transition is unlikely, and if it occurs, should have a visible signature such as the formation of a disoriented chiral condensate rajag . In the absence of such phenomena, if $`\lambda _s`$ does not correspond to the chemistry at $`T_c`$, then the second assumption above would be demonstrated to be false. Also, with increasing ion energy the system moves closer to zero baryon number, as seen at RHIC. This makes our third assumption better at the RHIC and LHC. The fourth assumption can be removed by a computation with full QCD. Such a computation is planned, where also the final assumption can be removed, and will be reported elsewhere.

## IV Summary

In summary, we have presented new and precise results on quark number susceptibilities over a wide range of temperatures and quark masses in the high temperature phase of QCD. The main aim of this study was to obtain the continuum extrapolation of these susceptibilities. As shown in Figure 2, while there is a significant change in going from $`N_t=4`$ to $`N_t=6`$ lattices, there is no statistically significant change in the susceptibilities beyond this. The results obtained for $`N_t=6`$ lattices can then be taken as the continuum limit. This limit is shown in Figure 3 and Table 1.

There is a strong discrepancy between the continuum extrapolated lattice results and HTL computations for these quantities— varying between about 12% at the highest $`T`$ to even larger at smaller $`T`$. The off-diagonal susceptibility, $`\chi _{ud}`$, is zero, in contrast to the HTL results. There is a somewhat larger discrepancy between the continuum extrapolated lattice results and the skeleton graph resummed results. The conjecture that there is a strong non-perturbative component to the quark number susceptibilities is supported by an observed strong correlation between the smallest quark bilinear screening mass, $`M_S`$, and the susceptibility $`\chi _3`$ (Figure 4).

It is interesting to note that the continuum extrapolated results for the strange and light quark susceptibilities can be used to give a surprisingly good description of the chemical composition of hadrons at freezeout in SpS and RHIC experiments (Figure 5). This has to be treated as a preliminary estimate due to the many caveats which we have listed, and some of which we plan to check in future lattice measurements.