text:
# Complete Reconstruction of the Neutralino System

## I Introduction

In the minimal supersymmetric extension of the Standard Model (MSSM), the spin-1/2 partners of the neutral gauge bosons, $`\stackrel{~}{B}`$ and $`\stackrel{~}{W}^3`$, and of the neutral Higgs bosons, $`\stackrel{~}{H}_1^0`$ and $`\stackrel{~}{H}_2^0`$, mix to form neutralino mass eigenstates $`\chi _i^0`$ ($`i`$=1,2,3,4). The neutralino mass matrix in the $`(\stackrel{~}{B},\stackrel{~}{W}^3,\stackrel{~}{H}_1^0,\stackrel{~}{H}_2^0)`$ basis

$`=\left(\begin{array}{cccc}M_1& 0& m_Zc_\beta s_W& m_Zs_\beta s_W\\ 0& M_2& m_Zc_\beta c_W& m_Zs_\beta c_W\\ m_Zc_\beta s_W& m_Zc_\beta c_W& 0& \mu \\ m_Zs_\beta s_W& m_Zs_\beta c_W& \mu & 0\end{array}\right)`$ (5)
is built up by the fundamental supersymmetry parameters: the U(1) and SU(2) gaugino masses $`M_1`$ and $`M_2`$, the higgsino mass parameter $`\mu `$, and the ratio $`\mathrm{tan}\beta =v_2/v_1`$ of the vacuum expectation values of the two neutral Higgs fields which break the electroweak symmetry. Here, $`s_\beta =\mathrm{sin}\beta `$, $`c_\beta =\mathrm{cos}\beta `$ and $`s_W,c_W`$ are the sine and cosine of the electroweak mixing angle $`\theta _W`$. In CP–noninvariant theories, the mass parameters are complex. Without loss of generality $`M_2`$ can be taken real and positive so that the two remaining non–trivial phases may be attributed to $`M_1`$ and $`\mu `$. Neutralinos are produced in $`e^+e^{}`$ collisions, either in diagonal or in mixed pairs: $`e^+e^{}\stackrel{~}{\chi }_i^0\stackrel{~}{\chi }_j^0`$ $`(i,j=1,2,3,4)`$. If the collider energy is sufficient to produce the four neutralino states in pairs, the underlying fundamental SUSY parameters $`\{|M_1|,M_2,|\mu |,\mathrm{\Phi }_1,\mathrm{\Phi }_\mu ;\mathrm{tan}\beta \}`$ can be extracted from the masses $`m_{\stackrel{~}{\chi }_i^0}`$ ($`i`$=1,2,3,4) and couplings. Partial information from the lowest $`m_{\stackrel{~}{\chi }_i^0}`$ ($`i`$=1,2) neutralino states is sufficient to extract $`\{|M_1|,\mathrm{\Phi }_1\}`$ if the other parameters have been pre–determined in the chargino sector kneur ; CKMZ .

## II Mixing formalism

In the MSSM, the four neutralinos $`\stackrel{~}{\chi }_{1,2,3,4}^0`$ are mixtures of the neutral U(1) and SU(2) gauginos and the SU(2) higgsinos. Since the matrix $``$ is symmetric, one unitary matrix $`N`$ is sufficient to rotate the gauge eigenstate basis $`(\stackrel{~}{B}^0,\stackrel{~}{W}^3,\stackrel{~}{H}_1^0,\stackrel{~}{H}_2^0)`$ to the mass eigenstate basis of the Majorana fields $`\stackrel{~}{\chi }_i^0`$: $`_{diag}=N^{}N^{}`$ with $`(\stackrel{~}{\chi }_1^0,\stackrel{~}{\chi }_2^0,\stackrel{~}{\chi }_3^0,\stackrel{~}{\chi }_4^0)=N(\stackrel{~}{B},\stackrel{~}{W}^3,\stackrel{~}{H}_1^0,\stackrel{~}{H}_2^0)`$. The mass eigenvalues $`m_i`$ in $`_{diag}`$ can be chosen positive by a suitable definition of the unitary matrix $`N`$, which can be parametrized by 6 angles and 10 phases. It is convenient to factorize the matrix $`N`$ into a diagonal Majorana $`𝖬`$ and a Dirac–type $`𝖣`$ component in the following way:

$$N=𝖬𝖣,\mathrm{with}𝖬=\mathrm{𝖽𝗂𝖺𝗀}\{\mathrm{e}^{\mathrm{i}\alpha _1},\mathrm{e}^{\mathrm{i}\alpha _2},\mathrm{e}^{\mathrm{i}\alpha _3},\mathrm{e}^{\mathrm{i}\alpha _4}\}$$
(6)
One overall Majorana phase is unphysical and may be chosen to vanish and the matrix $`𝖣`$ can be written as a sequence of 6 two-dimensional rotations. Complete analytical results are given in CKMZ .

## III The neutralino quadrangles

The unitarity constraints on the elements of the mixing matrix $`N`$ can be formulated by means of unitarity quadrangles which are built up by the links $`N_{ik}N_{jk}^{}`$ connecting two rows $`i`$ and $`j`$,

$$M_{ij}=N_{i1}N_{j1}^{}+N_{i2}N_{j2}^{}+N_{i3}N_{j3}^{}+N_{i4}N_{j4}^{}=0\mathrm{for}ij$$
(7)
and by the links $`N_{ki}N_{kj}^{}`$ connecting two columns $`i`$ and $`j`$

$$D_{ij}=N_{1i}N_{1j}^{}+N_{2i}N_{2j}^{}+N_{3i}N_{3j}^{}+N_{4i}N_{4j}^{}=0\mathrm{for}ij$$
(8)
of the mixing matrix. The areas of the six quadrangles $`M_{ij}`$ and $`D_{ij}`$ are given by

$$\mathrm{area}[M_{ij}]=\frac{1}{4}(|J_{ij}^{12}|+|J_{ij}^{23}|+|J_{ij}^{34}|+|J_{ij}^{41}|),\mathrm{area}[D_{ij}]=\frac{1}{4}(|J_{12}^{ij}|+|J_{23}^{ij}|+|J_{34}^{ij}|+|J_{41}^{ij}|),$$
(9)
where $`J_{ij}^{kl}`$ are the Jarlskog–type CP–odd “plaquettes”: $`J_{ij}^{kl}=\mathrm{}\mathrm{m}N_{ik}N_{jl}N_{jk}^{}N_{il}^{}`$. The matrix $`N`$ is CP violating, if either any one of the plaquettes is non–zero, or, if the plaquettes all vanish, at least one of the links is non–parallel to the real or to the imaginary axis. The phases of the neutralino fields are fixed (modulo a common phase) and the orientation of the neutralino quadrangles $`M_{ij}`$ and $`D_{ij}`$ in the complex plane is physically non–trivial. Since the neutralino mass matrix involves only two invariant phases $`\mathrm{\Phi }_1`$ and $`\mathrm{\Phi }_\mu `$, all the physical phases of $`N`$ are fully determined by these two phases in the mass matrix. Therefore the experimental reconstruction of the unitarity quadrangles overconstrains the neutralino system and numerous consistency relations can be exploited to scrutinize the validity of the underlying theory CKMZ .

## IV Threshold behavior of neutralino production

In CP–invariant theories the S-wave excitation giving rise to a steep rise $`\lambda ^{1/2}`$ of the cross section for the nondiagonal pair near threshold signals opposite CP parities of produced neutralinos. Not all nondiagonal pairs of neutralinos can be produced in the S–wave; if $`\{ij\}`$ and $`\{ik\}`$ have opposite CP parities, the pair $`\{jk\}`$ has the same and it will be excited in the P-wave characterized by the slow rise $`\lambda ^{3/2}`$ of the cross section. It is important to realize that CP-violation may change the threshold behaviour, cf. Fig.1a, and in particular may allow S-wave excitations in any nondiagonal pair, e.g. finding the $`\{ij\}`$, $`\{ik\}`$ and $`\{jk\}`$ pairs to be excited in the S–wave would uniquely signal CP violation.

## V Extracting the fundamental parameters

The measurements of the chargino–pair production processes $`e^+e^{}\stackrel{~}{\chi }_i^+\stackrel{~}{\chi }_j^{}`$ ($`i,j`$=1,2) carried out with polarized beams can be used for a complete determination of the basic SUSY parameters $`\{M_2,|\mu |,\mathrm{\Phi }_\mu ;\mathrm{tan}\beta \}`$ in the chargino sector with high precision snow\_c .

Each of the four characteristic equations for the neutralino mass squared can be cast into the form CKMZ

$`(\mathrm{}\mathrm{e}M_1)^2+(\mathrm{}\mathrm{m}M_1)^2+u_i\mathrm{}\mathrm{e}M_1+v_i\mathrm{}\mathrm{m}M_1=w_i(i=1,2,3,4)`$ (10)
The coefficients $`u_i`$, $`v_i`$ and $`w_i`$ are functions of the parameters $`\mathrm{tan}\beta `$, $`M_2`$, $`|\mu |`$, $`\mathrm{\Phi }_\mu `$ pre–determined in the chargino sector, and the mass $`m_{\stackrel{~}{\chi }_i^0}^2`$; the coefficient $`v_i`$ is necessarily proportional to $`\mathrm{sin}\mathrm{\Phi }_\mu `$ because physical masses are CP–even. Each neutralino mass defines a circle in the $`\{\mathrm{}\mathrm{e}M_1,\mathrm{}\mathrm{m}M_1\}`$ plane. Thus the measurement of three neutralino masses leads to an unambiguous determination $`M_1`$, cf. Fig. 2a. With only two light neutralino masses, the two–fold ambiguity can be resolved by exploiting the measured cross section $`\sigma \{12\}`$, as shown in Fig. 2b. However, if the phase $`\mathrm{sin}\mathrm{\Phi }_\mu `$ vanishes, there remains a two–fold discrete sign ambiguity CKMZ .

If all four masses are experimentally accessible the complete reconstruction of the mass and mixing parameters is easy. The four–state mixing of neutralinos in the MSSM is reflected in the sum rules for the neutralino couplings CKMZ . Therefore evaluating these sum rules experimentally, it can be tested whether the four–neutralino system $`\{\stackrel{~}{\chi }_1^0,\stackrel{~}{\chi }_2^0,\stackrel{~}{\chi }_3^0,\stackrel{~}{\chi }_4^0\}`$ forms a closed system, or whether additional states at high mass scales mix in, signaling the existence of an extended gaugino system.

## VI The supersymmetric Yukawa couplings

A linear collider with polarized beams TDR  offers the possibility to verify very accurately the fundamental SUSY assumption that the Yukawa couplings, $`g_{_{\stackrel{~}{W}}}`$ and $`g_{_{\stackrel{~}{B}}}`$ are indentical to the SU(2) and U(1) gauge couplings $`g`$ and $`g^{}`$. Varying the left–handed and right–handed Yukawa couplings leads to a significant change in the corresponding left–handed and right–handed production cross sections. Combining the measurements of $`\sigma _R`$ and $`\sigma _L`$ for the process $`e^+e^{}\stackrel{~}{\chi }_1^0\stackrel{~}{\chi }_2^0`$, the Yukawa couplings $`g_{_{\stackrel{~}{W}}}`$ and $`g_{_{\stackrel{~}{B}}}`$ can be determined to quite a high precision as demonstrated in Fig. 1b. The $`1\sigma `$ statistical errors have been derived for an integrated luminosity of $`𝑑t=100`$ and $`500`$ fb<sup>-1</sup> and for partially polarized beams.

## VII Summary

The measurement of the processes $`e^+e^{}\stackrel{~}{\chi }_i^0\stackrel{~}{\chi }_j^0`$ ($`i,j`$=1,2,3,4), carried out with polarized beams and combined with the analysis of the chargino system $`e^+e^{}\stackrel{~}{\chi }_i^+\stackrel{~}{\chi }_j^{}`$ ($`i,j`$=1,2), can be used to perform a complete and precise analysis of the basic SUSY parameters in the gaugino/higgsino sector $`\{M_1,M_2,\mu ;\mathrm{tan}\beta \}`$. The closure of the neutralino system can be verified by exploiting the sum rules for production cross sections.

JK was supported in part by the KBN Grant 5 P03B 119 20 (2001-2002) and the European Commision 5-th framework contract HPRN-CT-2000-00149. GMP was partially supported by the DPF/Snowmass Travel Fellowship from the Division of Particles and Fields of the American Physical Society, and of the Snowmass 2001 Organizing Committee. The authors would like to thank S.Y. Choi and P.M. Zerwas for many lively and interesting discussions.


format: text

text:
# Gas cooling in simulations of the formation of the galaxy population

## 1 Introduction

Gas-dynamical and radiative processes couple with gravitational instability to play a crucial role in the formation of galaxies. In hierarchically clustering universes galaxies are formed via the dissipative cooling and condensation of gas within dark matter halos. The analytic understanding of galaxy formation within this framework is well developed (e.g. White 1996 and references therein). It is most practically embodied in the so-called ‘semi-analytic’ (SA) models. These implement simple analytic treatments of baryonic processes within a model for the nonlinear development of the dark matter distribution to give Monte Carlo realisations of the evolution of the galaxy population (White & Frenk 1991; Lacey & Silk 1991; Kauffmann, White & Guiderdoni 1993; Cole et al. 1994; Somerville & Primack 1998).

A full description of galaxy formation requires treatment of many complex physical processes. The huge dynamic range between the scales on which individual stars form and those of a ‘typical’ region of the Universe makes it impossible to simulate all relevant processes together. Some kind of ‘sub-grid’ model must therefore be adopted to describe the unresolved physics of star-formation and feedback within a dynamical treatment of larger scale evolution. Semi-analytic modelling adopts this philosophy not only on small scales, but also on scales where the evolution can, at some computational cost, be simulated directly. In SA models the merging history of dark matter halos is taken as the principal dynamical input. The following additional key physical processes are then included: (1) shock heating and virialisation of gas within halos, (2) radiative cooling and condensation of this gas, (3) star formation in the cold dense gas, (4) energy feedback from supernovae and stellar winds, (5) metal enrichment, (6) stellar evolution, and (7) galaxy merging. Early work by White & Rees (1978) and White & Frenk (1991) showed that regulation of star formation by feedback is necessary to avoid substantial over-production of stars at early times and in small objects, i.e. at least the first four of the above processes must be taken into account to obtain a realistic model. More recent semi-analytic models have incorporated all these processes in a consistent manner into Monte Carlo halo merger trees constructed using the (extended) Press-Schechter theory (Kauffmann et al. 1993; Cole et al. 1994, 2000; Kauffmann, Nusser & Steinmetz 1997; Somerville & Primack 1998; Benson et al. 1999). The output from such models has been tested against a wide range of observational data, showing that many properties of the observed galaxy population can be successfully reproduced. A recent development has been to replace the Monte Carlo merger trees with trees constructed directly from the output of high resolution $`N`$-body simulations (Kauffmann et al. 1999, hereafter KCDW). This removes uncertainties related to the Press-Schechter halo merger model and allows the evolution of the galaxy distribution to be followed explicitly and in detail. With the high resolution of Springel et al. (2001b, hereafter SWTK) it is also possible to follow the merging of the brighter galaxies explicitly, thus removing item (7) from the list of processes to be treated by semi-analytic recipes.

Although the models underlying these recipes are physically plausible and the resulting galaxy populations are in reasonable agreement with observation, the validity and the limitations of each model remain unclear. Among the seven ingredients listed above, the first two are well understood and can be incorporated in hydrodynamic simulations in a direct and transparent manner. Recently, Benson et al. (2001) studied the evolution of cooling gas associated with galaxy formation both in smoothed particle hydrodynamics (SPH) simulations and in a semi-analytic model using Monte Carlo merger trees. Comparing the output of the two models, they found reasonable agreement in properties such as the global fractions of gas in cold, hot and uncollapsed phases, and the amount of cold gas predicted in a dark halo as a function of its mass. Their results suggest that the assumptions of semi-analytic models give an appropriate description of the cooling of gas onto ‘protogalaxies’. Since the Monte Carlo merger histories they assign to each halo in their simulation are unrelated to its actual formation history, it is unclear, however, whether gas behaves as predicted by the SA model on an object-by-object basis.

In this paper we address this issue using a set of numerical simulations of rich galaxy cluster formation. We carry out two SPH simulations including radiative cooling, but excluding star-formation and the related feedback. The simulations differ only in the specific implementation of SPH employed. We compare the ‘galaxy’ populations of these simulations with those produced by implementing a stripped-down SA model in which star-formation and feedback are also suppressed. This comparison tests for differences arising from the very different treatments of hydrodynamics and cooling, and is unaffected by differences in the recipes which the SA and SPH approaches typically use for star formation and feedback. Because both techniques follow the dark matter distribution using the same $`N`$-body particle representation, we are able to compare the galaxy populations on an object-by-object basis.

The remainder of the paper is organised as follows. Sections 2 and 3 describe our SPH simulations and present their results. In Section 4, we explain our semi-analytic modelling highlighting features included specifically for the present study. We study and compare the effective resolution of our SA and SPH techniques in Section 5. In Section 6, we describe how we match the galaxies formed by the two techniques, and we show how their masses compare. We discuss the implications of these results in our concluding Section 7.

## 2 The $`N`$-body/SPH simulations

We work with a flat $`\mathrm{\Lambda }`$-dominated Cold Dark Matter universe, with matter density $`\mathrm{\Omega }_\mathrm{m}=0.3`$, cosmological constant $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ and expansion rate $`H_0=70`$km s-1Mpc-1. We carry out further simulations of the cluster studied by SWTK and Yoshida et al. (2000a,b) which is the second most massive cluster in the GIF simulation of KCDW. As in these earlier papers, we resimulate the cluster with a multi-mass technique in which only the region immediately surrounding the object of interest is represented at high resolution; the rest of the simulated volume is followed at much lower resolution and serves merely to provide the correct tidal gravitational field on the high resolution region. Within this region we split each particle into a dark matter particle of mass $`1.2\times 10^{10}h^1M_{}`$ and a gas particle of mass $`2.1\times 10^9h^1M_{}`$, and we set the gravitational softening length to be $`30h^1`$kpc. The simulations are carried out with the parallel tree $`N`$-body/SPH code GADGET (Springel, Yoshida & White 2001a); further details are presented in Yoshida (2001). During the simulations, we dump 50 snapshots of the particle data spaced logarithmically in the cosmic expansion parameter from redshift $`z=20`$ to $`z=0`$. We use these outputs for identifying and tracing galaxies, and also for constructing dark matter halo merger trees for use in the semi-analytic simulation. For the latter we throw out the gas particles and increase the mass of each dark matter particle to $`1.4\times 10^{10}h^1M_{}`$.

We implement gas cooling in the SPH part of the code in a similar manner to that in Katz, Weinberg & Hernquist (1996, hereafter KWH). As in KWH we compute the abundances of ionic species by assuming collisional equilibrium for a gas of primordial composition consisting of 76% hydrogen and 24% helium. We call these simulations ‘S0-SPH’. They include a time-dependent uniform UV background radiation, as in Davé et al. (1999), but its effect, especially on the destruction of small objects, is very small at the resolution we are using (Weinberg, Katz, & Hernquist 1997). Hence the inclusion of the UV heating is not important for the analysis presented in this paper. We fix a minimum SPH smoothing length which prevents the gas smoothing length from dropping below a quarter of the gravitational softening length. We set the number of neighbouring particles in the SPH smoothing kernel, $`N_{\mathrm{NGB}}`$, to be 32. This smoothing in the SPH simulation poses an effective resolution limit for gas cooling; gas cooling can take place efficiently only within objects that become dense enough and which contain enough gas particles. We will discuss this resolution effect in more detail in Sections 4 and 5.

An important feature of our simulations is that we employ two different implementations of SPH. One is an energy- and momentum-conserving scheme based on taking geometric means of the pairwise hydrodynamic forces between neighboring particles; it is very similar to that employed by KWH, Weinberg et al. (1997), Davé et al. (1999) and White, Hernquist & Springel (2001) among others. The other is a novel formulation of SPH recently suggested by Springel & Hernquist (2001), which manifestly conserves momentum, energy and entropy in non-shock regions, even when adaptive smoothing lengths are employed. In a series of numerical tests, Springel & Hernquist showed that several of the commonly employed formulations of SPH can give rise to an over-cooling phenomenon which results from spurious losses of specific entropy in poorly resolved cooling flows. This effect is particularly strong when the geometric symmetrisation of the hydrodynamical forces is employed. Our two simulations highlight the extent to which such differences in the formulation of SPH can affect the final result. In the following, we concentrate on results from the simulation carried out with the new ‘entropy’ scheme, but we stress the major points at which they differ from those found using the geometric averaging technique. Note that while we will often refer to the latter technique as the ‘conventional’ SPH technique, it suffers from the overcooling problem more severely than other ‘conventional’ techniques which have also been used in studies of galaxy formation (e.g. Benson et al. 2001). Pearce et al. (2001) found that ignoring cold dense particles in the SPH density estimate for hot particles successfully reduces the over-cooling of gas. Although we do not employ this phase-decoupling technique in our SPH simulations, the problem of density overestimates for hot particles in the vicinity of the condensed cold phase is much alleviated in the new entropy scheme (see Figure 8 in Springel & Hernquist, 2002), so the effect of the decoupling is expected to be small in our simulations.

We begin by presenting the basic properties of our simulated cluster. The virial radius of the main cluster at the present epoch is $`1.48h^1`$Mpc, and the corresponding virial mass is $`7.7\times 10^{14}h^1M_{}`$. We define the virial radius $`R_{\mathrm{vir}}`$ as the radius of the sphere centred on the most bound particle of the FOF group having overdensity 200 with respect to the critical density. The virial mass $`M_{\mathrm{vir}}`$ is then the enclosed mass within $`R_{\mathrm{vir}}`$. Figure 1 shows the projected gas mass distribution around the main cluster. Bright clumps seen in Figure 1 are groups of cold gas particles, which we will later define as ‘galaxies’. The diffuse component in the figure consists mainly of a shock-heated hot gas in and around the cluster.

Figure 2 shows the distribution of gas in the thermodynamic phase plane. The cold, dense gas particles appear clearly as a narrow tail in the plot, whereas the shock-heated gas and the adiabatically cooled, unshocked gas appear in the upper and in the lower left portions, respectively. The rectangular box approximately separates out the gas particles in the ‘galaxy’ phase. In practice we identify galaxies by linking the cold gas particles using the standard FOF technique with a small linking length. We describe our galaxy identification in the next section.

## 3 Cold and hot gas in dark matter halos

In this section we show the distribution and the properties of gas in the dark halos in the simulations. We locate the dark halos by running a FOF group finder with linking parameter $`b=0.164`$. We assigned virial radii and virial masses to the dark halos in the manner described above. We found 833 dark halos with virial mass larger than $`10^{11}h^1M_{}`$ at the present epoch.

### 3.1 Gas mass fraction

In Figure 3 we plot the total gas mass and cold gas fraction against the host halo mass in the ‘entropy’ version of the S0-SPH simulation. Gas quantities are calculated within a sphere of radius $`R_{\mathrm{vir}}`$. We define the cold gas fraction as the total mass of all the cold ($`T<100,000`$ K) gas divided by the total gas mass within the halo. Three sets of plots for the outputs at z=0, 1, and 2, are shown in Figure 3. In the top panels the solid line indicates the global mean baryon fraction in the simulation, $`f_\mathrm{b}`$=0.15. A tight correlation is seen, showing that the baryon fractions of the halos are very close to the universal baryon fraction. At mass scales larger than $`10^{12}h^1M_{}`$ (corresponding about 80 dark matter particles), the scatter is very small. At early times, there is a clear tendency for the gas fraction to lie slightly above the global fraction. This is a consequence of the rapid loss of pressure support behind the accretion shock when cooling is efficient.

The bottom panels in Figure 3 show that the cold gas fraction is a decreasing function of halo mass and shows relatively little scatter above halo masses of $`10^{12}h^1M_{}`$. In this regime the cold gas mass fraction at given halo mass becomes progressively lower at lower redshift. This behaviour is expected as a result of the decrease in cooling efficiency at lower redshift produced by the lower gas density within halos. It confirms that the new SPH implementation substantially reduces over-cooling of gas. In the simulation carried out with the conventional SPH implementation, we found the opposite behaviour; the cold gas mass fraction gets progressively higher at lower redshift. This (mis-)behaviour agrees with the experiments of Springel & Hernquist (2001) who showed that an SPH formulation with geometric symmetrisation can substantially overestimate the cooling rates in the centres of halos, leading to the accumulation of too much cold gas. The cold gas fractions are near 80% for most halos at $`z=0`$ in our simulation using this scheme, much larger than the values plotted in Figure 3. Even for the main cluster halo itself, the cold gas fraction is overestimated by a factor of 2.

Although the small scatter and the systematic trends in cold gas fraction seen above a halo mass of $`10^{12}h^1M_{}`$ in Figure 3 are in good agreement with simple semi-analytic models (e.g. White & Frenk 1991 and Section 5 below), the increase in scatter below $`10^{12}h^1M_{}`$ is unexpected. This is a consequence of resolution problems in the SPH simulation. When a halo contains relatively few gas particles, the smoothing inherent in the SPH technique reduces the estimated gas density to the point where the cooling rate is seriously underestimated. As is clearly visible in the figure, the effect is worse at $`z=0`$ than at higher redshifts because of the lower mean cooling rates at later times. We will return to the effective resolution limit implied by this problem in Section 5.

At redshift $`z=1`$ and particularly at $`z=2`$, another interesting effect is visible at small masses in Figure 3. As the typical cold gas fraction approaches 100% the cold gas mass within halos rises systematically above the nominal maximum given by the halo mass times the global baryon fraction. At the lowest masses this effect is about a factor of two. Although this may, in part, be due to SPH artifacts in the strong cooling and small particle number limit, much of the effect is undoubtedly real. When cooling is efficient shocks are unable to raise infalling gas to the virial temperature of a halo, no significant hot gas atmosphere is formed and the accretion shock moves in close to the halo centre. The total gas mass within the nominal halo radius then increases since pressure effects no longer impede infall at larger radii. We note in passing that this effect is not correctly included in the semi-analytic model we discuss below, which assumes that the maximum mass which can cool in any halo is $`f_{\mathrm{baryon}}M_{\mathrm{halo}}`$.

We close this section by mentioning a different cooling problem. As noted above, it has been known for more than twenty years that feedback of some kind is needed in hierarchical cosmogonies to prevent excessive cooling and star-formation at early times. Katz & White (1993; see also Evrard, Summers & Davis 1995; Suginohara & Ostriker 1998) demonstrated that such over-cooling occurs in hydrodynamic simulations of CDM universes provided these are able to resolve the relevant early phases of evolution. Subsequent work by Lewis et al. (2000) and Pearce et al. (2001) confirmed that the cold gas fraction indeed increases with increasing resolution. In our simulation, the cold gas fraction of the main cluster is about 23% at the present epoch (see the rightmost point in Figure 3). Although this value is substantially smaller than the 55% which we found in the simulation with the conventional SPH implementation, it still appears too high to be compatible with observational constraints on the fraction of baryons locked up in stars in rich clusters (Balogh et al. 2001). This again confirms that efficient feedback mechanisms are required to regulate cooling and star formation in realistic models. Nevertheless, since our aim is to compare the treatment of gas cooling in SPH and SA simulations, this ‘cooling catastrophe’ is not directly relevant to our analysis.

### 3.2 ‘Galaxies’ in the SPH simulations

In our two versions of the S0-SPH simulation, we define cold, dense clumps of gas as ‘galaxies’. We locate these objects by running a FOF group finder over the gas particles with a very small linking parameter, $`b=0.016`$. We note that the choice of linking parameter is not particularly critical; adopting a slightly smaller or larger value makes only a minor difference to the number of galaxies and their masses. We also checked that for our chosen value most cold, dense particles in the more massive halos are linked into galaxies, and multiple galaxies within a single halo are successfully separated into discrete objects. In the following we require a clump to contain ten or more gas particles to be considered as a galaxy, corresponding to a minimum mass of $`2.1\times 10^{10}h^1M_{}`$. We refer to such objects as SPH galaxies, in contrast to the SA galaxies defined using the semi-analytic model. The mass of an SPH galaxy is simply the mass of its grouped particles. We note that it may not be appropriate to regard $`N_{\mathrm{min}}`$ SPH particles as a distinct group, where $`N_{\mathrm{min}}N_{\mathrm{NGB}}`$, because of the SPH smoothing with $`N_{\mathrm{NGB}}`$ particles. Nevertheless we include such small groups in our SPH galaxy catalogue, because, at early epochs, many of them are progenitor galaxies with which we can keep track of the formation history in detail. We give an extensive discussion on the resolution of our SPH simulation in section 5.

## 4 Semi-analytic simulations of galaxy formation

In the galaxy formation models of KCDW and of SWTK, outputs of high resolution N-body simulations are used to define populations of dark halos and to follow their formation and evolution. Although the merging of halos is followed explicitly by these simulations (as is also the evolution of halo substructure in SWTK), all aspects of the evolution of the baryonic component are included using simple, physically motivated, recipes taken from earlier semi-analytic models. In this paper we follow the procedures of KCDW, simplified to exclude processes other than gas cooling and ‘galaxy’ merging. This eliminates the free parameters normally adjusted in SA models, since these are associated with star formation and feedback. The standard cooling model contains no free parameters (see below) and the merger model of KCDW implies merger rates in good agreement with those seen in simulations where merging can be followed explicitly (SWKT).

For each dark halo in each simulation output, we assign the following related halo properties: the virial radius $`R_{\mathrm{vir}}`$, the virial mass $`M_{\mathrm{vir}}`$, and the circular velocity $`V_{\mathrm{vir}}^2=GM_{\mathrm{vir}}/R_{\mathrm{vir}}`$. These halos are linked to their progenitors in the previous output to build up merger history trees within which the growth of galaxies is followed. Notice that since we use the dark matter distribution from the $`N`$-body/SPH simulation to set up these trees, there is an exact correspondence between the halos and their merger histories in our SA and SPH simulations. Any differences in the resulting galaxy populations must therefore be produced purely by differences in how gas dynamics, cooling and galaxy merging are handled. We now turn to the assumptions used to model gas cooling.

### 4.1 Gas cooling

Gas cooling is treated as in White & Frenk (1991). We assume that the gas in a dark halo is shock-heated to the virial temperature and is initially distributed like a singular isothermal sphere with density profile $`\rho _\mathrm{g}r^2`$. We define the local cooling time of the gas, $`t_{\mathrm{cool}}`$, as the ratio of the thermal energy density of the gas to the cooling rate per unit volume:

$$t_{\mathrm{cool}}(r)=\frac{3}{2}\frac{kT\rho _\mathrm{g}}{\overline{\mu }m_\mathrm{p}n_\mathrm{e}^2(r)\mathrm{\Lambda }(T,Z)},$$
(1)
where $`\overline{\mu }`$ is the mean molecular weight, $`m_\mathrm{p}`$ is the proton mass, $`n_\mathrm{e}`$ is the electron number density, $`T`$ is the gas temperature, which we approximate with the virial temperature of the halo $`T=35.9(V_{\mathrm{vir}}/\mathrm{kms}^1)^2`$ K, and $`\mathrm{\Lambda }(T,Z)`$ is the cooling rate. We employ the tabulated cooling functions of Sutherland & Dopita (1993).

We define a cooling radius $`r_{\mathrm{cool}}`$ such that the cooling time at this radius is equal to the time for which the halo has been able to cool quasi-statically. We approximate this time with the halo’s dynamical time $`R_{\mathrm{vir}}/V_{\mathrm{vir}}`$. The cooling rate for the gas within this cooling radius is then given by

$$\frac{\mathrm{d}M_{\mathrm{cool}}}{\mathrm{d}t}=4\pi \rho _\mathrm{g}(r_{\mathrm{cool}})r_{\mathrm{cool}}^2\frac{\mathrm{d}r_{\mathrm{cool}}}{\mathrm{d}t}.$$
(2)
In reality, cooling rates depend strongly on the metallicity of the gas, but for simplicity (and consistency) we use zero metallicity cooling functions in both our SA and our SPH models.

We note that this cooling model is extremely simplified and corresponds quite poorly to the detailed dynamics seen in hydrodynamic simulations of the evolution of galaxy-sized objects in a CDM universe. The good agreement which we find below with the gas accumulation rates in our best SPH simulation is thus quite surprising.

### 4.2 The stripped-down semi-analytic model

Our SPH simulations do not include star formation or feedback processes. In order to make a direct comparison, we remove these processes also from our semi-analytic simulation. In this stripped-down model, star formation and feedback are switched off, and no cooling cut-off is implemented in massive halos (c.f. KCDW). As in the SPH simulations, galaxies in our stripped-down SA simulation do not possess properties such as luminosity or morphological type. They are defined purely by their mass of cold gas, their positions and velocities, and their formation histories. There are only two populations of galaxies in this model. Each dark halo carries exactly one central galaxy which is identified with the most-bound particle of the halo. Many halos have also one or more satellite galaxies. Satellites are galaxies that were central galaxies in earlier outputs, but whose host halos merged at some time into the larger halo where they now reside. Satellites are identified with the particle to which they were attached when they were last a central galaxy. They are assumed to merge with the central galaxy of their new halo on a dynamical friction timescale (see KCDW). In this implementation, only the central galaxy accretes new gas that cools within the halo; hence the mass of a satellite galaxy stays constant.

## 5 The resolution of the SA and SPH simulations

In this section we discuss the effective resolution limits of our two simulation techniques. The principal goal of our paper is to compare the predictions of these techniques for the cooling of gas onto protogalaxies, but before this is possible we must establish the range of protogalaxy masses for which the predictions of each technique are independent of the resolution parameters – particle mass, force softening, SPH smoothing, timestep – employed in the corresponding simulation. With this in mind, we define the effective resolution of a simulation to be the minimum galaxy mass for which the simulated galaxy population has identical properties to that in a much higher resolution simulation of the same system carried out with the same technique. Notice that with this definition there is no guarantee that the galaxy populations produced by simulations using our two techniques will agree above the larger of their two resolution limits, although we find below that they do. Notice also that for the physical processes included here – cooling, but no star formation or feedback – there is reason to worry that our definition may fail; as the resolution of a simulation is increased the masses of all galaxies might increase systematically towards their maximal possible values, corresponding to the total baryon content of their individual halos. We find below that a ‘cooling catastrophe’ of this type does not, in fact, occur.

We begin by studying the resolution limit of our SA simulation technique. This is relatively straightforward because SWTK have already carried out a series of higher resolution simulations of the physical system studied here. For their highest resolution model, S4, the mass resolution is about 250 times higher and the spatial resolution about 50 times higher than in our S0-SA simulation. We have run our ‘stripped-down’ SA modelling on SWTK’s S4 simulation to derive a galaxy population which can be compared directly with those of this paper. Figure 4 compares the cumulative mass functions of the two SA simulations. Above a mass of about $`1.5\times 10^{10}h^1M_{}`$ the agreement is good with the masses of galaxies in the high resolution simulation lying slightly below those of galaxies in S0-SA. At lower mass the mass function of S0-SA shows a clear change in slope. This is easily understood. Our halo merger trees keep track of all halos with 10 or more particles. The maximum mass of a galaxy which could form in an unresolved halo is thus given by the global baryon fraction times the mass of a 9 particle halo. This is $`1.9\times 10^{10}h^1M_{}`$, slightly larger than the mass at which we see the change in slope. The convergence tests of SWTK indicate similar resolution limits for their full semi-analytic model.

It is important to note the very steep mass functions in Figure 4 which are a consequence of the absence of feedback. The total amount of gas cooled in the SA simulations is a strong function of resolution, rising from 16% in S0-SA to 25% in S4-SA. (We have not attempted to include the effects of an ionizing background in the SA modelling, which would reduce the strength of this effect.) It turns out, however, that the additional gas which cools off in small objects does not reduce the reservoir of diffuse gas enough to substantially affect the masses of the larger objects which form later. Similarly merging is not efficient enough for the increased supply of low-mass merger candidates to cause a substantial increase in the masses of larger systems. Thus, although there is no formal convergence as resolution is increased, the properties of the higher mass galaxy population vary very little.

It is more difficult to study the mass limit for our SPH model, because we cannot carry out a similar convergence test to that we decribed above for the SA models. Our decision to impose a 10 particle lower limit on our SPH ‘galaxies’ clearly imposes a lower bound of $`2\times 10^{10}h^1M_{}`$ on the resolution limit, but Figure 3 shows that the effective resolution limit is substantially higher. Below a halo mass of about $`10^{12}h^1M_{}`$ many halos contain substantially less cold gas than expected and some contain none at all. This is a consequence of the smoothing associated with SPH density estimates which leads to an underestimate of the cooling in systems where the number of gas particles is comparable to that required within the SPH smoothing kernel. We plot the $`z=0`$ cumulative galaxy mass function for the ‘entropy’ version of our S0-SPH simulation in Figure 4. Above $`1.5\times 10^{11}h^1M_{}`$ this function agrees well with that of our two SA simulations, except for a slight shift to higher mass. Below $`1.5\times 10^{11}h^1M_{}`$ the mass function has a clear change in slope. We take this mass scale as a conservative estimate for the mass resolution of our SPH simulation. This effective resolution limit in our S0-SPH simulation is about an order of magnitude larger than the corresponding limit for our S0-SA simulation. In light of this a robust comparison of SA- and SPH-galaxies can only be made at masses larger than about $`10^{11}h^1M_{}`$.

## 6 Object-by-object comparison

We now turn to a direct comparison of the properties of the individual objects which form in our SA and SPH simulations. The dark matter distributions of the two simulations are identical at all times, by construction, and as a result it is possible to identify individual galaxies in the two models almost without ambiguity. We begin by comparing the masses of the central galaxies. Every halo with at least 10 particles is automatically assigned a central galaxy in our semi-analytic model (see Section 4.2) whereas defining the central galaxy in the SPH simulation is less clear. For simplicity, we define the most massive SPH-galaxy within a halo to be its central galaxy, and we tag other galaxies as satellites. Note that this definition corresponds closely to that used in the SA model for halos with multiple galaxies (see KCDW).

In Figure 5 we compare the distributions of SA- and SPH-galaxies more massive than $`6\times 10^{10}h^1M_{}`$ in a cubic region, 15$`h^1`$Mpc on a side, centred on our cluster. The two distributions are very similar, and an almost one-to-one correspondence can be found. Most of these galaxies are the central galaxies of their dark matter halos except in the main cluster itself where there are many satellites. Thus the galaxy distribution largely reflects that of the centres of the massive halos. Most differences arise when galaxies lie on opposite sides of our mass threshold in the two samples.

In Figure 6 we plot SA mass against SPH mass for the central galaxies of all halos with 10 or more dark matter particles. Halos with no SPH-galaxy are shown on the left of this diagram at an arbitrary SPH mass. Above our SPH effective resolution limit ($`10^{11}h^1M_{}`$) the central galaxy masses predicted by our two techniques agree extremely well over more than two orders of magnitude in mass. Below this limit the masses of the SPH-galaxies are systematically low and many halos contain no SPH-galaxy even though the corresponding SA-galaxy mass is well above $`2\times 10^{10}h^1M_{}`$, the limit imposed by our requirement that an SPH-galaxy contain at least 10 particles. We emphasize that the agreement above $`10^{11}h^1M_{}`$ limit is remarkable considering the crudity of the cooling model used in the SA treatment and the fact that our two implementations of SPH lead to galaxy masses differing systematically by a factor of two. (In this section we only show results from the ‘entropy’ implementation.)

It is also interesting to see the evolution of the masses in two models. In Figure 7 we compare the masses of the three most massive galaxies at $`z=0`$ and their progenitors at earlier epochs. It is clearly seen that even the time evolution is remarkably similar in the two models. Small deviations from the perfect linear relation at output instants are due to time offset of satellite merging onto the central galaxy, because of the simplified prescription of satellite merging in our SA model.

A similar comparison for satellite galaxy masses is shown in the bottom panels in Figure 6. In this case matching SA- and SPH-galaxies is much more complicated than for central galaxies. Once an object becomes a satellite its position will not be the same in the two kinds of simulation, and indeed it may merge with the central galaxy in one of the simulations but not in the other. We use the histories of the satellite galaxies in the two simulations in order to carry out this matching. Given a satellite galaxy identified in an output of the semi-analytic simulation, say at redshift $`z_1`$, we track it back through earlier outputs until we find the time when it was last a central galaxy. Let us denote this time $`z_2`$. (Recall that, in our semi-analytic scheme every satellite galaxy was a central galaxy at some earlier time.) We then search in the SPH simulation at $`z_2`$ and identify the SA progenitor of our satellite with the most massive SPH galaxy in its halo. In some cases the halo may contain no SPH galaxy and we then consider our original SA satellite to be unmatched. Otherwise we locate the particles of the matched SPH galaxy in the SPH simulation at $`z_1`$ and consider the SPH galaxy which contains the majority of them to be the counterpart of the original SA satellite. This counterpart may be an SPH satellite, the SPH central galaxy, or (rarely) an ‘orphan’ galaxy outside any halo. Occasionally it can fall below the 10 particle threshold to be considered a bona fide SPH galaxy and we then again consider the original SA satellite to be unmatched.

We also carry out the reverse procedure. Starting from an SPH satellite at $`z_1`$, we trace back its progenitors (defined at each time to be the SPH galaxy which donates the largest number of particles to the galaxy at the subsequent output) until we find one which was the most massive SPH galaxy in its halo. This one we then identify with the SA central galaxy of the same halo. We take particular care to investigate SPH satellite galaxies which appear to have no progenitors. In each output, some new SPH galaxies appear for the first time as satellites; in the previous output, their member particles are already cold and are just condensing at the center of a halo. Due to discrete sampling of the output times, there is always a slight offset between the true formation epoch of SPH galaxies and the time when they are first identified. Their host halo can merge with a larger system during this time. To take this into account, we tag the central SA galaxy in any halo that hosts such a forming SPH galaxy as its counterpart, because it is clear that such pairs indeed correspond closely. In this way we can partially compensate for incompleteness in our SPH–SA galaxy matching procedure. The $`z_1`$ descendent of the matched $`z_2`$ SA galaxy is identified as the counterpart of the original SPH satellite. Again, this counterpart may be a satellite, a central galaxy, or (rarely) an ‘orphan’. We find 40 out of 84 SPH satellite galaxies at $`z=0`$ matched to SA satellites. Another 34 are matched to SA central galaxies. The remaining 10 are unmatched, mostly because they were formed under peculiar conditions in large halos without being at the center of the halos. It is likely that many are associated with subhalos within the large halo, but our simulation does not have adequate resolution to trace their formation in detail.

There are 8 SA satellites which are matched to SPH central galaxies, and 74 SA satellites more massive than $`2\times 10^{10}h^1M_{}`$ which are unmatched. A few of the latter may, in fact, correspond to unmatched SPH satellites; most correspond to systems in which the SPH gas did not cool. The small number of SA satellites matched to SPH centrals in comparison with the number of SPH satellites matched to SA centrals is a consequence of the lower effective resolution of the SPH simulation. We have checked that the merging rate, estimated as the fraction of all satellites which merge between two output times, is actually very similar in the two simulations.

Satellite-satellite matches are plotted as filled circles in Figure 6, whereas SA satellite – SPH central and SA central – SPH satellite matches are plotted as open circles. Unmatched satellites are plotted at an arbitrary small value of the ‘missing’ mass and so appear along the left and bottom edges of the diagrams. The numbers of matches of each type are listed explicitly. The number given for the unmatched SA satellites corresponds to objects with masses above $`2\times 10^{10}h^1M_{}`$, the minimum mass for an SPH galaxy.

For satellite–satellite matches the correspondence between the masses found using our two techniques is still reasonably good although with considerably more scatter than was the case for central galaxies. There is a tendency for the SPH masses to be larger than the SA masses suggesting that gas is continuing to cool onto the galaxies after the time when they are last identified as central objects. It should be noted that a group of points appearing in the lower portion of the plot, indicating SPH satellites matched to small mass SA satellites, are mostly those which are matched by our particular treatment for forming SPH galaxies, as explained above. In such cases, the halos hosting the forming galaxies merge into larger halos before the next output. Their central SA galaxies then become satellites, and their masses stay constant, whereas the SPH galaxies grow in mass continually until (and even after) they actually become satellites. Unsurprisingly, for the satellite – central galaxy matches the central galaxy mass is almost always substantially larger than the satellite galaxy mass, regardless of which is SPH and which SA. It is noticeable that as one moves to higher redshift, the number of unmatched SPH satellites increases substantially relative to the number of unmatched SA satellites. This is probably a reflection of the limited time resolution of our merger trees. At early times structure is building up very rapidly and cooling is efficient. The spacing between our outputs may then be too large to resolve the formation and merging of significant numbers of halos. These can form an SPH galaxy and turn it into a satellite between two levels of our numerical merger tree. In such cases the SA scheme will miss the satellite altogether. This highlights a different resolution issue from those discussed above. The SA scheme must store data sufficiently frequently that no significant stages of halo formation are missed.

In general, we conclude from the overall level of agreement seen for individual central and satellite galaxies that the standard recipe of gas cooling in the semi-analytic model gives an adequate description of the cooling process in the SPH model.

## 7 Summary and discussion

We have studied the cooling and condensation of gas onto protogalaxies using two very different treatments of the gas physics within the same simulation of dark matter evolution. One treatment uses schematic recipes taken from semi-analytic galaxy formation models, while the other tries to follow the hydrodynamics ‘correctly’ and in detail using an SPH approach. We have also tested two different implementations of the SPH equations, confirming Springel & Hernquist’s (2001) conclusion that the outcome of a simulation can depend sensitively on such details. In our case the two SPH approaches lead to ‘galaxy’ masses which differ by a factor of 2, even for objects which are apparently well resolved. We believe the energy and entropy conserving implementation of Springel & Hernquist (2001) to be the most reliable and have concentrated on it when comparing SPH results to those of the semi-analytic simulation.

We have characterised the effective resolution limit of simulations using each of our gas modelling techniques by the lower limit for which the simulated galaxy mass function agrees with that of a much higher resolution simulation carried out using the same technique. For the SA modelling this galaxy mass limit is just below the total baryon content of the least massive halo included in construction of the merger trees (10 particles in the models used here). For the SPH modelling the effective resolution limit is significantly higher, particularly at low redshifts when cooling times exceed the age of the Universe in most resolved halos. For the simulations studied here the SPH mass resolution limit at $`z=0`$ is roughly an order of magnitude larger than the SA limit, corresponding to a galaxy of about 75 SPH particles.

Above the resolution limit of both techniques (i.e. for galaxy masses above $`10^{11}h^1M_{}`$) the $`z=0`$ galaxy populations found by the two methods correspond remarkably well on an object-by object basis. The agreement is particularly good for the most massive galaxies in each halo (the ‘central’ galaxies). For satellite galaxies the agreement is still reasonably good for objects which match between the two types of simulation, but two significant types of discrepancy occur. One concerns objects which merge with the central galaxy in one of the simulations but not in the other. This is inevitable given the crudity of the recipe used to determine merging probabilities in the SA model. The difficulty can be avoided, at least for the more massive galaxies, with $`N`$-body simulations of high enough resolution. The dark matter substructures associated with satellite galaxies can then be tracked explicitly and no SA recipe is needed to follow their merging (see SWTK).

A more interesting discrepancy is the significant number of SPH satellites, which are unmatched in the SA simulation. Some of these may reflect insufficient time resolution in the simulation outputs from which the SA merger trees are constructed; halos may collapse, form galaxies, and merge into a larger halo between two of the stored outputs. Other unmatched SPH satellites may reflect more fundamental failures of the SA scheme. A forming dark halo may be sufficiently irregular to cool gas onto two galaxies rather than one, or cooling shocks associated with large-scale structure may fragment to make galaxies without any associated halo. The number of unmatched SPH satellites is clearly greater at early times than at $`z=0`$, suggesting that the process producing these objects is more effective when cooling times are short and structure growth is more rapid. This issue clearly warrants further investigation. Although the unmatched SPH satellites are a small population above the effective resolution limit of the SPH simulation, they probably signal a real failure of the SA scheme.

It is instructive to compare the computational efficiency of our two schemes. The comparison of this paper is based on implementing each of the schemes on top of a dark matter simulation with given $`N`$. For the SA scheme the overhead in post-processing, although non-negligible, is not dominant. Thus the computational cost is essentially that of running the high resolution cosmological $`N`$-body simulation. The SPH scheme requires an additional $`N`$ SPH particles, the force calculation for the SPH particles is considerably more complex than for the dark matter particles, and substantially shorter timesteps are needed to get convergent results in an SPH simulation with cooling than in the corresponding pure dark matter simulation. As a result, somewhat more than twice the memory and roughly an order of magnitude more CPU time are needed to carry out the SPH simulation than the corresponding pure $`N`$-body simulation. (The actual CPU factor for the simulations of this paper was about 20, but in neither case were the simulation parameters optimally chosen for the particular configuration treated.) We have seen, however, that the effective galaxy mass resolution of the SPH simulation is roughly an order of magnitude worse than that of SA simulation. Thus an SPH simulation with ten times as many particles would be needed to get a galaxy population with the same effective mass resolution as the SA simulation. This would require 20 times as much memory and more than two orders of magnitude more CPU time than the SA simulation.

In practice these considerations imply that simulations can be carried out using the SA technique which are far beyond the current capabilities of the SPH technique. For example, the S4 cluster simulation of SWTK has 250 times better mass resolution than the SA simulation of this paper and so of order 2000 times better resolution than our SPH simulation. It will be a long time before we can consider an SPH cluster/galaxy formation simulation with 2000 times more particles than the one analysed above. A further advantage of the SA technique is that the assumptions about uncertain aspects of the galaxy formation modelling (e.g. the mode and efficiency of star-formation or feedback) can be changed without any need to rerun the simulation. An example of this was shown in Figure 4, where the predictions of the high resolution S4 simulation in the absence of star-formation and feedback could be calculated from the stored simulation output. Thus for simulations concentrating on the global properties of the galaxy population and their dependence on the ‘subgrid’ physics (star formation, feedback, etc.) the SA technique is clearly the method of choice.

On the other hand, it is obvious that the SA technique can only be as good as the recipes on which it is based. Our work has shown that the standard gas cooling prescription works surprisingly well, but has also uncovered some apparent discrepancies. Discrepancies of this kind can only be studied and understood through comparison with simulations (or analytic models) which attempt to follow the relevant physics from first principles. Such models or simulations are also needed, of course, to design the SA recipes in the first place. In addition there are many aspects of galaxy formation where SA modelling can only give crude indications. Confirmation through direct simulation is then indispensable. Examples include the internal structure of galaxies, the structure of the intergalactic medium, and the comparison of both of these with the wealth of available data. Galaxy formation and its interaction with the environment is a complex and highly nonlinear process involving many different aspects of astrophysics. Gaining some understanding will require concerted use of all the observational and theoretical weapons at our disposal. SA and SPH simulation techniques will both remain major parts of our theoretical arsenal.

We thank Simone Marri and Hugues Mathis for helpful comments.


text:
# Spiral Molecular Front in Galaxies: Quick Transition from Atomic to Molecular Hydrogen in Spiral Arms

## 1 Introduction

The neutral (HI) and molecular hydrogen (H2) gases are the two major components of the interstellar matter. Elmegreen (1993) has proposed a simple model to predict the molecular fraction by assuming a quasi-static equilibrium between the two phases in individual interstellar clouds. However, the dynamical transition between HI and H2 gas phases are not well understood in real galaxies: Are the GMCs (giant molecular clouds) long-lived clouds living for several galactic rotations crossing the spiral arms for many times, coexisting with the HI gas? Alternatively, are they formed quickly from HI gas when the inter-arm HI encounters the arms, and does molecular gas return to HI again, when going out of the arms?

This question can be answered by studying the distribution of molecular fraction, the ratio of molecular-gas density to the total gas density, in the gas disk of spiral galaxies. Global variations of the molecular fraction in disk galaxies have been obtained by our earlier works (Sofue et al 1995; Honma et al 1995), applying the Elemegreen’s theory of atomic-to-molecular gas transition. We showed that the inner gas disk is dominated by molecular gas, while the outer disk by HI, and the transition between HI and H2 occurs in a narrow annulus in the disk, across which the molecular fraction varies drastically. We call this transition region the “molecular front (MF)”. Such front phenomenon was found also in the direction perpendicular to the galactic plane, so that the H2 disk is sandwiched by a fat HI disk (Imamura and Sofue 1997). However, these analyses do not tell us about dynamics and time scale of the HI-to-H2 transition.

In the present paper, we investigate the MF phenomenon across spiral arms by obtaining a two-dimensional distribution of molecular fraction in a galactic disk. Azimuthal variation of molecular fraction across spiral arms will provide us with information about dynamical evolution of the HI-H2 phase transition and molecular gas formation in spiral arms.

## 2 Molecular Front

#### 2.0.1 Radial and Vertical Molecular Front

Many observations of the neutral hydrogen and the molecular gas have been carried out, and have revealed two important aspects. First, the neutral hydrogen gas is broadly distributed in the disk, and is often extended outside the optical disks, while it is deficient in the central regions. On the other hand, molecular gases are concentrated in the central a few kpc regions. The radial variation of the distributions of HI and H2 gases have been quantified by analyzing the molecular fraction as a function of the galactco-centric distance. The molecular fraction, $`f_{\mathrm{mol}}`$, is defined by

$$f_{\mathrm{mol}}=\frac{\rho (\mathrm{H}_2)}{\rho (\mathrm{HI})+\rho (\mathrm{H}_2)}=\frac{2\times n(\mathrm{H}_2)}{n(\mathrm{HI})+2\times \mathrm{n}(\mathrm{H}_2)}.$$
(1)
Sofue et al. (1995) have investigated radial variation of the molecular fraction in several spiral galaxies. They found that the molecular fraction is almost unity in the central a few kpc region, and decreases suddenly at a certain radius to almost zero, beyond which the gas is almost totally in the HI phase. Such a step-like variation of the molecular fraction is called the molecular front.

Honma et al. (1995) modeled the molecular front by applying the phase transition theory of Elmegreen (1993). According to this theory, the molecular fraction is determined by three parameters: interstellar pressure $`P`$, UV radiation field $`U`$ , and metallicity $`Z`$. Honma et al showed that, if these three parameters are approximated by exponential functions of galacto-centric radius, the observed molecular front in the radial direction can be well reproduced. They applied this model to several real galaxies, and obtained good coincidence with observed data.

Spiral galaxies have density waves and spiral arms. As the spiral pattern speed is much slower than the gas’s rotation, the gases collide with the spiral potential, resulting in galactic shock waves, where the gas density and pressure increase suddenly. The compression results in star formation, leading to increase of the radiation field intensity. Thus, the MF parameters, $`(P,U,Z)`$, will vary as a function of azimuthal angle, which may results in another type MF spatially correlated with the spiral arms. In order to investigate the existence of such spiral MF, we first examine the HI and CO-line data from the literature, and construct a 2D map of the molecular fraction in M51. Then, we perform numerical simulation using the transition model established by Elmegreen (1993), and compare with the observations.

### 2.1 Spiral Molecular Front in M51: Azimuthal Variation of Molecular Fraction

For obtaining two-dimensional map of molecular fraction, we make use of the data from the literature that provide sufficient angular resolution in both HI and CO. Here, we choose the ‘grand-design’ spiral galaxy NGC 5194 (M51), which was also studied in Honma et al. (1995). We use the 12CO ($`J=10`$) imaging data obtained with the Nobeyama 45-m telescope (Nakai et al 1994) with the beam width of $`16^{\prime \prime }`$ and grid spacing of $`15^{\prime \prime }`$, corresponding to an angular resolution of $`24^{\prime \prime }`$, and the HI data obtained with the VLA at $`34^{\prime \prime }\times 34^{\prime \prime }`$ resolution from Rots et al. (1990). Integrated-intensity maps of HI and CO are shown in Fig. 1 at the same resolution of $`34^{\prime \prime }`$. These two maps reveal the typical behaviors of atomic and molecular hydrogen gases: HI is extended broadly in the disk, having central deficient, while CO is concentrated in the central region.

—- Fig. 1 —-

We calculate the column density of HI and H2 from the integrated intensity, adopting the following relation between the column density and intensity:

$`n(\mathrm{HI})`$ $`=`$ $`C(\mathrm{HI})\times \mathrm{I}(\mathrm{HI}),`$ (2)
$`n(\mathrm{H}_2)`$ $`=`$ $`C(\mathrm{CO})\times \mathrm{I}(\mathrm{CO}),`$ (3)
where $`C(\mathrm{HI})`$ and $`C(\mathrm{CO})`$ are the conversion factors. The value of $`C(\mathrm{HI})`$ can be easily calculated as $`1.82\times 10^{18}`$ cm-2 K-1 km-1 sec. In order to estimate the amount of the total molecular gas mass from CO intensity, we adopt a conversion factor of $`C(\mathrm{CO})=1.1\times 10^{20}`$ cm-2 K-1 km-1 sec from Arimoto et al. (1996), and assume no radial dependence of the conversion factor. The difference of the angular resolution was corrected, so that both HI and CO maps have the same resolution of $`34^{\prime \prime }`$, and the grid separation was corrected by using bicubic interpolation. Using Eq. 1, we calculated the molecular fraction, and show a two dimensional distribution map of the molecular fraction in Fig. 2.

—- Fig. 2 —-

Fig. 2 shows that the molecular fraction is almost unity in the central $`r<60^{\prime \prime }`$ region, or at $`r<2.8`$ kpc. Toward the northwest and southeast, where the bisymmetric spiral structure is not seen clearly in both HI and CO distribution, the molecular fraction decreases drastically to zero at $`r>110^{\prime \prime }`$, or $`r>5.1`$ kpc. This sudden radial decrease of the molecular fraction confirms the radial MF at $`r5`$ kpc.

In the northeast and southwest, where the molecular spiral arms are prominent, a high-$`f_{\mathrm{mol}}`$ region extends until $`7`$ kpc. In the southeast direction, $`f_{\mathrm{mol}}`$ suddenly decreases to 0.6 at $`r60^{\prime \prime }`$, and increases again to a value greater than 0.9, and finally decreases to zero. This re-increase of the molecular fraction results from the spiral arms, as already mentioned by Kuno et al. (1997).

Fig. 2 shows that the molecular front exists along the spiral arms in addition to the global radial variation. This implies that the circularly flowing gas experiences quick transition from atom to molecule by encountering the spiral arms, and from molecule to atom when going out from the spiral arms. However, the present angular resolution is not high enough to investigate the scale length of the variation across the arms in more detail. We can estimate only the upper limit of the scale length, which is about 700 pc in the azimuthal direction. Since the rotation velocity is about 200 km/sec (Kuno & Nakai 1997) and pattern speed of the density waves will be $`20`$ km s-1 kpc-1, the upper limit of time scale of atom-to-molecule and molecule-to-atom phase transitions is estimated to be several Myr.

## 3 Simulation of the Molecular Front

### 3.1 Numerical Methods

In previous section, we have shown the existence of MF in azimuthal direction. However, because of poor angular resolution, we could not estimate the exact MF scale length in the arms. Here, we try to estimate the scale length by numerical simulations. The azimuthal variances of the gas density is calculated using the hydrodynamical code, not taking into account the self-gravity of the gas. The variation of metallicity and radiation field is assumed to have an axisymmetric exponential distribution, as adopted by Honma et al. (1995), but no azimuthal variation is assumed.

For calculating the gas dynamical behavior of the disk, we use the freely downloadable and usable hydrodynamical code, VH-1 (Blondin & Lufkin 1993). This is a multidimensional hydrodynamics code for an ideal compressible fluid written in FORTRAN, developed by the numerical astrophysics group at the University of Virginia on the bases of the Piecewise Parabolic Method Lagrangian Remap (PPMLR) scheme of Colella & Woodward (1984). The PPMLR has the advantage of maintaining contact discontinuities without the aid of a contact steepener, and is good to be applied to a galactic scale hydrodynamical simulation. The code does not include self-gravity, artificial viscosity, variable gamma equation of state, and radiative heating and/or cooling. Hence, we assume that the interstellar gas is ideal, inviscid, and compressible.

#### 3.1.1 Gravitational Potential

In order to simulate the evolution of spiral structure in the gas disk, we assume a given gravitational potential, which comprises the following two terms: (i) a static axisymmetric potential, and (ii) a nonaxisymmetric, rotating bar potential. The self-gravity of the gas is not taken into account. The potential is expressed by

$$\mathrm{\Phi }(R,\varphi )=\mathrm{\Phi }_0(R)+\mathrm{\Phi }_1(R,\varphi ).$$
(4)
We adopt a “Toomre disk” (Toomre 1981) potential for the axisymmetric component as given by:

$$\mathrm{\Phi }_0(R)=\frac{c^2}{a}\frac{1}{(R^2+a^2)^{1/2}},$$
(5)
where $`a`$ is the core radius, and $`c=v_{\mathrm{max}}(27/4){}_{}{}^{1/4}a`$ with $`v_{\mathrm{max}}`$ being the maximum rotation velocity. Through our numerical simulation, we fix the core radius and maximum circular velocity to be $`a=\sqrt{2}`$ kpc and $`v_{\mathrm{max}}=200`$ km s-1 .

The nonaxisymmetric potential is taken from Sanders (1977), assuming rigid rotation at a pattern speed $`\mathrm{\Omega }_p`$, which has the form

$$\mathrm{\Phi }_1(R,\varphi )=\epsilon \frac{aR^2}{(R^2+a^2)^{3/2}}\mathrm{\Phi }_0(R)\mathrm{cos}2\left(\varphi \mathrm{\Omega }_pt\right),$$
(6)
where $`\epsilon `$ is the strength of the bar of the order of $`\epsilon =0.15`$. Spiral shocked arms of gas are produced by this potential.

#### 3.1.2 Initial Condition

Initially, we set $`256\times 256`$ two-dimensional cells corresponding to $`12.8\times 12.8`$ kpc field, setting the field center at the coordinates origin. The initial number density of the gas is taken to be 5 cm-1 for the inner disk at $`R8`$ kpc, and 1 cm-1 at $`R>8`$ kpc. Initial rotation velocity of each gas cell is set so as for the centrifugal force to balance the gravitation. The bar pattern speed $`\mathrm{\Omega }_p`$ is taken to be 23 km s-1 kpc-1, and the strength of the bar $`\epsilon `$ is taken to be 0.10.

### 3.2 Elmegreen’s Parameters

We specify three parameters: interstellar pressure $`P`$, radiation field strength $`U`$, and metallicity $`Z`$. Following Honma et al. (1995), we represent $`U`$ and $`Z`$ by the following equations.

$`{\displaystyle \frac{U}{U_{}}}`$ $`=`$ $`\mathrm{exp}\left({\displaystyle \frac{RR_{U_s}}{R_U}}\right),`$ (7)
$`{\displaystyle \frac{Z}{Z_{}}}`$ $`=`$ $`\mathrm{exp}\left({\displaystyle \frac{RR_{Z_s}}{R_Z}}\right).`$ (8)
Here, $`R_U`$ and $`R_Z`$ are the scale radii of the radiation field and metallicity, respectively, and $`R_{U_s}`$ and $`R_{Z_s}`$ are the radii at which the quantities are normalized to the solar values.

It is possible that radiation field depends on the azimuthal direction, because star formation along the spiral arms is enhanced and newly born OB stars radiate strong UV emission. Thus, it seems that there is little UV radiation in the inter-arm regions. However, Greenawalt et al. (1998) presented deep H$`\alpha `$ emission-line images of several spiral galaxies including M51, concluding that a half of total H$`\alpha `$ emission is contributed by the diffuse ionized gas. This observation shows that the fluctuation of ionizing radiation field in arm and inter-arm regions is not more than twice at their resolution of a few arcsecs. In the present simulation, the resolution is much lower, and the arm-interarm UV fields may be assumed to have intensity fluctuation less than twice: We here assume that the UV field is uniform along an azimuthal circle, and depends only on the galactocentric distance. Although they are a scope beyond the present paper and will be subject for simulations in the future, variation of radiation fields by star-forming arms should be taken into account, and UV radiation transfer through the arms, particularly through dark lanes, should be solved, in order to obtain a more realistic, higher resolution behavior of the molecular-fraction across spiral arms.

The interstellar pressure is estimated as $`P\rho (\mathrm{gas})`$, in the same manner as Honma et al. (1995), and express it by the following relation:

$$\frac{P}{P_{}}\frac{\rho }{\rho _{}}.$$
(9)
We adopt the following normalizing parameters: $`\rho _{}=2\mathrm{H}\mathrm{c}\mathrm{m}^3`$, $`R_U=R_Z=5.2`$ kpc, and $`R_{U_s}=R_{Z_s}=11.5`$ kpc. The latter two parameters are as same as those assumed by Honma et al. (1995). Although they assumed constant metallicity for NGC 5194, we here assume an exponential function as above (see e.g. Arimoto et al 1996).

### 3.3 Simulated Spiral Molecular Front

Under these assumptions and initial conditions, we first calculated the time evolution of the distribution of total gas. Fig. 3a shows the density distribution after several galactic rotations. The density distribution in the initially uniform-density disk is strongly disturbed by the oval potential, and high-density arms are formed, which evolve into well-developed bisymmetric spiral arms. In the same time, the global radial density distribution is regulated to have a roughly exponentially-decreasing structure. Inside the bar, straight high-density arms are formed along the axis of the bar, as well as many faint spiral-like arms. We are here interested in the global spiral pattern in order to compare with the observed gas density distributions, and will not discuss the smaller and fainter features, which will be finally smeared out when we compare the result with the observation.

— Fig, 3a, b, c, d —

The simulated distribution of gas density is, then, used to calculated the distribution of gas pressure $`P`$, which is further combined with the assumed distributions of UV radiation field $`U`$ and the metallicity $`Z`$ as given above. Thus, we can calculated the molecular fraction at each grid point, and obtain a map of $`f_{\mathrm{mol}}`$, as well as maps of atomic and molecular gases. Fig. 3b, c and d shows snapshots of the thus calculated maps for HI, H2 and $`f_{\mathrm{mol}}`$. In order to compare with the observations, we smooth the map of molecular fraction by a Gaussian convolution, and compare it with the observation of M51 in the same resolution in Fig. 4.

— Fig. 4 —

## 4 Discussion

The simulated results can be summarized as follows. Detailed arm shapes are not well reproduced, because the simulation uses a fixed bar-potential. In order for detailed arm features to be reproduced, we need self-gravitating $`N`$-body simulation, including both stars and gas, is necessary, which is a subject for the future. However, the simulation in Fig. 3 and 4 shows global agreement with the observations, and reproduces the following characteristics of the distribution of molecular fraction in M51. Around the nucleus $`f_{\mathrm{mol}}`$ is almost unity, and maintains high values $`r4`$ kpc, where $`f_{\mathrm{mol}}>0.8`$. At $`r4`$ to 5 kpc, $`f_{\mathrm{mol}}`$ decreases drastically to 0.4 or lower, and becomes nearly zero in the outer region. It is clear that the transition region, where $`f_{\mathrm{mol}}`$ is $``$ 0.5 – 0.6, is distributed in a narrow annulus of radius of 4.5 kpc, which means that the HI-to-H2 transition occurs rather abruptly within a narrow range of radius. This confirms the radial molecular front (Sofue et al. 1995, Honma et al. 1995).

In addition to the radial MF in the gas disk, azimuthal fluctuation of molecular fraction is also significant. At $`r`$ 6 kpc, $`f_{\mathrm{mol}}`$ changes from nearly zero before encountering the arm to $``$ 0.8 within $``$ 200 pc in the azimuthal direction in the spiral arm. This narrow transition range corresponds to a time scale of $``$ 2 Myr. This indicates that the molecular gas observed along spiral arms is formed only within $``$ 2 Myr. This time scale is consistent with the H2 forming time scale $`>7\times 10^5`$ yr predicted by Takahashi & Williams (2000) using a molecular dynamics simulation.

On the other hand, the decaying time scale of molecular gas is longer, according to the more gradual decrease of gas density after the galactic shock passage. If we define a molecular arm an armed region with $`f_{\mathrm{mol}}>0.5`$, the width of the molecular arm is estimated to be $``$ 2.5 kpc in the azimuthal direction. This corresponds to a lifetime of molecular gas of about 25 Myr, which is consistent with the semi-analytically calculated GMCs lifetimes $``$ 40 Myr (Elmegreen 1991).

If we consider the radial variation of metallicity $`Z`$, these local (arm) $`f_{\mathrm{mol}}`$ structures will depend on the galactocentric distance, even for the same gas density. Fig. 3 shows that, for the same density of $`5`$ H cm-3, the molecular fraction is $``$ 0.5 – 0.6 in the inner inter-arm region, while it is 0.3 – 0.4 in the outer arm regions. This reflects the fact that the molecular fraction depends not only on the gas density, but also on the strength of radiation field and metallicity, particularly the latter.

We are grateful to Dr. M. Honma for helpful advice, and Mr. J. Koda for discussions on numerical simulations. The numerical simulations were performed using the VH-1 code developed by the numerical astrophysics group at the University of Virginia.

References

Arimoto, N., Sofue, Y., & Tsujimoto, T. 1996, PASJ 48, 275

Blondin, & Lufkin, 1993, ApJS, 88, 589

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

Elmegreen, B. G., 1991, in “The Physics of Star Formation and Early Stellar Evolution, ed. Lada. C. J., Kylafis, N. D. (Kluwer Academics), 35

Elmegreen, B. 1993, ApJ. 411, 170

Greenawalt, B., Walterbos, R. A. M., Thilker, D., Hoopes, C. G. 1998 ApJ, 506, 135.

Imamura, K, & Sofue, Y., 1997, AAp, 319, 1

Honma, M., Sofue, Y., & Arimoto, N., 1995 AAp, 304, 1

Kuno, N., Nakai, N., Handa, T., & Sofue, Y., 1994, PASJ, 47, 745

Kuno, N., & Nakai, N., 1997, PASJ, 49, 279

Nakai, N., Kuno, N., Handa, T., & Sofue, Y., 1994, PASJ, 46, 527

Rots, A. H., Crane, P. C., Bosma, A., Athanassoula, E., & van der Hulst, J. M., 1990, AJ, 100, 387

Sanders, R. H. 1977 ApJ 217 916

Sofue, Y. 1994 ApJ 423 207

Sofue, Y., Honma, M., & Arimoto, N., 1995, AAp, 296, 33.

Takahashi, J., & Williams, D. A., 2000, MNRAS, 314, 27

Toomre, A., 1981, in The Structure and Evolution of Normal Galaxies, ed S. M. Fall & D. Lynden-Bell (Cambridge: Cambridge Univ. Press), 111

Figure Captions

Fig. 1. (a) HI integrated-intensity map of NGC 5194 (M51) at $`34^{\prime \prime }`$ resolution (Rots et al. 1990)(left panel). The contour levels are 100, 300, 500, …, and 1900 mJy/beam and the intensity key is shown at the bottom.

(b) 12CO ($`J=10`$) integrated intensity map of M51 (Nakai et al. 1994), smoothed to the same spatial resolution ($`34^{\prime \prime }`$) as that of HI. The contour interval is 6.6 K km sec-1, and the intensity key is shown at the bottom. $`10^{\prime \prime }`$ in the maps correspond to linear scale of 466 pc.

Fig. 2. Molecular fraction ($`f_{\mathrm{mol}}`$) in NGC 5194 (M51). Contour levels are $`f_{\mathrm{mol}}=0.1,0.2,\mathrm{},0.9`$.

Fig. 3. (a) A snapshot of numerically simulated total density distribution of interstellar gas.

(b) The same for HI gas.

(c) The same for H2 molecular gas, simulating a CO-line intensity distribution.

(d) The same for molecular fraction $`f_{\mathrm{mol}}`$.

The value keys are shown at the bottom.

Fig. 4. Comparison of the simulated molecular fraction (left) smoothed to the same resolution as that of the observation of M51 (right).