format: text

text:
# Virial shocks in galactic haloes?

## 1 Introduction

The standard lore in the idealized picture of galaxy formation by spherical infall of gas inside dark-matter haloes is that the gas is first heated to the halo virial temperature behind an expanding virial shock. It is then supported by pressure in a quasi-static equilibrium while it is cooling radiatively and is slowly contracting to a disc where it can eventually form stars. The cooling process thus determines important galaxy properties such as the star-formation rate and the metal enrichment, so it is necessarily an important ingredient in the galaxy formation process.

However, it is not at all clear that a stable shock can persist in the halo gas away from the disc under the conditions valid in many galactic haloes. In the absence of a virial shock, the gas is not heated to the virial temperature until it falls all the way to the disc, where the collapse stops and the gas is heated in a thin layer. This may alter some of the assumed processes of disc formation and in particular the star formation rate in it. It may work against blowout by supernova-driven winds in dwarf galaxies. The result of heating near the disc instead of at the virial radius may result in weakening the soft x-ray emission from such haloes and producing a high flux of $`L_\alpha `$ instead. In this paper we evaluate the conditions for the existence of a virial shock in galactic haloes.

Initial density perturbations are assumed to grow by gravitational instability, reach maximum expansion, and collapse into virial equilibrium at roughly half the maximum-expansion radius. During the initial phase, and roughly until shells start crossing each other near the virial radius, the gas pressure is negligible compared to the gravitational force, so the shells of gas and dark matter move in a similar manner. Once interior to the virial radius, where shells tend to cross and the gas density becomes high enough, the gas pressure becomes an important player in the dynamics. Its hydrodynamic properties allow transfer of bulk kinetic energy into internal energy and the pressure prevents gas element from passing through other gas elements and from being compressed without limit. This makes the infall velocity vanish at the centre. Since in the cold infalling gas the typical velocity is higher than the speed of sound, the information about this inner boundary condition cannot propagate outwards in time, and these supersonic conditions create a shock. After the gas crosses the shock, it is heated up, the speed of sound increases, and the flow becomes subsonic.

The shock transfers the kinetic energy that has been built during the collapse into internal gas energy just behind the shock. A stable spherical shock would slowly propagate outwards through the infalling gas, leaving behind it hot, high-entropy gas that is almost at rest. The temperature of the post-shock gas roughly equals the virial temperature. The persistence of the shock depends on sufficient pressure by the post-shock gas, which supports it against being swept inwards due to the gravitational pull together with the infalling matter. Radiative gas cooling makes the gas lose entropy and pressure, which weakens the pressure support behind the shock front. Our approach here is to evaluate the existence of a virial shock by analyzing the gravitational stability of the supporting gas behind the shock in the presence of significant cooling.

In §2 we first summarize the standard analysis of an adiabatic shock and then generalize the gravitational stability criterion to the case where cooling is important. In §3 we describe our spherical hydrodynamic Lagrangian code, which includes gravitating dark-matter and gas shells, artificial viscosity, radiative cooling and centrifugal forces. We test the code in this section and in Appendix A. In §4 we apply the numerical code to simulations which demonstrate the shock formation and test the validity of the analytical model. In §5 we apply the shock stability criterion to realistic haloes forming in cosmological conditions. In §6 we summarize our results and discuss potential astrophysical implications.

## 2 Shock stability analysis

Our goal here is to derive a criterion for the existence of a virial shock in terms of the properties of the infalling gas just in front of the shock front. It is based on a gravitational stability analysis of the post-shock gas. We first remind ourselves of the standard stability analysis in the simple adiabatic case, and then derive a more general criterion for stability in the radiative case, under certain assumptions and using a perturbation analysis.

### 2.1 The standard adiabatic case

Throughout this paper, we treat the baryons as an ideal monoatomic gas. Their equation of state could therefore be written as

$$P=(\gamma 1)e\rho ,$$
(1)
where $`P`$ is the pressure, $`e`$ is the specific internal energy, $`\rho `$ is the density of the gas and $`\gamma `$ is the adiabatic index. Along an isentrope (an adiabatic process of constant entropy) the pressure and density are related via $`P\rho ^\gamma `$, so the adiabatic index is defined by

$$\gamma =\left(\frac{\mathrm{ln}P}{\mathrm{ln}\rho }\right)_\mathrm{s}.$$
(2)
For a monoatomic gas $`\gamma =5/3`$.111As the temperature exceeds the binding energy of the hydrogen and helium atoms, electrons become detached from the nuclei and $`\gamma `$ becomes smaller. Once the gas becomes fully ionized, the original value of $`5/3`$ is restored, but with a different effective density. This should have only a marginal effect on our results, and is ignored in this paper.

The virial shock is assumed to be a spherical accretion shock which propagates outwards slowly while infalling gas crosses it inwards.

The kinetic energy of the infalling gas is transformed at the shock front into thermal energy — the post-shock gas is thus heated to a temperature close to the virial temperature of the system of dark-matter halo and gas, $`V_{\mathrm{infall}}^2k_\mathrm{B}T_{\mathrm{vir}}`$. Because the original temperature of the infalling gas is negligible compared to the virial temperature, the system obeys the strong-shock limit. When we denote the pre-shock and post-shock quantities by subscripts 0 and 1 respectively, the jump conditions across the shock are in this case (Zel’dovich & Raiser, 1966):

$$\rho _0=\frac{\gamma 1}{\gamma +1}\rho _1,$$
(3)
$$(u_0u_\mathrm{s})=\frac{\gamma +1}{\gamma 1}(u_1u_\mathrm{s}),$$
(4)
$$P_1=\frac{2\rho _0u_0^2}{\gamma +1},$$
(5)
$$T_1=\frac{\mu }{k_bN_a}\frac{P_1}{\rho _1}=\frac{\mu }{k_bN_a}\frac{2\gamma 1}{(\gamma +1)^2}u_0^2,$$
(6)
where $`u`$ stands for radial velocity, $`u_\mathrm{s}`$ is the shock velocity, $`N_a`$ is Avogadro’s number, $`N_a/\mu `$ is the average number of molecules per unit mass, and $`k_B`$ is Boltzman’s constant.

According to standard shock theory, the post-shock gas is always sub-sonic (in the frame of reference of the moving shock) because of the increase of the sound velocity behind the shock. This gas is thus capable of providing the necessary pressure to support the shock against the gravitational pull inwards applied by the self-gravity of the gas and the dark-matter halo as well as the pressure applied by the infalling matter at the shock front.

The criterion for gravitational stability of this post-shock gas in the adiabatic case is the standard Jeans stability criterion: $`\gamma >4/3`$ (e.g., Cox, 1980, Chapter 8)

If the post-shock gas is gravitationally unstable, it falls into the galaxy centre on a dynamical timescale and can no longer support the shock. As a result, the shock weakens and it is swept inwards.

The Jeans criterion can be qualitatively understood in terms of the following heuristic derivation. For a shell of radius $`r`$, we compare the gravitational pull inwards, $`a_\mathrm{g}=GM/r^2`$ (where $`M`$ is the mass interior to $`r`$), to the pressure pushing outwards, $`a_\mathrm{p}=\rho ^1P`$. We assume an isentrope, $`P\rho ^\gamma `$. We also assume homology, such that the local density scales like the mean density in the sphere interior to $`r`$, $`\rho M/r^3`$. Then $`P`$ can be replaced by $`P/r`$ and we obtain

$$\frac{a_\mathrm{p}}{a_\mathrm{g}}\rho ^{\gamma 4/3}.$$
(7)
If $`\gamma <4/3`$, we have an unstable configuration. Starting in hydrostatic equilibrium, $`a_\mathrm{p}/a_\mathrm{g}=1`$, a perturbation involving contraction is associated with a larger $`\rho `$, and therefore $`a_\mathrm{p}/a_\mathrm{g}<1`$ by eq. (7), implying that the pressure cannot prevent collapse. If $`\gamma >4/3`$, the pressure force increases until it balances the increased gravitational pull. We note that even this simple derivation of the Jeans criterion had to assume homology — an assumption that we will have to adopt also in our analysis of the radiative case below.

### 2.2 Shock stability under radiative cooling

We wish to replace the adiabatic Jeans criterion by a more general stability condition that will be valid also in the radiative case. This criterion must depend on the cooling rate and should therefore be naturally expressed in terms of time derivatives. We generalize the adiabatic $`\gamma `$ of eq. (2) by an effective $`\gamma `$ following a comoving volume element along its Lagrangian path:

$$\gamma _{\mathrm{eff}}\frac{d\mathrm{ln}P/dt}{d\mathrm{ln}\rho /dt}.$$
(8)
We expect that the system would be stable when $`\gamma _{\mathrm{eff}}`$ is larger than a certain critical value, the analog of the requirement $`\gamma >4/3`$ in the adiabatic case.

In our Lagrangian analysis all the quantities ($`r`$, $`u`$, $`\rho `$, $`P`$, etc.) refer to comoving shells; they are all functions of the gas mass $`m`$ interior to radius $`r`$ and time $`t`$. Derivatives with respect to time following a comoving volume element will be denoted by an upper dot, and derivatives with respect to m will be denoted by a prime.

The effective gamma can be related to its adiabatic analog given the cooling rate and other post-shock gas quantities. The time derivative of eq. (1) yields:

$$\dot{P}=(\gamma 1)(\dot{e}\rho +e\dot{\rho }).$$
(9)
Energy conservation in the presence of radiative losses can be expressed by

$$\dot{e}=P\dot{V}q=\frac{P\dot{\rho }}{\rho ^2}q,$$
(10)
where $`q`$ is the radiative cooling rate \[to be discussed below, e.g., eq. (20)\] and $`V=\rho ^1`$ is the specific volume. Substituting $`\dot{e}`$ from eq. (10) into eq. (9), and using it in eq. (8), we obtain

$$\gamma _{\mathrm{eff}}=\gamma \frac{\rho }{\dot{\rho }}\frac{q}{e}.$$
(11)
Note that in the limit $`q/e\dot{\rho }/\rho `$ we reproduce the adiabatic case; the process is nearly adiabatic when the cooling timescale is long compared to the contraction timescale.

We assume that in the region close behind the shock the pattern of the velocity field is homologous. By this we mean that at any given time the (radial) velocity is proportional to the radius (as in a Hubble flow), namely

$$u/r=u_1/r_\mathrm{s},$$
(12)
thus providing a boundary condition for the post-shock gas. The homology is shown to be a valid approximation in the simulations discussed below, where the post-shock shell trajectories are roughly parallel to each other in the $`\mathrm{log}rt`$ plane, at any given time close enough to shock crossing. The time evolution of the density can then be evaluated via the continuity equation in Lagrangian form for the spheri-symmetric case,

$$\frac{\dot{\rho }}{\rho }=𝐮=\frac{1}{r^2}\frac{}{r}(r^2u)=\frac{3u_1}{r_\mathrm{s}},$$
(13)
where the last equality results from the assumed homology, eq. (12). The homology thus implies that $`\dot{\rho }/\rho `$ at a given $`t`$ is a constant in $`m`$ throughout the post-shock region. Eq. (11) can then be simplified:

$$\gamma _{\mathrm{eff}}=\gamma +\frac{r_\mathrm{s}}{3u_1}\frac{q}{e}.$$
(14)
We start with a hypothetical unperturbed state for the post-shock gas, where we assume that the net force vanishes, $`\ddot{r}=0`$. The system adjusts itself to this state on a timescale associated with the speed of sound $`c_\mathrm{s}`$, provided that it is much higher than the infall velocity $`u`$. This is expected to be the case in the sub-sonic post-shock medium, where $`c_\mathrm{s}`$ becomes high and $`u`$ becomes low. The unperturbed equation of motion in Lagrangian form is then

$$\ddot{r}=4\pi r^2P^{}\frac{GM}{r^2}=0,$$
(15)
where $`M`$ is the total mass interior to radius $`r`$.

We then introduce a perturbation due to a homologous infall velocity $`u`$. Over a short time interval $`\delta t`$, it introduces a small displacement inwards, $`\delta r=u\delta t`$. In order to distinguish between stability and instability we wish to determine whether the induced acceleration, $`\ddot{\delta r}`$, is positive or negative, tending to decrease or increase the velocity respectively. Note that under homology, eq. (12), the relative displacement is

$$\delta r/r=u_1\delta t/r_\mathrm{s}.$$
(16)
Writing the equation of motion, eq. (15), but for the perturbed quantities $`P+\delta P`$ and $`r+\delta r`$, and subtracting the unperturbed eq. (15), we obtain to first order

$$\ddot{\delta r}=4\pi r^2(\delta P)^{}+\frac{4GM\delta r}{r^3}.$$
(17)
We next manipulate the right-hand side of eq. (17) to obtain a simple expression involving $`\gamma _{\mathrm{eff}}`$.

In the second term we use the homology, eq. (16), and then the unperturbed equation of motion, eq. (15), to obtain

$$\frac{4GM\delta r}{r^3}=\frac{16\pi r^2u_1\delta t}{r_\mathrm{s}}P^{}.$$
(18)
The manipulation of the first term is somewhat more elaborate. We use the definition of $`\gamma _{\mathrm{eff}}`$, eq. (8), to write

$$\delta P=(\dot{\rho }/\rho )P\gamma _{\mathrm{eff}}\delta t=(3u_1/r_\mathrm{s})P\gamma _{\mathrm{eff}}\delta t,$$
(19)
where the second equality is due to eq. (13). Note that the $`m`$ dependence in this term is only in the product $`P\gamma _{\mathrm{eff}}`$. We now express $`\gamma _{\mathrm{eff}}`$ in terms of the cooling rate $`q`$ as in eq. (14), and need to take the derivative $`(Pq/e)^{}`$. We make here the standard assumption that the radiative cooling rate is proportional to density,

$$q=\rho \mathrm{\Lambda }(T),$$
(20)
with $`\mathrm{\Lambda }(T)`$ the macroscopic cooling function and $`T`$ the post-shock temperature. The immediate post-shock medium is assumed to be isothermal, reflecting via the jump conditions an assumed approximate uniformity of the pre-shock gas over a short time interval. Using eq. (1) we have $`P/e=(\gamma 1)\rho `$, and together with eq. (20) it becomes

$$Pq/e=(\gamma 1)\mathrm{\Lambda }\rho ^2.$$
(21)
In the computation of $`(Pq/e)^{}`$, we first replace $`\rho ^{}=(d\rho /dP)P^{}`$, then use eq. (8) to write $`d\rho /dP=\gamma _{\mathrm{eff}}^1\rho /P`$, use eq. (21) backwards to replace $`(\gamma 1)\mathrm{\Lambda }\rho ^2/P`$ by $`q/e`$, and finally use eq. (14) to obtain $`(Pq/e)^{}=3u_1/r_s[2\gamma _{\mathrm{eff}}^1(\gamma \gamma _{\mathrm{eff}})P^{}]`$. We thus have in the first term of the rhs of eq. (17)

$$(\delta P)^{}=\frac{3u_1\delta t}{r_\mathrm{s}}P^{}[\gamma 2\gamma _{\mathrm{eff}}^1(\gamma \gamma _{\mathrm{eff}})].$$
(22)
With the right-hand side of eq. (17) given by eq. (22) and eq. (18), the first-order equation finally becomes,

$$\ddot{\delta r}=\frac{12\pi r^2u_1\delta tP^{}}{r_\mathrm{s}}\left[\gamma 2\gamma _{\mathrm{eff}}^1(\gamma \gamma _{\mathrm{eff}})\frac{4}{3}\right].$$
(23)
Since $`u_1`$ and $`P^{}`$ are both always negative, the desired sign of $`\ddot{\delta r}`$ is determined by the sign of the expression inside the square brackets. Note that in the adiabatic case, $`q=0`$, we have $`\gamma _{\mathrm{eff}}=\gamma `$, so we recover the standard stability criterion, $`\gamma >4/3`$. In the radiative case, $`\gamma _{\mathrm{eff}}\gamma `$, we finally obtain the generalized stability criterion:

$$\gamma _{\mathrm{eff}}>\frac{2\gamma }{\gamma +2/3}\gamma _{\mathrm{crit}}$$
(24)
For a monoatomic gas, where the adiabatic value is $`\gamma =5/3`$, the threshold for stability is $`\gamma _{\mathrm{crit}}=10/7=1.43`$, which is close but not identical to the adiabatic threshold $`4/3`$.

### 2.3 Stability in terms of pre-shock quantities

Next, we wish to express $`\gamma _{\mathrm{eff}}`$ and the stability criterion in terms of the properties of the pre-shock gas; the infall velocity $`u_0`$ and the gas density $`\rho _0`$ at $`r_\mathrm{s}`$. We use the jump conditions, eq. (3) through eq. (6), in eq. (14). In eq. (4) we assume $`u_\mathrm{s}=0`$, namely that the shock is temporarily at rest, which should be valid when the shock is marginally stable (or unstable). This is because a stable shock is pushed outwards by the post-shock gas, while cooling reduces the pressure, slows the outward motion, and eventually causes it to halt and then be swept inwards by the infalling matter and gravitational pull. The transition from stability to instability can thus be associated with a transition from expansion to contraction of the shocked volume.

According to eq. (1) and eq. (5) we have

$$e_1=\frac{1}{(\gamma 1)}\frac{P_1}{\rho _1}=\frac{2u_0^2}{(\gamma +1)^2}.$$
(25)
According to eq. (13) and eq. (4) with $`u_\mathrm{s}=0`$ we have

$$\frac{\dot{\rho }}{\rho }=\frac{3u_1}{r_\mathrm{s}}=3\frac{(\gamma 1)}{(\gamma +1)}\frac{u_0}{r_\mathrm{s}}.$$
(26)
With these and eq. (20) we obtain the desired expression for the effective $`\gamma `$ of the post-shock gas in terms of the pre-shock conditions:

$$\gamma _{\mathrm{eff}}=\gamma \frac{(\gamma +1)^4}{6(\gamma 1)^2}\frac{\rho _0r_\mathrm{s}\mathrm{\Lambda }(T_1)}{|u_0|^3}.$$
(27)
For a monoatomic gas, $`\gamma =5/3`$, we obtain

$$\gamma _{\mathrm{eff}}=5/318.96\frac{\rho _0r_\mathrm{s}\mathrm{\Lambda }(T_1)}{|u_0|^3}.$$
(28)
Based on eq. (24), the criterion for stability of a $`\gamma =5/3`$ gas finally becomes

$$\frac{\rho _0r_\mathrm{s}\mathrm{\Lambda }(T_1)}{|u_0|^3}<0.0126.$$
(29)
The post-shock temperature is related to the pre-shock infall velocity using the jump condition, eq. (6), which for $`\gamma =5/3`$ gives

$$T_1=\frac{3}{16}\frac{\mu }{k_bN_a}u_0^2.$$
(30)
For a given cooling function $`\mathrm{\Lambda }(T)`$, eq. (29) is a simple criterion for determining whether a stable shock can form at some radius $`r_\mathrm{s}`$ of the halo. It is in a form that can be directly tested against hydrodynamic simulations (§3), and can serve for evaluating shock stability under realistic conditions in cosmological haloes (§5).

Under the simplifying assumption that the gas is unclumped, the cooling rate is given by eq. (20). The macroscopic cooling function $`\mathrm{\Lambda }(T)`$ is related to the microscopic $`\mathrm{\Lambda }_{\mathrm{mic}}(T)`$, the energy-loss rate of a particle, via $`\mathrm{\Lambda }(T)=(N_a^2\chi ^2/\mu ^2)\mathrm{\Lambda }_{\mathrm{mic}}(T)`$, where $`\chi `$ is the number of electrons per particle. We assume a Helium atomic fraction of 0.1 for $`\mu `$ and $`\chi `$. The microscopic cooling function is shown in Fig. 1 for three different values of mean metallicity $`Z`$. The cooling at temperatures below $`10^4`$K is very slow because the main available cooling agent is molecular hydrogen, which is very inefficient. At temperatures slightly above $`10^4`$K the cooling function peaks due to $`L_\alpha `$ emission from atomic hydrogen. At very low metallicities, a second peak arises near $`10^5`$K due to recombination of atomic Helium. Metals give rise to a higher peak at $`10^5`$K and slightly above, due to line emission from the heavier atoms. At $`10^6`$K and above, the cooling is dominated by brehmstralung, and the cooling function increases slowly. We use the cooling function as derived by Sutherland & Dopita (1993), and presented in their table, in the manner described in Somerville & Primack (1999).

## 3 The spherical hydro code

We test the validity of the shock stability criterion using numerical simulations based on a spherical hydrodynamics code which follows the evolution of shells of dark matter and gas. Since the problem we intend to examine is of global spherical symmetry, and since we need to follow the cooling and the shock with high precision, we use a one-dimensional code. Most of the simulations presented here were run using $`2000`$ gas shells and $`10,000`$ dark-matter shells. A comparable resolution in a three-dimensional code would require on the order of $`10^{10}`$ and $`10^{12}`$ particles respectively, which is impractical. We use no smoothing in the dark-matter shell-crossing scheme, we introduce small-scale smoothing at the halo centre to avoid an artificial singularity there, and we include small artificial viscosity in the hydrodynamics. Tests of the code performance are described in Appendix A.

### 3.1 Dark matter

The dark-matter particles are represented by infinitely thin spherical shells of constant mass and of radii $`r`$ that vary in time. The shell of current radius $`r`$ obeys the equation of motion

$$\frac{d^2r}{dt^2}=\frac{G(M+m)}{(r+a)^2}+\frac{j^2}{r^3},$$
(31)
where $`M`$ and $`m`$ refer respectively to the mass of dark matter and gas within the sphere of radius $`r`$. The last term is a centrifugal acceleration, determined by the the specific angular momentum $`j`$ of the particle represented by the shell. This $`j`$ is assigned to each shell at the initial conditions and is assumed to be preserved during the simulation. The parameter $`a`$ is the smoothing length that becomes effective only near the centre; it has been set to be $`50`$pc throughout this work. The dark-matter shells are allowed to cross each other (and the gas shells). The dark-matter mass is evaluated by

$$M(r)=\underset{r_i<r}{}\mathrm{\Delta }M_i+\frac{1}{2}\underset{j}{}\delta (r=r_j)\mathrm{\Delta }M_j,$$
(32)
where the shell radii and masses are denoted by $`r_i`$ and $`\mathrm{\Delta }M_i`$ respectively, $`i=1,\mathrm{},n_\mathrm{d}`$. The second term adds half the mass of a shell when $`r`$ coincides with one of the shells. Generally, this summation requires $`n_\mathrm{d}^2`$ calculations for the dark matter alone. The particles are kept sorted by radius. When two shells cross each other, we re-sort the array by exchanging pairs which violate the order. This kind of sorting algorithm, termed ‘Shell’s Method’ in Numerical Recipes (Press, 1997), is natural in cases where only a few shells cross each other in each timestep. When two shells cross, they exchange an energy of $`G\mathrm{\Delta }M_i\mathrm{\Delta }M_j/r`$. In order to conserve energy, the radius at which the shells cross must be known with great precision. We therefore reduce the timestep to a small value, $`t_{\mathrm{sc}}`$, when two shells are about to cross each other (see below).

### 3.2 Gas

The hydrodynamic part of the code is based on Lagrangian finite elements in the form of spherical shells. The basic equations governing the dynamics of each shell are

$$\frac{d^2r}{dt^2}=\frac{1}{\rho }(P+\sigma )\frac{G(M+m)}{(r+a)^2}+\frac{j^2}{r^3},$$
(33)
$$\frac{de}{dt}=\rho ^2(P+\sigma )\frac{d\rho }{dt}q,$$
(34)
$$\rho =\frac{dm}{4\pi r^2dr}.$$
(35)
$$P=(\gamma 1)e\rho ,$$
(36)
An artificial viscosity term, $`\sigma `$, is added to the pressure for numerical purposes, as explained below. The smoothing length effective at the center, $`a`$, is the same as for the dark matter, eq. (31). As in the model described before, the loss of internal energy due to radiative cooling is represented by the cooling rate $`q`$.

The gas is divided into discrete shells. The mass enclosed within a shell, $`\mathrm{\Delta }m`$, is assumed constant in time, while the inner and outer shell boundaries move independently in time. Each boundary is characterized by a temporal position $`r`$, velocity $`v`$ and specific angular momentum $`j`$. The acceleration \[eq. (33)\] is evaluated at the boundary position. The variables $`\rho `$, $`P`$, $`e`$, $`q`$, and $`T`$ for each shell are evaluated within the shell between the boundaries.

In particular, the pressure term in eq. (33) is evaluated at the outer boundary $`r_i`$ of shell $`i`$ using eq. (35):

$$\frac{1}{\rho }(P+\sigma )=\frac{4\pi r_i^2}{\mathrm{\Delta }m}[(P+\sigma )_{i+1}(P+\sigma )_i],$$
(37)
where $`\mathrm{\Delta }m=(\mathrm{\Delta }m_i+\mathrm{\Delta }m_{i+1})/2`$.

The boundary conditions for the outer boundary of the system are $`P=\sigma =0`$, and zero mass beyond the outer boundary.

Since gas shells cannot cross each other, the gas mass in the sphere interior to each gas shell is constant throughout the simulation:

$$m(r_i)=\underset{j=1}{\overset{i}{}}\mathrm{\Delta }m_j.$$
(38)
For the evolution of the dark-matter shell at $`r`$, we evaluate the gas mass that appears in eq. (31) using

$$m(r)=\underset{j=1}{\overset{i1}{}}\mathrm{\Delta }m_j+\frac{r^3r_{i1}^3}{r_i^3r_{i1}^3}\mathrm{\Delta }m_i,$$
(39)
where $`i`$ refers to the gas shell for which $`r_{i1}r<r_i`$.

### 3.3 Integration and Timestep

The discrete integration of $`r`$ and $`v`$ is performed by a Runge Kutta fourth-order scheme (Press, 1997). The state of the system at the beginning of each timestep is kept in memory until the timestep is completed, such that it is possible to return to the beginning of the timestep and retry with a smaller timestep if the convergence criteria are not met. The timesteps are set such that the position $`r`$ and velocity $`v`$ do not change by too much during a single timestep. For a given accuracy parameter $`ϵ_{\mathrm{rk}}`$, we demand that the difference between the forth-order displacement $`\mathrm{\Delta }r_4`$ and the analogous first-order displacement $`\mathrm{\Delta }r_1`$ obeys $`|\mathrm{\Delta }r_4\mathrm{\Delta }r_1|/r<ϵ_{\mathrm{rk}}`$, both for the dark-matter and the gas. The similar requirement is applied to the change in velocity over a timestep. If this condition is not fulfilled, we reduce the timestep by a certain factor and repeat the calculation over this timestep. We use here as our default $`ϵ_{\mathrm{rk}}=0.1`$.

In addition, we make sure the timestep for each shell does not violate the Courant condition, for an accuracy parameter $`ϵ_\mathrm{c}`$. This implies $`c_\mathrm{s}\mathrm{\Delta }t/\mathrm{\Delta }r<ϵ_\mathrm{c}`$, where $`c_\mathrm{s}^2=(dP/d\rho )_s=\gamma P/\rho `$ is the speed of sound. We use here as our default $`ϵ_\mathrm{c}=0.3`$.

A third limitation on the timestep comes from the desire to conserve energy when shells cross. When two shells are about to cross each other within the current timestep $`dt`$, we set the timestep to $`\mathrm{min}(dt,t_{\mathrm{sc}})`$, and keep it small until they actually cross. We use here as our default $`t_{sc}=10^4`$Gyr.

The values for $`ϵ_\mathrm{c},ϵ_{\mathrm{rk}}`$ and $`t_{sc}`$ were chosen empirically such that energy is conserved and the dynamics converges to our satisfaction, in the sense that it does not change by much when smaller parameters are used. We demonstrate in Appendix A how well these requirements are met.

Once we have computed the new radii and velocities of the shells at the end of the timestep, we correct the energy of the gas for the $`PdV`$ work term using the states of the system at the beginning and at the end of the timestep. The cooling is explicitly subtracted from the internal energy after the hydrodynamic timestep is completed. Once the final state of the system is ready, it is copied onto the memory array of the initial state, and the simulation is ready to execute a new timestep.

### 3.4 Initial conditions

The simulation starts at high redshift, $`z=100`$, with a small spherical density perturbation. The initial density fluctuation profile is set to be proportional to the linear correlation function of the assumed cosmological model, representing the typical perturbation under the assumption that the random fluctuation field is Gaussian (see Dekel, 1981, and Appendix C). The amplitude of the density fluctuation at the initial time, averaged over a given mass, determines the time of collapse, as desired. The initial velocity field is assumed to follow a quiet Hubble flow and the radial peculiar velocities build up in time. We assume the standard $`\mathrm{\Lambda }`$CDM cosmology with $`\mathrm{\Omega }_\mathrm{m}=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`h=0.7`$ and $`\sigma _8=1`$.

### 3.5 Angular momentum

We assume that in a real system the orbits of dark-matter particles, and the initial orbits of the gas particles, are quite elongated. Cosmological N-body simulations show that the velocity distribution tends to be more radial than tangential (Ghigna et al. , 1998; Safran & Dekel, 2003) and already for an isotropic distribution the eccentricities are about 1:6. The processes we study in this paper occur away from the galactic disc at a radius on the order of the virial radius, namely in a regime where the centrifugal force can be expected to be negligible compared to the gravitational force and the gas pressure force. The prescribed angular momentum for the shells is thus mainly for numerical purposes, to avoid divergent densities of gas or dark matter shells when they pass through the halo centre. Our results concerning the virial shock are insensitive to the actual way by which we assign angular momentum to each shells.

In the current study we practically assume that the dark-matter particles are almost on radial orbits. The angular momentum of the gas is prescribed such that the shells, once they lost their energy by radiation, would settle into an exponential disk with pure circular motions and a characteristic radius of a few kpc, smaller than the inner characteristic radius of the halo. Our spherical ‘disc’ thus contains gas that is cold and dense compared to the shocked gas.

### 3.6 Artificial viscosity

It is impossible to follow the discontinuous behavior across the shock using the conventional continuity equation for the density and standard conservation of energy and momentum. The jump conditions can be calculated explicitly, as in eq. (3) to (6) (termed ‘the characteristic method’ or ‘Godonov’s method’). Alternatively, as proposed by Von-Newman, one can slightly smear the discontinuities and then solve them within the framework of the standard hydrodynamic equations. By adding an artificial pressure term in a few shells around the shock, the differential equations become solvable and one can continue the calculation without affecting the energy and the dynamics of the shock (while its internal structure naturally changes).

Artificial viscosity is applied when the inner and outer shell boundaries at $`r_1`$ and $`r_2`$ approach each other, $`\mathrm{\Delta }v=v_2v_1<0`$, and when the volume of the shell decreases, $`dV/dt=4\pi (r_2^2v_2r_1^2v_1)<0`$. The artificial viscosity then takes the form

$$\sigma =a_2\rho (\mathrm{}v)^2+a_1\rho c_s|\mathrm{\Delta }v|.$$
(40)
The quadratic, common form of artificial viscosity smears discontinuities over about 3 shells. The Linear discontinuity affects a slightly larger range, and is usually added with a smaller coefficient $`a_1`$. The coefficients $`a_1`$ and $`a_2`$ are varied for different shells in the course of the simulation in order to overcome a specific numerical problem in the cold ‘disc’, where the gravitational and centrifugal forces balance each other and the pressure force is negligible. In this case the gas is not a standard hydrodynamic gas because the pressure does not regulate large discontinuities, and information is not transported because of the low sound speed. When a ‘disc’ shell vibrates, it is artificially heated by the artificial viscosity in every contraction until its pressure grows and stops the process. If we are not careful to properly tune the artificial viscosity we may end up with one ‘disc’ shell that has been heated to $`10^7`$K while the rest of the ‘disc’ is at $`10^4`$K. This imposes an undesired drastic decrease in the corresponding timestep. In order to overcome this numerical problem, we gradually turn off the quadratic term ($`a_2`$) of the artificial viscosity inside the ‘disc’. We define a ‘disc’ radius $`R_{\mathrm{disc}}`$ to be the largest radius for which the difference between the gravitational and centrifugal forces is less than 1/4 of the gravitational force. Once at $`r<0.6R_{\mathrm{disc}}`$, we continuously decrease the parameter $`a_2`$ in eq. (40) according to $`a_2=(r0.3R_{\mathrm{disc}})/(0.3R_{\mathrm{disc}})`$ and make it completely vanish at $`r<0.3R_{\mathrm{disc}}`$. This prescription was found by trial-and-error to properly solve the numerical problem in most cases. The linear term of the viscosity, being proportional to the speed of sound, is anyway very small in the cold ‘disc’, so effectively no artificial viscosity is applied in the inner ‘disc’.

Appendix A provides tests and examples of the hydrodynamic simulations in some detail.

## 4 Virial shock in the simulations

### 4.1 Existence of a virial shock

We now investigate the formation of a virial shock using the spherical hydrodynamical simulations described above. We wish to test in particular the validity of the analytic stability criterion developed in §2.

In order to mimic a typical perturbation in a random Gaussian field (Dekel 1981, Appendix B), the initial density-fluctuation profile was set to be proportional to the correlation function, normalized such that the mean density fluctuation in a sphere enclosing $`M_\mathrm{i}=10^{11}M_{}`$ was $`\delta _\mathrm{i}=0.09`$ at $`z=100`$. For example, the shell encompassing $`M3\times 10^{10}M_{}`$ is expected to collapse at $`z=3`$, and $`M10^{12}M_{}`$ is expected to collapse at $`z=0`$.

Fig. 2 shows the time evolution of the radii of Lagrangian shells in a simulation of the adiabatic case, with the cooling turned off. We find that a shock exists at all times. It appears as a sharp break in the flow lines, associated with a discontinuous decrease in infall velocity \[eq. (4)\]. Shown in the figure is the shock radius, defined by the outermost shell for which the inner and outer shell boundaries approach each other and the volume of the shell decreases \[the same conditions that have been used for turning on the artificial viscosity in eq. (40)\]. The shock gradually propagates outwards, encompassing more gas mass and dark matter in time. The gas below the shock is pressure supported and at quasi-static equilibrium. Not shown here are the dark-matter shells, which collapse, oscillate and tend to increase the gravitational attraction exerted on the gas shells.

Shown in comparison is the evolution of the virial radius, computed from the simulation density as the radius within which the mean overdensity is $`\mathrm{\Delta }_\mathrm{v}`$ times the mean cosmological background density. The virial overdensity $`\mathrm{\Delta }_\mathrm{v}`$ is provided by the dissipationless spherical top-hat collapse model; it is a function of the cosmological model, and it may vary with time. For the Einstein-deSitter cosmology, the familiar value is $`\mathrm{\Delta }_\mathrm{v}176`$ at all times.222This can be derived from the top-hat formalism of Appendix B, once the final radius is assumed to be fixed at half the maximum-expansion radius but the overdensity is evaluated at the time when the top-hat sphere would have collapsed to a singularity. For the family of flat cosmologies ($`\mathrm{\Omega }_\mathrm{m}+\mathrm{\Omega }_\mathrm{\Lambda }=1`$), the value of $`\mathrm{\Delta }_\mathrm{v}`$ can be approximated by Bryan & Norman (1998)

$$\mathrm{\Delta }_\mathrm{v}(18\pi ^2+82x39x^2)/(1+x),$$
(41)
where $`x\mathrm{\Omega }_\mathrm{m}(z)1`$, and $`\mathrm{\Omega }_\mathrm{m}(z)`$ is the ratio of mean matter density to critical density at redshift $`z`$. For example, in the $`\mathrm{\Lambda }`$CDM cosmological model that serves as the basis for our analysis in this paper ($`\mathrm{\Omega }_\mathrm{m}=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$), the value at $`z=0`$ is $`\mathrm{\Delta }_\mathrm{v}340`$. We see in Fig. 2 that the shock radius almost coincides with the virial radius at all times. This is hardly surprising, as the shock is likely to appear at the outermost radius at which shell crossing first occurs, which is near the virial radius (to be demonstrated in Fig. 8 below).

Fig. 3 is the result of a similar simulation, but now with realistic radiative cooling for $`Z=0`$. We see that a stable shock does not exist in this case before $`t=3.9`$Gyr. During this period, the cooling makes the gas lose its pressure support and lets it collapse freely under gravity into the halo centre. The collapse stops by the assumed angular momentum, in a ‘disc’ whose marked radius can be identified at the bottom of the plot by the abrupt change of the infalling flow lines into horizontal lines. The matter in the ’disc’ is angular-momentum supported. As is visible in the figure, a shock, in the sense of a discontinuity in velocity and density, is present at the edge of the disc. Once the stability criterion is met, a shock forms and propagates outward abruptly. The propagation of the shock causes it to re-enter a regime for which $`\gamma _{\mathrm{eff}}`$ is below the critical value. Consequently, the post shock gas becomes non-supportive again and falls. The oscillatory behavior of the shock continues with an increasing period until it stabilizes at the largest radius for which the stability criterion is met. The shock never expands beyond the virial radius because shells do not tend to cross there (Fig. 8 below).

Fig. 4 shows the evolution of the total mass interior to the characteristic radii in the adiabatic and radiative simulations of Fig. 2 and Fig. 3. In the adiabatic case, the shock mass practically coincides with the virial mass, and the gas never forms a disc. In the radiative case, the shock is initially at the disc radius, and only after 3.9Gyr does it start propagating outwards.

### 4.2 Testing the stability criterion

We first test the assumption made in §2.2 that the post-shock gas has a homologous velocity profile in the vicinity of the shock. In Figures 2 and 3 we see that the flow lines beneath the shock tend to be nearly parallel in the $`\mathrm{log}rt`$ plane. This means that $`d\mathrm{log}r/dt=u/r`$const., namely homology, eq. (12). A similar behavior has been found in all the simulations that we have performed.

Next, we use the simulations to test the validity of the stability criterion derived in §2, eq. (29), where $`\gamma _{\mathrm{eff}}`$ is given by eq. (27). Recall that the $`\gamma _{\mathrm{eff}}`$ is expressed in terms of the pre-shock quantities. In order to map the value of $`\gamma _{\mathrm{eff}}`$ in the different regions of the free falling gas in the forming halo, we ran the same simulation as in the previous section except that the cooling rate was set to be very high, such that the virial shock never develops. Fig. 5, top panel shows the flow lines in this case, on which overlaid are 4 contours of equal $`\gamma _{\mathrm{eff}}`$ values, evaluated via eq. (27) with eq. (30) and the cooling rate from Sutherland & Dopita (1993). As shells are falling into the halo, their $`\gamma _{\mathrm{eff}}`$ is gradually increasing. Also, as time progresses, the value of $`\gamma _{\mathrm{eff}}`$ at the same radius is increasing. By following the value of $`\gamma _{\mathrm{eff}}`$ just above the “disc” radius (shown as the break at the bottom of the plot), in comparison with the critical value for stability $`\gamma _{\mathrm{crit}}`$, we can therefore use our model to predict when we expect the virial shock to form. This is shown in the middle panel of Fig. 5. We see that at early times (and smaller masses) we have $`\gamma _{\mathrm{eff}}<\gamma _{\mathrm{crit}}`$, predicting no stable shock. The system is predicted to enter the stable-shock regime at about $`t=3.9`$Gyr, where $`\gamma _{\mathrm{eff}}`$ becomes larger than $`\gamma _{\mathrm{crit}}`$. A comparison with the realistic radiative simulation described in Fig. 3 yields that this model prediction is very accurate: the shock indeed starts forming at $`t=3.9`$Gyr, and is globally stable thereafter.

Fig. 6 shows the actual evolution of $`\gamma _{\mathrm{eff}}`$ at either the ‘disc’ radius or the shock radius, whichever is larger, as computed directly from the pre-shock quantities in the simulation with realistic cooling. Shown at the same times are the characteristic masses, already shown in Fig. 4. The virial shock is first generated at time $`t=3.9`$Gyr, when the virial mass is $`10^{11}M_{}`$. Starting at this time, the shock is propagating outwards very rapidly. As a result of this fast expansion, the $`\gamma _{\mathrm{eff}}`$ of the pre-shock infalling matter at the shock, which is decreasing with $`r`$, drops below the threshold. This makes the shock lose its pressure support, it becomes temporarily unstable and its expansion slows down until it is eventually swept back on a dynamic time scale. The associated drop in total mass behind the shock, seen around $`t=5.8`$Gyr, is due to the fact that the dark matter is not swept back with the gas. Once the shock is shrunk to a low enough radius, $`\gamma _{\mathrm{eff}}`$ rises again to above $`\gamma _{\mathrm{crit}}`$; the shock becomes stable again and it resumes its associated expansion towards the virial radius. After the conditions for the shock stability are first met, the shock is visible most of the time. In the rest of this paper, we treat haloes at this state as ones containing a virial shock.

Fig. 7 presents results similar to Fig. 6 for two other simulations with different initial overdensities, and therefore different masses collapsing at different times. The small difference seen in one case between the $`\gamma _{\mathrm{eff}}`$ at which the shock actually forms and the predicted $`\gamma _{\mathrm{crit}}`$ may be totally due to numerical inaccuracies in the simulation. Such inaccuracies may occur when a dark-matter shell crosses a gaseous shell, which, near the threshold, may lead to a slightly premature shock formation. This is seen in the convergence test of our code described in Appendix A, when we compare the shock formation times in table 2. Thus, the model stability criterion is found to be valid within the accuracy of the simulations in all the cases studied, indicating that the model is not limited to a special range of masses and collapse times.

## 5 Shock stability in cosmology

The analysis of §2 thus provides a successful criterion for shock stability, eq. (29), as a function of the pre-shock properties of the infalling gas at radius $`r`$: the density, velocity and metallicity. In order to apply this criterion to a given protogalaxy in a cosmological background, we wish to evaluate the gas density and velocity just before it hits the disc, for a gas shell initially encompassing a total mass $`M`$ that virializes at redshift $`z_\mathrm{v}`$. In this calculation we assume an Einstein-deSitter cosmology, as a sensible approximation at $`z>2`$ (where $`\mathrm{\Omega }_\mathrm{m}>0.9`$).

We assume a given universal baryonic fraction $`f_\mathrm{b}`$, and a global spin parameter $`\lambda `$ which determines the ratio of disc to virial radius. The initial mean density perturbation profile, $`\overline{\delta }_\mathrm{i}(M)`$, is given at some fiducial time in the linear regime; it is the average profile derived from the power spectrum of initial density fluctuations, as described in Appendix C. In the cosmological toy model used here we approximate the power spectrum as a power law, $`P_kk^n`$, where $`n2.4`$ to mimic the $`\mathrm{\Lambda }`$CDM power spectrum on galactic scales.

We follow gas shells from the initial perturbation till they approach the disc using a two-stage model. During the expansion, turn-around, and until an assumed virialization at half the maximum-expansion radius, we assume no shell crossing, the total mass interior to the shell remains constant in time, and we follow the radii, density and velocity of the shell via the spherical top-hat model (see Appendix B). From the virial radius inwards we assume that the gas shells, which do not cross each other, contract inside the fixed potential well of an isothermal dark-matter halo. This idealized model involves several crude approximations, such as the instantaneous transition at the virial radius, and neglecting the effect of the angular momenta of the individual gas particles at small radii, but we show using spherical simulations that this model predicts the minimum halo mass for which a stable shock first appears to an accuracy better than 25%. This allows us to use the model for exploring the critical mass as a function of cosmological parameters such as galaxy formation time, metallicity, spin parameter, fluctuation power spectrum, and baryonic fraction.

### 5.1 Toy model until virialization

For a given shell $`M`$ and initial mean perturbation profile $`\delta (M)`$ \[standing for the $`\overline{\delta }_\mathrm{i}(M)`$ of Appendix B\], the top-hat model \[eq. (70) and eq. (64)\] yields the implicit solution

$$r(M,\eta )=r_\mathrm{v}(M)(1\mathrm{cos}\eta ),$$
(42)
$$t(M,\eta )=t_\mathrm{v}(M)(\eta \mathrm{sin}\eta ),$$
(43)
where the mass dependence enters via the virial quantities

$$r_\mathrm{v}=C_\mathrm{i}M^{1/3}\delta (M)^1,$$
(44)
$$t_\mathrm{v}=r_\mathrm{v}/v_\mathrm{v}=C_\mathrm{i}^{3/2}G^{1/2}\delta (M)^{3/2},$$
(45)
with $`v_\mathrm{v}^2=GM/r_\mathrm{v}`$. The coefficient $`C_\mathrm{i}=(6/\pi )^{1/3}(0.15/\rho _{\mathrm{ui}})`$ is determined by $`\rho _{\mathrm{ui}}`$, the cosmological density at the initial time when $`\delta (M)`$ is given, independent of $`M`$.

The velocity of the shell $`M`$ is

$$u=\frac{r/\eta }{t/\eta }=v_\mathrm{v}\frac{\mathrm{sin}\eta }{(1\mathrm{cos}\eta )}.$$
(46)
At virialization, $`\eta =3\pi /2`$, it is simply $`u=v_\mathrm{v}`$.

In order to evaluate the local density, we follow the radii of two adjacent shells, encompassing masses $`M`$ and $`M+dM`$ respectively, at a given time $`t`$, e.g., the time when shell $`M`$ virializes (at half its maximum expansion radius). Let $`\eta `$ correspond to shell $`M`$ at that time, and $`\eta +d\eta `$ to shell $`M+dM`$. In order to express $`d\eta `$ in terms of $`dM`$ we use the fact that the time $`t`$ is the same for the two shells: $`0=dt=(t/M)dM+(t/\eta )d\eta `$. Using eq. (43) this gives

$$d\eta =\frac{3}{2}\frac{(\eta \mathrm{sin}\eta )}{(1\mathrm{cos}\eta )}\frac{\delta ^{}}{\delta }dM,$$
(47)
where we denote $`\delta ^{}d\delta /dM`$. Expressing $`dr`$ in terms of $`dM`$ and $`d\eta `$ based on eq. (42), we obtain using eq. (47) and after some algebra

$$\frac{dr}{r}=\frac{1}{3}\frac{dM}{M}\left(1\frac{3M\delta ^{}}{\delta }\left[1\frac{3}{2}\frac{\mathrm{sin}\eta (\eta \mathrm{sin}\eta )}{(1\mathrm{cos}\eta )^2}\right]\right).$$
(48)
At virialization of shell $`M`$, $`\eta =3\pi /2`$, the quantity in square brackets equals $`(10+9\pi )/4`$. Not surprisingly, if the initial perturbation is of uniform density, $`\delta ^{}=0`$, we are left with $`dr/r=(1/3)(dM/M)`$, the straightforward result of $`Mr^3`$. Recall that the virial radius $`r_\mathrm{v}`$ of shell $`M`$ can be obtained either from the universal density at the time of virialization using eq. (66), or from the initial perturbation using eq. (70). The desired local density $`\rho `$ can be obtained from eq. (48) via $`dM=4\pi r^2\rho dr`$.

If the initial perturbation profile is a power law, $`\delta (r)r^{(n+3)}`$, using $`Mr^3`$ we have $`3M\delta ^{}/\delta =(n+3)`$. So finally

$$\frac{dr}{r_\mathrm{v}}=\frac{1}{3}\frac{dM}{M}\left(1+(n+3)\frac{(10+9\pi )}{4}\right).$$
(49)
eq. (48) \[or eq. (49) in the power-law case\] allows us to compute the desired radii of the two adjacent shells at the time of virialization of shell $`M`$.

### 5.2 Toy model after virialization

Given the radii of the two adjacent shells at $`r_\mathrm{v}`$, we enter the shell-crossing regime and continue to follow the shells down to the disc radius by numerical integration. The shell radius $`r`$ and velocity $`u`$ are related via energy conservation. Assuming that the gas shells contract without crossing each other inside a dark-matter halo that is a fixed isothermal sphere, the total mass interior to the shell that originally encompassed a total mass $`M`$ is

$$M(r)=f_\mathrm{b}M+(1f_\mathrm{b})\frac{M}{r_\mathrm{v}}r,$$
(50)
and the gravitational potential at $`r`$ is

$$\varphi (r)=\frac{GM(r)}{r}(1f_\mathrm{b})\frac{GM}{r_\mathrm{v}}\mathrm{ln}\left(\frac{r_\mathrm{v}}{r}\right).$$
(51)
The integration is performed by advancing $`r`$ according to the velocity $`u`$ and then recalculating $`u`$ according to energy conservation:

$$(1/2)u^2+\varphi (r)=(1/2)v_\mathrm{v}^2+\varphi (r_\mathrm{v})=\mathrm{const}.$$
(52)
We follow shell $`M`$ for the time it falls from $`r=r_\mathrm{v}`$ to the disc radius $`r=\lambda r_\mathrm{v}`$, and shell $`M+dM`$ for the same time interval. Denoting the separation between the shells at the end of this time interval by $`dr`$, we compute the desired gas density by

$$\rho =\frac{f_\mathrm{b}}{4\pi r^2}\frac{dM}{dr}.$$
(53)
The resultant values of $`r`$, $`u`$ and $`\rho `$ are inserted into eq. (28) in order to obtain an approximation for $`\gamma _{\mathrm{eff}}`$ and then to evaluate stability by eq. (24). This allows us to check stability for the case where mass $`M`$ virializes at redshift $`z_\mathrm{v}`$, with metallicity $`Z`$, spin parameter $`\lambda `$, baryonic fraction $`f_\mathrm{b}`$, and a given power spectrum.

### 5.3 Model versus simulations

In Fig. 8 and Fig. 9 we compare the evolution of the quantities of a given gas shell according to the toy model described in the previous subsections and according to the spherical hydro simulation described in the earlier sections. We follow a specific shell that hits the disc at about $`t3.8`$Gyr, just before the shock starts propagating into the halo \[see Fig. 3\]. The quantities shown as a function of radius $`r`$ are total mass $`M`$ interior to $`r`$, radial velocity $`u`$, gas density $`\rho `$, and the corresponding value of $`\gamma _{\mathrm{eff}}`$. For $`M`$ and $`u`$ the evolution starts at the top-left corner and ends at the bottom left, while for $`\rho `$ the upper part of the curve corresponds to the expansion phase and the lower part to the contraction phase. The evolution of $`\gamma `$ is followed only during part of the contraction phase.

In Fig. 8 we calibrate the toy model to match the simulation at the maximum expansion radius. We see that while the mass interior to the shell is reproduced by the model only to a limited accuracy in the last stages of the collapse, the velocity, density and the resulting value of $`\gamma _{\mathrm{eff}}`$ are recovered very well by the model. This allows us to predict quite accurately the point where $`\gamma _{\mathrm{eff}}`$ exceeds $`\gamma _{\mathrm{crit}}`$.

Since we wish to use the toy model without an exact knowledge of the conditions at maximum expansion, we normalize the model evolution in Fig. 9 based on $`M`$ and $`z_\mathrm{v}`$.

The slight deviations in $`u`$ and in $`\rho `$ now translate into a larger error in $`\gamma _{\mathrm{eff}}`$. The error in the toy model originates mostly from the slight ambiguity in the definition of the virial radius. On one hand we assume it to equal half the maximum-expansion radius, and on the other hand we assume it to represent an overdensity of $`\mathrm{\Delta }_\mathrm{v}`$ as in eq. (41). These two assumptions are not fully consistent with the actual behavior of the virializing system in the simulation. Nevertheless, we see below that our approximate model allows us to estimate the critical halo mass below which the shock does not form to an accuracy of better than 25%, which is quite satisfactory for our purpose here.

### 5.4 Critical mass for shock formation

Fig. 10 shows for several different cases the critical halo mass, below which a shock does not propagate into the halo, versus the redshift at which this critical mass virializes. For each case we compare the model prediction to the shock formation as actually seen in the simulation. The cases differ by the mean metallicity, $`Z=0`$ and $`Z=0.05`$ for the lower and upper sets of points respectively, and by the amplitude of the initial perturbation, corresponding to a range of shock-formation redshifts at every given $`Z`$. The assumed baryonic fraction is always $`f_\mathrm{b}=0.13`$, but the assumed spin parameter may be different for the different shells in a given simulation because we set it for each shell such that the final disc has an exponential surface density profile. However, the $`\lambda `$ values vary in the range $`0.02`$ to $`0.05`$, compatible with the distribution of spin parameter in cosmological simulations (Bullock et al. 2001). We see that the model predicts the critical mass with an accuracy better than 25%, such that we can use it for mapping the parameter space in more detail.

Fig. 11 shows for several different choices of parameters the critical halo mass for shock formation versus the halo virialization redshift as predicted by the model. A virial shock does not form in haloes of masses below the line. The lines are not always monotonic due to the non-monotonic features in the cooling curves (Fig. 1). Shown in comparison is $`M_{}`$, the characteristic mass for haloes forming at $`z`$ according to the Press-Schechter approximation (Lacey & Cole, 1993). The default values of the parameters, used unless specified otherwise, are $`Z=0`$, $`\lambda =0.05`$, $`f_\mathrm{b}=0.13`$ and $`n=2.4`$.

The upper panel has the metallicity varying from $`Z=0`$ to $`Z=0.3`$. The critical mass tends to be higher at higher redshifts (especially for $`Z=0`$) because the higher density implies more efficient cooling. It is striking that even for the case of zero metallicity, for which the cooling is not at its maximum efficiency, an $`M_{}`$ halo cannot produce a shock until a relative late redshift, $`z2.1`$. The addition of a small amount of metals, $`Z=0.05`$, increases the cooling rate significantly (see Fig. 1) such that $`M_{}`$ haloes start producing virial shocks only after $`z1.6`$.

The second panel has $`\lambda `$ varying as marked. The shock forms slightly earlier if the disc is smaller (lower $`\lambda `$), because the conditions become more favorable for shock formation closer to the centre. At high redshifts the increase in infall velocity happens to balance out the increase in density, temperature and cooling rate. The post shock temperature there is a few $`10^6`$K. The bottom line is that the critical mass is not too sensitive to $`\lambda `$.

The third panel has $`f_\mathrm{b}`$ varying as marked. The critical mass is monotonic with the baryonic fraction because the cooling rate is monotonic with gas density. The parameter $`f_\mathrm{b}`$ can be interpreted as the fraction of the baryons that actually take part in the shock formation. This can be smaller than the universal baryonic fraction if some of the gas falls into the halo in the form of dense clumps. Even with $`f_\mathrm{b}`$ as low as $`0.05`$, meaning that most of the gas is not participating in the cooling, an $`M_{}`$ halo would not produce a shock until $`z2.4`$. The conclusion is that the critical mass is not too sensitive to $`f_\mathrm{b}`$ either.

The bottom panel explores three values for the initial power index $`n`$ approximating the power spectrum of $`\mathrm{\Lambda }`$CDM on galactic scales. The dependence of the critical mass on $`n`$ in this regime is weak.

## 6 Discussion

The heating of the gas behind a virial shock in haloes has been a basic component in galaxy formation theory (Rees & Ostriker, 1977). We studied the conditions for the existence of such a virial shock in spherical haloes. We first pursued an analytic stability analysis in the presence of cooling, and then demonstrated its validity using high-resolution spherical hydrodynamical simulations. The obtained criterion for shock stability in terms of the post-shock quantities is

$$\gamma _{\mathrm{eff}}\frac{d\mathrm{ln}P/dt}{d\mathrm{ln}\rho /dt}>10/7.$$
(54)
In terms of the pre-shock gas properties, this condition reads

$$\frac{\rho _0r_\mathrm{s}\mathrm{\Lambda }(T_1)}{|u_0|^3}<0.0126,$$
(55)
where $`\rho _0`$ and $`u_0`$ are the gas density and infall velocity in front of the shock, $`r_\mathrm{s}`$ is the shock radius, $`\mathrm{\Lambda }(T)`$ is the cooling function which depends on the metallicity $`Z`$, and $`T_1u_0^2`$ is the post-shock temperature as a function of the pre-shock infall velocity.

Based on this criterion, we find that a virial shock forms only in big haloes forming at late redshifts. A virial shock does not form in smaller haloes forming early where the cooling behind the shock efficiently removes its pressure support. For example, we find that most galactic haloes that have collapsed and virialized by $`z2`$ did not produce a virial shock. Haloes less massive than $`10^{11}M_{}`$ never produce a shock even if the gas is of zero metallicity. If the metallicity is non-negligible (e.g. $`Z0.05`$), this lower bound to shock formation rises to $`710^{11}M_{}`$. When a shock does not exist, the gas is not heated to the halo virial temperature until it falls all the way to the disc at the inner halo.

Forcada-Miro & White (1997), in an unpublished work, have pursued independently a numerical analysis along similar lines, involving a more detailed treatment of the cooling processes involved. They also find that the virial shock radius is significantly reduced due to the cooling in haloes of small masses, $`M<10^{11}M_{}`$. In their case the shock never completely disappears because of a different feature in their numerical scheme; they put all the cooled post-shock gas in one central “shell” to avoid numerical difficulties at the centre. This makes the inner boundary of the system follow the shock quite closely in cases where there is efficient cooling behind the shock, and allows the presence of a small-radius shock even in such cases. Overall, our numerical results are in encouraging agreement, and our analytic model provides a natural explanation for their numerical results as well.

The most severe uncertainty when attempting to apply our results to real galaxies arises from the assumed spherical symmetry in both the model and the numerical simulations. The validity of this approximation for the asymmetric halo configurations in the hierarchical clustering scenario is an open question to be addressed in future work. Nevertheless, we notice that Katz et al.  (2001) and Fardal et al.  (2001) observe in their cosmological simulations that a large fraction of the mass accreted onto haloes indeed remains cold and is never heated to the virial temperature. Toft et al.  (2002) find in their simulations, using a similar treeSPH code as the one used by Katz et al.  (2001), that the soft X-ray radiation is mainly emitted from within the innermost $`20kpc`$ of their haloes, well inside the virial radius, in encouraging agreement with our results. On the other hand, it is not obvious that the resolution in these simulations is adequate for studying the shock physics involved; our estimates indicate that three-dimensional simulations with proper resolution are not practical at present (Appendix A).

Another complication may arise from radiative effects. Even when there is no virial shock, the kinetic energy of the gas eventually turns into radiation when the gas infall motion is brought to a halt at the disc. At such densities, the width of the shock front is much smaller than the width of the cooling front behind the shock, $`10^2`$pc versus $`10^2`$pc. Thus, the gas in a thin shell behind the shock is heated to a temperature corresponding to its kinetic energy, and it cools by radiating soft X-rays. The X-ray radiation is expected to generate an ionized $`H_{\mathrm{𝐼𝐼}}`$ bubble, in which the ionization rate balances the recombination rate. The Strömgren radius of this bubble is relatively small, on the order of a few kiloparsecs, because the high gas density implies a high recombination rate. The recombination process then generates a flux of $`L_\alpha `$ radiation, emitted at the inner few kpc of the halo.

A naive inspection of cross sections might indicate that the $`L_\alpha `$ radiation would be trapped inside the halo. This could in principle affect the shock stability in three different ways: by increasing the radiation pressure, by heating up the infalling matter, and by slowing down the radiative cooling responsible for the shock instability. It has been argued by Rees & Ostriker (1977) that the radiation pressure at these low temperatures must be insignificant compared to the gas pressure even if all the internal energy was drained from the baryons into the radiation field. One might add that since the radiation pressure behaves like a $`\gamma =4/3`$ gas, it could at most make the system marginally stable. When work is done on the radiation field, any leakage of radiation out would turn the energy into cooling rather than $`PdV`$, and will thus reduce the effective gamma, making the system unstable.

Partial heating of the infalling gas should not affect our analysis as long as the temperature of the infalling matter is significantly below the virial (post-shock) temperature such that the strong shock approximation remains valid. The effect of the reduced cooling rate is yet to be investigated. In practice, we do not expect the radiation trapping to be very efficient, because the effective opacity is reduced by thermal broadening and by the systematic blue shift due to the gas infall motion. When the opacity is high, the radiation heats up the gas, which enhances the thermal broadening and the collisional ionization rate. This reduces the opacity and allows for radiation escape. The system is likely to reach a steady state in which it gradually cools. This process is under current investigation.

Feedback effects may further complicate the picture and affect shock stability. The energy fed back to the gas from stars, supernovae and AGNs may heat the halo gas and expel part of it. Merging substructures may have additional complicated effects. These effects cannot be captured by our idealized spherical analysis, and a proper study would probably require high-resolution three-dimensional hydrodynamical simulations. While observations and certain theoretical considerations indicate that feedback effects are likely to be important in galaxies as large as $`10^{11}M_{}`$ and may thus affect the shock stability (e.g. Dekel & Silk 1986; Dekel & Woo 2003, and references therein), it has proven difficult for numerical simulations to reproduce such effects in any but small dwarf galaxies, indicating that feedback effects may not be so crucial for the understanding of shock stability. Until the dust settles on the role of feedback effects, our preliminary conclusions based on the spherical analysis should be taken with a grain of salt.

The general absence of a virial shock might have three direct implications, which we study in associated papers.

First, as explained above, when the gas is heated at the disc rather than near the halo virial radius, the generated X-ray radiation serves to ionize the gas and is not emitted outwards. The result would be a suppression of the X-ray emission in the range $`5\times 10^5`$ to $`2\times 10^6`$K. This may help explaining the missing X-ray problem pointed out by Pen (1999) and Benson et al.  (2000). Pen (1999) argue that there is an order-of-magnitude discrepancy between the soft X-ray flux as observed by Cui et al.  (1996), after subtracting the contribution of quasars, and the predicted flux from haloes constructed by a Press-Schechter hierarchical model under the assumption of shock heating to the virial temperature.

Second, the infall energy, via the ionizing X-ray, is efficiently transformed into $`L_\alpha `$ radiation at the inner few kpc of the halo. A related increase in the $`L_\alpha `$ flux has indeed been seen in the cosmological simulations of Fardal et al.  (2001). This may explain the observed high flux of $`L_\alpha `$ emitters at high redshift (e.g. Pentericci et al. , 2000, 2001; Breuck et al. , 2000). Based on the high observed flux and the assumption that the $`L_\alpha `$ is emitted from stars, Pentericci et al.  (2000) estimate large masses for the $`L_\alpha `$ emitters, but the much higher flux per unit mass predicted by our model may lead to significantly lower mass estimates. Based on our analysis, most of the $`L_\alpha `$ flux is expected to be emitted from the inner few kpc of the halo, where the gas is at $`10^4`$K. Neglecting line shifts and broadening, the halo might be opaque to $`L_\alpha `$ , thus eventually emitting its energy from an outer photosphere where the halo becomes transparent. However, a careful study of the thermal broadening and the systematic redshifts within the halo is required in order to determine whether the system is opaque or transparent to the $`L_\alpha `$ photons. This is a subject of an ongoing investigation.

Finally, the direct collapse of cold gas into the disc may have interesting theoretical consequences to be worked out. It may induce an efficient star burst in analogy to the burst originating in the shock between two colliding gas clouds. In turn, the strong inwards flow of gas may prevent an efficient gas removal by supernova-driven winds. In particular, current cosmological semi-analytic models (SAM) of galaxy formation (Kauffmann et al. , 1993; Kaffmann et al. , 1999; Cole et al. , 1994; Somerville & Primack, 1999; Maller et al. , 2001, and related works) use the standard picture of heating behind a virial shock in their modeling. This has strong effects on the disk formation rate, star formation rate, feedback etc. Other semi-analytic models (Efstathiou, 2000; White & Frenk, 1991) also appeal to the slow gas infall rate as a mechanism that regulates the gas input into the disc. Since the cooling time for a $`10^{11}M_{}`$ halo is relatively short, the SAM predictions for such haloes may be only slightly affected by the inhibition of heating. However, given some metal enrichment, no heating is expected for haloes as massive as $`7\times 10^{11}M_{}`$, for which the cooling time is longer, and the effect on the SAM predictions may be more severe. Shocks, when present, are also expected to alter the properties of the gas, for example - extinct dust particles. These effects can change SAMs that incorporate dust extinction.

## Acknowledgments

We acknowledge advice from Z. Barkat and E. Livne, J. Ostriker and stimulating discussions with S. Balberg, E. Bertschinger, T. Broadhurst, D. Gazit, Y. Hoffman, W. Mathews, A. Nusser, N. Shaviv, and S.D.M. White. This research has been supported by the Israel Science Foundation grant 213/02, by the German-Israel Science Foundation grant I-629-62.14/1999. and by NASA ATP grant NAG5-8218.

## Appendix A Testing the hydro code

The numerical code, Hydra, has been developed specifically for simulating the evolution of a single spherical halo through collapse and feedback processes. A proper computation of the cooling and shock formation requires high precision. In this appendix we describe a few of the tests performed in order to verify that the code works properly. In the following three subsections we test for energy conservation, spatial convergence, and the performance of the code in a self-similar case.

### A.1 Energy Conservation

Our numerical scheme does not use the total energy equation in the integration of the partial differential equations. Furthermore, the total energy of the system is not a straightforward sum of other variables that are involved in the calculation. The requirement of energy conservation is therefore an independent test for the accuracy of the numerical scheme. Energy conservation is harder to achieve than spatial convergence for several reasons. First, the error in total energy is systematic, in the sense that when dark-matter shells cross each other the energy tends to increase. Second, since our system is only marginally bound, the total energy is a small difference between two large quantities. We notice that energy conservation is simpler to achieve when there is no cooling, or when dark matter is absent (and thus there is no shell crossing).

The total energy of the system at time $`t`$ is the sum of terms:

$$E=K_\mathrm{d}+T_\mathrm{d}+K_\mathrm{g}+T_\mathrm{g}+U+Q,$$
(56)
where subscripts d and g refer to dark matter and gas respectively, $`K`$ stands for kinetic energy, $`T`$ stands for potential energy, $`U`$ is the gas internal energy, and $`Q`$ is the thermal energy lost to radiation by time $`t`$.

For the dark matter, these are straightforward sums over the discrete dark-matter shells:

$$K_\mathrm{d}=\mathrm{\Delta }M_i\frac{1}{2}\underset{i=1}{\overset{n_\mathrm{d}}{}}\left(v_i^2+\frac{j_i^2}{r_i^2}\right),$$
(57)
$$T_\mathrm{d}=\underset{i=1}{\overset{n_\mathrm{d}}{}}\frac{G\mathrm{\Delta }M_i}{r_i+a},$$
(58)
where $`_i`$ is the total mass interior to dark-matter shell $`i`$, as defined in §3.

For the gas shells, recall that the quantities $`r`$, $`v`$ and $`j`$ are given at the inner and outer shell boundaries, $`i1`$ and $`i`$ respectively, so we compute the shell energies by averaging over the two boundary values:

$`K_\mathrm{g}=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{n_\mathrm{g}}{}}}\mathrm{\Delta }m_i\left({\displaystyle \frac{v_ir_i^3+v_{i1}r_{i1}^3}{r_i^3+r_{i1}^3}}\right)^2`$
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{n_\mathrm{g}}{}}}\mathrm{\Delta }m_i\left({\displaystyle \frac{j_ir_i^2+j_{i1}r_{i1}^2}{r_i^3+r_{i1}^3}}\right)^2,`$
$$T_\mathrm{g}=\underset{i=1}{\overset{n_\mathrm{g}}{}}\frac{G\mathrm{\Delta }m_i(_i+_{i1})/2}{[(r_i^3+r_{i1}^3)/2]^{1/3}+a}.$$
(60)
The internal energy is a straightforward sum

$$U=e_i\mathrm{\Delta }m_i.$$
(61)
The energy radiated away, $`Q=𝑑tq𝑑m`$, is computed by

$$Q_=\underset{j=1}{\overset{n_\mathrm{t}}{}}\mathrm{\Delta }t^j\underset{i=1}{\overset{n_\mathrm{g}}{}}q_i^j\mathrm{\Delta }m_i,$$
(62)
where $`\mathrm{\Delta }t^j`$ is the length of timestep $`j`$, and $`q_i^j`$ is the cooling rate in shell $`i`$ at timestep $`j`$ (in units of $`\mathrm{erg}\text{ }\mathrm{g}^1\mathrm{s}^1`$).

In a run with 10,000 dark-matter shells and 2,000 gas shells, we require and obtain energy conservation at the level of 1% in a Hubble-time, using a typical Runge Kutta timestep of about $`5\times 10^6\text{Gyr}`$. (Such a run takes about 10 hours on an Alpha-6 DEC processor).

We check the conservation first by varying the accuracy parameters presented in §3.3, and then by varying the number of shells. The three cases presented in Table 1 demonstrate that the results converge when the accuracy is increased. The simulations shown in this table are of the standard case with realistic cooling shown in Fig. 3. When cooling is shut off, energy conservation is much better. With the nominal choice of accuracy parameters the final energy is $`0.9999`$ of the initial energy.

### A.2 Spatial Convergence: 3D versus 1D

A proper treatment of the competition between the pressure increase due to contraction and the pressure decrease due to cooling requires high temporal and spatial resolution. In particular, when the spatial resolution is increased, the shock appears earlier. Table 2 shows results from simulations of the case with realistic cooling (Fig. 3), all with the same accuracy criteria ($`ϵ_\mathrm{c}`$, $`ϵ_{\mathrm{rk}}`$, and $`t_{sc}`$), but with different spatial resolutions. The average distance between gas shells near the center ranges from about $`80`$pc to $`2`$kpc. With the poorest resolution of 125 gas shells the virial shock appears almost immediately after the virialization of the first shells of the simulation. The energy changed by about 75% during this simulation. Even if we assume that the precision of a 3D calculation is as good as that of an analogous 1D calculation (actually SPH codes converge slower than finite element schemes for problems involving shocks), we still need to cube the number of particles or grid points in order to achieve the same resolution. A three-dimensional simulation with $`2\times 10^6`$ gas particles and $`2.5\times 10^8`$ dark-matter particles, which is close to the limit of what is computationally feasible today, would correspond to the unsatisfactory case with the lowest spatial resolution in table 2.

### A.3 A self-similar case

When the initial conditions are scale free (unlike the initial conditions assumed in the body of this paper, motivated by $`\mathrm{\Lambda }`$CDM), and when the cooling function is also scale free (unlike the realistic cooling function used above), the results should be self similar. This can provide a test for the accuracy of our numerical code. We follow Bertschinger (1985a, b) in using an initial perturbation consisting of a point-mass embedded in a uniform-density background. Far from the point mass, the system should be self similar. We ran a simulation of such a case using our code with gas only ($`f_\mathrm{b}=1`$) and no cooling, starting at $`z=200`$ with an overdensity of 10 inside the innermost $`2`$kpc.

The upper panel of Fig. 12 shows the density profile at different times. As expected, a shock appears at every time as a density jump by a factor of 4 \[$`=(\gamma +1)/(\gamma 1)`$ for $`\gamma =5/3`$\], and the post-shock gas settles to a complete rest after it is shocked. (The slope of the post-shock density profile is somewhat different from Bertschinger (1985b), because our calculation assumes a $`\mathrm{\Lambda }`$CDM cosmology rather than the Einstein-deSitter assumed by Bertschinger.) The lower panel shows the same profiles after they were scaled to the same time (5Gyr) according to the scaling relation of Bertschinger (1985b): $`rt^{8/9}`$. We see that our simulations recover the expected scaling relation almost perfectly.

Fig. 13 shows an analogous test for the case where both gas and dark matter are present, with $`f_\mathrm{b}=0.13`$. The results are similar except for the somewhat higher noise level caused by the dark-matter component.

## Appendix B Top-Hat Model

Consider a bound spherical perturbation encompassing mass $`M`$, whose mean density fluctuation profile at some fiducial initial time in the linear regime is $`\overline{\delta }_\mathrm{i}(M)`$, embedded in an Einstein-deSitter (EdS) cosmological background when the universal expansion factor is $`a_\mathrm{i}`$. We wish to express the shell radius $`r`$ as a function of time in terms of $`M`$.

The implicit solution for a closed “mini-universe”, via a conformal time parameter $`\eta `$ specific to this perturbation, is

$$r=r_\mathrm{v}(1\mathrm{cos}\eta ),$$
(63)
$$t=\frac{r_\mathrm{v}}{v_\mathrm{v}}(\eta \mathrm{sin}\eta ).$$
(64)
Maximum expansion occurs when $`\eta =\pi `$, and then a collapse to half this radius, which we identify with virialization, is obtained at $`\eta =3\pi /2`$, with virial radius $`r_\mathrm{v}`$ and corresponding virial velocity $`v_\mathrm{v}^2=GM/r_\mathrm{v}`$. We normalize the universal expansion factor $`a`$ by identifying it at the initial time with the shell radius $`r`$. Assuming $`\eta _\mathrm{i}1`$, we have $`r_\mathrm{i}(1/2)r_\mathrm{v}\eta _\mathrm{i}^2`$ and $`t_\mathrm{i}(1/6)(r_\mathrm{v}/v_\mathrm{v})\eta _\mathrm{i}^3`$, so $`a_\mathrm{i}=r_\mathrm{i}`$ yields $`a_\mathrm{i}=(9r_\mathrm{v}v_\mathrm{v}^2/2)^{1/3}t_\mathrm{i}^{2/3}`$. The EdS expansion factor, $`at^{2/3}`$, can now be related using eq. (64) to the perturbation’s $`\eta `$ at any time:

$$a=(9r_\mathrm{v}v_\mathrm{v}^2/2)^{1/3}t^{2/3}=(9/2)^{1/3}r_\mathrm{v}(\eta \mathrm{sin}\eta )^{2/3}.$$
(65)
The mean density within the perturbation relative to the universal density at the same time becomes a straightforward function of $`\eta `$:

$$\frac{\overline{\rho }}{\rho _\mathrm{u}}=\frac{a^3}{r^3}=\frac{9}{2}\frac{(\eta \mathrm{sin}\eta )^2}{(1\mathrm{cos}\eta )^3}.$$
(66)
This is a standard result of the top-hat model.

In order to relate the density to the small initial perturbation at $`\eta _\mathrm{i}1`$, we obtain from eq. (66) by a proper Taylor expansion to the first non-vanishing order:

$$\overline{\delta }_\mathrm{i}0.15\eta _\mathrm{i},$$
(67)
where $`\overline{\delta }\overline{\rho }/\rho _\mathrm{u}1`$ is the mean fluctuation. Using this in the linear term of eq. (63) we obtain

$$a_\mathrm{i}=(1/2)r_\mathrm{v}\eta _\mathrm{i}^2=(1/0.3)r_\mathrm{v}\overline{\delta }_\mathrm{i}.$$
(68)
This allows us to write the mean density at any $`\eta `$, using eq. (63), as

$$\frac{\overline{\rho }}{\rho _{\mathrm{ui}}}=\frac{a_\mathrm{i}^3}{r^3}=\frac{\overline{\delta }_\mathrm{i}^3}{0.3(1\mathrm{cos}\eta )^3}.$$
(69)
Recalling that $`\overline{\rho }=M/[(4\pi /3)r^3]`$ we finally obtain at any $`\eta `$

$$r=C_\mathrm{i}\frac{M^{1/3}}{\overline{\delta }_\mathrm{i}(M)}(1\mathrm{cos}\eta ),$$
(70)
where

$$C_\mathrm{i}\left(\frac{6}{\pi }\right)^{1/3}\frac{0.15}{\rho _{\mathrm{ui}}}.$$
(71)
In particular, at $`\eta =3\pi /2`$, we obtain for the virial radius $`r_\mathrm{v}=C_\mathrm{i}M^{1/3}/\overline{\delta }_\mathrm{i}(M)`$. The constant $`C_\mathrm{i}`$ is independent of $`M`$; the universal density $`\rho _{\mathrm{ui}}`$ \[$`=(1+z_\mathrm{i})^3\rho _{\mathrm{u0}}`$\] is determined by the choice of the fiducial redshift $`z_\mathrm{i}`$ at which $`\overline{\delta }_\mathrm{i}(M)`$ is given.

## Appendix C Initial Profile

We adopt in the linear regime the typical density fluctuation profile for the assumed power spectrum of fluctuations. For a Gaussian random field, this profile is proportional to the two-point correlation function (Dekel, 1981):

$$\delta _\mathrm{i}(r)=\delta _{0\mathrm{i}}\frac{\xi (r)}{\xi (0)},$$
(72)
where $`\delta _{0\mathrm{i}}`$ specifies the amplitude normalization. For a given power spectrum $`P(k)`$, the correlation function is given by (Peebles, 1993, eq. 21.40)

$$\xi (r)=4\pi _0^{\mathrm{}}k^2𝑑kP(k)\frac{\mathrm{sin}kr}{kr},$$
(73)
and the local variance is

$$\xi (0)=4\pi _0^{\mathrm{}}k^2𝑑kP(k).$$
(74)
The mean density fluctuation interior to radius $`r`$, containing mass $`M=(4\pi /3)\rho _{\mathrm{ui}}r^3`$ when the fluctuation is small, is

$$\overline{\delta }_\mathrm{i}(r)=\frac{3}{r^3}_0^rr^2𝑑r\delta _\mathrm{i}(r).$$
(75)
This involves the integral (Peebles, 1993, eq. 21.62)

$$J_3(r)_0^rr^2𝑑r\xi (r)=\frac{4\pi r^3}{3}_0^{\mathrm{}}k^2𝑑kP(k)\stackrel{~}{W}_s(kr),$$
(76)
where $`\stackrel{~}{W}_s(kr)`$ is the Fourier transform of the top-hat window,

$$\stackrel{~}{W}_s(kr)=3\left[\frac{\mathrm{sin}kr}{(kr)^3}\frac{\mathrm{cos}kr}{(kr)^2}\right]=\frac{3}{kr}j_1(kr),$$
(77)
with $`j_1`$ the spherical Bessel function.

We use in the simulations of this paper the $`\mathrm{\Lambda }`$CDM power spectrum as from the fitting formula of Bardeen et al.  (1986)

$$P(k)=AkT^2(k),$$
(78)
$$T(k)=\left\{1+\left[ak/\mathrm{\Gamma }+(bk/\mathrm{\Gamma })^{3/2}+(ck/\mathrm{\Gamma })^2\right]^\nu \right\}^{1/\nu }$$
(79)
with $`a=6.4h^1`$Mpc, $`b=3.0h^1`$Mpc, $`c=1.7h^1`$Mpc, $`\nu =1.13`$ and $`\mathrm{\Gamma }=0.21`$ (the $`\tau `$CDM model of Efstathiou et al. , 1992). The normalization is such that $`\sigma _8=1.`$

In the cosmological toy model we approximate the $`\mathrm{\Lambda }`$CDM by a power-law power spectrum $`P_kk^n`$, for which the two-point correlation function is also a power law, $`\xi (r)r^{(n+3)}`$, and then

$$\overline{\delta }_\mathrm{i}(r)=\left(\frac{r}{r_1}\right)^{(n+3)},$$
(80)
where $`r_1`$ provides the normalization. In terms of mass we obtain

$$\overline{\delta }_\mathrm{i}(M)=\left(\frac{M}{(4\pi /3)\rho _{\mathrm{ui}}r_1^3}\right)^{(n+3)/3}.$$
(81)
This serves as the input to eq. (70), or eq. (48).

We normalize the initial perturbation such that a specific mass $`M_1`$ reaches virialization at some cosmological epoch $`a_\mathrm{v}=1/(1+z_\mathrm{v})`$. Using eq. (65) and eq. (68) we obtain the linear analog to the nonlinear fluctuation growth rate:

$$\frac{\overline{\delta }(\eta )}{\overline{\delta }_\mathrm{i}}\frac{a(\eta )}{a_\mathrm{i}}=0.3\left(\frac{9}{2}\right)^{1/3}(\eta \mathrm{sin}\eta )^{2/3}.$$
(82)
At virialization, this gives $`\delta _\mathrm{v}\overline{\delta }(\eta =3\pi /2)1.58`$. Then:

$$\overline{\delta }_\mathrm{i}(M_1)=\delta _\mathrm{v}\left(\frac{a_\mathrm{i}}{a_\mathrm{v}}\right).$$
(83)
The normalization parameter $`\delta _{0\mathrm{i}}`$ (or $`r_1`$) at $`a_\mathrm{i}`$ is obtained by equating this with eq. (75) \[or eq. (81)\] at $`M=M_1`$.