text:
# FORMATION OF PROTO-GLOBULAR CLUSTER CLOUDS BY THERMAL INSTABILITY

## I INTRODUCTION

Globular clusters (GCs) are the oldest relics in our galaxy, presenting a cosmological challenge, as their ages often burst the allowed bounds of the expansion age. Since there are at least three distinct populations of GCs, the old halo and the middle-aged disk populations of GCs in our Galaxy (Zinn 1985), and the young GCs in interacting galaxies (Ashman & Zepf 1992), we may well need more than one model to explain their formation and evolution. In the present study, we focus on the formation of the halo GCs. Half of the total mass of our galaxy appears to be contained within $`50\mathrm{kpc}`$, whereas roughly half of the luminosity of the disk as well as half of the total number of halo GCs are distributed within the solar galactocentric distance ($`8\mathrm{kpc}`$). So halo GCs clearly know about the dissipation that followed galaxy formation and we must place their formation in the protogalactic context. This encourages models that form halo clusters during the formation of the protogalaxy, such as thermal instabilities in the protogalactic halo (Fall & Rees 1985) or pregalactic cloud collisions (Gunn 1980; Lake 1987; Kang et al. 1990; Kumai, Basu, & Fujimoto 1993).

Many models of globular cluster formation assume the presence of proto-globular cluster clouds (PGCCs) in pressure equilibrium with hot halo gas (e.g., Gunn 1980; Dopita & Smith 1986; Brown, Burkert & Truran 1991; Kumai, Batsu, & Fujimoto 1993). The model by Fall & Rees (1985, FR85 hereafter) has been most widely adopted to explain the formation of such PGCCs: the thermal instability drives catastrophic cooling of over-dense regions in the hot halos of protogalaxies and condenses them to dense clouds. In the present study, adopting FR85 model, we follow the thermal and dynamical evolution of over-dense clouds by numerical simulations in order to examine more quantitatively the physical properties of resulting PGCCs.

Although the thermal instability was extensively investigated in several areas of astrophysics, in recent years it has received much attention especially for X-ray cluster cooling flows: (cf. Balbus 1986; Malagoli, Rosner & Bodo 1987; Balbus & Soker 1989; Malagoli, Rosner, & Fryxell 1990). According to the results from numerical simulations in one-dimensional (1D) geometries (David, Bregman & Seab 1988 \[DBS88 hereafter\]; Brinkmann, Massaglia & Müller 1990 \[BMM90\]), the over-density undergoes a quasi-static compression in near pressure equilibrium when a cloud is small enough to adjust to pressure change faster than it cools, while larger clouds may undergo a “supersonic compression stage”, leading to a density increase much higher than what is expected from the isobaric compression. The plane-parallel calculations by DBS88 showed that the density contrast increases to the “isobaric” ratio of the temperature of the hot background medium ($`T_h10^710^8\mathrm{K}`$) and the minimum temperature ($`T_c10^4\mathrm{K}`$) where the cooling becomes inefficient. On the other hand, the spherically symmetric calculations of BMM90 showed that the density contrasts can be more than three orders of magnitude higher than those found in the plane-parallel calculations of DBS88. The difference between two geometries is mainly due to the fact that the gravitational acceleration increases as $`GM/R^2`$ when the spherical cloud collapses, while it remains the same (independent of $`R`$) when the plane sheet compresses. Also isotropic spherical compression (focusing) leads to much higher density accelerating the cooling and rise in density. However, 2D simulations of the collapse of an elongated blob by BMM90 showed a shape instability to a pancake-shape leading to evolution best represented by the 1D plane-parallel case. Thus geometry is clearly important in the nonlinear evolution. On the other hand, 2D simulations by Vázquez-Semadeni, Gazol, & Scalo (2000) showed that the condensation of fluctuating density field via thermal instability leads to a network of transient filaments which later accrete onto isolated round blobs. These blobs form at intersections between filaments just like galaxy clusters in large-scale structure formation in cosmological simulations. This implies that the 1D spherical simulation can actually represent a generic collapse of the highest density peaks in the randomly fluctuating density field in 3D.

Numerical simulations in a 2D stratified background (Yoshida, Hattori, & Habe 1991; Reale et al. 1991 and references therein), on the other hand, have shown that the dynamic instabilities, such as Rayleigh-Taylor and Kevin-Helmholtz instabilities due to buoyancy-driven oscillations can disrupt the thermally unstable bubbles if the timescale for buoyancy oscillation (Brunt-Väsälä period) is shorter than the local radiative cooling time. Thus it is believed that the classic theory of thermal instability can not explain the effective mass depletion from the cooling flows of galactic clusters. We will examine the role of such dynamic instabilities for the formation of PGCCs in the protogalactic halo environment in the next section. The details of our models are given also in this section. The main results are in §3. Conclusions are given in §4.

## II MODEL CALCULATIONS

### (a) Protogalactic Halo Model

We adopt an idealized model for the protogalactic halo gas as in FR85; an isothermal sphere based on the Milky Way’s properties, with a circular velocity $`V_c=220\mathrm{km}\mathrm{s}^1`$ and corresponding gas temperature $`T_h=1.7\times 10^6\mathrm{K}`$.

Many previous studies used a gas distribution specified by the condition that the free fall time to the center of the protogalaxy equals the cooling time, yielding the halo gas density at a distance $`R_g`$ from the center $`\rho _11.7\times 10^{25}\mathrm{g}\mathrm{cm}^3(R_g/10\mathrm{kpc})^1`$. The use of $`\rho _1`$ has lead previous workers to consider halo densities of $`n_h=0.1\mathrm{cm}^3`$ (FR85; Kang et al. 1990) or even $`n_h=1\mathrm{cm}^3`$ (Vietri & Pesce 1995). For a primordial gas of $`\mathrm{H}`$ and $`\mathrm{He}`$ with an assumed ratio $`n(\mathrm{He})/n(\mathrm{H})=0.1`$, the halo gas mass density is $`\rho _h=(2.34\times 10^{24}\mathrm{g}\mathrm{cm}^3)n_h`$. For the isothermal sphere in hydrostatic equilibrium, however, the gas density is given by $`\rho _2=5.8\times 10^{25}\mathrm{g}\mathrm{cm}^3f_B(R_g/10\mathrm{kpc})^2`$, where $`f_B`$ is the ratio of the baryonic to total matter density. Realistic values should be $`0.02<n_h<0.07\mathrm{cm}^3`$ for $`R_g10\mathrm{kpc}`$ (assuming $`f_B0.1`$) with some degree of central concentration. A value as high as $`n_h=1\mathrm{cm}^3`$ is greater than the total density of baryonic and dark matter. Thus we take $`0.05\mathrm{cm}^3`$ as a fiducial value and express the halo gas density as

$$n_h(0.05\mathrm{cm}^3)R_{10}^q,$$
(1)
where $`R_{10}`$ is $`R_g`$ in units of $`10\mathrm{kpc}`$ and $`1<q<2`$.

Although the density irregularities in protogalaxies likely result from a superposed spectrum of perturbations, we consider, for simplicity, isolated spherically symmetric clouds resulted from sinusoidal perturbations. So the initial gas density of the cloud $`n_{cloud}`$ decreases gradually with radius $`r`$ from its center to its edge $`R_c`$ as follows:

$$n_{cloud}=n_h\left[1.0+\delta \mathrm{cos}\left(\frac{\pi }{2}\frac{r}{R_c}\right)\right]\mathrm{for}r<R_c$$
(2)
where $`\delta `$ is the amplitude of the initial density perturbation and $`R_c`$ is the cloud radius.

The initial temperature throughout the cloud is set by the condition of hydrostatic equilibrium. To start in the nearly linear regime, we always begin with $`\delta 1`$, unless otherwise noted. Initial parameters of the model clouds are summarized in Table 1.

### (b) Time Scales

The free fall time $`t_{ff}^{pg}`$ is the time it takes a particle dropped from $`R_g`$ to reach the center of the protogalaxy,

$$t_{ff}^{pg}=\left(\frac{\pi }{2}\right)\left(\frac{R_g}{V_c}\right)=5.6\times 10^7\mathrm{yrs}\left(\frac{R_g}{10\mathrm{kpc}}\right)_.$$
(3)
The free fall time of a uniform gas cloud including only its self-gravity is

$$t_{ff}^{sg}=\left(\frac{3\pi }{32G\rho }\right)^{1/2}=4.4\times 10^6\mathrm{yrs}\left(\frac{n}{100\mathrm{cm}^3}\right)_,^{1/2}$$
(4)
where $`n`$ is the cloud density and any dark matter contribution is ignored.

The cooling time for the halo gas of $`T_h`$ is given by

$$t_{cl,h}=\frac{1.5n_hkT_h}{\mathrm{\Lambda }}(4.0\times 10^7\mathrm{yrs})R_{10}^q,$$
(5)
where the cooling rate of $`\mathrm{\Lambda }/n_h^2=5.5\times 10^{24}\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1`$ at $`T_h=1.7\times 10^6`$ and the halo density given in Eq. (1) are used to calculate the latter relation. For an isobaric perturbation of the density contrast $`\delta `$, the cooling timescale of the perturbed gas becomes $`t_{cl,c}t_{cl,h}(1+\delta )^2`$, considering that the cooling rate is almost constant for temperatures within factor of two of $`T10^6\mathrm{K}`$. We define the cooling distance as the distance over which the sound wave of the hot gas can travel within the cooling timescale of the perturbed gas,

$$l_{cool}c_ht_{cl,c}(8.1\times 10^3\mathrm{pc})R_{10}^q(1+\delta )_,^2$$
(6)
where $`c_h=198\mathrm{km}\mathrm{s}^1`$ is the sound speed of the hot gas of $`T_h`$.

The ratio of the cooling time to the self-gravity time for initial clouds is

$$\eta =6.2\times 10^2R_{10}^{q/2}\left(\frac{1+\delta }{2}\right)_.^{3/2}$$
(7)
For the range of density considered here, $`0.01<n_h<1\mathrm{cm}^3`$, cooling takes place too fast for self-gravity to be important initially. However, the cooling becomes inefficient once the gas cools to $`10^4`$K, and the gravitational time scale decreases as the density increases. The density in the condensed clouds increases by a factor of $`\chi _{isob}=(T_h/T_c)(\mu _c/\mu _h)350`$ for an isobaric case, where $`\mu _c`$ and $`\mu _h`$ are the molecular weights of the cold and hot gas, respectively. The density enhancement is even higher for supersonic collapses, resulting in the cloud pressure higher than the background pressure. The typical additional enhancement in supersonic collapses is $`10^210^{3.5}`$ (see the next section). Thus, self-gravity becomes dynamically important in late stages. Also note that self-gravity would be more important for lower density environments, since $`\eta n_h^{1/2}`$.

The growth of thermal instabilities in a stratified medium is affected by the buoyancy-driven dynamical instabilities (see below), so the local buoyancy oscillation period is also important. For our halo model it can be estimated by the local scale height $`H`$ divided by the sound speed $`c_h`$,

$$t_b\frac{H}{c_h}(5.0\times 10^7\mathrm{yrs})R_{10},$$
(8)
since $`HR_g`$ for the power-law density distribution assumed here. For linear initial perturbation (e.g., $`\delta <<1`$), the cooling time is about the same as the buoyancy oscillation period ($`t_{cl,c}t_b`$) and the cooling distance is about the size of scale height of the halo ($`l_{cool}R_g`$). This means the global stability needs to be considered (Balbus & Soker 1989) for such perturbations.

### (c) Effects of Dynamical Instabilities

When a bubble-shaped perturbation condenses due to radiative cooling, it falls along the direction of gravity and so it is subject to the Rayleigh-Taylor and Kevin-Helmholtz instabilities. Nonlinear evolution of such perturbations in a stratified halo of X-ray clusters was followed via numerical simulations (Hattori & Habe 1990; Yoshida, Hattori, & Habe 1991; Reale et al. 1991). The bubble is disrupted by these instabilities on a buoyancy oscillation period if it does not cool and condense significantly. Thus the cooling timescale must be shorter than the buoyancy oscillation period in order to prevent the disruption of the bubble. That depends on three parameters in order of importance: the ratio of the local buoyancy oscillation period to the local cooling time, the initial density contrast of the perturbation, and the ratio of the perturbation radius to the local scale-height (Reale et al. 1991). We note here the first and second parameters are not completely independent, since the cooling time of the perturbed gas is dependent on the initial density contrast.

Following Reale et al. (1991), we define for our protogalactic halo model,

$$\xi =\frac{t_b}{t_{cl,c}}1.25(1+\delta )^2R_{10}^{(1q)},$$
(9)
using Eq. (6). According to their 2D numerical simulations, for $`\xi 1`$ the perturbations lead to fragmentation due to the dynamical instabilities, but the fragments can be thermally unstable. For $`\xi >5`$, the perturbation cools and condenses so rapidly that the dynamical instabilities do not have time to disrupt the bubble. These numerical simulations find that thermal instabilities cannot be an effective mechanism for mass depletion from the X-ray cluster cooling flows, since typically $`\xi 0.01`$ there. According to Eq. (9), on the other hand, $`\xi 5`$ for $`\delta 1`$ in a protogalactic halo, so the galactic halo provides more favorable conditions for thermal instabilities.

Thus the initial density contrast of the perturbations is the next important physical parameter. Here, we will take the generic view that the protogalaxy was formed from the hierarchical clustering of smaller substructures and contains some clumpiness and turbulences. We suggest only nonlinear (i.e. $`\delta \mathit{}\mathit{>}1`$) inhomogeneities in the protogalactic halo can condense due to thermal instabilities without being disrupted by the buoyancy-driven dynamical instabilities.

### (d) Numerical Method

The gas dynamical equations including self-gravity and cooling are:

$$\frac{d\rho }{dt}+\rho u=0,$$
(10)
$$\frac{du}{dt}=\frac{1}{\rho }P+g,$$
(11)
$$\frac{de}{dt}=\frac{1}{\rho }(Pu)+gu\frac{\mathrm{\Lambda }}{\rho }_,$$
(12)
where $`e=(\gamma 1)p/\rho +(1/2)u^2`$ is the total energy of the gas per unit mass and $`g`$ is the gravitational acceleration. We use the Piecewise Parabolic Method (PPM) (Colella & Woodward 1984) in one-dimensional spherical geometry and treat gravity and cooling as source terms. The gravitational acceleration is softened near the center according to

$$g=\frac{GM(<r)}{[\mathrm{max}(r,s)]^2}_,$$
(13)
where $`G`$ is the gravitational constant, $`M(<r)`$ is the mass within the radius $`r`$, and $`s0.04R_c`$ is the softening length. In most previous numerical simulations of the thermal instability (e.g., , the references in §1), self-gravity was ignored because the gravitational time scale was longer than the cooling time scale and/or the emphasis was on the effects of the thermal instability alone. In our models, the gravitational time scale is longer than the cooling time scale initially, but not at later times. Thus, we must follow the transition of the inward flow from being driven by background pressure at early times to self-gravity in the last stages of collapse.

The simulations start at $`t=0`$ with the clouds at rest in pressure equilibrium and cease when the background gas has cooled to $`10^4\mathrm{K}`$ (see below). The standard mirror condition for the reflecting boundary at the center is used, while the flow is assumed to be continuous across the outer boundary.

### (e) Radiative Cooling Rate

Our radiative cooling rate $`f(T)`$ applies to a primordial gas in ionization equilibrium (cf. Shapiro & Kang 1987). (The cooling rate per unit volume $`\mathrm{\Lambda }=n^2f(T)`$ where $`f(T)`$ is a function of temperature only.) The success of FR85 model is hinged on imprinting of the gravitationally unstable mass scale of PGCCs at $`10^4\mathrm{K}`$ which is similar to the characteristic mass scale of GCs. This requires that formation of $`\mathrm{H}_2`$molecules should be delayed longer than a free fall time of PGCCs. Kang et al. (1990) suggested that the globular cluster might have formed only during the phase when the protogalaxies are bright in UV background radiation. Hence we assume a model in which UV background radiation prevents formation of $`\mathrm{H}_2`$, so that the cooling rate is zero below $`10^4\mathrm{K}`$ and the minimum temperature is set at $`10^4\mathrm{K}`$.

The criterion for the thermal instability is that the over-density continues to cool faster than its surroundings, requiring $`d\mathrm{ln}f(T)/d\mathrm{ln}T<1`$ to destabilize isochoric condensations and $`d\mathrm{ln}f(T)/d\mathrm{ln}T<2`$ to destabilize isobaric condensations (Balbus 1986). Fig. 1 shows $`f(T)`$ and includes dotted lines where $`d\mathrm{ln}f(T)/d\mathrm{ln}T=1`$ and $`d\mathrm{ln}f(T)/d\mathrm{ln}T=2`$.

We ignore the poorly understood process of thermal conduction. With the classical Spitzer (1979) conductivity, thermal conduction may dominate over the cooling. If we assume that the heat flux is saturated when the particle mean free path is comparable to the length scale of the temperature variation, the conductivity is one or two orders of magnitude smaller than the Spitzer value (see Gray & Kilkenny 1980). Brinkmann, Massaglia & Müller (1990) showed that flow quantities differ $`<`$3% when they added the reduced thermal conduction to their simulations of hot gas in clusters of galaxies.

## III Simulation Results

### (a) Classification by Cloud Size

We calculated spherical collapses for a wide range of halo densities ($`0.01\mathrm{cm}^3n_h1\mathrm{cm}^3`$) and cloud sizes ($`0.04<R_c/l_{cool}<2.5`$). Initial parameters of the model clouds are summarized in Table 1, where $`M_c`$ is the mass of the cloud contained within the initial perturbation.

Radiative cooling of the cloud leads to three different types of flow motions depending on its size relative to the cooling distance, as recognized in previous studies (FR85, DBS88). According to our simulations, they can be summarized as follows: 1) $`(R_c/l_{cool})<0.2`$, small isobaric compression regime, 2) $`0.2<(R_c/l_{cool})<11.5`$, intermediate supersonic compression regime, and 3) $`(R_c/l_{cool})>11.5`$, large isochoric cooling regime. While in small clouds the collapse remains isobaric and quasi-static, in large clouds the central gas cools almost isochorically and the thermal instability fails. Examples of these limiting cases are plotted in Fig. 2. Spherical collapses of clouds with $`R_c/l_{cool}0.04`$ and $`R_c/l_{cool}2`$ in the background halo with $`n_h=0.1\mathrm{cm}^3`$ and $`T_h=1.7\times 10^6\mathrm{K}`$ are shown. Although the transitions between types are obviously continuous and gradual, we classified our models into “isobaric”, “supersonic”, and “isochoric” types according to the characteristics of flow motions in the simulations. The last column of Table 1 shows the types.

The intermediate size clouds belong to the most interesting regime. In order to see the evolution of an intermediate size cloud and the geometrical effects, spherical collapse as well as plane-parallel collapse were calculated for one model: $`T_h=1.7\times 10^6\mathrm{K}`$, $`n_h=0.1\mathrm{cm}^3`$, $`l_{cool}10^3\mathrm{pc}`$, and the scale height of the plane slab, $`H`$, and the radius of the cloud, $`R_c`$, are $`H=R_c=410\mathrm{pc}0.4l_{cool}`$. Fig. 3 shows the evolution at $`t/t_{cl,c}=1.3,1.5,1.7,1.9,`$ and 2.2, where $`t_{cl,c}=5.0\times 10^6\mathrm{yrs}`$. In both geometries, the gas cools isobarically until $`T2.5\times 10^5\mathrm{K}`$, below which the cooling becomes exponentially rapid and out of pressure equilibrium. As the pressure drops in the middle, a compression wave moves in. The infalling flow accelerates up to $`0.5c_h`$, moves to the center and then bounces back as an accretion shock. In the plane-parallel case, the accretion shock smoothes out the pressure gradient leading to near equilibrium with the background pressure. In the spherical case, on the other hand, the pressure of the shocked gas increases beyond the initial value owing to the central focusing. Thus the density enhancement factor can be as high as $`10^6`$ ($`n_h10^4\mathrm{cm}^3`$), which is $`10^{3.5}`$ higher than the isobaric compression ratio. Although there is a large positive pressure gradient between the cloud core and the background gas, there exists a halo of supersonic accretion flow instead of the out-flow, because the high density and high pressure core is gravitationally bound.

Fig. 4 shows the density profiles of the clouds of intermediate size ($`R_c0.42l_{cool}`$) with different halo densities, SM25 ($`n_h=0.01\mathrm{cm}^3`$) and SM4 ($`n_h=0.1\mathrm{cm}^3`$). This suggests that the density profiles of PGCCs can be represented by a core with a halo of $`\rho r^2`$. Since the gravitational time scale decreases with the density as $`n^{1/2}`$ while the cooling time scale decreases as $`n^1`$, the relative importance of self-gravity increases in lower density models. For example, the model SM25 stops expanding after $`t_6215`$, while the cloud in SM4 expand outward slowly. But the expansion speed is about $`0.7\mathrm{km}\mathrm{s}^1`$, much less than the sound speed of $`10^4\mathrm{K}`$ gas.

Vietri & Pesce (1995) proposed that, instead of the quasi-static evolution into the two-phase medium, over-dense regions create an implosion shock that leads to $`\mathrm{H}_2`$formation and cooling of the post-shock gas to $`10^2\mathrm{K}`$ (assuming no external UV radiation). However, we found that the implosion shock develops only if the initial density perturbation is strongly nonlinear ($`\delta >>1`$) or has a discontinuous top-hat distribution, both of which are unlikely in real density perturbations. Even in these cases, the implosion shock is nearly isothermal due to the rapid cooling of the post-shock gas. When the implosion shock rebounds at the center, it becomes an accretion shock. The post accretion shock gas still cools rapidly and condenses, so Vietri & Pesce’s “small dispersing clouds” are an artifact of their methods.

### (b) Physical Properties of PGCCs

According to our numerical calculations, only the intermediate size clouds in the supersonic compression regime, that is, $`0.4\mathrm{kpc}<R_c(n_h/0.05\mathrm{cm}^3)<2\mathrm{kpc}`$ and $`10^{5.6}\mathrm{M}_{}<M_c(n_h/0.05\mathrm{cm}^3)^2<10^8\mathrm{M}_{}`$, become PGCCs via thermal instability. Smaller clouds in the isobaric limit are not gravitationally bound (see in next subsection), while bigger clouds in the isochoric limit are too large to outrace the evolution of the background gas to become distinct structures. The relevant cloud sizes increase for lower halo gas densities (so for greater $`R_g`$), since the cooling distance is inversely proportional to the gas density. But the formation time becomes longer at lower densities and the length scale increases to a scale on which other dynamical effects may suppress the density growth.

In the transition from small to intermediate clouds, we see larger compressions owing to the supersonic infall with the compression factor increasing with the mass of the cloud (see Figs. 1 and 2). Thus the final radii of PGCCs, $`R_{pgcc}`$, lie in the narrow range: $`R_{pgcc}(n_h/0.05\mathrm{cm}^3)3060\mathrm{pc}`$. This owes to the dynamics of the supersonic infall whereby the shock in a larger cloud has a greater opportunity to accelerate and contributes more to its overall compression. This is extremely interesting, as the half mass radii of globular clusters varies by only about a factor of $`3`$ over two decades in mass (c.f. Fall & Rees 1977).

PGCCs form within $`2t_{cl,c}(2\times 10^7\mathrm{yrs})`$ $`(0.05\mathrm{cm}^3/n_h)`$, while the background halo gas cools in a time scale of $`2t_{cl,h}(8\times 10^7\mathrm{yrs})(0.05\mathrm{cm}^3/n_h)`$. So PGCCs would lose the background pressure support in the same time scale, unless the halo gas is continuously heated by stellar winds and supernova explosion from Pop III stars or accretion shocks or shocks generated by merging sub-clumps. Thus only self-binding intermediate clouds keep their identity, while clouds disperse.

### (c) Fragmentation of PGCCs

Although our numerical simulations are limited to one dimension and cannot follow the thermal history of PGGCs below $`10^4\mathrm{K}`$, we expect the dense core of PGCCs may fragment further owing to self-gravity. If the shape of PGCCs is more like a pancake, then the cloud would break up first by other dynamical instabilities into more or less spherical fragments and then only gravitationally unstable fragments undergo further collapses. Thus it is beyond our study to make the detail predictions on fragmentation processes. Here we attempt to make a rough estimate of the minimum mass scale of unstable PGCCs.

For small isobaric clouds with $`M_c<10^{5.6}\mathrm{M}_{}`$ $`(n_h/0.05\mathrm{cm}^3)^2`$, the minimum unstable mass can be the “critical mass”, $`M_{cr}`$ adopted by FR85, which was derived for an isothermal sphere confined by an external pressure (McCrea 1957),

$$M_{cr}=1.18\left(\frac{kT_c}{\mu m_H}\right)^2G^{3/2}p_h^{1/2}$$
(14)
For $`T_c=10^4\mathrm{K}`$ and $`p_h=1.17\times 10^{11}\mathrm{dyne}\mathrm{cm}^2`$, $`M_{cr}4.7\times 10^6\mathrm{M}_{}>M_{small}`$, so the small clouds are too small to become unstable gravitationally.

For intermediate size clouds, as shown in Figs. 3-4 (spherical cases), the density distribution of a typical PGCC can be approximated by a dense core with a halo of $`\rho r^2`$. The flow in the core is almost static while the outer region is infalling supersonically. Thus the core of PGCCs can be approximated by an almost uniform static sphere. This is certainly different from the singular isothermal distribution ($`\rho r^2`$ all the way to $`r=0`$), which is often assumed in the previous studies mentioned in the introduction. In this case, PGCCs are confined by self-gravity and the background pressure is negligible, so $`M_{cr}`$ given in Eq. (14) is not relevant. In a simple approach where the virial equilibrium is used (McCrea 1957) (e.g., $`2K+\mathrm{\Omega }0`$, where $`2K=3MkT/(\mu m_H)`$ is the thermal energy and $`\mathrm{\Omega }=(3/5)GM^2/R`$ is the gravitational potential energy for a uniform sphere), the minimum unstable mass is given by

$$M_{vir}=5.46\left(\frac{kT_c}{\mu m_HG}\right)^{3/2}\rho _{c}^{}{}_{}{}^{1/2},$$
(15)
where $`T_c`$ and $`\rho _c`$ are the temperature and density of condensed PGCCs. In fact this mass scale is almost the same as Jeans mass of the isothermal uniform sphere (Spitzer 1979) defined as

$$M_J=\rho \lambda _J^3=5.57\left(\frac{kT_c}{\mu m_HG}\right)^{3/2}\rho _{c}^{}{}_{}{}^{1/2}.$$
(16)
Thus we take this definition of Jeans mass as the minimum mass of unstable PGCCs as $`M_J(5.9\times 10^6\mathrm{M}_{})T_4^{3/2}n_2^{1/2}`$ where $`T_4=T_c/10^4\mathrm{K}`$ and $`n_2=n_c/(100\mathrm{cm}^3)`$. Thus supersonic compression lowers the Jeans mass by the square root of the additional compression factor. Owing to the spherical convergence, the highest compression factor we have seen in our numerical simulation is $`10^6`$, up to $`10^{3.5}`$ times higher than the isobaric compression. For PGCCs considered here, $`0.1<n_2<100`$, so the Jeans mass in the intermediate clouds is $`10^{5.5}\mathrm{M}_{}<M_J(n_h/0.05\mathrm{cm}^3)^{1/2}`$ $`<10^{7.2}\mathrm{M}_{}`$. Considering that the critical star formation efficiency may range between 0.1 and 0.5 in order for the resulting cluster to be self-bounded (Elmegreen et al. 1999), this Jean mass scale can be consistent with the observed GC mass distribution. If the evolution of more realistic clouds leads to irregular pancake-like shapes (§2.2, BMM90), however, the density enhancement for real clouds would be smaller than what is found in the spherical simulations, but still larger than the isobaric compression.

The use of shocks to lower the Jeans mass to a value that is closer to the characteristic mass of GCs was considered by Gunn (1980) in the context of chaotic protogalactic collapse and by Lake (1987) who invoked shocks created by cloud-cloud collisions. The cases considered by Gunn and Lake both have so much shear that they would be unlikely to lead to bound structures, but shocks generic to the collapse of the intermediate mass clouds don’t have this problem and still serve to reduce the Jeans mass to a value that is characteristic of GCs.

## IV CONCLUSION

We re-examined Fall & Rees (1985) model for the formation of PGCCs via thermal instabilities by numerical simulations of overdense clouds in a protogalactic halo. The key idea of this model is that characteristic mass scale of GCs, $`M_c10^6\mathrm{M}_{}`$, can be explained by imprinting of critical (Jeans) mass of gas cloud at $`10^4\mathrm{K}`$ which has cooled from a hot halo gas in pressure equilibrium via thermal instability. If there were significant cooling below $`10^4\mathrm{K}`$ due to $`\mathrm{H}_2`$and metals, the cloud will continue to cool and the imprinting is not possible (Kang et al 1990). Then the model fails. So we considered the cases where the radiative cooling is ineffective below $`10^4\mathrm{K}`$, because the formation of $`\mathrm{H}_2`$molecules is delayed due to UV radiation from a central AGN or diffuse radiation from halo gas.

The main conclusions are summarized as follows:

1. Unlike X-ray cluster cooling flows, a protogalactic halo provides much more favorable conditions for thermal instability to operate. So PGCCs can form from marginally nonlinear (i.e. $`\delta \mathit{>}\mathit{}1`$) density inhomogeneities in a protogalactic halo via thermal instability. Linear perturbations are likely to be disrupted by the buoyancy-driven instabilities (Balbus & Soker 1989; Reale et al. 1991).
2. The size of thermally unstable perturbations is determined by the cooling distance, $`l_{cool}`$, over which a sound wave can travel in a cooling time. According to our numerical simulations in one-dimensional spherical geometry, the clouds of $`0.2l_{cool}<R_c<l_{cool}`$ can condense to form PGCCs via thermal instability. Smaller clouds ($`R_c<0.2l_{cool}`$) cool only isobarically and become gravitationally unbound clouds, while larger cloud ($`R_c>l_{cool}`$) cool isochorically and so do not condense to form any distinct structures.
3. The thermally unstable clouds, which have mass $`10^{5.5}\mathrm{M}_{}<M_{PGCC}(n_h/0.05\mathrm{cm}^3)^2<10^8\mathrm{M}_{}`$, are compressed supersonically, and the density enhancement of $`10^310^6`$ results. This is much higher than what is expected from an isobaric compression. Considering the larger compression factor, we estimate the Jeans mass for PGCCs as $`10^{5.5}\mathrm{M}_{}<M_J(n_h/0.05\mathrm{cm}^3)^{1/2}<10^{7.2}\mathrm{M}_{}`$. The density distribution of simulated PGGCs can be approximated by an isothermal distribution with a constant core.

Although we focus on the formation of old halo GCs in the Milky Way in this work, we note globular clusters exist in some extremely different systems. The largest of the dwarf spheroidals, Fornax, Sagittarius, NGC 147, NGC 185 and NGC 205 all possess GCs. The specific frequency of clusters in these dwarf systems is as large as in bright ellipticals and larger than in spirals (Hodge 1988). If gas were in hydrostatic equilibrium in these parent galaxies, the temperature would be just $`10^410^5`$K. Thus a different model other than thermal instability model for the formation of PGCCs or the formation of GCs themselves is necessary for these systems.

Many models for GC formation have been suggested so far, but none of them seem to be complete or can explain all the observed properties of GCs. One of the aspects of FR85 model that has not been touched upon in the present work is the effect of the turbulent motions in the halo. Turbulent energy equivalent to the virial temperature should be injected at the largest scale and cascade down to smaller scales. According to 2D simulations by Vázquez-Semadeni, Gazol, & Scalo (2000) the turbulence cannot suppress the thermal instability, because the time scale for the thermal instability remains shorter than the turbulent crossing time scale ($`R_{\mathrm{cloud}}/v_{\mathrm{turb}}`$) at smaller scales than the virialization-scale. This is because the turbulent velocity decreases faster than the scale, so the turbulent crossing time increases at smaller scales. On the other hand, the turbulence can provide some heating to the halo gas. Since the halo gas cools radiatively in about $`10^8\mathrm{yrs}`$, additional heating is necessary to maintain the halo temperature at $`T_h`$ for much longer than the cooling time scale in order to explain the observed age spread of GCs, typically a few Gyr.

As discussed in §2, it is necessary to consider more realistic multi-dimensional simulations of thermally unstable clouds in a protogalactic halo environment, which is currently under study. The present work, based on one-dimensional considerations, however, provides useful guidance for the up-coming multi-dimensional works.

ACKNOWLEDGEMENTS The work by HK was supported in part by Pusan National University Research Grant, 1999. The work by GL was supported in part by NASA HPCC/ESS, ATP and LTSA at the University of Washington. The work by HK and DR was supported in part by grant 1999-2-113-001-5 from the interdisciplinary research program of the KOSEF. HK and DR acknowledge generous hospitality by the University of Washington, while this work had been carried out.

### REFERENCES

Ashman, K. M. & Zepf, S. E. 1992, ApJ, 384, 50.

Balbus, S. A. 1986, ApJL, 303, L79.

Balbus, S. A. & Soker, N. 1989, ApJ, 341, 611.

Brinkmann, W., Massaglia, S. & Müller, E. 1990, A&A, 237, 536 (BMM90).

Brown, J. H., Burkert, A. & Truran, J. W. 1991, ApJ, 376, 115.

Colella, P. & Woodward, P. R. 1984, J. Comp. Phys., 54, 174.

David, L. P., Bregman, J. N. & Seab, C. G. 1988, ApJ, 329, 66 (DBS88).

Dopita, M. A. & Smith, G. H. 1986, ApJ, 304, 283 .

Elmegreen, B. G., Efremov, Y. N., Pudritz, R. E., & Zinnecker, H., 1999, astro-ph/9903136

Fall, S. M. & Rees, M. J. 1977, MNRAS, 181, 37.

Fall, S. M. & Rees, M. J. 1985, ApJ, 298, 18.

Gray, D. R. & Kilkenny, J. D., 1980, Plasma Phys., 22, 81.

Gunn, J.E. 1980, in Globular Clusters, ed. D. Hanes & G. Madore, (Cambridge: Cambridge University Press), p. 301.

Hattori, M. & Habe, A. 1990, MNRAS, 242, 399.

Hodge, P. 1988, PASP, 100, 568.

Kang, H., Shapiro, P. R., Fall, S. M. & Rees, M. J. 1990, ApJ, 363, 488.

Kumai, Y., Baku, B. & Fujimoto, M. 1993, ApJ, 404, 144.

Lake, G. 1987, in IAU Symposium 127, Structure and Dynamics of Elliptical Galaxies, ed. T. de Zeeuw, (Dordrecht: Reidel), p. 331.

Malagoli, A., Rosner, R. & Bodo, G. 1987, ApJ, 319, 632.

Malagoli, A., Rosner, R. & Fryxell, B. 1990, MNRAS, 247, 367.

McCrea, W. H. 1957, MNRAS, 117, 562.

Reale, F., Rosner, R., Malagoli, A., Peres, G. & Serio, S. 1991, MNRAS, 251, 379.

Shapiro, P. R. & Kang, H. 1987, ApJ, 318, 32.

Spitzer, L. 1979, Physical Processes in the Interstellar Medium, (New York: Wiley-Interscience).

Vázquez-Semadeni, E., Gazol, A., & Scalo, J., 2000, astro-ph/0001027

Vietri, M. & Pesce, E. 1995, ApJ, 442, 618.

Yoshida, T., Hattori, M. & Habe, A. 1991, MNRAS, 248, 630.

Zinn, R. 1985, ApJ, 392, 424.


format: text

text:
# XMM-Newton observations of the cluster of galaxies Sérsic 159-03

## 1 Introduction

The visible mass in clusters of galaxies is dominated by hot diffuse gas. Due to the high temperature of the gas it is predominantly visible in the X-ray band. X-ray spectroscopy is the key tool to understand the physics of this gas and its role in the structure and evolution of the cluster. Since the gas is extended, spatially resolved spectroscopy is required. Until recently medium-resolution spectroscopy using Gas Scintillation Proportional Counter (GSPC) or Charge-Coupled Device (CCD) technology was limited to a spatial resolution of a few arcminutes (ASCA and BeppoSAX). Therefore only in a handfull of nearby clusters the core could be resolved spatially. While other instruments aboard the Einstein and Rosat satellites had better spatial resolving power, they either lacked bandwidth or spectral resolving power. High spectral resolution data of clusters have only been obtained for the bright core of the Virgo cluster with the Einstein FPCS detector (Canizares et al. canizares (1979)). With XMM-Newton it is now possible to combine high-resolution spectroscopy of moderately extended sources using the Reflection Grating Spectrometer (RGS) with high-sensitivity, medium (CCD-type) spectral resolution imaging with the European Photon Imaging Camera (EPIC) on spatial scales down to a few arcseconds.

Here we report the XMM-Newton observation of the rich cluster of galaxies Sérsic~159$-$03, discovered by Sérsic (sersic (1974)), also named Abell S 1101. The cluster shows within the central 2′ a cooling flow of 230 M/year (Allen & Fabian allen (1997)), centered on the dominant cD galaxy ESO 291$-$9. The only redshift measurement of the cluster ($`z`$=0.0564) is from this central galaxy (Maia et al. maia (1987)). Using $`H_0`$=50 km s-1 Mpc-1 and $`q_0=0.5`$ we then have a luminosity distance of 343 Mpc and an angular size distance of 307 Mpc. We adopted a galactic column density of $`1.79\times 10^{24}`$ m-2 (Dickey & Lockman dickey (1990), using NASA’s w3nH tool).

## 2 Observations

The observations were obtained on May 11, 2000. Data processing was done using the development version of the Science Analysis System (SAS) of XMM-Newton. Due to enhanced and variable background, in particular at the end of the observation we only used 32000 s for EPIC and 36000 s for RGS. Background subtraction for both EPIC and RGS was done using an exposure of the Lockman hole, with similar data selections and the same extraction regions as for the cluster, scaling according to the exposure time.

In the EPIC data we excluded the 7 strongest point sources both in the Sérsic 159$``$03 and Lockman hole field. The thin filter was used for both MOS cameras. In both RGS instruments, we ignored data below $`10`$Å due to problematic background subtraction.

## 3 EPIC analysis

We first fitted the spectrum using isothermal models. We determined the abundances of the elements with respect to iron in the 1–6′ range, where the spectrum has good statistics and the effects of a temperature drop in the center and outer regions appear to be small (see later). We determine abundances relative to solar, where we take the solar abundances of Anders & Grevesse (anders (1989)), and then express all abundances relative to the iron abundance, the best determined abundance with a relative error of only 7.5 % (Table 1).

The well-determined abundances (O, Mg, Si, S) are all consistent with solar ratios; the others have in general larger error bars but are not inconsistent with solar values. Therefore we will adopt solar ratios in this paper, varying only the global metallicity (iron abundance).

As a next step we fitted spectra in 12 logarithmically spaced annuli with outer radii ranging between 8–800″. We also created a set of spectra corrected for the projection on the sky (assuming spherical symmetry). The radial temperature and abundance profiles of both sets of spectra agreed within their error bars, but since the observed (projected) spectra are less noisy we focus our analysis upon them. However, for the radial density profile we used the deprojected spectra.

The radial density profile is shown in Fig. 1. The best-fit $`\beta `$-model $`n_\mathrm{H}=n_0[1+(r/a)^2]^{3\beta /2}`$yields a core radius $`a`$ of 0.53$`\pm `$0.06′ (47$`\pm `$5 kpc), a value for $`\beta `$ of 0.57$`\pm `$0.02 and a central hydrogen density $`n_0`$ of $`(2.40\pm 0.22)\times 10^4`$ m-3. These values are consistent within the error bars with $`\beta `$ and $`a`$ as derived from Rosat HRI observations (Neumann & Arnaud neumann (1999)).

In Fig. 2 we show the radial distribution of iron and in Fig. 3 the temperature distribution. Empirically, the iron abundance is modelled well by an exponential decay with a central abundance of 0.51$`\pm `$0.04 and a scale height of 2.6$`\pm `$0.4′, although this representation is not unique. Within the central core radius (0.5′) the abundance is more or less constant.

Sérsic 159$``$03 was reported to have a cooling flow, however Fig. 3 shows only a very modest temperature drop: from $``$2.6 keV between 1–6′ to about 2.3 keV in the center. The strong iron abundance gradient and hence the strong Fe-L line complex around 1 keV may have mimicked partly a cooling flow. Furthermore, a strong temperature decrease beyond 5′ is clearly visible.

We investigated the temperature drop in the center further by fitting the spectrum within 1′ by a hot isothermal model plus an isobaric cooling flow model (Johnstone et al. johnstone (1992)). The temperature of the hot component was frozen to the 1–2′ temperature of 2.69 keV. The best-fit parameters are a mass-deposition rate of 540$`\pm `$220 M/year, metallicity 0.40$`\pm `$0.02, and a low-temperature cut-off of 1.37$`\pm `$0.17 keV. The cooling flow component comprises about 2/3 of the 2–10 keV flux in this region. The $`\chi ^2`$ of the fit was 351 for 308 degrees of freedom (a single temperature fit has only a slightly worse $`\chi ^2`$ of 373). The low-temperature cut-off of the cooling flow model is very significant; we could not obtain a satisfactory fit for very small cut-off temperatures, because these models produce too much Fe-L emission below 0.9 keV.

We also studied the azimuthal distribution of the abundances and temperature by using hardness ratios in selected energy bands. While the temperature map (2$`4`$ keV / 0.3$`0.5`$ keV ratio; not shown) shows no obvious deviations from the almost isothermal behaviour in the cluster core, the abundance distribution is not spherical at all. We produced a map of the equivalent width of the Fe-L blend between 0.9–1.2 keV as compared to the underlying 0.7–1.3 keV continuum (Fig. 4).

For the temperature range of interest and accounting for the only minor variations in temperature based upon the hardness ratio map in the core, the variations in equivalent width as discussed below translate almost directly into iron abundance variations.

The map of the equivalent width shows several features. First, the largest iron abundance is found about 20″ to the north of the brightness peak (compare Fig. 5). Further, there is an extension out to 2′ to the NE of the center, coinciding more or less with the orientation of the central cD galaxy. Finally, there is an elongated structure in the EW direction extending several arcminutes to the east and west, with at least two times more iron than in other regions at the same distance from the core.

## 4 RGS analysis

Firstly, we fitted the spectra with the isothermal collisional ionization equilibrium model present in the SPEX package (kmn (1996)). As in the MOS fits we took abundances relative to solar (Anders & Grevesse anders (1989)). The free parameters are the normalization, the temperature and the abundances of oxygen, neon and iron. In RGS2, the O viii $`\mathrm{Ly}\alpha `$ line falls on CCD 4, which failed, and therefore the oxygen abundance is not well constrained. The abundances of all the other relevant elements (C, N, Mg, Si, S, Ar, Ca, Ni) have been coupled to the global metallicity (Fe) consistenly with the MOS result. To this model we applied the cosmolgical redshift and galactic absorption like we did for MOS. The fit was already satisfactory in both instruments, but residuals in the region between 13 and 17 Å suggested the presence of a cooler component. We added it with the same abundances as the hotter component and the improvement is $`\mathrm{\Delta }\chi ^210`$. In Fig. 6 we show the spectrum obtained with RGS1, in Table 2 we reproduce the fitting parameters. The emission from the cooler component is about 2% of the emission from the hotter component, respectively $`5.9\times 10^{35}`$ W and $`2.5\times 10^{37}`$ W.

The two instruments give slightly different results but consistent within $`1\sigma `$, RGS1 gives systematically higher abundances and temperature, this could be due to residuals problems with the calibration.

The RGS fit is consistent with a model composed by two components: the dominant one has $`T=2.6\pm 0.5`$ keV with iron abundance of $`0.55\pm 0.19`$ times solar, a Ne abundance compatible with the Fe one and an oxygen abundance of $`0.19\pm 0.07`$ times solar which is small compared with the others. The cooler component has a temperature of $`0.8\pm 0.2`$ keV and its contribution is about 2 % of the total X-ray luminosity.

## 5 Discussion

Our results show that the central part of Sérsic 159$``$03 shows very little temperature structure: the average temperature in the innermost part is only 10 % lower than the temperature of the hot gas outside the core. From imaging data with little spectral information a large mass deposition rate has been deduced in the past (e.g. 288 M/year based upon Einstein IPC observations (White et al. white (1997)); 230 M/year with the Rosat PSPC (Allen & Fabian allen (1997)). These values will be partly biased by the strong abundance gradient in the core (Fig. 2). Our results with both EPIC and RGS are consistent with such a large mass deposition rate, but only if the emission measure distribution has a low-temperature cut-off around 1.4 keV. The RGS spectra give upper limits of 1 % of the emission measure of the hot gas at any temperature below 1.5 keV. This upper limit translates into a 0.1–10 keV luminosity of any cool gas that is less than 20 % of the predicted luminosity of a 230 M/year cooling flow without a low-temperature cut-off. The absence of cool gas is similar to what has been found in hotter clusters (A 1835: Peterson et al. peterson (2001); A 1795: Tamura et al. tamura (2001)). We refer to those papers for a further discussion of this effect.

Another interesting feature found from our high-resolution RGS spectrum is the small O/Fe ratio in the core of the cluster: 0.34$`\pm `$0.17. This is significantly smaller than the O/Fe ratio of 1.0$`\pm `$0.3 derived from EPIC spectral fitting in the outer parts of the cluster. The oxygen abundance is mainly determined through the strong O VIII Ly$`\alpha `$ line. This line is subject to resonance scattering. We estimate the optical depth of this line to be $``$0.5, which yields a relative line flux depression of no more than 15 % in the core. We conclude that there is a real decrease in the O/Fe ratio towards the center of the cluster. However the absolute oxygen abundance does not differ significantly between the core (0.19$`\pm `$0.07, as derived from the RGS) and the outer parts (0.19$`\pm `$0.06, as derived from EPIC), it is merely the iron abundance that increases towards the center. Since oxygen is almost totally produced by type II supernovae and iron mostly by type Ia supernovae, we conclude that the core of the cluster has been relatively rich in type Ia supernovae, while the type II supernovae are more or less uniformly distributed over the cluster. The oxygen from type II supernovae may have been ejected in protogalactic winds (Larson & Dinerstein larson (1975)). The dominance of type Ia supernovae in the center may be due to the dominance of elliptical galaxies with an old stellar population in the center of the cluster. Due to ram pressure stripping these galaxies can loose their metals to the intracluster medium (Gunn & Gott gunn (1972)).

The structure seen in the Fe-L equivalent width map (Fig. 4) in the central cD galaxy (within 30″ from the core) indicates that the iron is distributed inhomogeneously in the core. A possible explanation might be a recent merger of a gas-rich galaxy in the NE part of the cD galaxy. Iron is not expected to diffuse over large distances during the life-time of a cluster. Therefore the elongated iron-rich structure seen in the EW direction extending out to at least 3′ (0.3 Mpc) towards the east probably reflects the initial star formation in the early evolution phase of the cluster. The additional amount of iron in this structure is consistent with the total amount of interstellar iron from a solar-composition galaxy with $`10^{10}`$ M of gas. Rosat PSPC images show indications for an extended, elongated structure out to at least half a degree to the east of the core, with a width of a few arcminutes. The characteristic density in this structure is of the order of 100–200 m-3, similar to the average cluster density in our outermost annulus (Fig. 1). Perhaps these structures are associated with a supercluster; unfortunately the region around Sérsic 159$``$03 is poorly studied, so we cannot confirm this.

A possible cause of the strong temperature drop near the outer part of the cluster might also be associated with the transition from cluster to supercluster. Temperature drops by a factor of 2 or more and an associated metallicity gradient, on slightly larger spatial scales than in Sérsic 159$``$03 have been found e.g. in the A 3562/Shapley supercluster (Kull & Böhringer kull (1999); Ettori et al. ettori (2000)). Numerical models of cluster mergers (Ricker ricker (1998)) show that after the merger there can be large-scale temperature gradients in the outer parts of a cluster. This is due to a post-merger accretion shock caused by gas falling back after the merging of the cores. The temperature can drop by an order of magnitude at 10 core radii (about 5′ in Sérsic 159$``$03), qualitatively similar to what we find.

###### Acknowledgements.

This work is based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and the USA (NASA). The Laboratory for Space Research Utrecht is supported financially by NWO, the Netherlands Organization for Scientific Research.