text:
# Quantum entanglement in spinor Bose-Einstein condensates

## Abstract

We propose a scheme to generate and detect various kinds of quantum entanglement in a spin-1 Bose-Einstein condensate. It is shown that substantial many-particle entanglement can be generated directly in the spin-1 condensate by free dynamical evolution with a properly prepared initial state. The scheme also provides a simple method to generate three-mode entanglement in the second-quantization picture and to detect the continuous variable type of entanglement between two effective modes in the spin-1 condensate.

The achievement of atomic Bose-Einstein condensate (BEC)  has raised extensive experimental and theoretical studies in this area, and these studies have opened up some possibilities for important applications. An impressive example of these applications is to use BECs with internal degrees of freedom for generation of quantum entanglement \[2-7\], which is the essential ingredient of many quantum information protocols . To generate entanglement, one requires an experimental system which can be prepared in a pure state, and there should be significant interactions between the relevant particles which can be coherently controlled. BECs fulfil these requirements, and therefore they could provide an ideal experimental system for studying quantum entanglement. It has been shown that with coherent collision interactions in BECs, it is possible to generate substantial many-particle entanglement in a two-component condensate , or to generate entangled atomic beams from the outputs of a single condensate prepared in a specific internal level which enables spin-exchange collisions .

A good experimental system for demonstration of these proposals is the Sodium condensate confined in an optical trap \[9-11\]. Recently, Stamper-Kurn et al. have realized an optically trapped BEC in which all the three ground Zeeman levels of sodium atoms are involved . In such a three-component (so-called spin-1) condensate, there are spin-exchange collisions, which constantly mix the different spin components, and other kinds of internal spin dynamics, which result in complex nonlinear collision phase shifts \[12-17\]. In this paper, we will investigate whether one can directly prepare and detect quantum entanglement in such a spin-1 condensate. The motivation of this work is two-folds: firstly, direct detection of quantum entanglement in the realized spin-1 condensate eliminates the requirements of the level shifting and the output techniques in previous works \[2-4\], and thus allows for a easier experiment; and secondly, the spin-1 condensate with the presence of both the spin-exchange collisions and the nonlinear phase shifts is a more complicated system, and opens up possibilities for richer physical behaviors. We will show that various kinds of entanglements are available in this system. We can get substantial many-particle entanglement (the spin squeezing type of entanglement ) as well as three-mode and two-mode entanglement (the continuous-variable type of entanglement ). All these entanglements are useful for different purposes of applications.

The spin-1 BEC has been studied both theoretically \[12-17\] and experimentally \[9-11\], with the emphasis on the ground-state structure, the spin-mixing dynamics, and the response to the external magnetic field. Here, we will focus on its entanglement property. In Ref. , it was shown as a by-product that the spin ground state of the sodium condensate with the anti-ferromagnetic interaction should be entangled. Unfortunately, this ground-state entanglement was not confirmed by the experiment  because it is very hard to cool the system to the spin ground state (the energy gap between the spin ground and the first excited states is many orders smaller than the chemical potential), and the spin ground state is very fragile under the influence of the external magnetic field . Here, we describe a different method to generate and detect quantum entanglement by free dynamical evolution with a proper initial state. As we will see, this provides a feasible method of detecting entanglement in the spin-1 BEC with the current technology.

Before we enter the discussion of the spin-1 condensate, let us first clarify different kinds of entanglements available in a bosonic many-particle system. Consider $`N`$ bosonic identical particles with each particle having three internal levels $`|+1`$, $`|1`$, and $`|0`$. First, the $`N`$ particles are said to be entangled in the first quantization picture if their internal density matrix $`\rho `$ cannot be decomposed into

$$\rho =\underset{k}{}p_k\rho _1^k\rho _2^k\mathrm{}\rho _N^k\text{,}$$
(1)
where the coefficients $`p_k`$ are positive real numbers satisfying $`_kp_k=1`$, and $`\rho _i^k`$ is a $`3\times 3`$ density matrix of the $`i`$th particle. To check whether an $`N`$-particle system is entangled, a good experimentally detectable criterion is the squeezing parameter. For particles with only two internal levels $`|a`$ and $`|b`$, one can construct a collective spin operator $`𝐉_{ab}=_{i=1}^N𝐣_i`$, with the individual spin-$`1/2`$ operators $`j_i^x=\left(|a_ib|+|b_ia|\right)/2,`$ $`j_i^y=i\left(|a_ib||b_ia|\right)/2,`$ and $`j_i^z=\left(|b_ib||a_ia|\right)/2`$. The squeezing parameter $`\xi _{ab}^𝐧`$ is defined directly from the variance and mean values of the $`𝐉_{ab}`$ operator

$$\xi _{ab}^𝐧=\frac{N\left(\mathrm{\Delta }J_{ab}^{n_1}\right)^2}{J_{ab}^{n_2}^2+J_{ab}^{n_3}^2},$$
(2)
where $`J_{ab}^n=𝐧𝐉_{ab}`$, and the $`𝐧`$s are mutually orthogonal unit vectors. For two-level particles, the spin squeezing parameter (2) with a value less than $`1`$ provides a sufficient condition for many-particle entanglement  (such a case is also referred as spin squeezing ). For three-level particles, we can construct two orthonormal states $`|a`$ and $`|b`$ from arbitrary superpositions of the three levels $`|+1`$, $`|1`$, and $`|0`$, and define the $`𝐉_{ab}`$ and $`\xi _{ab}^𝐧`$ in the same way. Then, $`\xi _{ab}^𝐧<1`$ assures that the projected $`N`$-particle density matrix $`P\rho P`$ onto the subspace spanned by the states $`|a_i`$ and $`|b_i`$ ($`i=1,2,\mathrm{},N`$) is entangled. Since the projection $`P`$ is a local operator which cannot increase entanglement, the original density matrix $`\rho `$ for the $`N`$ three-level particles should also be entangled if $`\xi _{ab}^𝐧<1`$. For our later applications, we define the spin operators $`𝐉_\pm `$ and the squeezing parameters $`\xi _\pm ^\theta `$ respectively by choosing $`|a=|\pm =\left(|+1\pm |1\right)/\sqrt{2}`$, $`|b=|0`$, and $`n_1=[\mathrm{cos}\theta ,\mathrm{sin}\theta ,0]`$ in Eq. (2). Either of $`\xi _\pm ^\theta `$ less than $`1`$ for some $`\theta `$ is a signature of many particle entanglement. Note that the spin squeezing type of entanglement represented by a small $`\xi _{ab}^𝐧`$ is useful in practice since the inverse of $`\xi _{ab}^𝐧`$ gives exactly the improvement of the single-to-noise ratio for atomic clocks relative to the standard quantum limit with the Ramsey scheme involving the levels $`|a`$ and $`|b`$ .

The above system can also be described in the second quantization picture. We have three bosonic modes $`a_\alpha `$ $`\left(\alpha =0,\pm 1\right)`$ corresponding to the three internal levels, with the commutation relations $`[a_\alpha ,a_\beta ^{}]=\delta _{\alpha \beta }`$. These modes can also be entangled like the individual particles. Similar to the definition in Eq. (1), we say that two modes $`a_\alpha ,a_\beta `$ are entangled if their density matrix cannot be expressed as $`_kp_k\rho _\alpha ^k\rho _\beta ^k`$, and all the three modes $`a_\alpha ,a_\beta ,a_\gamma `$ are intrinsically entangled if their density matrix cannot be expressed as $`_kp_k\rho _\alpha ^k\rho _{\beta \gamma }^k`$, where $`\alpha ,\beta ,\gamma =0,\pm 1`$, $`p_k`$ are positive real numbers, and $`\rho _{\alpha (\beta ,\beta \gamma )}^k`$ is the a density matrix of the corresponding mode(s). Note that the mode entanglement in the second quantization and the particle entanglement in the first quantization are two different concepts, and they are useful for different purposes of applications. One can easily find example states which are mode entangled but not particle entangled, or particle entangled but not mode entangled . Mode entanglement of bosonic atoms could be useful in various continuous variable or qubit quantum information processing schemes currently based on entangled optical modes . For two modes $`a_{\pm 1}`$, there is an easy experimentally detectable entanglement criterion . One can define quadrature phase amplitudes $`X_{\pm 1}^\theta =\left(e^{i\theta }a_{\pm 1}^{}+e^{i\theta }a_{\pm 1}\right)/\sqrt{2}`$ and $`X_\pm ^\theta =\left(X_{+1}^\theta \pm X_1^\theta \right)/\sqrt{2}`$ for the two modes $`\pm 1`$. It has been shown in Ref.  that $`\left(\mathrm{\Delta }X_+^\theta \right)^2+\left(\mathrm{\Delta }X_{}^{\theta +\pi /2}\right)^2<1`$ for some angle $`\theta `$ provides a criterion for entanglement between the modes $`a_{+1}`$ and $`a_1`$ (a state with this property is also referred as a two-mode squeezed state, or a continuous variable entangled state ). For later applications in the spin-$`1`$ condensate, we are more interested in two effective bosonic modes $`a_{\pm 1}^{}`$, defined as $`a_{\pm 1}^{}=a_{\pm 1}a_0^{}/\sqrt{a_0^{}a_0}`$. The modes $`a_{\pm 1}^{}`$ are approximately bosonic if nearly all the atomic populations remain in the component $`a_0`$, i.e., $`a_0^{}a_0N`$. The definition of the new modes $`a_{\pm 1}^{}`$, which are proportional to $`a_{\pm 1}`$, is to remove the phase randomness in the condensate by phase locking both of the modes $`a_{\pm 1}`$ to the condensate mode $`a_0`$. For the modes $`a_{\pm 1}^{}`$, the quadrature phase amplitudes defined above reduce to $`\left(\mathrm{\Delta }X_\pm ^\theta \right)^2\xi _\pm ^\theta /2`$. The inequality $`\xi _+^\theta +\xi _{}^\theta <2`$ from the spin squeezing parameters thus becomes a sufficient criterion for entanglement between the modes $`a_{\pm 1}^{}`$. As a result, a measurement of the spin squeezing parameters $`\xi _\pm ^\theta `$ provides criteria for both many-particle entanglement and two-mode entanglement.

We will see that three-mode entanglement is also produced in a spin-1 condensate. A general criterion for the three-mode entanglement is still absent, but our system to a good approximation will be in a pure three-mode state $`|\varphi _{+1,1,0}`$, and in this case the three-mode entanglement can be easily characterized by the quantity $`E_3=\sqrt[3]{E\left(\rho _{+1}\right)E\left(\rho _1\right)E\left(\rho _0\right)}`$, where $`E\left(\rho _\alpha \right)=tr\left(\rho _\alpha \mathrm{log}_2\rho _\alpha \right)`$ denotes the von Neuman entropy of $`\rho _\alpha `$, and $`\rho _\alpha `$ is the reduced density operator of the mode $`\alpha `$ $`\left(\alpha =+1,1,0\right)`$. It is obvious that for pure states $`E_3`$ is non-increasing under local operations and positive if and only if the system state is intrinsically three-mode entangled. It will be shown that in a spin-1 condensate $`E_3`$ can be significantly larger than $`0`$. Note that different to the spin squeezing parameters it is not easy to directly measure $`E_3`$.

Now we turn to the physical system, which is a dilute gas of trapped bosonic atoms with hyperfine spin $`f=1`$. Assume that the system is very cold so that the collision interaction between atoms is effectively described by a pair-wise pseudo $`\delta `$-potential, which preserves the hyperfine spin of individual atoms and is rotationally invariant in the hyperfine spin space . With this symmetry assumption, the most general form of the Hamiltonian can be dived into two parts $`H=H_s+H_a`$, with the symmetric and asymmetric Hamiltonians $`H_s`$ and $`H_a`$ respectively given by  (setting $`\mathrm{}=1`$)

$`H_s`$ $`=`$ $`{\displaystyle \underset{\alpha }{}}{\displaystyle 𝑑𝐫\mathrm{\Psi }_\alpha ^{}\left(\frac{^2}{2m}+V\right)\mathrm{\Psi }_\alpha }`$ (4)
$`+{\displaystyle \frac{\lambda _s}{2}}{\displaystyle \underset{\alpha ,\beta }{}}{\displaystyle 𝑑𝐫\mathrm{\Psi }_\alpha ^{}\mathrm{\Psi }_\beta ^{}\mathrm{\Psi }_\alpha \mathrm{\Psi }_\beta },`$
$$H_a=\frac{\lambda _a}{2}\underset{\alpha ,\beta }{}𝑑𝐫\mathrm{\Psi }_\alpha ^{}\mathrm{\Psi }_\alpha ^{}^{}𝐒_{\alpha \beta }𝐒_{\alpha ^{}\beta ^{}}\mathrm{\Psi }_\beta \mathrm{\Psi }_\beta ^{},$$
(5)
where $`\mathrm{\Psi }_\alpha `$ $`\left(\alpha =1,0,+1\right)`$ is the atomic field operator associated with atoms in the spin state $`|f=1,m_f=\alpha `$. The first and the second terms of the Hamiltonian $`H_s`$ denote the kinetic energy and the potential energy, with $`m`$, the mass of the atom, and $`V`$, the trap potential, which has been assumed to be the same for all the three components. The $`3\times 3`$ spin matrices $`𝐒`$ in $`H_a`$ denote the conventional $`3`$-dimensional representation (corresponding to the spin $`f=1`$) of the angular momentum operator, with $`S_{\alpha \beta }^x=\left(\delta _{\alpha ,\beta 1}+\delta _{\alpha ,\beta +1}\right)/\sqrt{2}`$, $`S_{\alpha \beta }^y=i\left(\delta _{\alpha ,\beta 1}\delta _{\alpha ,\beta +1}\right)/\sqrt{2}`$, $`S_{\alpha \beta }^z=\alpha \delta _{\alpha \beta }`$. The interaction energy has been divided into the symmetric and the asymmetric parts, represented respectively by the last term of $`H_s`$ and the Hamiltonian $`H_a`$. The symmetric interaction remains the same if we exchange two arbitrary spin components, whereas the asymmetric one is spin-dependent. The symmetric and the asymmetric collision coefficients $`\lambda _s`$ and $`\lambda _a`$ are determined respectively by the average and the difference of the scattering lengths in two different collision channels. For sodium atoms, $`\lambda _a`$ is positive and about $`25`$ times smaller than $`\lambda _s`$ . Thus, the Hamiltonian $`H_s`$ is the dominant one for determining the dynamics of the motional state. However, the internal spin dynamics which we are interested in is completely determined by $`H_a`$ since the symmetric part $`H_s`$ is spin-independent.

For generating entanglement, we start with a BEC prepared in the internal level $`|m_f=0`$. To quantify the amount of various kinds of squeezing and entanglement obtainable in this system, we need to calculate the time evolution of the entanglement criteria $`\xi _\pm ^\theta `$ and $`E_3`$ defined before. For this purpose, we make use of the single motional mode approximation to the Hamiltonian $`H_a`$ following Refs. \[15-17\]. Note that the Hamiltonian $`H_s`$ is the dominant one for determining the motional state. When the system is at a very low temperature, the motional state is frozen to be approximately the ground state $`\phi \left(𝐫\right)`$ (normalized as $`𝑑𝐫\left|\phi \left(𝐫\right)\right|^2=1`$) of the spin-independent Hamiltonian $`H_s`$, and we can factorize the field operators as $`\mathrm{\Psi }_\alpha a_\alpha \phi \left(𝐫\right)`$ $`\left(\alpha =0,\pm 1\right)`$, where $`a_\alpha `$ is the usual bosonic annihilation operator. With this approximation, one can introduce a spin operator $`𝐋`$ with $`𝐋=_{\alpha ,\beta }a_\alpha ^{}𝐒_{\alpha \beta }a_\beta `$, and the spin-dependent Hamiltonian $`H_a`$ is simplified to $`H_a=\lambda _a^{}\left(𝐋^22N\right)`$, where $`\lambda _a^{}=\lambda _a𝑑𝐫\left|\phi \left(𝐫\right)\right|^4/2`$. This Hamiltonian has been used for studying the spin ground state structure . We use it here for calculating squeezing and entanglement in this system. For this purpose, it is convenient to express the Hamiltonian $`H_a`$ by the more relevant spin operators $`𝐉_\pm `$ defined before. For the specific system configuration considered here, $`L^z=0`$ for the initial state, and it will remain zero since $`L^z`$ is conserved. On the other hand, by definition we have $`L^x=2J_+^x`$ and $`L^y=2J_{}^y`$. Therefore, the Hamiltonian $`H_a`$ expressed by $`𝐉_\pm `$ has the form

$$H_a=4\lambda _a^{}\left(J_+^{x2}+J_{}^{y2}\right),$$
(6)
where we have dropped the irrelevant constant $`2\lambda _a^{}N`$. The only non-zero commutation relation between $`𝐉_+`$ and $`𝐉_{}`$ is $`[J_+^x,J_{}^x]=[J_+^y,J_{}^y]=\left(a_1^{}a_{+1}a_{+1}^{}a_1\right)/4`$, which remains very small if the population of the atoms is still dominantly in the level $`|m_f=0`$, as it is in our scheme. In this case, as an approximation we can consider separately the dynamics of $`𝐉_+`$ and $`𝐉_{}`$, respectively with the effective Hamiltonian $`J_+^{x2}`$ or $`J_{}^{y2}`$. It is well known that the Hamiltonian $`J_+^{x2}`$ or $`J_{}^{y2}`$ will produce a spin squeezing state, with the minimum squeezing parameter about $`\left(3/N\right)^{2/3}/2`$ after an evolution time proportional to $`N^{2/3}/\lambda _a^{}`$ . This is the intuitive reason why we can get substantial spin squeezing and entanglement in this system. To be more quantitative, we have calculated exactly the time evolution of the squeezing parameter $`\xi _\pm ^\theta `$ under the Hamiltonian (5) using the numerical approach. The system spin state $`|\varphi `$ under the Hamiltonian (5) can be expanded as $`|\varphi =_nc_n|N2n,n,n_{0,+1,1}`$ in the number basis, where the subscript $`\alpha `$ $`\left(\alpha =0,\pm 1\right)`$ represents the corresponding internal mode $`a_\alpha `$. One can solve numerically the time evolution of the expansion coefficients $`c_n`$, from which the spin squeezing parameters $`\xi _\pm ^\theta `$ are exactly obtainable. We find that the squeezing magnitudes $`\xi _\pm =\mathrm{min}_\theta \xi _\pm ^\theta `$ evolve in the same way for $`𝐉_+`$ and $`𝐉_{}`$, while their squeezing directions corresponding to the optimal angle $`\theta `$ are always orthogonal to each other. This means that the criterion $`\mathrm{min}_\theta \left(\xi _+^\theta +\xi _{}^{\theta +\pi /2}\right)/2`$ for entanglement between the two effective modes $`a_{\pm 1}^{}`$ evolves in the same way as the criterion $`\xi _\pm `$ for many-particle entanglement . The resultant time evolution of $`\xi _\pm `$ is shown in Fig. 1a, together with a prediction from the simple Hamiltonian $`J_+^{x2}`$ (or $`J_{}^{y2}`$). One can see that the results from the two approaches agree quite well, with only slight difference in the maximum obtainable squeezing and the corresponding evolution time. For $`10^5`$ atoms, about three orders reduction is obtainable in $`\xi _\pm `$, which is a clear demonstration that substantial many-particle entanglement and two-mode entanglement have been generated. Fig. 1b shows the time evolution of the three-mode entanglement $`E_3`$. One can see that significant three-mode entanglement can be generated with only $`10^3`$ atoms.

Finally, we would like to address some imperfections and approximations used in this scheme. First, in a practical case there would be some remaining external magnetic field due to the special way of preparing the condensate. It is remarkable that the linear Zeeman shift from a known homogeneous magnetic field has no influence on our scheme since it commutes with the spin-dependent Hamiltonian $`H_a`$ and will be eliminated if one transforms to the interaction picture. As a result, only quadratic Zeeman shift and linear Zeeman shift from an inhomogeneous field have influence. These high-order effects can be easily controlled very small in the relevant time scale. For instance, one can see from Fig. 1 that the time scale in this scheme is roughly set by the inverse of $`N\lambda _a^{}`$, which is about $`2\pi \times 10^2`$Hz for a typical density $`3\times 10^{14}`$cm<sup>-3</sup> of sodium atoms . Within this time, the quadratic and the inhomogeneous Zeeman shifts have negligible effects if the magnetic field and its gradient are respectively smaller than $`1`$mG and $`1`$mG/cm. Secondly, we have assumed that the three components with $`m_f=0,`$ $`\pm 1`$ have the same spatial wave function. This is a good approximation as long as there is no demixing instability leading to phase separation of different components . For the two-component condensate with components $`\alpha `$ and $`\beta `$ having equal macroscopic populations, a mean-field theory has shown that there is no demixing instability if the inter-component collision coefficient $`\lambda _{\alpha \beta }`$ is smaller than the geometric average of the intra-component collision coefficients $`\sqrt{\lambda _{\alpha \alpha }\lambda _{\beta \beta }}`$ . For the relevant dynamics in our scheme with only the level $`|m_f=0`$ macroscopically populated, Fig. 1 shows that this system is well approximated by two independent sets of two-component condensates. For each set of two-component condensate, we can apply the above result, and find there is no demixing instability in this scheme if the asymmetric collision coefficient $`\lambda _a`$ is positive, as it is for sodium atoms. This is different to the case where both of the levels with $`m_f=0,`$ $`+1`$ are equally populated, which leads to demixing instabilities . Finally, we assumed the single-mode approximation for the atomic motional state. This approximation is supported by the following arguments: first, in the two-component condensate more complete numerical simulations gives a squeezing which agrees quite well with the calculation from the single-mode approximation , and the present spin-1 system is well approximated by two independent sets of two-component condensates; and second, a direct energy analyses in Ref.  also shows the validity of the single-mode approximation for the spin-1 condensate with only the level $`|m_f=0`$ initially populated.

In summary, we have shown how to directly generate and detect various kinds of squeezing and entanglement in the experimentally realized spin-1 condensate. This helps for the final observation of quantum entanglement in Bose-Einstein condensates.

This work was supported by the Austrian Science Foundation, the EQUIP, the ESF, the European TMR network Quantum Information, and the Institute for Quantum Information Ges.m.b.H. L.M.D. thanks in addition the support from the Chinese Science Foundation.


format: text

text:
# Multi-Stage Three Dimensional Sweeping and Annealing of Disc Galaxies in Clusters

## 1 INTRODUCTION: HIGHLIGHTS OF RAM PRESSURE STRIPPING THEORY

In this paper we present new three dimensional numerical models of stripping from gas-rich disc galaxies in clusters, with radiative cooling included. We begin with a brief review the literature to put this work in context.

Gunn and Gott (1972) initiated the study of ram pressure stripping of disc gas in cluster galaxies by the intracluster medium (ICM) in their seminal study of the formation and evolution of galaxy clusters. Specifically, they estimated how the ICM ram pressure compared to the local gravitational binding of the gas at a representative point in a typical galaxy disc (e.g., the solar radius in the Milky Way), and found that stripping was likely to be very effective, and by implication prompt. The Gunn and Gott gravitational binding estimate included only the stellar and gas disc components, but no halo.

The first numerical hydrodynamical simulations were published within a few years by Lea and De Young (1976), and Gisler (1976). Lea and De Young found that 80-90% of the gas in the model galaxy was stripped within one crossing of a cluster like Coma, or within about $`10^9`$ yrs. However, calculations were two-dimensional, used a beam scheme algorithm, and an adiabatic equation of state. The gravity of the early type model galaxy was assumed to be dominated by a spherical, isotropic stellar potential, with only a single hot gas component, initially in hydrostatic equilibrium. Gisler also considered the gas mass evolution in cluster ellipticals, including the effects of: gas input by stellar mass loss, gas loss in galactic winds, and ICM stripping. His two-dimensional, fluid-in-cell simulations also produced substantial stripping.

In the late 1970s and early 1980s, satellite X-ray observations of the hot gas in clusters, and of galaxies within the nearest clusters were published (see review of Forman and Jones 1982). The X-ray emission characteristics were commonly interpreted in terms of stripping models, and provided support for those models as well as a more precise picture of the ICM.

At the same time, evidence mounted that all types of spirals in clusters are gas poor relative to their neighbors in the field, specifically that they are HI deficient (e.g., Bothun 1982, Giovanelli and Haynes 1983). Systematic dependences on galaxy position and velocity relative to cluster core values supported stripping as the cause of the deficits, rather than formation processes. More recent work has confirmed, and greatly extended, these results, e.g., in the Virgo cluster (Cayatte et al. 1990), and Coma (Bravo-Alfaro et al. 2000, also see Van Gorkom’s 1996 review of HI in cluster galaxies).

Kenney and Young (1986, 1988) noted that molecular gas, typically found at smaller radii in galaxy discs than the HI, is not stripped. This showed that stripping is not strong enough to remove the denser gas in the inner discs of spirals, and is not as complete as we might have expected from the earliest models.

Farouki and Shapiro (1980) carried out three-dimensional particle simulations to study the thickening of stellar discs following stripping. While their primary goal was to investigate the SO origin question, and their models had only 100 “gas” particles, their simulations were the first on stripping from discs. They found: that only the outer parts of discs would be swept in some cases, that the amount of stripping depended on ICM parameters, and that stripping might not always occur promptly. Their statement that - “The action of a strong restoring force on gas clouds venturing from the disc following ejection gives rise to complicated motions in the simulations” - is prescient and relevant to results we present below.

The simulations and analytic work of Takeda, Nulsen & Fabian (1984), advanced the theory in several areas. They also carried out two-dimensional simulations of spherical galaxies with a spherical gas distribution using a fluid-in-cell algorithm. In their simulations, they assumed the galaxy traveled on a nearly radial orbit, and used a cluster model with a hydrostatic, isothermal gas distribution, to derive density and velocity variations along the galaxy orbit. At the outermost excursion, the galaxy velocity and ICM density are low, stripping is unimportant, and the galaxy rebuilds its gas supply via stellar mass loss.

Takeda et al. identified three stripping regimes: 1) prompt, in which ram pressure immediately overwhelms gravitational binding and removes the gas, 2) continuous, in which it takes the ICM wind some time to provide the momentum to remove the gas (and the gas may be partially replenished on a comparable timescale, or even accrete with low relative velocities), 3) cyclical, in which gas is stripped in the cluster core, but replenished in the outer cluster.

Nulsen (1982), in a review of both hydrodynamic and kinetic processes in stripping concluded that turbulent shear viscosity would play an important role. However, Takeda et al. state that the shear viscosity included in their models did not affect their basic results, such as the existence of different stripping regimes. Recently, Mori & Burkert (2000) have resolved Kelvin Helmholtz turbulence in models of stripping from dwarf galaxies. While they find that it significantly enhances mass loss from dwarfs, the characteristic estimated timescale for the process is very long in normal, disc galaxies.

The work of Gaetz, Salpeter & Shaviv (1987, also see Shaviv and Salpeter, 1982) introduced more physical processes into models of stripping from spherical galaxies. These authors used a second order Eulerian hydrodynamics code to carry out two-dimensional simulations, with a steady-state cooling curve, and algorithms for star formation, gas consumption and replenishment. The various processes were characterized by 5 timescales and their ratios. They found that with these competing processes included it was possible to achieve a steady state, for galaxies moving through a uniform medium. The degree of stripping depended on the parameter $`n_{icm}v_{gal}^{2.4}/\dot{M}_{rep},`$ where the variables are the density of the ICM, the velocity of the galaxy through the medium, and the constant replenishment rate, respectively. Portnoy, Pistinner, & Shaviv (1993) updated and extended the models of Gaetz et al., and provided fitting formulae for the effects of various parameters.

More recent high resolution models of stripping from spherical galaxies have been presented by Balsara, Livio, & O’Dea (1994), and Stevens, Acreman, and Ponman (1999). The models of both groups were calculated with two dimensional, PPM codes, which included the effects of star formation and replenishment. Like Takeda et al., Stevens et al. found a number of stripping modes, including a cyclic one. Both works had high spatial resolution in the wake, and they argued that the numerical viscosity of previous works limited their ability to resolve shocks. However, some contradictory results, and questions about kinetic effects and the need for a two-fluid prescription (see Portnoy, et al. 193), leave the theory of stripping of hot gas from spherical galaxies in a rather ambiguous state.

Kritsuk (1983) examined the stripping of different types of interstellar clouds, located at different galactic radii, and with ICM flows at a range of angles relative to the disc plane. Kritsuk found that diffuse HI clouds would be readily swept in most cases, while dense molecular clouds would have a low stripping probability. His approach was to derive an equation of motion for clouds, subject to gravity from the stellar disc and ram pressure forces, including mass accretion and ablation. The latter effects, and cloud structural evolution were treated very approximately.

From the cloud equation of motion he derived an azimuthally dependent stripping radius in the galaxy disc. This work provided information about three-dimensional stripping, and except for one run of Farouki and Shapiro, is the first time the effects of inclination where examined. For example, Kritsuk’s Figure 2 shows a strong coupling between disc rotation and the ICM wind in highly inclined cases, a result that will be emphasized in the new models below.

A number of recent observations of asymmetric, extra-disc gas distributions, or wakes, in cluster galaxies give a direct view of the stripping process. E.g., van Driel and van Woerden (1989) mapped the HI emission from NGC 4694, a Virgo cluster lenticular, and located most of the emission in an HI tail, 36 kpc. in length. They favor a stripping explanation for the origin of the tail, though the presence of a small dwarf galaxy makes an interaction origin possible.

M96 (NGC 4406), a Virgo elliptical with X-ray, HI and infrared emission from an asymmetric plume (White et al. 1991) has long been considered a possible example of stripping. The giant Virgo elliptical NGC 4472 seems to be in a very similar circumstance. X-ray observations show an asymmetic gas plume of comparable size to the optical galaxy, and an apparent bow shock on the opposite side (Irwin and Sarazin 1996).

These examples confirm the general morphology of stripping suggested by the early hydrodynamic models. The presence of gas within cluster ellipticals also supports the conclusions that either gas is not entirely stripped from early-type galaxies, or replenishment is important, and perhaps leads to a steady or slowing changing state.

The discovery in the last decade that many clusters continue to accrete groups of galaxies, or merge with other clusters, has changed the conceptual framework of the field. Continuing infall of spirals into large clusters can account for diverse disc types in clusters, and provides the opportunity for more precisely constraining the stripping process. For example, Phookun and Mundy (1995) presented a well resolved HI map of the Virgo spiral NGC 4654, which showed a very asymmetric HI distribution, with the NW side compacted, as expected for a bow shock, and a long tail on the SE side. They argue that these asymmetries are the result of the combined effects of ram pressure and rotation if the cluster ’wind’ is not face-on to the disc. Indeed, we find a very similar morphology in inclinded disc models we present below. Vollmer et al. (1999) find similar, asymmetries in the Virgo disc galaxy NGC 4548. Both objects have kinematic asymmetries as well.

Kenney and Koopman (1999) also found morphological distortions in their H$`\alpha `$ imagery of the Virgo spiral NGC 4522. This is also true of the HI and radio continuum emission, according to van Gorkom (2000, private comm.)). Here the cluster wind (with a projected velocity of 1300 km/s) seems to be nearly face-on to the disc, which is observed nearly edge-on. The SE side appears compressed, while extraplanar filaments of gas emission are found on the NW side. (Also see Gavazzi et al. 1995, and Chaname,́ Infante, & Reisenegger 2000, for more objects with similar morphologies.)

Asymmetries are also common in the cluster galaxy HI surveys cited above, and Rubin, Waterman, & Kenney (1999) find, from optical spectroscopy, that kinematic disturbances are common in Virgo disc galaxies.

These observational results stimulated new (three dimensional) modeling studies of sweeping from disc galaxies, including: tree-SPH models of Abadi, Moore, & Bower (1999), sticky particle models of Vollmer (2000), and simulations made with a high resolution Eulerian finite difference code by Quilis et al. (2000). With an extensive, pre-existing gas disc, the continuum, hydrodynamic limit should be a good approximation for calculating the formation and structure of the bow shock. Abadi et al. and Quilis et al. both argue that stripping from disc galaxies is very prompt. Abadi et al. find that up to 80% of the gas of a spiral galaxy can be lost in about $`10^7`$ yr., in the core of a dense cluster like Coma.

Vollmer also found that stripping occurred in prompt events. Vollmer’s models were run with a time-varying wind to represent passage through different parts of the cluster. They found that for galaxies on radial orbits, a great deal of material could be pushed out into long tails in the wake, but a significant amount of this material falls back as the galaxy moved out of the cluster core (especially for higher angles between the wind and the disc rotation axis). A corollary of this is that the details of stripping depend strongly on the type of orbit a galaxy pursues through the cluster (as Takeda et al. concluded for spherical galaxies). Vollmer use a rather ad hoc prescription for the wind pressure, which only approximates pressure effects, so detailed results should be treated with caution.

Abadi et al. carried out a grid of runs in which constant ICM winds impact multicomponent model galaxies (bulge/disc/halo) and compared them to semi-analytic stripping criteria like those of Gunn and Gott, but generalized to multicomponent galaxies. Quilis et al. use a finite difference code to achieve higher spatial resolution, but present only a couple simulations of sweeping from galaxies with central gas holes.

Our simulations were generally run for a time of $`6\times 10^8`$ yrs., comparable to those of Vollmer, but those of Abadi et al. were run for only slightly more than $`10^8`$. Although they emphasize that most of the stripping has already occurred by this time, their figures show some disc gas remains nearby. They also show slightly increasing gas mass in the disc at the end of some runs. Thus, some fallback or removal must be completed on a longer timescale, as in Vollmer’s and our models.

More generally, we find that the dynamics of stripped gas is quite complex on timescales of a few times $`10^8`$ yrs. (like Farouki & Shapiro, and Kritsuk). We will also show that there are important evolutionary processes within the disc itself on this longer timescale. In particular, while the initial stripping of the outer disc occurs rapidly, we find that the internal disc instability facilitates continued ablative sweeping. Radiative cooling, or at least the absence of compression heating is an important part of this instability. We present new analytic models to help understand this multi-step stripping process. We also examine the evolution of the angular momentum of the gas disc, as another tool for understanding the complex dynamics.

## 2 STRIPPING THEORY AND A SIMPLE ANALYTIC MODEL

Gunn and Gott (1972) assumed that gas in the bulk of the disc was bound primarily to the local disc potential, and that in the case of a face-on wind, removal of gas to a distance of one vertical scale height was equivalent to stripping. Abadi et al. updated this criterion by adding bulge and halo gravity to the disc gravity, but like Gunn and Gott they considered only the local balance between gravity and ram pressure within the disc to derive a stripping radius.

In the case of a face-on wind impacting a disc bound by a massive halo, the gravitational force of the halo does not oppose the wind pressure initially, because the force vectors are perpendicular. As a gas element is swept vertically out of the disc (henceforth, the z-direction), to a distance of several scale heights, the disc gravity diminishes. On the other hand, the projection of the radial halo gravity vector into the z-direction, and opposite the ram pressure increases. This halo component can grow quite large, so that the greatest resistance to sweeping can be found at some distance from the disc.

In fact, because of this simple geometric effect, the halo gravity can exceed the ram pressure force, and bind the gas to a region substantially offset from the initial disc. The region is like a potential well that has been displaced by the addition of the ram pressure force. This effect, together with possible shadowing from the ram pressure, allows gas elements to hang up at some distance behind the original disc, rather than being completely removed.

The trajectories of some clouds in Kritsuk’s (1983) models show that they experienced this effect. Indeed, the dashed curve in Figure 1 of Kritsuk’s paper would seem, by its definition, to be part of the boundary of the displaced binding region.

The following simple analytic model illustrates these effects. Assume the galaxy disc is face-on into the wind, with the wind flowing from positive z values to negative. Then the ram pressure force on a gas element or cloud in the disc is,

$$F_w=\rho _wv_w^2A_{cl},$$
(1)
where $`\rho _w`$ is the mass density of the wind, $`v_w`$ is the velocity of the wind in the z-direction, and $`A_{cl}`$ is the cross-sectional area of a gas cloud in the galaxy.

For simplicitiy, assume that the halo is the only important source of gravity. Then the component of gravitational force on the cloud in the z-direction is,

$$F_h=\frac{GM(r)m_{cl}}{r^2}\frac{z}{r}=\frac{m_{cl}v_c^2z}{r^2},$$
(2)
where M(r) is the mass contained in a region of radius r, $`m_{cl}`$ is the mass of the cloud, z is the vertical distance, and $`v_c`$ is the local circular velocity.

In the second equality of equation (2) the halo mass interior to the cloud radius has been formally eliminated in favor of the local circular velocity, which we assume is constant in the region of interest. We equate the two forces to derive not a “removal radius”, but rather a curve of r versus z, which when rotated around the axis of symmetry gives the boundary surface of the displaced binding region, as follows:

$$\frac{z}{r}=\frac{\rho _wrv_w^2A_{cl}}{m_{cl}v_c^2}.$$
(3)
Generally, this an ovaloid region.

In terms of the initial orbital radius of the cloud in the disc, R, this becomes

$$\frac{zR}{r^2}=\frac{\rho _wRv_w^2A_{cl}}{m_{cl}v_c^2}=W(R),$$
(4)
where $`r^2=z^2+R^2`$. Let $`r^{}=r/R=1/sin\theta `$, $`z^{}=z/R=tan\theta `$. Then,

$$\frac{z^{}}{1+z^2}=W.$$
(5)
Solving for $`z^{}`$, we get,

$$z^{}=\frac{1}{2W}\pm [\frac{1}{4W^2}1]^{1/2}.$$
(6)
we need $`W<1/2`$ for a real solution, if $`W>1/2`$, for a given R, there is no binding region. If $`W<1/2`$ for some R, the cloud will move to the bound region, and probably pursue a complex trajectory determined by local pressure and shadowing forces.

Note that the $`R/v_c`$ is approximately equal to the local free-fall time. The remaining terms in the middle of eq. (4) can be written as $`(\rho _wv_w^2A_{cl})/(m_{cl}v_c)`$, which is the ram pressure force on the cloud divided by the cloud momentum, i.e., the inverse momentum change timescale. Thus, W is the ratio of the gravitational restoring time, $`\tau _{ff}`$, to the momentum change time due to wind pressure, $`\tau _{mom}`$.

## 3 SIMULATION CODE AND INITIAL CONDITIONS

### 3.1 Simulation Code

Our simulations were produced with the serial code Hydra 3.0, which has been made publicly available by H. Couchman, P. Thomas, and F. Pearce. The Hydra program implements smoothed particle hydrodynamics (SPH) and calculates gravity with an adaptive particle-particle (PP), particle-mesh ($`AP^3M`$) algorithm (for details see Couchman et al. 1995, Pearce & Couchman 1997). For a typical timestep in our models, adaptive refinements were carried out on about 1000 gas particles in the grid center, with a scaling factor of 2.67.

The simulations were all run using an adiabatic equation of state. Optically thin radiative cooling was calculated via the tables of Sutherland and Dopita (1993), which were supplied with the Hydra code. Cooling times were not used to limit the size of the computational timestep, since the dynamical time is usually longer than the cooling time. Particles evolve adiabatically and are only cooled at the end of a given timestep, at constant density.

We note for completeness that the Sutherland and Dopita cooling curves include atomic and ionic line and continuum processes for $`T10^4K`$. These curves were also calculated at several metallicities, and these alternate curves are included as options in the Hydra code. However, all of our calculations were carried out at solar metallicity. In the range $`T=10^410^6`$ the cooling is dominated by H and He line and recombination cooling, and so, is relatively insensitive to the metallicity.

Nonetheless, this is a very approximate treatment of the complex thermal physics involved in the interaction of hot ICM gas and the multi-phase disc gas. Processes such as X-ray photoheating, molecular cooling and the thermal effects of grain destruction are neglected. However, these effects are most important at $`T<10^4K`$. At higher temperatures it is likely that radiative cooling (e.g., near the peak of the atomic cooling curve) dominates compressive and shock heating in cool-to-warm gas elements originating in the galaxy disc, even when these elements are swept out of the disc (see Section 5.1). Our treatment does capture this basic effect, and confirms simple rate estimates. Several of the observational works cited in the introduction have discovered extradisc HI gas, including long filaments, apparently stripped from the cluster galaxy. Abadi et al. (1999) and Quilis et al. (2000) did not include cooling in their models, arguing that it is unnecessary since their codes cannot resolve molecular clouds. Vollmer’s (1999, 2000) sticky particle code is intrinsically isothermal. We agree that dense clouds cannot at present be resolved, but our disc gas particles can represent diffuse atomic gas and large cloud envelopes. Several of our results on the evolution and kinematics of stripped gas depend on cooling, as we discuss below.

### 3.2 Initial Conditions

In each of the computational runs, we start with a model disc that contains 29100 particles, divided up as follows: 9550 gas particles, 9550 disc star particles, and 10000 dark matter particles. The masses for each of the particle types, in code units, are: $`3.0\times 10^5`$ for gas and star particles, and $`6.0\times 10^4`$ for dark matter particles. The mass unit is $`1.0\times 10^{10}M_{}`$, so each star or gas particle represents $`3.0\times 10^5M_{}`$. This gives a total galaxy mass of $`6.57\times 10^{10}M_{}`$, a relatively small galaxy.

The model galaxy was built first by initializing the halo particles in an isothermal distribution, with an outer radius of about 0.3 grid units, and then run for several thousand steps to allow relaxation. Stellar and gas disc particles were initialized on circular orbits with a small, three-dimensional random velocity component added. The model was relaxed with no ICM wind after the addition of each component.

The model galaxy is placed in a grid ranging from 0.0 to 1.0 in units for each direction, and centred at (0.5,0.5,0.5) in (x,y,z). The code unit of length is about equal to 100 kpc. The gas disc has a radius of 15 kpc. The radius of the stellar disc is about 10 kpc. in x and y, and the dark halo is spherical with a radius of about 10 kpc.

When this initial galaxy is read into Hydra, it is surrounded with 80000 gas (ICM) particles. These particles were randomly dispersed in 50 layers of 1600 particles each. The layers were spaced 0.0196 units apart (1.96 kpc.) in the z-direction, so as to create a uniform ICM. This ICM particle spacing is also comparable to the cell size of the $`64^3`$ large-scale grid. As will be seen below, these initial conditions do allow us to (minimally) resolve the large-scale bow shock around the galaxy. Each of the ICM particles was given a mass 10% of that of a gas particle belonging to the galactic disc. In using this reduced mass we follow the recommendations of Abadi et al. (1999), who point out that if the mass is too high, the ICM particles will punch through the galaxy disc.

In these models the temperature of the initial ICM was set to 5.0 units, of $`4.6\times 10^5`$ K, so $`T=2.3\times 10^6`$ K. In addition, the density of the wind is approximately $`7.3\times 10^5cm^3`$. These values are typical for a moderate galaxy cluster.

In order to simulate the galaxy traveling through the ICM, we kept the galaxy stationary, while having the intracluster plasma flow past the galaxy as an effective wind. This flow was directed from +z to -z. A control simulation was also run with no wind, for a time equal to that of the other runs. No discernable change was seen in the model galaxy in this run, including no spiral or bar instabilities.

### 3.3 Boundary Conditions

The following boundary conditions were employed. First, in the x, y, and z directions, the particles are required to remain between 0.01 and 0.99 in code units. These conditions were checked every few timesteps. This is done to maintain efficient computing performance while still ensuring that few particles were lost from the grid (in a typical run, less than 0.3% of the intracluster gas particles). In addition, when particles reach the lower bound of $`z=0.01`$, they are randomly placed back into the top of the ICM flow, and given a mass 1/10 that of the gas particles in the galactic disc.

Thus, some of the particles put back into the ICM were formerly gas particles from the disc. When they are brought back into the grid, they are given the smaller ICM particle mass. So, in determinations of which gas disc particles are still on the grid, these particles are not included (for further details see Schulz (2000)).

## 4 MODEL RESULTS

In this section, we present the results of our new simulations, with further analysis in the following sections. We begin with a face-on ’fiducial’ model. Next we consider models in which the galaxy is tilted relative to the ICM wind, then a model with a larger wind flow velocity, and finally models in which the temperature of the disc gas was initially set to a high temperature, but allowed to cool, to crudely model the effects of a strong global starburst. The parameters of each of the data runs are shown in Table 1.

### 4.1 Fiducial Model

The fiducial model is a face-on case with a wind velocity of -977 km/s with the adopted scaling. We ran the code for 18000 timesteps, equivalent to a time span of $`5.64\times 10^8`$ yrs. with this scaling. At -977 km/s, a particle starting at the top crosses the entire grid over 5 times by the end of the run.

In Figures 1 and 2, we show the particles from the gas disc (only) at various timesteps over the evolution of the model. These figures bear a qualitative resemblance to Figure 12 of Farouki & Shapiro (1980).

During the early evolution of the model (Figure 1), the outer disc deforms rapidly. Mass loss begins as particles move towards the bottom of the grid. By the time of the second snapshot spiral waves appear in the outer disc, and a ring of compressed gas separates the inner and outer discs. This ring also divides the planar and nonplanar regions.

A clear bow shock also forms early in the ICM. This is shown in Figure 3 after about $`10^8`$ yrs. For clarity, wind particles above 0.6 code units in z have been removed, and only a section ranging from 0.425 to 0.575 units in y is shown in Figure 3. The opening angle of the shock is quite broad (see discussion below), so that the shock extends across most of the x-axis at the bottom of the grid. Very few particles are visible behind the bow shock, signifying that few, if any, wind particles are able to pass through the galaxy unimpeded. The bow shock visibly compresses the disc in the vertical direction.

The temperature contours for the bow shock are also shown in Figure 3. The numbers are in the code units described in Section 3.2. Postshock heating, as well as the heating around the galaxy disc gas by the reflected shock, is visible. Radiative cooling keeps gas within the disc from heating to these million degree temperatures.

Figure 2 shows that mass loss tapers off at the end of the run. Low-density material is removed first. This consists of gas found between spirals or in low-density spirals. Higher density spirals are stripped later, and then tend to “hang up,” or in some cases, even fall back onto the disc. These phenomena are predicted by the analytic model, and can be seen in Figure 4. The filamentary appearance of the stripped spirals is reminescent of HI filaments observed in several galaxies, as described above.

In particular, the last x-z plot in Figure 2 also shows the existence of a tail of gas extending out beyond the galaxy a distance of at least 20kpc. The development of this tail from a disc can be discerned in the preceding snapshots.

The shock angle, as measured relative to the bottom of the grid, varies over time (ranging from 60 to $`75^{}`$). These changes are due not only to stripping but to vertical and radial oscillations of the disc following the appearance of the bow shock. These changes can most easily be seen in the position of the shock at the bottom of the grid in images like Figures 1, 2, and 4.

In order to better understand the evolution of stripped particles, we examined small, arbitrarily defined gas clouds (groups of particles) at various times. Each of the clouds that we looked at was identified at 3159 timesteps and followed thereafter. The motions of individual clouds are quite complex and dynamic. An example of such a cloud is shown in Figure 5. Initially, it consisted of 129 gas particles. By 18000 timesteps, only 42 were left in the grid. Of the particles that were left, most appear to have been accreted back onto the disc. This is, however, a rather exceptional example. Most clouds in the outer disc are eventually entirely removed.

We examined both the mass entirely removed from the computational grid, and the mass lost from a fixed volume enclosing the original disc, as a function of time (see Figure 6). This volume is defined as a rectangle that goes from 0.3 to 0.7 grid units in x and y, and 0.425 and 0.575 units in the z-direction. The figure shows that the disc lost the most material within the first $`10^8`$ yrs. This result is roughly consistent with that of Abadi et al. (1999).

After about 375 million years only a few particles are still being stripped off the disc. On the other hand, after 100 million years, the rate of mass loss from the grid as a whole remains fairly constant. The difference between the two curves lies in the fact that some particles linger in the grid for long timescales, while other particles are slowly accelerating away from the effective potential minimum.

Figure 7 is a plot of total angular momentum versus time for particles that originated in the gas disc. The falloff in angular momentum over time appears quite similar to that of the gas mass. The ratio of mass loss to angular momentum loss rises rapidly within the first 100 million years but then settles down to a value of approximately 0.5, meaning that about twice as much angular momentum is being lost as gas mass. This is because the material with high specific angular momentum spreads the most, and thus, is easiest to strip. This high rate of angular momentum loss accounts for most of the compression we see in the disc with time.

We found, somewhat to our surprise, that the stellar disc, too, appears to be affected by the intracluster medium. While the stellar disc was not altered in size or shape, by the end of the run it had been displaced roughly 0.02 length units in the negative z-direction, which translates to $`6.17\times 10^{16}`$ km, or 2 kpc. We believe that this displacement is an indirect result of the bulk displacement of the gas disc which is not stripped. The mutual gravity of the gas and star discs, which are of equal mass in our models, prevents them from separating, and thus, pulls the stellar disc down relative to the fixed grid.

We estimate that the net kinetic energy that is transferred from the wind to the gas disc exceeds the disc self-gravity. Thus, it is adequate to account for the displacement.

The effect provides an interesting, if rough, check on the numerical model. Since the outer disc will be most affected by the wind force, we consider the vertical wind force acting an annular region in the outer part of the disc. That is, we essentially view the outer disc as analogous to a solar sail which is dragging its payload, where the payload is the stellar disc. Vertical force balance is achieved when

$$\rho _wv_w^2A_{ann}=\frac{GM(r)m_{ann}}{r^2}\frac{z}{r}.$$
(7)
Here, $`\rho _w`$ is the wind mass density, $`v_w`$ is the wind velocity in the z-direction, $`A_{ann}`$ is the area of the annular region in x and y, M(r) is the disc mass contained in r, $`m_{ann}`$ is the mass of particles in the annulus, r is the distance from a particle in the annulus to the potential center, and z is the vertical displacement of the disc.

We can substitute in $`v_c^2/r`$ for $`GM(r)/r^2`$, where $`v_c`$ is the effective centripetal velocity. We find that,

$$z=\frac{\rho _wv_w^2A_{ann}r^2}{m_{ann}v_c^2}=\frac{\rho _wv_w^2}{\rho _{ann}v_c^2}\left(\frac{r}{2h}\right)^22h,$$
(8)
where to obtain the last form we use the relation $`\rho _{ann}=A_{ann}(2h)`$, with disc thickness, 2h. With values from the fiducial model, we calculate that $`zh`$ in code units. Since $`h0.02`$, this estimate agrees well with the numerical result. The displacement of the disc is greater in a higher flow velocity model, as the above equation predicts.

### 4.2 Tilted Models

Next we describe two simulations run to study the effect of inclination on gas sweeping. In the first one, the galaxy was tilted by $`40^{}`$ from the z-axis. Six timesteps for this run are displayed in Figures 8 and 9. One noteworthy aspect of these plots is that, after 18000 timesteps ($`5.6\times 10^8`$ yrs.), almost all of the gas from the disc is either still in the disc or has been removed from the grid. Much less remains between the gas disc and the bottom of the grid than in the fiducial model. Another interesting feature is the appearance of a large “tail” of material at timestep 9200 ($`2.9\times 10^8`$ yr.). It includes particles ranging from 0.35 to 0.7 units in x, 0.25 to 0.5 in y, and 0.2 to 0.42 units in z. From step 9200 on, it can be seen moving down and to the right in the x-z plots. At step 18000, only a small portion of the original tail remains in the region from 0.6 to 0.8 in x, and from 0.0 to 0.15 in z.

This tail can first be seen at step 6000 as a strong spiral arm in the x-y view of Figure 8. Subsequent views in the two projections reveal that as the arm is rotated into a position where it is moving against the wind, it is stripped. This seems to be an example of the kind of rotationally aided stripping postulated by Phookun & Mundy (1995), and mentioned in the introduction.

Figure 10 shows the temperature contours for the bow shock. Although the temperatures are comparable to that of the fiducial model, the contours show that the heating is much less symmetric than before. Note, however, that the contouring was performed with a coarser mesh resolution than the simulations themselves. Granted this, and the limited particle numbers immediately downstream from the disc, the contours there are uncertain.

Figure 11 shows the z-velocity contours on the x-y projection of the disc at $`10^8`$ yrs. into this run. In the inner disc, within the ring, we see a more or less normal pattern of projected rotation, with only small asymmetries. Outside the ring, the asymmetries are large. At the top of the figure the disc rotates with the wind. There the velocities reach high amplitudes, and the contours are compressed. At the bottom the rotational flow opposes the wind, and velocity contours are stretched out. There is a general similarity to the kinematic observations of NGC 4848 presented in Vollmer (2000). The kinematic continuity between disc and stripped gas, which has been noted by a number of observers is also apparent in Figure 11. This is also apparent in Figure 12, which shows contours of y velocities in the x-z plane.

Overall, roughly the same number of particles were removed in this model as in the fiducial model (about 1000 particles over about 7600 timesteps). Figure 6 shows that the initial mass loss is not as rapid as the fiducial model, but by the end of the run, it is comparable.

The total angular momentum loss over time is shown in Figure 7). It is interesting that after about $`3.75\times 10^6`$ yrs., the quantity around the disc is actually larger than that on the whole grid! The explanation for this is that a large cloud of particles (the tail noted above), was counter-rotated, and then stripped, as a result of an interaction with the ICM.

In the second tilted simulation, the galaxy was inclined $`60^{}`$ about the z-axis. Generally this model is similar to the previous one. However, the increased tilt narrows the bow shock. By the end of the run, a large part of the gas disc has settled into the ring-like structure that we have seen in the previous models. A one-arm spiral develops, and is stripped much as in the previous run. About 500 fewer particles have been removed than when the galaxy was inclined $`40^{}`$.

### 4.3 Increased Flow Velocity Model

Next we describe a simulation with twice the flow velocity of the previous ones (-1953 km/s). Like the others, this model was run for 18000 timesteps. With the usual Courant-Friedrichs-Lewy (CFL) condition (see e.g., LeVeque et al. 1998, 278) the numerical timestep was reduced to about half its previous value. Thus, the total run time was about 288 million years.

Figures 13 and 14) show snapshots of the gas particles from the disc. Although the plots are quite similar in appearance to that of the fiducial model, the sweeping process occurs on a much shorter timescale.

About 1000 more particles were lost in this case than in the fiducial model, and in about half the time, see Figure 6. Here, the mass loss around the disc is much more rapid within the first 100 million years, but beyond that it essentially levels off. The ICM parameters of this run are close to those of Run C of Abadi et al. (1999, though our model galaxy is about 4 times less massive), and the resulting mass loss curve is also very similar. The total angular momentum loss is plotted in Figure 7.

Higher wind velocity produces more shock heating, and thus, higher temperatures than before. In this model they range from 14.7 million K to 29.5 million K.

We also observed the motions of individual clumps of particles over time in this model. For nearly every one of the clouds we looked at, all of the particles were removed from the grid by the end of the simulation. This is consistent with the sweeping criterion, in terms of the W parameter, discussed above. W is proportional to the square of the wind velocity, so by doubling the velocity for this model, we have effectively multiplied W values by 4, increasing the likelihood of stripping. More precisely, the volume of the “trapping region” in the halo, with $`W<1/2`$, is very small.

### 4.4 Models with Secondary Heating

We also tested the effect of impulsive internal heating in the galaxy disc, for example, heating as a result of a global burst of star formation. Specifically, we carried out two simulations in which the gas disc of the fiducial model was impulsively heated after 6000 timesteps. All other parameter values were kept the same as in the fiducial model.

In the first case, the gas disc was heated to approximately 1.4 million K. However, this temperature was still low enough that the galaxy was able to cool rapidly, and there were no noticeable differences in gas removal from the fiducial model.

In the second case, the gas disc was heated to about 3.5 million K. This heat was not immediately radiated, and the gas particles in the galaxy puffed out vertically over a timescale of less than 94 million years. By a time of $`2.9\times 10^8`$ yrs., 1100 more particles (11%) were removed than in the case without heating. After that, the effects were the same overall as in the fiducial model. We conclude that enhanced star formation may lead to only moderately enhanced stripping on cluster crossing timescales, and substantial (or prolonged) heating is required.

## 5 DISCUSSION

### 5.1 The Deterministic, Multistage Sweeping Process

There is a great deal of consistency in the results of the different runs described above. It appears that the stripping/sweeping process proceeds through the same distinct stages regardless of galaxy tilt, or ICM velocity and density, at least over a range of an order of magnitude in these parameters. Based on the model results above and previous work, we suggest that the process has the following characteristic stages.

(1) Prompt stripping of the outer disc gas from the disc plane. The effectiveness of this Gunn/Gott process has been emphasized in the studies of Abadi et al. (1999) and Quilis et al. (2000). Promptly stripped material is lightly bound to the disc plane, but it remains bound to the dark halo initially.

(2) Specifically, the stripped gas hangs up near the minimum of the effective potential downstream from the disc, as described in Section 2. Thus, stripping from the disc is not always equivalent to sweeping out of the galaxy.

Cooling plays an important role in the hang-up phenomenon. Warm clouds of a given mass, in a constant pressure environment, will have a smaller cross section than an equal mass of hot gas, $`A_{cl}\rho _{cl}^{2/3}T_{cl}^{2/3}`$. In our models most of the hung-up gas has a temperature of about $`10^4K`$. While it is heated by interaction with the ICM, this heating is not sufficient to push it over the peak of the cooling curve.

(3) Simultaneous with the first two stages is the formation of a bow shock in the ICM, and also of a reflected shock around the gas disc. Until recently models have not provided much information about the three dimensional structure of the shocks, although high resolution two dimensional models, like those of Stevens et al. (1999), are very suggestive.

The figures in Abadi et al. (1999) show the presence of bow shocks, the model galaxies of Quilis et al. (2000) have a large central hole, and as a result an annular bow shock forms with wind flow through the centre. These latter authors also contend that the wind flows through kiloparsec scale holes in the disc ISM, and that this enhances prompt viscous stripping. We find a similar effect in our models, in the efficient stripping of gas between dense spirals in the outer disc. As described in Section 4, we also find that the bow shock expands and contracts on large scales, following changes in the cross section of the gas disc (which are due to oscillations, stripping, and angular momentum evolution).

Despite holes and gaps, the global bow shock still forms in all the simulations above, and compresses the disc. Given the fact that the clumpy, multiphase nature of the ISM is not fully resolved, we might wonder if this result is correct. E.g., if all the gas is promptly swept except the cores of the largest clouds, then a global shock would probably not form. However, there are several reasons to think that this will not occur. Firstly, Crosthwaite, Turner & Ho’s (2000) recent high-resolution (about 350 pc) HI map of the nearby spiral IC 342 shows that almost all of the disc is covered with diffuse, resolved atomic gas with a column density of $`3\times 10^{20}cm^2`$, which is substantial. Studies of the diffuse ISM in galaxies suggest that this component is also enveloped in hotter ionized or partially ionized phases of significant column density (e.g., the review of Reynolds 1996 on this component in the Milky Way). In addition the whole is threaded and connected by magnetic fields. Thus, it appears that the column density of the diffuse ISM is generally large enough to prevent the ICM wind from promptly sweeping it, and leaving bare cores.

It also seems unlikely that the ICM could sweep through many small, as yet unresolved holes. Such holes would be connected to their surroundings by magnetic fields, and as pointed out by Crosthwaithe et al. they would be closed very quickly by the effects of shear. Indeed, the shear timescale is comparable to the timescale for the wind to flush the holes. Therefore, except for large central holes and gaps between strong spirals, it seems likely that the disc will respond as a cohesive whole, whether it is promptly swept, or only compressed.

(4) Gravitational instability is strongly suggested by the prompt formation of flocculent spirals. To confirm this, we have computed the value of the Toomre Q factor in annular rings in the fiducial model at various times. We find that $`Q3`$ at most radii in the initial model, but after the onset of the wind, the Q value typically drops to less than or about 1.0 at these radii.

Interestingly, the cause of this instability is the displacement of the disc relative to the halo potential center, which was described in Section 4.1. We can understand the triggering of the instability in terms of the simple gas disc plus halo model used in Section 2. The gravitational restoring force balances ram pressure in the midplane of the displaced disc. The side of the disc closest to the halo center of mass (facing the wind) feels an excess ram pressure. The opposite side of the disc is farther from the halo center and feels an excess restoring force. As a result the disc is squeezed from both sides, in what amounts to a tidal compression. These forces are equivalent to the weight of added mass in the disc, so their effect can be included in the usual gravitational instability criterion as an ’effective’ gas surface density that is larger than the true surface density.

Assuming that the tidal compression is balanced by a pressure increase within the disc, we can derive the following expression for the fractional pressure increase,

$$\frac{\mathrm{\Delta }P}{P}2\frac{\rho _wv_w^2}{\rho _gc^2},$$
(9)
where c is the sound speed of the disc gas, and $`\rho _g`$ is its mean density. The factor of 2 comes from the usual tidal force equation, and we have also assumed that the disc displacement distance z is comparable to the disc thickness h. If $`z>>h`$, then a factor of h/z should be inserted in the right hand side of the equation.

If we further assume that the disc gas is isothermal, then the fractional increase in the effective surface density is also given by the above expression. With the following representative values: $`\rho _w/\rho _g=10^4,v_w1000km/s,c5km/s`$, we get a fractional increase of order unity. The gas surface density in most galaxy discs is close to the critical value for gravitational instability (i.e., Toomre $`Q1`$), so we expect that doubling the effective surface density will halve Q, and induce a strong instability, as seen in the simulations.

This instability could not be discovered in cylindrically symmetric models. However, it has also not been noted in previous three dimensional simulations, which did not include radiative cooling. Heating by the wind would increase the sound speed and the equation above, reduce disc compression, and prevent gravitational instability. This prevents annealing (see below) and aids prompt stripping.

The cooling/heating balance is sensitive to both galactic and ICM parameters. We can demonstrate this with a simple estimate. We take the kinetic energy flux of the wind ($`\rho _{ICM}v^3`$) as the maximum heating per unit area. We equate this to the cooling rate at the peak of the standard ISM (atomic) cooling curve, multiplied by the disc thickness (i.e., the cooling per unit surface area of the disc, $`n_{ISM}^2\mathrm{\Lambda }_{max}(2h)`$). Using a disc thickness of $`2h1.0kpc.`$, the wind density and velocity of our fiducial model, and $`\mathrm{\Lambda }_{max}5.0\times 10^{22}ergscm^3s^1`$, we find that cooling dominates whenever $`n_{ISM}>0.02cm^3`$. Within our fiducial model this critical density is exceeded by factors of a few to ten in the gas disc. Thus, energy input from the wind is never able to heat the disc through the cooling peak. Moreover, once a bow shock forms, the wind does not directly impact the disc, and so the heating is likely to be less. Even stripped filaments typically exceed the critical density, and they remain at $`T=10^4K`$ in our models.

Note that the critical density is especially sensitive to the wind velocity. E.g., the critical density would increase to $`n_{ISM}1.8cm^3`$ with the ’Coma’ wind parameters of Abadi et al., while their ’Virgo’ wind is close to our high wind speed model. In the case of such a high critical density, heating would generally dominate, except perhaps, in the densest inner disc. When heating dominates we would expect that it would facilitate the stripping of the diffuse ISM, as described below.

(5) When gravitational instability occurs, the spiral waves enhance the radial tranport of angular momentum in the disc. The additional angular momentum deposited at large radii drives expansion, the resulting column density reduction facilitates stripping. The loss of angular momentum at smaller radii leads to compression. The ’annealed’ inner disc within the ring is resistant to stripping. The timescale for the combined gravitational instability and annealing is about $`(13)\times 10^8yrs.`$, i.e., a few disc dynamical times.

(6) The last stage observable in our models is the continued mass loss as flocculent spirals in the outer disc shear and expand until their column density is reduced to the point where they are stripped, as described in the previous section. Here we merely emphasize the fact that since the column density is just low enough for the material to be lifted, these filaments are very likely to become hung up in the halo. Most of this material will eventually be swept, but some falls back onto the disc. As shown by Vollmer (2000), this is especially true when the lingering time is comparable to the cluster crossing time for galaxies on radial orbits. In that case, the wind dimishes before the filaments are swept.

In addition to these sweeping phases, which are clearly discernable in our simulations, we can speculate about the role of several additional processes that were not included in our models. The first of these is enhanced star formation in the annealed disc. The annealed disc has a somewhat higher gas density throughout, but especially in the ring, so star formation is a natural consequence. Our models with impulsive heating show that even a burst of star formation may not lead to much more stripping. However, the energy input and subsequent enhanced dissipation may drive some disc spreading.

Since we find annealing to be most important at ram pressures of order or somewhat lower than those of the Virgo cluster, annealing induced star formation may be most important in poor clusters, and young, high redshift clusters. At the same time, annealing requires some minimum ram pressure to induce gravitational instability, so annealing induced star formation is probably less important in galaxy groups than tidally induced star formation.

One of the best candidates for an annealed disc that we are aware of is the Virgo spiral NGC 4580. Its image in the atlas of Sandage & Bedke (1994) is remarkably similar to snapshots of the fiducial model on comparable scales, at late times. Sandage & Bedke classify it as type Sc(s)/Sa (!), and state that it is “peculiar enough to be outside the classification system.” It appears to have knots of young stars thoughout, but concentrated in the broken ring. Incidently, expanded views of the fiducial model at late times, show that it’s ring is also broken, and is in fact, a spiral that wraps tightly at that radius. Other examples of possible annealed morphologies in Virgo can be seen in the paper of Koopmann, Kenney and Young (2001).

A most impressive set of candidate objects is seen in Figure 11 of Oemler, Dressler, and Butcher’s (1997) study of the four rich clusters at z=4. Apparently recent star formation reveals these blue optical rings.

Star formation and nuclear activity may also be induced by the fallback of gas clouds stripped from the disc, but not swept out of the halo. Because the mass involved is generally not large, we would not expect this process to dominate the global SFR of the galaxy disc, but the mass could be more than sufficient to feed an AGN.

### 5.2 Annealing Versus Other Instabilities

As explained in the previous subsection, the annealing process depends crucially on an induced gravitational instability, and on efficient angular momentum transport mediated by the spiral waves formed in the instability. In this subsection we address two important questions about these conclusions. The first is, are we sure that this gravitational instability is the result of ram compression as stated above? The second is, are there other dynamical instabilities that might mimick the characteristic gas morphologies formed by annealing. That is, are those morphologies unique?

In Section 3.2 we noted that a run of the model galaxy with no cluster wind, confirmed that the model disc remains stable in isolation for times longer than the characteristic dynamical time (e.g., a free-fall time). There is no evidence for intrinsic gravitational instability. There was also no significant change in the global disc angular momentum in this control run.

This result, together with the consistency of the annealing phenomenon in the different model runs described in the previous section, and the analytic theory, provide us a good deal of confidence that we understand the annealing process. However, there remains the possibility that the disc is near the gravitational stability threshold, so that any significant disturbance, including modest transients, could trigger the onset of instability. In this case the initial wind impact, could be the accidental instability trigger, rather than the effective weight of steady compression.

We are skeptical of this possibility because we have used this code with the same initial galaxies in other applications (including models of galaxy collisions) without triggering internal instabilities. However, we have also carried out a test simulation designed to search for internal instability triggered by transient finite amplitude disturbances. Specifically, we increased the mass of all particles by 5%, and introduced a pseudo-tidal velocity disturbance, also with a 5% amplitude. The latter consisted of an outward velocity impulse added to the x-component of each velocity vector, and an inward part added to the y,z velocities, with the amplitude of each added velocity modulated in latitude as in a pure tidal perturbation.

As expected this disturbance generated spiral and ring waves in the disc, and these propagating waves persisted through the run. Generally, these wave morphologies look very different from those of annealed galaxies. However, at the timestep show in Figure 15, the morphology appears quite similar to those in Figure 1.

In summary, by introducing both a modest tidal and gravitational perturbation we generated a significant response, and found that it can mimick the annealing morphology, though only for short times. The waves in the test galaxy will eventually decay, leaving a disc that is little changed from the initial one. In the annealed galaxies the changes are permanent. This is clear from the mass and angular momentum loss diagrams of the previous section. In the test simulation, no mass is lost, and the angular momentum changes by less than 1% during the run.

The radial distribution of angular momentum does change slightly in the test run. Since we have argued that wind compression is like added weight on the disc, it is likely that by adding enough mass to gas particles we could replicate the annealing instability without a wind. However, the estimates of the previous section show that we would have to increase the disc mass by about 100%, not the much smaller amount used in the test. In any case, it is clear that the gravitational instability in the annealing models is not the result of numerical inaccuracy, or inherent instability of the disc to a modest perturbation.

The test run also suggests that while the morphology of annealed galaxies is not completely unique, good mimicks are rare. It is evidently very difficult to produce this morphology by pure tidal disturbances. Mass transfer or accretion events may be able to do it, but again the appearance of such galaxies will generally be dominated by transients with a very different appearance.

### 5.3 Scaling in Sweeping

Disc galaxies on radial orbits in clusters will experience most of their stripping in crossing the cluster core where the ICM densities and their orbital speeds are highest. In light of the results above and previous work, it appears that in a core traversal stripping can be terminated in one of three ways: 1) in high ram pressure cases most of the gas is promptly stripped, 2) gravitational instability and angular momentum transfer cause annealing, 3) rapid core traversal ends stripping on a timescale shorter than the annealing time. In large, dense clusters the annealing and core crossing time are comparable, but the ram pressure is generally large, so prompt stripping is dominant. In clusters comparable to Virgo or smaller, annealing may be important. Annealing will also be important for galaxies on nearly circular orbits, though in such cases star formation heating may drive secular mass loss.

In most cases, either prompt stripping or annealing dominate, so stripping goes to saturation and ends on a timescale of order a few times $`10^8`$ yrs. Thus, we can speak of discrete mass loss events, and ask whether the mass loss in such events is a regular function of the parameters? Figure 16 shows the percentage of gas mass retained in our simulations (including some additional runs not described in Section 4) as a function of the ram pressure $`\rho _wv_{w}^{}{}_{}{}^{2}`$ in the wind, normalized to that of the fiducial model. This log-log plot shows that the mass retained is in fact a power-law function of the ram pressure, with an index of about 0.21. The plots includes runs in which the either the wind density and velocity were changed, but there is no dependence on these parameters individually, only the ram pressure.

Suppose that the cutoff radius for stripping/sweeping is determined by a universal critical value of the parameter W(R) of Section 2, i.e., $`W(R_{cr})=W_{cr}`$. Since, as noted in Section 2, W is the ratio of the free-fall time to the momentum change time, this seems to be a reasonable assumption. Using the definition of W, we can write,

$$W_{cr}=\frac{\rho _wv_w^2R_{cr}}{\rho _gv_c^2(2h)},$$
(10)
where $`\rho _g2h`$ is the gas column density written as the radius-dependent mean disc density times the disc thickness. In late type spirals, the gas column density generally falls off as $`1/R`$, but we will make the slightly more general assumption that the dependence is of the form, $`R^p`$. We will also assume that the circular velocity in the disc goes as $`R^{m/2}`$, where $`m=0,1`$ for flat and Keplerian rotation curves, respectively.

If we compare identical galaxies in two different cluster environments, with these assumptions, the ratio of the disc cutoff radii in the two cases is given by,

$$\frac{R_2}{R_1}=\left(\frac{\rho _{w1}v_{w1}^2}{\rho _{w2}v_{w2}^2}\right)^{1/(1+p+m)}.$$
(11)
The gas mass out to the critical radius goes as $`M_g(R)R^{2p}`$, so that the ratio of retained gas masses is,

$$\frac{M_g(R_2)}{M_g(R_1)}=\left(\frac{R_2}{R_1}\right)^{2p}=\left(\frac{\rho _{w1}v_{w1}^2}{\rho _{w2}v_{w2}^2}\right)^{\frac{2p}{1+p+m}}.$$
(12)
Thus, this equation predicts that the index of the mass retention versus ram pressure power-law is $`s=(2p)/(1+p+m)`$. For a flat-rotation curve galaxy, with a gas surface density that falls as 1/R with radius, m=0, p=1, and the index is $`s=0.5`$. This is a much steeper falloff than we find in our simulations, yet this is more or less what we would expect for prompt stripping from a typical late-type disc.

In addition, the loss of angular momentum in the annealing process drives radial compression of the gas disc, and reduces the mass loss. This can be viewed as making the effective critical radius in the initial disc larger. The loss of angular momentum increases as the ram pressure increases. Suppose, for example, that the total angular momentum loss increases linearly with the ram pressure. In a m=0, p=1 disc with fixed central surface density, the total angular momentum scales with the total disc area, or the square of the effective critical radius. In this case, the effective critical radius increases as the square root of ram pressure, and the index s=0, i.e., mass loss is independent of ram pressure.

However, the models show a more modest dependence of angular momentum loss on ram pressure, more like a square root dependence than a linear. This leads to $`s1/4`$, which is quite close to the value $`s0.21`$ derived directly from the models. Therefore, we attribute the weak dependence of mass loss on wind parameters to the annealing process.

In Figure 16 we have also plotted results from Abadi et al. (1999), who used ICM densities ranging from $`5.64\times 10^{24}g/cm^3`$, to an order of magnitude higher, and relative velocities of $`10003000km/s`$ (see their Table 2). The model galaxy used in their simulations was not identical to ours, so to normalize to our fiducial model we must multiply by the ratio of the factors $`\rho _gv_c^2(2h)`$ for the two model galaxies (see eq. (9)). We have estimated this ratio in a very simple way. First it appears that the mean gas column densities $`\rho _g(2h)`$ in the two models are comparable, so we have taken them as equal. The quantities $`v_c^2`$ are assumed to scale with the total mass of the model galaxies, which have a ratio of about 4. Serendiptously, this rough normalization gives almost overlapping curves for the two sets of models; we might have expected a greater offset.

The overlap may be due to a partial cancellation of two effects. Prompt mass loss may be aided by compressional heating in the Abadi et al. models, and so, is greater. However, gravitational instability and angular momentum transport drive additional mass loss on longer timescales in our models. Apparently, the end result is not greatly different. The fact that the form of the Abadi et al. mass retained versus ram pressure relation is very similar to ours does not depend on the normalization. The combined results show that this relation extends over two orders of magnitude in dimensionless ram pressure.

More generally, the figure indicates that the mass loss from typical late-type disc galaxies, in a single, saturated, face-on, stripping event, is determined by the dimensionless parameter,

$$S=\left(\frac{\rho _wv_w^2}{\rho _dv_c^2h}\right)/\left(\frac{\rho _wv_w^2}{\rho _dv_c^2h}\right)_{fiducial}=\frac{\tau _{mom}}{\tau _{mom,fid}}$$
(13)
i.e., the ratio of the momentum change timescale of Section 2 relative to that of the fiducial model.

## 6 CONCLUSIONS

In summary, we have presented a small grid of three-dimensional, self-consistent, thermohydrodynamical models of stripping of gas from disc galaxies traversing the ICM of a galaxy cluster. We confirm many of the general results of other recent simulations. However, we have focussed on cooling effects and on models with lower ram pressures than those presented by Abadi et al. (1999) and Quilis et al. (2000). As in those works, we maintained a constant ICM wind on the target disc. We also ran our simulations for a relatively long time, up to $`6\times 10^8yrs.`$, primarily to follow dynamical processes to completion, but also beause this is appropriate for some applications.

We find that stripping is largely deterministic. If the ram pressure is great enough to strip most of the disc promptly, then it is also likely strong enough to sweep the gas directly out of the galaxy halo, and into the general ICM. If the ram pressure is not so great, then a good deal of the gas may become hung-up in the effective potential minimum described in Section 2, for a time of order $`10^8yrs.,`$ before being swept. A small amount of this material falls back into the galaxy. As Vollmer (2000) emphasized, a good deal more would fall back when the wind duration is short, as in galaxies on radial orbits.

The hang-up effect is one type of longer term mass loss, and may account for observations of HI clouds and filaments in the wake of some cluster disc galaxies. In our models this material is usually concentrated within a couple galaxy diameters downstream from the disc. The location of the effective potential mimimum depends on the structure of the dark halo of the galaxy. The dark halos of our model galaxies were quite concentrated.

Another important process in prolonging the gas removal is the gravitational instability and inner disc annealing process described above. Cooling, and a substantial tidal compression are required for this instability. After flocculent spirals form, angular momentum transport and the shearing of these waves occur on comparable timescales. The result is removal of the flocculent spirals from the outer disc; they blow away as gas filaments. The inner disc is radially compressed, and with higher mean column densities it is ’annealed’ against further disruption. The boundary of this inner disc is often marked by a ringlike structure, which may in fact, consist of tightly wound spirals. The compression may induce subsequent star formation.

The primary reason these processes were not discovered in previous simulations is that they are fully three-dimensional, and require cooling. The role of cooling in retaining gas is somewhat counter-intuitive. Stripping criteria do not generally include gas temperature amoung their various dependencies. However, the role of cooling is indirect, it allows the onset of gravitational instability, which promptly generates an annealed or sweeping-resistant inner disc, as described above.

Despite all the complexities, we find that the net amount of sweeping at the end of our runs, when the sweeping has indeed terminated, is a simple power-law function of one dimensionless parameter. This ’S’ parameter, described in the previous section, includes contributions from the ICM ram pressure, the galaxy halo (via the maximum circular velocity), and the gas disc, as might be expected.

These results are modified in cases where ICM wind impacts the gas disc at a significant angle to the disc symmetry axis. It is still true in such cases that there is prompt stripping and downstream hang-up. However, in our models we find more material gets hung-up, and the long term mass loss takes longer than in corresponding face-on cases. Of course, the ram pressure is reduced relative to the face-on case with the same ICM density and relative velocity. However, the differences are not simply explained by a projection cosine. This is illustrated by the coupling between rotation and mass loss discussed in Section 3. This coupling is evidently responsible for part of the time-delay in the mass loss in the angled cases. In addition, the loss in angular momentum for a given amount of mass removal was higher in the tilted models than in face-on models with the same ICM parameters. This explains why higher inclination models appear to have more compact ring structures as well. Nonetheless, the mass loss process in the inclined cases is still proceeds in regular stages.

We conclude with a couple more speculations. The first is that star formation, and the formation of star clusters or small dwarf galaxies, may occur in the gas that gets hung-up in the halo, or that is swept out in relatively dense filaments. It is tempting to interpret the dwarf companion (“North Object”) of NGC 1427A, studied by Chaname et al. 2000, and the dwarf in the wake of NGC 4694 (van Driel & van Woerden 1989), as examples. In the former case, both galaxy and companion have essentially the same velocity as determined by optical spectroscopy. In the latter, there is a peak in the HI wake superimposed on the optical dwarf.

The fate of star clusters, and dwarf companions formed in this fashion is not clear. Small groups or clusters of stars formed in swept filaments may remain gravitationally attached to the filament and pulled out of the galaxy, adding to the population that is not bound to any single galaxy. Dwarf galaxies formed in the effective potential minimum, with little angular momentum relative to the parent galaxy, may separate from their placental gas and fall back into the galaxy.

Our final speculations concern the evolutionary effects of the stripping power-law derived above. It is clear from Figure 16 that as the ram pressure increases from that of our fiducial model to that of Abadi et al.’s Coma model, the mass loss in a stripping event goes from small (25%) to nearly complete. If galaxy clusters complete their growth on a relatively long timescale, via subcluster mergers, then the ICM density and free-fall velocity increase on a comparable timescale. Except in clusters that grow quickly to Coma proportions, we might expect that sweeping was much less efficient in clusters at high redshift. (Note, however, that annealing also requires a minimum wind pressure to generate the compression needed for gravitational instability.)

If ICMs do develop on long timescales, galaxies on radial orbits can experience repeated cycles of sweeping, annealing, and compression induced starbursts. This process could contribute to the Butcher-Oemler effect. There is evidence that the dominant processes driving the Butcher-Oemler effect are galaxy collisions and effects related to the merger of subclusters (e.g., Couch et al. 1994). The annealing and compression induced starbursts may well be the relevant process in subcluster mergers.

Successive cycles of annealing may also be an effective means of building up bulges, and converting spirals to S0s on a timescale of several cluster crossing times. The power-law functions for mass and angular momentum loss in discrete events could be used in simulations of cluster evolution to investigate these possibilities.

## 7 ACKNOWLEDGEMENTS

We are grateful to D. C. Smith for helpful conversations and computer tips, to J. van Gorkom and J. Kenney for advice and results in advance of publication, and to an anonymous referee for many helpful suggestions. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This research has also made use of NASA’s Astrophysics Data System Bibliographic Services.