Ram Pressure Stripping of the Multiphase ISM: A Detailed View from TIGRESS Simulations

We use the TIGRESS framework developed by Kim & Ostriker (2017) to evolve the multiphase, turbulent, magnetized ISM with which the ICM interacts. We refer the reader to Kim & Ostriker (2017) for full details of the methods and tests. TIGRESS solves the ideal MHD equations in a local shearing box representing a ∼kiloparsec patch of differentially rotating galactic disks using the Athena code (Stone et al. 2008; Stone & Gardiner 2009). Local Cartesian coordinates x and y correspond to the local radial and azimuthal directions of global galactocentric coordinates such that (x, y) = (R − R₀, R₀[ϕ − Ω₀ t]), while z is the vertical coordinate. The simulation domain is corotating with galactic rotation speed at the center of the simulation domain, Ω₀ ≡ Ω(R₀), leading to inertial forces, including the Coriolis force and the tidal potential in the momentum equation. The flat rotation curve is assumed, $d\mathrm{ln}{\rm{\Omega }}/d\mathrm{ln}R=-1$. We adopt shearing-periodic boundary conditions in the horizontal directions (Stone & Gardiner 2010) and outflow boundary conditions in the vertical directions. The bottom vertical boundary conditions are modified for the ICM inflows (see Section 2.3).

We solve Poisson's equation to obtain gravitational potential from gas and newly formed young stars using the fast Fourier transform method with horizontally shearing-periodic and vertically open boundary conditions (Gammie 2001; Koyama & Ostriker 2009). The gravitational potential of old stellar disks and dark matter halos is held fixed and only exerts vertical gravity. We introduce a sink particle when a gas cell is experiencing unresolved self-gravitating collapse as indicated by the Larson–Penston density threshold at ${\rho }_{\mathrm{LP}}\equiv (8.86/\pi ){c}_{s}^{2}/G{\rm{\Delta }}{x}^{2}$, where c_s ≡ (P/ρ)^1/2 is the local sound speed and Δx is the side length of a cubic grid cell used in the simulation (Gong & Ostriker 2013). We adopt additional criteria for the sink particle creation, including a converging flow check (in all three directions) and a local potential minimum check. Typically, ρ_LP ∼ 100 cm⁻³ for 8 pc resolution and ∼ 300 cm⁻³ for 4 pc resolution. Note that the typical mass of sink particles is in a range between a few × 10³ M_⊙ and 10⁵ M_⊙, representing star clusters rather than individual stars. We treat each particle as a population of stars with the fully sampled initial mass function (IMF) of Kroupa (2001).

We use the STARBURST99 stellar population synthesis model to obtain the SN rate and FUV luminosity for each star cluster (Leitherer et al. 1999). In addition to clustered SNe occurring at the position of the sink particle, we produce a massless particle with a probability of 50% for each SN event to model a runaway star ejected from a binary OB star (Eldridge et al. 2011). The total SN rate still matches that of the original STARBURST99 model.

For each SN event, we first identify the cells with distances from the explosion center smaller than R_SNR = 3Δx and calculate the total mass M_SNR and volume V_SNR of the feedback region (or the SN remnant (SNR)). If M_SNR/M_sf < 1, where ${M}_{\mathrm{sf}}=1540\,{M}_{\odot }{({n}_{\mathrm{amb}}/{\mathrm{cm}}^{-3})}^{-0.33}$ is the shell formation mass at a given ambient medium density n_amb = M_SNR/V_SNR (Kim & Ostriker 2015), we inject 10⁵¹ erg, dividing into thermal and kinetic energy with the Sedov stage energy ratio of 0.72:0.28. Otherwise, we inject the terminal momentum of the SNR ${p}_{\mathrm{SNR}}=2.8\times {10}^{5}\,{M}_{\odot }\,\mathrm{km}\,{{\rm{s}}}^{-1}{({n}_{\mathrm{amb}}/{\mathrm{cm}}^{-3})}^{-0.17}$ as calibrated in Kim & Ostriker (2015). The total and metal masses of SN ejecta, M_ej = 10 M_⊙ and ${Z}_{\mathrm{SN}}{M}_{\mathrm{ej}}=2\,{M}_{\odot }$ with ${Z}_{\mathrm{SN}}=10\,{Z}_{\odot }$, are traced using passive scalars. See Section 2.4 for details.

We use the total FUV luminosity from star clusters in the simulation domain to set the instantaneous photoelectric heating rate by interstellar radiation field (Bakes & Tielens 1994; Weingartner & Draine 2001). We apply the mean attenuation factor using the plane-parallel approximation as in Kim et al. (2020b). As a result, the heating rate varies in time self-consistently but is spatially constant.

Optically thin cooling is included in the energy equation using a tabulated cooling rate coefficient Λ(T) from Koyama & Inutsuka (2002) at T < 10^4.2 K and from Sutherland & Dopita (1993) at T > 10^4.2 K (collisional ionization equilibrium at solar metallicity is adopted) depending only on temperature. Although we follow the metallicity of gas in each cell (see Section 2.4), we note that we do not use the metallicity information to set the cooling rate.

More self-consistent treatment of radiation and chemistry and hence cooling and heating rates is being developed for the TIGRESS framework (J.-G. Kim et al. 2022, in preparation), which will enable further study of RPS with more realistic ISM. This extension is particularly important to pursue the extraplanar molecular gas in RPS galaxies.

In this work, we make use of the solar neighborhood model of the TIGRESS simulation suite, setting the parameters for the gravitational potential of stars and dark matter (see below; Equation (4)). We adopt the angular velocity of galactic rotation Ω₀ = 28 km s⁻¹ kpc⁻¹, giving rise to the orbit time t_orb = 2π/Ω₀ = 224 Myr. We use a vertically elongated rectangular box with the outer dimensions of (L_x , L_y , L_z ) = (1024, 1024, 7168) pc. A uniform, cubic grid cell is used with the side length of Δx = 8 pc at which we achieve convergence of overall properties of the ISM, SFR, and outflows (see Kim & Ostriker 2017, 2018; Kim et al. 2020b). This model is referred to as the noICM model throughout the paper (identical to the solar neighborhood model, R8, presented in other works). Additional details of the solar neighborhood model can be found in Kim & Ostriker (2017) for initial conditions, overall evolution, numerical convergence, and technical details; in Kim & Ostriker (2018) and Vijayan et al. (2020) for galactic fountains and winds; and in Mao et al. (2020) for the properties of gravitationally bound clouds and their connection with SFRs.

The simulation starts from an idealized initial condition with horizontally uniform, vertically stratified gas profiles with the initial gas surface density of Σ_gas = 13 M_⊙ pc⁻². We introduce the initial velocity perturbation and set thermal pressure to ensure that the disk is in rough hydrostatic equilibrium. Soon after the simulation begins, the initially imposed velocity perturbation dissipates, and the gas cools. The overall vertical contraction occurs owing to the reduction of turbulent and thermal pressure, leading to a burst of star formation. SNe and FUV heating from newly formed massive stars respectively offset turbulence dissipation and gas cooling, recovering vertical support against gravity. The disk expands vertically, reducing SFRs and hence feedback. The reduction of feedback causes another disk contraction, and the cycle repeats. Each cycle has a period similar to the vertical oscillation timescales of 40–50 Myr (see Kim et al. 2020b). Although the first burst is a consequence of the idealized initial setups, our simulations soon enter a self-consistently regulated state after a few star formation–feedback cycles (t > 100 Myr in this model).

To overview the evolution in a quasi-steady state far from the initial burst, Figure 1 shows the horizontally averaged physical quantities in the space (vertical coordinate z) and time (t) plane as defined by

$$\begin{eqnarray}&&\left\langle q(z;t)\right\rangle \equiv \displaystyle \frac{\int q(x,y,z;t){dxdy}}{{L}_{x}{L}_{y}}\end{eqnarray} \tag{ 1 }$$

for a physical quantity q of interest. We only show a self-regulated state over t ∼ 250–500 Myr. From top to bottom, we show (a) hydrogen number density n_H, (b) outgoing mass flux ρ v_out, (c) turbulent pressure (or, equivalently, vertical momentum flux) $\rho {v}_{z}^{2}$, (d) thermal pressure P, and (e) magnetic pressure P_B ≡ B²/(8π). Here the outward vertical velocity is defined by ${v}_{\mathrm{out}}\equiv {v}_{z}\mathrm{sgn}(z)$ such that v_out is positive (red) for outflow and negative (blue) for inflow about the midplane. Note that the midplane (z = 0) in our simulation defines the symmetric plane of the fixed gravitational potential of stars and dark matter. Even without ICM inflows, the gas distribution can be largely asymmetric, as stochastic SN explosion cannot be perfectly symmetric. As a result, the total gravitational potential, including the gravitational potential of gas and young star clusters, can be asymmetric.

Zoom In Zoom Out Reset image size

Within the simulation duration shown in Figure 1, we can visually identify four strong outflow launching epochs (t ∼ 250, 320, 380, and 420 Myr; see panel (b)). These epochs are associated with strong star formation events. The outflows in our simulations show clear multiphase nature, consisting of fast, hot winds that escape the simulation domain and slow, warm fountains that fall back to the midplane (Kim & Ostriker 2018). The hot winds (T ≳ 10⁵−10⁶ K) can be easily identified by the high thermal and turbulent pressure gas in the extraplanar region z > 1 kpc with a steeper slope in the z-t plane (panels (c) and (d)). The outgoing mass flux in panel (b) associated with the hot winds is always red (only outflows). However, the warm fountains (T ≲ 10⁴ K) are evident with alternating colors in panel (b), implying that the warm outflows are always followed by inflows (see also panel (a)). Only a small fraction of the warm outflow has reached high velocity to escape the simulation domain (see Kim & Ostriker 2018; Vijayan et al. 2020). The magnetic pressure in panel (e) is overall subdominant, especially in low-density gas far from the midplane. The magnetic field strength grows over time via galactic dynamo and ranges from a few to 10 μG in the warm and cold medium with comparable turbulent and mean field strengths (see Kim et al. 2019). This is consistent with the observed magnetic field strength of neutral hydrogen in the solar vicinity (Heiles & Troland 2005). The full complexity of star formation/feedback and multiphase outflow/inflow cycles in the TIGRESS simulation suite is extensively discussed in Kim et al. (2020b).

We take the first snapshot shown in Figure 1 (t ∼ 245 Myr) as initial conditions and restart the simulation with an ICM inflow. Here and hereafter, we exclusively use the term ISM in simulations to denote the gas that was in the simulation domain before injecting the ICM. At this time, gas surface density in the noICM model is reduced to Σ_gas = 9.5 M_⊙ pc⁻² as gas turns into stars and leaves the simulation domain as outflows through the vertical boundaries. We model the ICM inflow as a constant, unmagnetized, vertical inflow through the bottom boundary (i.e., face-on interaction). We set the ICM metallicity to Z_ICM = 0.1 Z_⊙ (e.g., Urban et al. 2011), which serves as a tracer of the gas origin together with other passive scalars (see Section 2.4).

The ICM inflows are characterized by two parameters: hydrogen number density of the ICM n_ICM = ρ_ICM/(1.4271m_H) and inflow velocity v_ICM. We adopt the ICM sound speed c_s,ICM = 300 km s⁻¹. The total pressure of the ICM at injection is

$$\begin{eqnarray}\begin{array}{rcl}{P}_{\mathrm{ICM}} & \equiv & {\rho }_{\mathrm{ICM}}({v}_{\mathrm{ICM}}^{2}+{c}_{s,\mathrm{ICM}}^{2})={\rho }_{\mathrm{ICM}}{v}_{\mathrm{ICM}}^{2}(1+{{ \mathcal M }}_{\mathrm{ICM}}^{-2})\\ & \approx & 1.73\times {10}^{4}(1+{{ \mathcal M }}_{\mathrm{ICM}}^{-2}){k}_{{\rm{B}}}\,{\mathrm{cm}}^{-3}\,{\rm{K}}\\ & & \times \left(\frac{{n}_{\mathrm{ICM}}}{{10}^{-4}\,{\mathrm{cm}}^{-3}}\right){\left(\frac{{v}_{\mathrm{ICM}}}{{10}^{3}\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right)}^{2},\end{array}\end{eqnarray} \tag{ 2 }$$

which is dominated by the ram pressure for our chosen v_ICM ≥ 10³ km s⁻¹ (or Mach number of the ICM ${{ \mathcal M }}_{\mathrm{ICM}}\,\gt 3.3$). ³ In the simulations, however, as soon as the ICM sweeps up the ISM, a reverse shock thermalizes the inflowing ICM, and it is the hot ICM with the total pressure P_ICM dominated by the thermal term that interacts with the ISM.

While the ISM is pushed away from the galactic disk owing to the interaction with the ICM, the stellar and dark matter components are not immediately disturbed in RPS galaxies. This is particularly true for our simulations because we use a fixed analytic potential for stellar and dark matter gravity. The gas weight under the external gravity ⁴ is

$$\begin{eqnarray}&&{{ \mathcal W }}_{\mathrm{ext}}\equiv {\int }_{0}^{\infty }\rho \left|\displaystyle \frac{d{{\rm{\Phi }}}_{\mathrm{ext}}}{{dz}}\right|\,{dz},\end{eqnarray} \tag{ 3 }$$

where the functional form of the external gravitational potential is

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Phi }}}_{\mathrm{ext}}(z)\equiv 2\pi G{{\rm{\Sigma }}}_{* }{z}_{* }\left[{\left(1+\displaystyle \frac{{z}^{2}}{{z}_{* }^{2}}\right)}^{1/2}-1\right]\\ \,+\,2\pi G{\rho }_{\mathrm{dm}}{R}_{0}^{2}\mathrm{ln}\left(1+\displaystyle \frac{{z}^{2}}{{R}_{0}^{2}}\right).\end{array}\end{eqnarray} \tag{ 4 }$$

We adopt the parameters representing solar neighborhood conditions: galactocentric distance R₀ = 8 kpc, stellar surface density Σ_* = 42 M_⊙ pc⁻², stellar scale height z_* = 245 pc, and midplane dark matter density ρ_dm = 6.4 × 10⁻³ M_⊙ pc⁻³ (Zhang et al. 2013; McKee et al. 2015). When the gas is stripped far away from the stellar disk (i.e., the mean gas position is much larger than the stellar disk scale height, z ≫ z_*), the stellar gravity (the first term in Equation (4)) becomes nearly constant such that ∣dΦ_*/dz∣ = 2π GΣ_*. The ISM weight can then be well approximated by ${{ \mathcal W }}_{\mathrm{ext}}\approx {{ \mathcal W }}_{\mathrm{GG}}$, where

$$\begin{eqnarray}\begin{array}{c}\begin{array}{cc}{{ \mathcal W }}_{{\rm{G}}{\rm{G}}} & \,\equiv \,2{\pi }G{{\rm{\Sigma }}}_{{\rm{g}}{\rm{a}}{\rm{s}}}{{\rm{\Sigma }}}_{\ast }\\ & =\,5.27\times {10}^{4}{k}_{{\rm{B}}}\,{{\rm{c}}{\rm{m}}}^{-3}\,{\rm{K}}\,\left(\frac{{{\rm{\Sigma }}}_{{\rm{g}}{\rm{a}}{\rm{s}}}}{9.5\,{M}_{\odot }\,{{\rm{p}}{\rm{c}}}^{-2}}\right)\left(\frac{{{\rm{\Sigma }}}_{\ast }}{42\,{M}_{\odot }\,{{\rm{p}}{\rm{c}}}^{-2}}\right).\end{array}\end{array}\end{eqnarray} \tag{ 5 }$$

This "restoring" force per area (often called the "anchoring" pressure), originally presented in Gunn & Gott (1972), has been conveniently compared with the ICM ram pressure to determine the stripping condition (e.g., Kenney et al. 2004; Vollmer et al. 2006; Chung et al. 2007; Jaffé et al. 2018; Köppen et al. 2018). Note that, for our adopted gravitational potential, the dark matter contribution in the vertical gravity keeps increasing with z and becomes comparable to that of stars at the vertical boundaries z ∼ 3.5 kpc. Therefore, ${{ \mathcal W }}_{\mathrm{GG}}$ slightly underestimates the maximum ${{ \mathcal W }}_{\mathrm{ext}}$ in our simulations by 25%.

Table 1 lists the ICM models. Column (1) is the model name; we adopt a nomenclature including the strength of the ICM ram pressure presented in Column (4). The higher-resolution models (Δx = 4 pc) have their name ending with "h" (see Column (6)). For the higher-resolution models, we refine the original data cube from the noICM model using a zero gradient prolongation (i.e., volume- and area-averaged quantities in finer grid cells are the same as their parent cell values). Therefore, the initial conditions are identical across all models. Columns (2) and (3) are the number density and inflow velocity of the ICM, respectively, which set the ICM pressure (Column (4); Equation (2)). Column (5) shows the ratio of the ICM pressure and the maximum ISM weight under the stellar gravity (Equation (5)), which is a rough estimate for the relative strength of the ICM pressure to the maximum ISM weight. Finally, we list the spatial resolution in Column (6).

Table 1. ICM Model Parameters

Model	nICM	vICM	PICM/kB	${P}_{\mathrm{ICM}}/{{ \mathcal W }}_{\mathrm{GG}}$	Δx
	(10−4 cm−3)	(103 km s−1)	(104 cm−3 K)		(pc)
(1)	(2)	(3)	(4)	(5)	(6)
ICM-P1	0.5	1	0.94	0.18	8
ICM-P3(h)	1	1.4	3.6	0.69	8(4)
ICM-P7(h)	2	1.4	7.2	1.4	8(4)
ICM-P14	2	2	14	2.7	8

Note. Column (1): model name. Column (2): hydrogen number density of the ICM. Column (3): relative velocity of the ICM and the ISM disk. Column (4): ICM pressure. Column (5): ratio of the ICM pressure to the ISM weight. ${{ \mathcal W }}_{\mathrm{GG}}=5.27\times {10}^{4}{k}_{{\rm{B}}}\ \,{\mathrm{cm}}^{-3}\,{\rm{K}}$ is an approximate ISM weight estimated by Equation (5). Column (6): spatial resolution.

Download table as:  ASCIITypeset image

We consider four different ICM conditions (with two additional higher-resolution runs), covering P_ICM/k_B ∼(1 − 14) × 10⁴ cm⁻³ K. Since the ISM condition is fixed, the relative strength of the ICM−ISM interaction simply increases as P_ICM increases, with ICM-P1 and ICM-P3(h) having ${P}_{\mathrm{ICM}}/{{ \mathcal W }}_{\mathrm{GG}}\lt 1$ and ICM-P7(h) and ICM-P14 having ${P}_{\mathrm{ICM}}/{{ \mathcal W }}_{\mathrm{GG}}\gt 1$. Throughout the paper, the former two models are referred to as the weak ICM models, and the latter two models are referred to as the strong ICM models. Our parameter choice brackets the relative strength of the ISM−ICM interaction seen in NGC 4522, a prototypical galaxy undergoing RPS in the Virgo Cluster (Kenney et al. 2004; Vollmer et al. 2006). In addition, the anchoring pressure of our simulation (Equation (5)) is comparable to that near the truncation radius of NGC 4522 (Kenney et al. 2004; Chung et al. 2009; Lee et al. 2017; Lee & Chung 2018). Thus, our weak/strong ICM models can represent the evolution of the inner/outer part of the truncation radius of NGC 4522.

In the TIGRESS framework, the gas is divided into five thermal phases based on its temperature, corresponding to typical discriminators of the three-phase ISM (but including thermally unstable phases; McKee & Ostriker 1977). Each cell is exclusively assigned a cold, unstable, warm, intermediate (warm-hot ionized medium), or hot phase following the temperature criteria in Table 2. We often combine the cold, unstable, and warm phases and call them the cool phase.

The hot gas in the noICM model is created by SN shocks, while the warm and cold phases are maintained via radiative heating due to FUV radiation. With ICM inflows, a significant fraction of the hot ICM is directly added. The hot gas can accelerate the cooler gas directly through its pressure gradient (both ram and thermal pressure), but another significant (likely dominant) acceleration mechanism, as we shall show, is by the momentum transfer through the mixing of the hot gas into the cooler gas (Section 4; see also Fielding & Bryan 2022). It is therefore critical to separate the origin of gas and trace the fraction of different origin gas in different thermal phases. We utilize passive scalars to track the mass fractions of the initial ISM, SN ejecta, and ICM in each cell. Here we use the term "passive scalar" to denote the multiplication of gas density and tracer field (or specific scalar).

Practically, we follow the total metallicity Z, SN ejecta mass fraction ${f}_{\mathrm{SN}}^{}$, and ICM mass fraction ${f}_{\mathrm{ICM}}^{}$. The metallicity tracer field Z is initialized with Z_⊙ = 0.02 at the beginning of the noICM simulation, while the other two tracer fields are initialized to zero everywhere. In the code, additional continuity equations for ρ Z, $\rho {f}_{\mathrm{SN}}^{}$, and $\rho {f}_{\mathrm{ICM}}^{}$ are solved with the velocity field of gas. For each SN event, we add the total and metal densities of SN ejecta, ρ_ej ≡ M_ej/V_SNR and ${Z}_{\mathrm{SN}}{\rho }_{\mathrm{ej}}$, respectively, to passive scalars in the feedback region (of course, the SN ejecta density is added to the gas density as well). As the noICM model has evolved for ∼250 Myr before the restart with the ICM, the ISM disk's metallicity has been enriched by SN ejecta. When we restart simulations with ICM inflows, we adopt metallicity and SN ejecta fraction inherited from the noICM model. The ICM inflow with the ICM tracer field ${f}_{\mathrm{ICM}}^{}=1$ is then added and followed by another passive scalar. In addition, the ICM metallicity is set to Z_ICM = 0.1 Z_⊙ and mixed into the metal passive scalar ρ Z as the ICM interacts with the existing gas. When sink particles are formed and accrete gas, passive scalars are also locked into particles, which represent the metallicity of star-forming gas.

As mentioned above, we adopt three distinct metallicities for different origin gases:

• 1.  
  Genuine ISM—initial gas from the beginning of the noICM model (Z_ISM,0 = Z_⊙);
• 2.  
  SN ejecta—gas added in the feedback region by SNe (${Z}_{\mathrm{SN}}=10\,{Z}_{\odot }$); and
• 3.  
  ICM—gas added from the bottom boundaries as the ICM inflow (Z_ICM = 0.1 Z_⊙).

The metallicity is a good proxy for the composition of the gas in simulations, providing potential observational imprints. In each cell, the metallicity is connected to the SN ejecta and ICM mass fractions as

$$\begin{eqnarray}&&Z={Z}_{\mathrm{ISM},0}(1-{f}_{\mathrm{SN}}^{}-{f}_{\mathrm{ICM}}^{})+{Z}_{\mathrm{SN}}{f}_{\mathrm{SN}}^{}+{Z}_{\mathrm{ICM}}{f}_{\mathrm{ICM}}^{}.\end{eqnarray} \tag{ 6 }$$

As presented in Schneider et al. (2020), the dominance of the mixing-driven momentum transfer from hot to cool phase can be simply evidenced by the linear relation between the source tracer field and velocity of the cool phase. The predicted outflow velocity of the cool phase in the case of mixing-driven momentum transfer is

$$\begin{eqnarray}&&{v}_{{\rm{z}}}^{\mathrm{cool}}={v}_{\mathrm{ICM}}{f}_{\mathrm{ICM}}^{\mathrm{cool}}\end{eqnarray} \tag{ 7 }$$

for ICM-accelerated outflows. In our simulations, both SNe and ICM create the hot gas, so that both SN ejecta and ICM tracer fields can leave imprints on the accelerated cool outflows. However, as the relation only holds for the fresh tracer field (the tracer field that is first transferred from hot to cool phase), the ICM mass fraction provides an ideal tracer for this purpose, especially for the first acceleration of the ISM. In contrast, as we restarted the simulations from the noICM model after many feedback−star formation cycles, a lot of SN ejecta is already mixed into the cool phase that has been accelerated, fallen back, and reaccelerated many times. The total SN ejecta mass fraction in the cool phase is no longer representative of the amount of the hot gas that is currently mixed. Still, we can see signs of SN-accelerated gas from the relatively metal-enriched gas that is moving faster (Section 4.2).

To summarize the general response of the ISM as a whole to the ICM inflows, Figure 2 plots the vertical profiles of the ICM mass fraction in the hot phase (left) and vertical momentum density (=upward mass flux; right) as a function of time for all ICM models. The former is defined by the mass-weighted horizontal average of ${f}_{\mathrm{ICM}}^{}$ of the hot gas, $\overline{{f}_{\mathrm{ICM}}^{\mathrm{hot}}}\equiv {\left\langle \rho {f}_{\mathrm{ICM}}^{}\right\rangle }^{\mathrm{hot}}/{\left\langle \rho \right\rangle }^{\mathrm{hot}}$. From top to bottom, we show all models in ascending order of the ICM ram pressure, including two high-resolution models shown in the third and fifth rows. We plot the reference line of $\overline{{f}_{\mathrm{ICM}}^{\mathrm{hot}}}=0.5$ that defines the mean boundary of the ICM–ISM. Owing to the multiphase structure of the ISM, actual boundaries between the ICM and ISM are much more complex (see Section 3.3). The positive (upward) mass flux below the interface simply represents the mass flux of the ICM inflows. As soon as the simulations restarted with the ICM inflows, the ISM in the bottom half is quickly pushed up, and the ICM−ISM interface approaches the midplane in 10−50 Myr depending on the ICM inflow strength. Then, the interface either remains near the midplane with clear separation in the ICM fraction (weak ICM models) or continuously marches upward (strong ICM models). This dichotomy is in excellent agreement with the expectation based on the simple stripping condition listed in Column (5) of Table 1.

Zoom In Zoom Out Reset image size

Using the position of the ICM−ISM interface, we divide overall evolution into three stages: compression stage, early stripping stage, and active stripping stage. The earliest time at which the ICM−ISM interface reaches z = − 500 and 500 pc, respectively, defines the beginning of the early and active stripping stages (marked by left- and right-pointing triangles in Figure 2).

Because the multiphase ISM is porous and has low-density channels through which the ICM can penetrate, the ICM gradually pollutes the ISM in the upper disk even when the interface stays near the disk midplane. The only exception is the ICM-P1 model, where the penetration of the ICM is not effective, and the ISM in z > 0 remains unpolluted over the entire simulation duration (${f}_{\mathrm{ICM}}^{}\lt 1 \%$). As a result, the mass flux evolution in the upper disk of the ICM-P1 model is qualitatively similar to that in the noICM model (Figure 1(b)), indicating that the outflows are still driven by SN feedback. In the ICM-P3 model, the ICM mass fraction in the upper half increases quickly and becomes larger than 10%. The mass flux in this model is overall enhanced, while the fountain component (alternating positive and negative signs) still exists, implying insufficient acceleration of the cool phase with the marginally weak ICM. The high-resolution model (ICM-P3h) behaves essentially the same, while the late-time evolution shows a strong outflowing epoch. Given the inherently stochastic nature of the evolution, this difference should not be interpreted as systematic resolution dependence. Rather, the qualitative similarity of the early evolution (t < 350 Myr) means the convergence of the overall evolution.

The ICM penetration in both ICM-P7 and ICM-P14 models is highly efficient and leads to the immediate enhancement of the ICM mass fraction in the upper disk. The ICM−ISM interface continuously moves upward in the ICM-P14 model, but the ICM-P7 model spends quite a long time (∼70 Myr) in the early stripping stage with a significantly larger ICM mass fraction (>10%) in the upper disk than that in the ICM-P3 model. The net mass flux in the upper disk is always positive in the strong ICM models, demonstrating the dominant role of the ICM in driving outflows and implying RPS in action. Again, the high-resolution model (ICM-P7h) shows a very similar evolution to its low-resolution counterpart.

We emphasize that the multiphase RPS occurs continuously in both early and active stripping stages for the strong ICM models. In the early stripping stage, the ICM finds low-density channels in the porous ISM to penetrate. In doing so, the ICM begins to shred the ISM and transfer mass, momentum, and energy while mixing occurs. In the active stripping stage, which only exists in the strong ICM models, the ICM fills the volume in a wide range of the disk. The ICM mixing and momentum transfer occur almost all around the simulation volume, and the ISM is effectively accelerated and removed from the simulation domain. We will delineate the mixing-driven stripping in Section 4.

Figure 3 shows the time evolution of (a) total (ICM+ISM) gas surface density, (b) surface density of new stars Σ_new−star, and (c) surface density of gas passed through the vertical boundaries Σ_out. Σ_new−star is defined by summing up the total mass of stars formed since we restart the simulations, and Σ_out is calculated by integrating the net mass flux at the vertical boundaries (both escaped to the top and injected from the bottom) over time. A clear dichotomy between the weak and strong ICM models is visible. The strong ICM models lose gas and stop forming stars after ∼100 Myr. The weak ICM models retain (or even gain) gas and form more stars than the noICM model. Despite the highly complex interaction between the multiphase ISM and ICM revealed in our simulations (as discussed in Section 3.3), the simple stripping condition estimated by Equation (5) provides a reliable prediction for the fate of the gas disk. This is in part because we only model the face-on, plane-parallel interaction, ideal for the simple criteria to work best.

Zoom In Zoom Out Reset image size

Even without the ICM, the noICM model also loses its gas through star formation and outflows powered by SN feedback (Kim & Ostriker 2018; Kim et al. 2020b). The mean SFR and mass outflow rate over the simulation duration of 250–500 Myr are Σ_SFR = 3.1 × 10⁻³ M_⊙ kpc⁻² yr⁻¹ and ${\dot{{\rm{\Sigma }}}}_{\mathrm{gas},\mathrm{out}}=7.7\,\times {10}^{-4}\,{M}_{\odot }\,{\mathrm{kpc}}^{-2}\,{\mathrm{yr}}^{-1}$, respectively. With the ICM inflows, the total gas mass within the domain can increase as the ICM is added to the system unless the outflow and SFRs are greatly increased. The ICM-P1 model closely follows the evolution curve of the noICM model as the enhancement in SFRs is compensated by the decrease in outflow rates (panel (a)). In the ICM-P3 model, Σ_out becomes negative, implying that there is a net inflow through the boundaries (panel (c)). Overall gas mass still decreases when taking into account the loss due to star formation. In the ICM-P7 and ICM-P14 models, the gas mass decreases quickly, as the outward mass fluxes are greatly enhanced after the compression stage (see also Figure 2). The gas compression by the ICM causes enhancement of the early star formation in all ICM models (see Section 5 for in-depth analysis). The half-mass stripping time defined by the time interval between ${{\rm{\Sigma }}}_{\mathrm{gas}}(t)={{\rm{\Sigma }}}_{\mathrm{gas},\max }$ and ${{\rm{\Sigma }}}_{\mathrm{gas},\max }/2$ is ∼ 130 and 60 Myr for the ICM-P7 and ICM-P14 models, respectively. At a similar timescale, star formation in the ICM-P7 and ICM-P14 models is completely quenched, showing a flattening in Σ_new−star (panel (b)). Σ_out flattens later after complete stripping (the ICM flows freely; panel (c)).

Overall qualitative behaviors are converged with the resolution, but the later-time evolution shows differences mainly due to the stochasticity of the evolution.

To delineate the interaction between the ICM inflows and the multiphase ISM, Figure 4 shows snapshots at t = 275 Myr for two representative models, ICM-P3h (top) and ICM-P7h (bottom). At this epoch, the ISM in the bottom half has already pushed up the midplane in both models, so we only show the upper disk z > −0.5 kpc. Visual impressions between projected density and slices are substantially different. The projected density shows an overall density distribution with mild fluctuations, with similar sharp cutoffs at around z = 0 for both models. The immediate impression might be that there is a well-defined ICM−ISM interface near the midplane. However, the slices unveil a highly porous, multiphase structure with a large density and temperature contrast. In addition, significant penetration of the ICM (see panel (f)) depicts a significant difference between the two models. We emphasize that, compared to the projection maps, slices (or thin projections) of gas physical quantities are often more useful to deliver visual insights. This is particularly true when looking into the interaction between the different phases with large contrasts in physical properties.

The density (panel (b)) and temperature (panel (c)) slices of both models show a large, cool gas structure that begins to face the ICM near the midplane. This structure looks like a continuous, single structure in density and temperature, but v_z (panel (d)) and hence $\rho {v}_{z}^{2}$ (panel (g)) show a sharp change between the left and right sides. The enhanced ICM mass fraction ${f}_{\mathrm{ICM}}^{}$ (panel (f)) and the relatively large outward velocity (panel (d)) on the right side of this structure imply that this is the gas originally in the bottom half accelerated by the ICM. The accelerated cool gas from the lower disk remains intact and appears as a high ram pressure chunk at around z = 0.5 kpc in the ICM-P3h model, but it is already substantially shredded and fragmented in the upper disk in the ICM-P7h model. Another interesting cool gas structure that shows a difference between the two models is one located at z ∼ 1 kpc near the left edge of the slice. In the ICM-P3h model, this structure is still falling (panel (d)), while the same structure has already significantly shredded and accelerated by the interaction with the ICM in the ICM-P7h model. There are plenty of similar falling gases (fountain flows) in the ICM-P3h model across a wide range of z. In the ICM-P7h model, such infalling fountain flows are no longer prevalent, but they are generally more compressed and even outflowing.

Generally, the ICM in the ICM-P7h model manages to intrude almost the entire regions of the upper disk. It is visually evident from the enhanced ${f}_{\mathrm{ICM}}^{}$ seen in most regions, which correspond to the enhanced ram and thermal pressure and outward velocity. The metallicity (panel (e)) contains more complex information due to the additional contribution of the high-metallicity SN ejecta. As noted in Section 2.4, the ISM has been enriched in the noICM model over ∼250 Myr. The mean ISM metallicity at the beginning of the ICM models (shown in bright orange in panel (e)) is thus larger than the initial metallicity Z_⊙. Without the ICM, it is the high-metallicity gas injected by SNe that is filling low-density regions (it is still the case for z > 2 kpc in the ICM-P3h model). With the ICM, now the low-metallicity gas is filling the low-density regions (more evident in the ICM-P7h model), which also show high pressures (panels (g) and (h)). As the high-pressure ICM compresses the ISM, the magnetic pressure is enhanced in cooler, denser gas (panel (i)).

In Figure 5, we select snapshots at t = 340 Myr for the ICM-P3h (top) and ICM-P7h (bottom) models. The late-time evolution differs in the two standards more substantially and is even evident in the density projection (panel (a)). Since the ICM inflow alone cannot keep pushing the ISM away in the ICM-P3h model, the bulk ISM is falling back as star formation has suppressed. In contrast, the strong ICM inflow alone can continue to strip the ISM in the ICM-P7h model (see also Figure 2). While disturbed significantly, in the ICM-P3h model, the ICM fraction in the cool phase is still less than a few percent, and the overall visual impression is not very different from what is shown in Figure 4. The hot ICM keeps penetrating through the low-density channels, creating shearing interfaces between the cool and hot phases in which the majority of mixing occurs. The volume fractions of the hot and cool phases are comparable all over the upper disk. We note that when the hot gas is only created by SNe in the noICM model, the hot gas volume fraction is ∼20%–30% near the midplane and increases to ∼50% at z ∼ 1 kpc and to ∼100% at z > 2 kpc (Kado-Fong et al. 2020). Star formation continues in this model at slightly higher rates than the noICM model. Once stars form, they fall faster than gas, as stars do not feel the ICM pressure. This causes SN feedback in the ICM-dominated regions below the midplane, sometimes creating metal-enriched hot bubbles between the hot ICM and cool ISM (panels (e) and (f)).

Zoom In Zoom Out Reset image size

At this time, the ICM-P7h model already loses ∼30% of its total mass, and almost all ISM is pushed above z ∼ 0.5 kpc (panel (a)). The cool ISM is highly fragmented and confined in smaller volumes (panels (a) and (b)). The dense gas structure facing the ICM at z ∼ 0.5 kpc is the leftover from the first major stripping of the main ISM disk and falling back to the midplane. Stars were just born in this strongly compressed structure, which is the final major star formation event before complete quenching (there will be additional extraplanar star formation in the structure far above the midplane later). Star clusters formed at high z during the previous evolutionary stage have fallen below the ICM−ISM interface, some of which still host SNe and create metal-enriched bubbles at z ∼ 0 kpc. Even in the cool gas above z > 2 kpc, the ICM fraction is quite high, ∼10%, as the ICM is continuously mixed into transfer momentum and maintains momentum flux against the weight at that height (Section 4.2). The acceleration was not sufficient to blow away this structure though. Except for these distinctive large structures, smaller clouds have already been ablated as evidenced by many tadpole structures whose density and temperature are respectively higher and lower than those of the typical hot ICM. The intermediate-phase gas populated in the wakes of the front clouds in part escapes the domain and in part condenses back to the cool phase, especially when the wakes meet the large cool structure in the back. From panels (d) to (g), it is evident that such acceleration is most efficient in the envelope of the large cool phase structure, where both ICM fraction and vertical velocity are relatively high.

In this subsection, we define physical quantities averaged over a range of volume and time to understand how different gas phases exchange their mass, momentum, and energy at different regions and times. We begin by integrating the conserved form of the MHD equations over the entire horizontal area A = L_x L_y and chosen vertical range ($z\in ({z}_{\min },{z}_{\max })$) to obtain a set of conservation equations

$$\begin{eqnarray}&&\dot{q}+[{{ \mathcal F }}_{q}({z}_{\max })-{{ \mathcal F }}_{q}({z}_{\min })]A=+{\dot{q}}_{\mathrm{source}}-{\dot{q}}_{\mathrm{sink}},\end{eqnarray} \tag{ 8 }$$

where q = M, p, and E for mass, vertical momentum, and total energy, which are respectively defined as

$$\begin{eqnarray}&&M\equiv {\int }_{{z}_{\min }}^{{z}_{\max }}\rho {dV},{\rm{and}}\quad p\equiv {\int }_{{z}_{\min }}^{{z}_{\max }}\rho {v}_{z}{dV},\end{eqnarray} \tag{ 9 }$$

and

$$\begin{eqnarray}&&E\equiv {\int }_{{z}_{\min }}^{{z}_{\max }}\left(\displaystyle \frac{1}{2}\rho {v}^{2}+\displaystyle \frac{P}{\gamma -1}+{P}_{B}\right){dV}.\end{eqnarray} \tag{ 10 }$$

The horizontally averaged fluxes (the square bracket term on the left-hand side of Equation (8)) are respectively defined for mass, vertical momentum, and total energy as

$$\begin{eqnarray}&&{{ \mathcal F }}_{M}(z)\equiv \displaystyle \frac{1}{A}\int \rho {v}_{z}{dxdy},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{{ \mathcal F }}_{p}(z)\equiv \displaystyle \frac{1}{A}\int \left(\rho {v}_{z}^{2}+P+{P}_{B}-\displaystyle \frac{{B}_{z}^{2}}{4\pi }\right){dxdy},\end{eqnarray} \tag{ 12 }$$

and

$$\begin{eqnarray}&&{{ \mathcal F }}_{E}(z)\equiv \displaystyle \frac{1}{A}\int \left(\rho {v}_{z}\left(\displaystyle \frac{{v}^{2}}{2}+\displaystyle \frac{\gamma }{\gamma -1}\displaystyle \frac{P}{\rho }\right)+{{ \mathcal S }}_{z}\right){dxdy},\end{eqnarray} \tag{ 13 }$$

where the adiabatic index γ = 5/3 and the vertical component of the Poynting flux is ${{ \mathcal S }}_{z}\equiv ({v}_{z}{B}^{2}-{B}_{z}v\cdot {\boldsymbol{B}})/(4\pi )$. Note that we do not include the gravitational potential term in the energy flux such that the work done by gravity appears as a sink term on the right-hand side of Equation (8), similarly to the momentum sink term by weight.

In the mass conservation equation, we have mass sink due to star formation (${\dot{M}}_{* }$) and mass source due to SN ejecta (${\dot{M}}_{\mathrm{SN}}={\dot{N}}_{\mathrm{SN}}{M}_{\mathrm{ej}}$) if stars are born and SNe are exploded in the particular volume of interest. For our chosen stellar population synthesis model (STARBURST99 with the Kroupa IMF; Leitherer et al. 1999), 1 SN ejects 10 M_⊙ for ∼100 M_⊙ of new star formed, implying ${\dot{M}}_{\mathrm{SN}}\sim 0.1{\dot{M}}_{* }$ on average. We note that we did not subtract the mass of stars exploding as SNe in the simulations. The gravitational weight acts as a sink of vertical momentum, ${\dot{p}}_{\mathrm{sink}}={ \mathcal W }A$ and

$$\begin{eqnarray}&&{ \mathcal W }\equiv \displaystyle \frac{1}{A}{\int }_{{z}_{\min }}^{{z}_{\max }}\rho \displaystyle \frac{d{\rm{\Phi }}}{{dz}}{dV},\end{eqnarray} \tag{ 14 }$$

where Φ = Φ_ext + Φ_sg + Φ_tidal includes the external potential (Equation (4)), the self-gravity potential returned from the solution of Poisson's equation for both gas and star cluster particles, and the tidal potential arising from the local rotating frame Φ_tidal = − qΩ² x². Finally, there are energy sources from SN energy injection and shearing-box stress ${\dot{E}}_{\mathrm{source}}={\dot{N}}_{\mathrm{SN}}{E}_{\mathrm{SN}}\,+{w}_{\mathrm{xy}}$, where w_xy = (qΩL_x )∫[ρ v_x δ v_y − B_x B_y /(4π)]dydz is integrated over either of the x boundaries (Hawley et al. 1995) and sinks from net radiative cooling ${\dot{E}}_{\mathrm{cool}}=\int ({n}_{{\rm{H}}}^{2}{\rm{\Lambda }}(T)-{n}_{{\rm{H}}}{\rm{\Gamma }}){dV}$ and the work done by gravity

$$\begin{eqnarray}&&{\dot{E}}_{{\rm{\Phi }}}\equiv {\int }_{{z}_{\min }}^{{z}_{\max }}\rho v\cdot {\rm{\nabla }}{\rm{\Phi }}{dV}.\end{eqnarray} \tag{ 15 }$$

As we exclusively assign the gas in different thermal phases, everything is separable by thermal phases. Since stars are formed from the cool (more precisely, cold) gas, mass sink by star formation can be accounted for in the cool phase. We include SN mass and energy source terms in the hot phase net mass and energy flux calculations, respectively; only a small fraction (≲10%) of unresolved SNe inject mass and energy into the form of the cooler phases. By integrating Equation (8) over a time interval Δt, we obtain

$$\begin{eqnarray}&&\sum _{\mathrm{ph}}{\dot{q}}_{\mathrm{net}}^{\mathrm{ph}}=0,\end{eqnarray} \tag{ 16 }$$

where

$$\begin{eqnarray}&&{\dot{q}}_{\mathrm{net}}^{\mathrm{ph}}\equiv \displaystyle \frac{{\rm{\Delta }}{q}^{\mathrm{ph}}-{\rm{\Delta }}{q}_{\mathrm{source}}^{\mathrm{ph}}+{\rm{\Delta }}{q}_{\mathrm{sink}}^{\mathrm{ph}}}{{\rm{\Delta }}t}+\left[{F}_{q,u}^{\mathrm{ph}}-{F}_{q,l}^{\mathrm{ph}}\right]A,\end{eqnarray} \tag{ 17 }$$

where ph = cool, int, and hot (see Table 2). Here Δq means the temporal difference of mass, vertical momentum, and total energy. Similarly, the temporal difference of cumulative mass and energy injection by SNe and stellar mass formed defines SN source and star formation sink. The weight and gravity work terms are defined by time integration. The net cooling term is calculated by time integration of the instantaneous net cooling rate from simulation outputs. Similarly, the time-averaged fluxes at upper and lower faces are defined by

$$\begin{eqnarray}&&{F}_{q,u/l}\equiv \displaystyle \frac{1}{{\rm{\Delta }}t}{\int }_{t}^{t+{\rm{\Delta }}t}{{ \mathcal F }}_{q}({z}_{\max /\min }){dt}.\end{eqnarray} \tag{ 18 }$$

By fully accounting for the temporal changes of each quantity, fluxes through vertical faces, and source/sink in each phase, ${\dot{q}}_{\mathrm{net}}^{\mathrm{ph}}$ represents any loss/gain of mass, momentum, and energy through phase transitions within spacetime bins of interest.

In practice, we measure each term using output snapshots dumped every ∼1 Myr. The cadence of snapshots was not fine enough to satisfy the conservation perfectly (Equation (16)), as this post-processing assumes that each variable involved in the time integration (e.g., cooling and flux terms) is constant over the snapshot interval of 1 Myr. ⁵ Bearing this caveat in mind, we analyze mass, momentum, and energy transfers between phases in space and time bins and only consider them to be reliable when the net changes are distinctive over the level of nonconservation errors.

4.1.1. Mass Transfer

We first consider mass conservation to understand phase transition between thermal phases. Figure 6 plots the net mass gain (red) and loss (blue) rates of cool, intermediate, and hot phases from left to right within the spacetime bins. Here we consider Δz = 400 pc thickness slabs centered at z = −3.6, −3.2, ⋯ , −0.4, 0, 0.4 ⋯ , 3.2, 3.6 kpc and a time interval of Δt = 10 Myr. The mass sink by star formation is added in the cool phase (left column) using new star particles formed in given spacetime bins. Similarly, the mass source by SN ejecta is subtracted in the hot phase (right column). The positive (red) and negative (blue) values in Figure 6 are solely due to phase transition within the spacetime bins.

Zoom In Zoom Out Reset image size

Without the ICM inflows (noICM; row (a)), the net gain of the hot and intermediate phases at the midplane stands out as thin red strips. In this case, it is clustered SNe that create the hot phase via shock heating of the ambient medium, showing the net loss (blue strip) in the cool phase. In our simulations, superbubbles expand into an inhomogeneous ambient medium, also creating a lot of the intermediate phase through the ablation of the cool phase shown as the net gain in the intermediate phase.

Both hot and intermediate phases subsequently cool above and below the midplane slab. Sometimes, this "cooling" region is further extended to z ∼ 1 kpc. It can be either direct cooling of the hot shocked gas (shell formation; e.g., Kim & Ostriker 2015) or mixing of the hot and cool gas (interface mixing; e.g., Kim et al. 2017; El-Badry et al. 2019). If the hot bubbles have completely cooled within the thickness of the midplane slab we consider (z = ±200 pc), the net gain/loss in different phases will not be visible. In other words, most superbubbles in our simulations can expand with their radii larger than 200 pc over 10 Myr before they cool.

The upper half of the ICM-P1 model shows overall similar results to the noICM model. As the ICM ram pressure gets stronger, noticeable differences begin to appear. In the ICM-P3h model, the layer in which the hot and intermediate phases gain the mass remains thin, while the cooling region is more extended toward higher z and becomes more prominent. The ICM−ISM interface stays near the midplane, and the effect of the ICM appears as the enhanced gain of the hot (and intermediate) phase in the midplane slab as the hot ICM shocks the cool ISM. At t ∼ 400 Myr, the net gain in the hot phase is visible beyond the midplane, which is due to the successful penetration of the hot ICM through the entire upper disk (see Figure 2). This breakout causes cool-to-hot phase transition followed by hot-to-cool phase transition at t ∼ 420–450 Myr.

With even stronger ICM inflows, the mass-gaining layer for the hot and intermediate phases gets thicker (slightly thicker for the intermediate phase). Now, the cool-to-hot phase transition is dominated by the interaction between the hot ICM and the cool ISM rather than SN feedback, especially at late times when star formation is nearly quenched. The large energy flux carried by the hot ICM can heat and ablate a significant amount of the cool phase, converting the cool ISM into the hot phase while populating the intermediate phase in mixing layers. Above the cool-to-hot phase transition layer, there is an extended region where the hot-to-cool phase transition occurs. In this region, best viewed in ICM-P7h, the hot and intermediate phases cool back to the cool phase. What is happening here is more like precipitation of the hot (with intermediate) phase into the volume filling cool phase that has been pushed by earlier interaction. This is somewhat different from the cooling of the mixed gas itself that drives phase transition from hot to cool in the cloud wakes as seen in radiative cloud-crushing simulations with large sizes (e.g., Armillotta et al. 2016; Gronke & Oh 2018, 2020; Sparre et al. 2020; Li et al. 2020; Kanjilal et al. 2021; Abruzzo et al. 2022).

As RPS is much more efficient in ICM-P14, the hot-to-cool phase transition layer moves quickly outside the simulation domain. At t > 350 Myr, only the cool-to-hot phase transition occurs within our simulation domain, and the gas escapes mostly in the form of the hot phase (see also Figure 11). While not followed in our simulations, the hot-to-cool phase transition can still occur far above the disk midplane.

4.1.2. Momentum Transfer

Figure 7 shows net vertical momentum gain/loss rates per unit area for each phase. The weight of each phase is included as a sink, but only the weight of the cool phase is significant. To understand the plot, it is important to keep in mind that the vertical momentum flux ρ v_z has a sign and the positive value (red) means the gain of the upward momentum flux and the loss of the downward momentum flux. The noICM model shows the change of sign in each phase at about z = 0. The hot phase loses its upward momentum flux in the upper disk (blue in z > 0) and downward momentum flux in the lower disk (red in z < 0)—in short, the hot phase loses the outward momentum flux and the cool phase gains it. In other words, as SN-driven superbubbles expand, the cool phase is accelerated outward by transferring the outward momentum flux of the hot phase. Note that the weight term of the cool phase dominates the net gain (see also Figure 8). This means that the continuous momentum transfer from the hot phase (or SN momentum injection) enables the cool phase (mass-dominating component) to remain vertically extended (more extended than by the thermal and magnetic support).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the ICM models, as soon as the ICM−ISM interface reaches the midplane z = 0, Figure 7 becomes just upward momentum gain/loss in the upper disk. As the ICM pressure gets stronger, the hot ICM dominates the momentum transfer, occurring in a larger region of the upper disk. In the weak ICM models, the momentum transfer is still limited within z < 1–2 kpc. In these models, the vertical momentum flux gain in the cool phase is simply counterbalanced by the increased weight, confining the disk within the simulation domain (no significant stripping).

In the ICM-P7h model, significant hot-to-cool momentum transfers occur all over the upper disk as the ICM fills up a larger volume in this region. In contrast to the mass transfer, which shows a clear dichotomy of the cool-to-hot and hot-to-cool phase transition layers at different heights (Figure 6(d)), the vertical momentum is always transferred from the hot phase to the cool phase. However, the momentum transfer to the intermediate phase closely follows the mass transfer; it gains momentum when it gains mass. On the one hand, near the ISM–ICM interface, where the mass and momentum transfers have opposite signs in both cool and hot phases, the cool phase is accelerated and gains vertical momentum, while it is shredded and losing mass by hydrodynamic instabilities. Here, the momentum transfer is mainly due to the drag force from the hot gas. The intermediate phase is newly populated by the accelerated and shredded cool phase, gaining both mass and momentum. On the other hand, in the upper region of the disk, the hot phase, together with the intermediate phase, is continuously mixed and cooled back to the cool phase, delivering both mass and momentum from hot (and intermediate) to cool phase. There, the momentum transfer is more mixing dominated.

To help us further understand vertical momentum transfer in detail, Figure 8 shows a decomposition of each term in Equation (17) in a different phase. From top to bottom, each row shows (a) net gain/loss (identical to Figure 7(d)), (b) time-dependent term, (c) kinetic, (d) thermal, and (e) magnetic flux terms, and (f) sink term due to weight. In the third column (hot), the ICM shock front marching upward is visible at late times. The hot ICM is thermalized at the shock—gaining thermal flux and losing kinetic flux. The thermalized hot ICM loses its thermal flux (pressure) owing to the interaction with the ISM. This turns into the kinetic flux gains of all phases. The hot kinetic flux gain is due to the mass loading from the cool phase (Figure 6), and the cool kinetic flux gain represents acceleration by drag force. Note that the thermal flux gains in the cool and intermediate phases are minimal, as the majority of the thermal flux transferred to these phases is radiated away (see Section 4.1.3). The magnetic flux term is subdominant, while this can be as important as the thermal term in weaker or no ICM models at the midplane. The weight term is dominated by the cool phase as shown in panel (f), which is larger than the kinetic flux gain of the cool phase at early times. At late times, the weight term becomes comparable to the kinetic flux gain of the cool phase, implying effective stripping. In short, the hot ICM's momentum flux transferred to the cool phase has been used to extend the ISM vertically and support the increased weight. The kinetic momentum flux of the cool phase shows consistent gains at later times (active stripping phase).

In the ICM-P14 model, the ICM pressure is so strong that the momentum gain in the cool phase is now dominated by the kinetic flux gain, resulting in actual acceleration of the entire cool phase at a velocity larger than the escape velocity of the simulation domain. The majority of the cool ISM is ablated, while it is accelerated and quickly stripped away from the simulation domain.

4.1.3. Energy Transfer

In the simulations, in addition to the energy added by SNe and ICM inflows, there is radiative heating by FUV radiation in the cool phase. However, this radiative heating is balanced by cooling within the same phase. The remaining energy transfer across thermal phases is simple: always from hot to cooler phases regardless of the source of energy in the hot phase. As an example, Figure 9 shows the decomposition of the energy transfer terms for the ICM-P7h model. The sink term due to radiative cooling in panel (e) is much larger than the residual energy flux, which is then shared by the actual kinetic flux gain in panel (c) and work done against gravity in panel (f) of cooler phases. As seen in Figure 8, there is little energy that ends up arriving in the thermal energy/pressure. In this particular model, the source term due to SNe is only significant at early times. Immediately after the ICM−ISM interface reaches the midplane, the overall energy loss is much larger than the explicit source term by SNe (panels (a) and (d)), implying that it is the energy mostly delivered by the ICM inflows (panel (c)). Compared to the noICM model, the additional energy inputs from the ICM enhance radiative cooling in all phases, including, while limited, cooling in the hot phase. We discuss the enhancement of X-ray luminosity as a function of RPS strengths (Section 6.1.2).

Zoom In Zoom Out Reset image size

Regarding the stripping of the ISM, the main question is how the cool ISM gets accelerated. If the cool ISM is accelerated as a whole by the drag force (or ram pressure force) from the ICM, the cool ISM simply gains outward momentum as is. In reality, the hot ICM penetrates through low-density channels in the multiphase ISM. Large relative velocities between the hot ICM and the cool ISM cause hydrodynamical instabilities, shredding the cool ISM and creating turbulent mixing layers. At the same time, strong radiative cooling due to the high cooling rate of the mixed gas in the mixing layers results in mass (and associated momentum and enthalpy) influx from the hot to cool gas. The competition between shredding and hot gas cooling leads to either net gain or loss of the cool gas mass. As the velocity associated with the mass gain (hot gas velocity) is in general significantly larger than the velocity associated with the mass loss (cool gas velocity), the cool ISM is accelerated by the momentum transfer associated with the mixing, while the direct acceleration is still in play.

As we show in Section 4.1.1, there is a significant phase transition that occurs in the large regions of the ISM disk, especially in the strong ICM models, with corresponding momentum transfers seen in Section 4.1.2. It is thus clear that the mixing must play a role. In this subsection, we seek evidence that the mixing-driven momentum transfer is actually a dominant mechanism for the cool ISM acceleration and stripping.

The mixing-driven momentum transfer simply means that the more the hot phase mixes in, the faster the cool phase moves. If this is the dominant mechanism, there must be a linear correlation between the velocity of the cool ISM and the mass fraction of the hot ICM (Equation (7); see also Schneider et al. 2020; Tonnesen & Bryan 2021). The same is true for the SN-driven outflows. In our simulations, however, the SN ejecta tracer field is a less sensitive probe of the mixing-driven momentum transfer since the SN tracer field has been accumulated in all gas phases over many star formation–feedback cycles. On the other hand, the ICM mass fraction provides a telltale sign of the acceleration of the cool phase by the mixing of the hot ICM. Given the large difference in the metallicity of the ICM and ISM, this imprints a noticeable difference in the metallicity of the fast-moving cool ISM.

We choose two epochs (t = 260–280 Myr and t = 360–380 Myr) and select the cool phase within z = 1–2 kpc. Figure 10 plots the mass distribution in the metallicity Z and vertical velocity ${v}_{z}^{\mathrm{cool}}$ plane (first and third columns) and the ICM mass fraction ${f}_{\mathrm{ICM}}^{\mathrm{cool}}$ and vertical velocity ${v}_{z}^{\mathrm{cool}}$ plane (second and fourth columns). The solid lines in the second and fourth columns are corresponding predictions from Equation (7). The corresponding predictions for metallicity assuming ${f}_{\mathrm{SN}}^{}=0$ using Equation (6) are shown in the first and third columns. For the fourth column, we apply an offset using the mean ${f}_{\mathrm{ICM}}^{\mathrm{cool}}$, as the baseline ICM fraction in the cool phase is nonzero and gradually increases. Except for the ICM-P1 model, all show tight correlations between the ICM mass fraction in the cool phase and the outward vertical velocity of the cool phase in the early epoch. The high-velocity component accelerated by the ICM appears as the low-metallicity component (anticorrelated component) in the Z–${v}_{z}^{\mathrm{cool}}$ panels. Similarly, the mixing of the hot gas created by SNe produces a correlation between the metallicity excess and vertical velocity, but this correlation is less clear. The mixing-driven acceleration by the SN-origin hot gas is only visible at the early epoch of the ICM-P3h model—the high-velocity component correlated with metallicity Z/Z_⊙ ≥ 1.2. The contribution of the ICM mixing increases as the ICM pressure increases, dominating the SN-origin mixing component. Generally, the outflow velocity is above the simple linear prediction, indicating that the acceleration by ram pressure still contributes on top of the mixing-driven momentum transfer. The direct ram pressure drag is more important in the earlier evolution and the stronger ICM model when the ISM reacts to the ICM as a whole. At late times, cool clouds are more fragmented so that the cross section to the ICM inflows gets smaller while the surface area of the mixing layers gets larger.

Zoom In Zoom Out Reset image size

In the late epoch (we omit the ICM-P14 model since there is almost no cool phase gas left at this epoch), a significant fraction of the cool phase gas falls back in the weak ICM models. Since this gas has been pushed upward mostly by the ICM previously, ${f}_{\mathrm{ICM}}^{\mathrm{cool}}$ is generally enhanced. The correlation between ${f}_{\mathrm{ICM}}^{\mathrm{cool}}$ and ${v}_{z}^{\mathrm{cool}}$ becomes less clear in the ICM-P3h model owing to the lack of continuous acceleration by the ICM but remains tight in the ICM-P7h model. At this epoch, the mean metallicity of the cool phase in both ICM-P3h and ICM-P7h models is greatly reduced compared to that of the noICM model shown in contours. The metallicity of the cool phase decreases at least 0.1 dex in both models owing to the ICM mixing.

We now ask, when the ISM is stripped by mixing with phase transition, which phase dominates in the outflows, hot (by shredding and escape before cooling) or cool (by cooling of stripped gas)? In Section 4.1.1, we show that the mass-dominating cool phase is shredded near the ICM−ISM interface and first stripped in the form of the intermediate and hot phases. Then, significant cooling occurs within the simulation domain (z ∼ 2−3 kpc) for the marginally strong ICM model (ICM-P7), but the majority of the stripped gas escapes the simulation domain before cooling in the strongest pressure model (ICM-P14).

To provide a more quantitative view, in Figure 11 we measure the outgoing fluxes (F_q,out) normalized by the injected fluxes in the strong ICM models at z = 3 kpc. The injected fluxes include the ICM inflow fluxes (F_q,in) measured at z = − 3 kpc and those injected by SNe (${F}_{q,\mathrm{SN}}$). ⁶ Throughout the simulation, the contributions of SN feedback to total injected mass, momentum, and energy fluxes are 5%, 28%, and 17% for the ICM-P7h model and 3%, 14%, and 6% for the ICM-P14 model, respectively. In the top row (ICM-P7h), there are significant and continuous outflowing mass, momentum, and energy fluxes from the cool phase. This is in stark contrast to the weak ICM models (and the noICM model), where all outflow fluxes at z = 3 kpc are dominated by the hot phase at a level of a few percent of the injected fluxes (e.g., Kim et al. 2020b) and frequently truncated by inflows. The cool outflows in the ICM-P7h model consist of both directly accelerated cool phase (mostly at early times) and additional cool phase created in the hot-to-cool phase transition layer (see Figure 6). Roughly 10%–20% of the injected momentum flux is transferred to the outflowing cool phase. The energy flux in the hot phase at this distance is about 10%–20% of the injected flux, as significant thermal energy is transferred to the cooler phases and then radiated away. As a consequence, the energy flux in the cool and intermediate phases is much lower (a few percent) again because of large cooling losses (see Figure 9).

Zoom In Zoom Out Reset image size

In the bottom row, strong outgoing fluxes in the cool phase exist only at very early times. The hot phase soon dominates all fluxes as the cool gas is shredded and evaporated to the hot phase. This mass loading (or entrainment) to the hot outflowing gas increases the hot gas mass flux, at maximum, by a factor of two compared to the injected mass flux. Similar behavior is seen in the ICM-P7h model at t > 400 Myr. The maximum momentum transfer efficiency to the cool phase reaches up to 50%. To summarize, both ICM-P7h and ICM-P14 models show the mixing-driven ISM stripping soon after the direct cool phase acceleration, and the cool phase dominates the outflows in the early epoch. Then in ICM-P14, the hot phase takes over the outflows, whereas the cool phase outflow is maintained to some extent in ICM-P7.

Finally, we also make use of the metallicity of the outflowing gas to find the contribution of the ICM in accelerating the gas. Figure 12 plots the metallicity of the outflowing gas (v_z > 0) at z = 3 kpc for the (a) weak and (b) strong ICM models, along with the noICM model in both panels as a reference. In the ICM-P1 model, the metallicities of the outflow in all three phases are essentially unchanged from those in the noICM model. This is expected, as the ICM cannot penetrate directly to the upper disk. The ICM is mixed into the ISM near the midplane, but the mass flux is insignificant. The outflowing gas is mostly driven out by SNe with enhanced metallicity.

Zoom In Zoom Out Reset image size

The ICM makes a noticeable difference in the ICM-P3h model. At t ∼ 300–330 Myr, the cool outflow metallicity is clearly reduced, while the hotter phases still show metallicities similar to those in the noICM and ICM-P1 models. The reduced metallicity in the cool outflow means the ICM mixing-driven acceleration as evidenced in Figure 10. At later times (t > 350 Myr), the outflow metallicity of all phases is significantly reduced, equally in all phases, and increases again. The decrease of metallicity indicates that the mixing of the ICM into the ISM is the main driver of the outflows, while the later increase of the metallicity signals that the SN feedback plays a major role in driving outflows.

In the strong ICM models shown in Figure 12(b), the metallicity of outflows is reduced at all times and keeps decreasing. This makes it plain that SNe in the strong ICM models are not a major driver of outflows, except very early in time in the ICM-P7h model as seen in Figure 10. In Figure 12(b), we also find that the outflow metallicity in the cool and intermediate phases is very similar for both strong ICM models, while the hot outflow metallicity is more reduced with stronger ICM pressure. The distribution of ${f}_{\mathrm{ICM}}^{\mathrm{cool}}$ shown in Figure 10 shows ${f}_{\mathrm{ICM}}^{\mathrm{cool}}\lt 0.1-0.2$. Having demonstrated that the mixing is the main mechanism to drive outflows (or stripping), the limited range of the ${f}_{\mathrm{ICM}}^{\mathrm{cool}}$ implies that the outflowing cool gas would have been ablated and evaporated before mixing in too much hot ICM. The maximum ICM fraction then sets the maximum velocity and minimum metallicity difference of the cool gas accelerated by the ICM.

The general expectation is that strong ICM ram pressure that can strip the gas in galaxies will reduce SFRs in galaxies. At the same time, mild ICM ram pressure that compresses the gas in galaxies may enhance SFRs. Indeed, in our simulations, we find both effects depending on the ICM strengths and evolutionary stages. In short, SFRs are enhanced locally (inside the truncation radii) and temporarily (before active stripping), but the gas stripping eventually quenches star formation.

Figure 13(a) plots the time evolution of SFR surface density defined by the total mass of young stars formed in the last t_bin = 10 Myr:

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{SFR}}({\rm{\Delta }}t={t}_{\mathrm{bin}})\equiv \frac{{\rm{\Sigma }}{m}_{\mathrm{sp}}({t}_{m}\lt {t}_{\mathrm{bin}})}{{L}_{x}{L}_{y}{t}_{\mathrm{bin}}},\end{eqnarray} \tag{ 19 }$$

where m_sp and t _m are the mass and mass-weighted mean age of the sink particle representing star clusters, respectively. This roughly corresponds to SFRs traced by Hα (e.g., Kennicutt & Evans 2012). Figure 13(b) shows the box-and-whisker plots, presenting the distributions of Σ_SFR over two periods separated by t = 330 Myr, before and after quenching of star formation in the strong ICM models. The enhancement of SFRs compared to the noICM model in the early epoch is common in all models with the ICM. At later times (t > 330 Myr), such enhancement of SFRs persists in the weak ICM models, while the gas stripping quenches star formation in the strong ICM models. The enhancement levels in Σ_SFR are ∼30%–50% in the weak ICM models for more than 200 Myr explored in this paper. In addition, the temporal modulation of Σ_SFR in these models gets stronger with higher peaks.

Zoom In Zoom Out Reset image size

The enhancement of SFRs in the early epoch is mainly due to the compression of the overall ISM disk in the vertical direction. The introduction of the ICM inflows simply pushes the ISM from the lower disk to the midplane, effectively supplying more gas for star formation. In the weak ICM models, this additional gas remains near the midplane, where the majority of star formation takes place. However, strong ICM inflows can blow away the ISM altogether in ∼100 Myr. Figures 13(c) and (d) show the time evolution and the box-and-whisker plots of the dense gas surface density Σ_dense selected by n_H > 10 cm⁻³. The first compression increases the peak dense gas mass at ∼270 Myr by about a factor of two in all ICM models. The corresponding enhancement of SFRs is delayed by ∼10 Myr, a freefall time of gas at n_H = 10 cm⁻³ (Mao et al. 2020). The enhancement of Σ_dense persists in the weak ICM models. In the strong ICM models, however, the dense gas mass quickly decreases over time owing to shredding and stripping by the ICM. In the ICM-P7 model, the dense gas still exists for a longer time than the ICM-P14 model. Some of the dense gas pushed far above the disk manages to form stars at late times (see Section 5.2).

One of the intriguing properties commonly found among the RPS galaxies is star-forming patches outside the stellar disk that remains intact. In Figure 14(a), we show the vertical distance of the newly formed star clusters (sink particles) from the midplane z_sf over time. The size of symbols represents the mass of star clusters. The black star symbols are for the noICM model.

Zoom In Zoom Out Reset image size

On the one hand, for the strong ICM models, the bulk ISM keeps moving away from the midplane. As a consequence, z_sf increases over time. This continues for the ICM-P14 model, while the ICM-P7 and ICM-P7h models show turnover. Although one may get an impression that the ISM continuously moves upward in these models (see Figure 2), the main gas reservoir is fragmented, and a large chunk of dense gas falls back (Figure 5(b)). As a result, two star-forming sites near and far from the midplane are visible at late times of the ICM-P7 model. On the other hand, for the weak ICM models, as more and more gas moved the upper disk, the ISM weight shortly dominates the ICM pressure. The entire ISM disk falls back, and so does the star formation location. This introduces larger-amplitude vertical oscillations of z_sf in the ICM-P3 and ICM-P1 models than the noICM model, in which a small-amplitude vertical oscillation is naturally introduced by the asymmetry (Figure 1).

Figure 14(b) plots the metallicity of sink particles. Each symbol is now color-coded by z_sf shown in Figure 14(a). In the noICM model, the metallicity of new star clusters increases over time, as the star-forming gas is continuously metal-enriched by mixing of the high-metallicity SN ejecta. The injected SN ejecta first goes into the hot phase and then quickly cools and mixes into the cool phase (see the top row of Figure 6). We find that the metallicity within the cool phase is nearly homogeneous, implying the efficient mixing of SN ejecta to the cold, star-forming gas. The metallicity of new stars born within the main ISM disk follows a common enrichment trend even with the ICM inflows, implying that they are born in the genuine ISM. However, in the strong ICM models, star clusters formed in the stripped gas far from the midplane at late times (marked by black circles) show lower metallicities compared to the enrichment trend. These star clusters are born in the gas that has experienced significant mixing with the low-metallicity ICM. The gradual mixing of the ICM in the ICM-P3 and ICM-P3h models also reduces the metallicity of new stars at late times (t > 400 Myr), while higher SFRs with insignificant mass contribution from the ICM in the ICM-P1 model result in an even higher metallicity of new stars.

While reduced, the metallicity is still much higher than the ICM metallicity, implying that the composition of the star-forming gas in the extraplanar region is dominated by the genuine ISM. In our simulations, there is no sign of the complete shredding and recondensation of the star-forming cold gas in the stripped tails within the simulation domain z < 3.5 kpc, which will be generally true for the extraplanar star formation within a few kiloparsecs away from the disks of RPS galaxies.

The multiphase nature of the ICM−ISM interaction is often neglected when developing theoretical understandings based on simple analytic models, although multiwavelength observations have revealed the multiphase gas involved in RPS galaxies such as cold molecular gas via CO (e.g., Vollmer et al. 2008; Verdugo et al. 2015; Moretti et al. 2018b; Jáchym et al. 2017, 2019), cold and warm neutral gas via H i (e.g., Cayatte et al. 1990; Chung et al. 2009; Scott et al. 2010), warm ionized gas via Hα (e.g., Fumagalli et al. 2014; Boselli et al. 2016), and hot gas via X-rays (e.g., Sun et al. 2010; Poggianti et al. 2019). We show that the mass, momentum, and energy transfer from the hot ICM to the ISM via gas mixing is likely the dominant mechanism for stripping in our multiphase RPS simulations. This is a qualitatively different process from that of a simple acceleration due to ram pressure without phase transition and mixing. A wealth of observational signatures will be imprinted on different gas phases.

6.1.1. Imprints in Metallicity of Stripped Tails

The main observational imprint of the mixing-driven acceleration model is on the anticorrelation of metallicity and cool gas velocity in the stripped tails (Figure 10). If the ICM has distinctively lower metallicity than the ISM as we assumed, the fast-moving part of the stripped gas should have lower metallicity than the genuine ISM. ⁷ Assuming ${f}_{\mathrm{SN}}^{}=0$, Equations (6) and (7) give the slope in the metallicity and velocity correlation dZ/dv = (Z_ICM − Z_ISM)/v_ICM. For Z_ICM/Z_ISM = 0.1 and v_ICM = 1000 km s⁻¹, this implies that the metallicity reduction in the mixed cool gas is

$$\begin{eqnarray}&&\frac{{Z}_{\mathrm{mix}}}{{Z}_{\mathrm{ISM}}}=1-0.09\left(\frac{{\rm{\Delta }}v}{100\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right),\end{eqnarray} \tag{ 20 }$$

resulting in roughly 10% reduction for 100 km s⁻¹ difference in outflow velocity. The metallicity reduction decreases by a factor of two if Z_ICM/Z_ISM = 0.3 and v_ICM = 1400 km s⁻¹.

It is also noteworthy that the mixed ICM fraction in the cool gas cannot be arbitrarily high. Although efficient cooling helps to keep the mixed gas cool, the cool ISM can evaporate if the energy flux from the hot ICM is too large to be radiated away in the mixing layer. In our models, the difference of ${f}_{\mathrm{ICM}}^{}$ between high- and low-velocity cool gas is typically less than 0.1. This limits the dynamic range of outflow velocity less than < 200 km s⁻¹ even in ICM-P14 with v_ICM = 2000 km s⁻¹. Similar results are also seen in numerical simulations of an isolated galaxy that is experiencing the hot ICM inflows (Tonnesen & Bryan 2021), where they follow accelerated clouds in the stripped tails of scales of tens to hundreds of kiloparsecs. Although the stripped tails that traveled very far from the disk can have much larger ${f}_{\mathrm{ICM}}^{}$ of clouds up to 0.8, the range of ${f}_{\mathrm{ICM}}^{}$ at a particular position is limited to ∼0.1–0.2, translating into the velocity difference ≲ 200–400 km s⁻¹. In addition to the maximum ICM fraction, there must be a significant fraction of total gas across a wide range of outflow velocities for such mixed gas to be visible. Since the gas mass fraction is sharply decreasing at high velocities in general (Kim et al. 2020a), observable velocity ranges and hence the metallicity differences can be further limited. Finally, it can possess a cleaner signal only if the stripping occurs quicker than the enrichment by SNe. Such a condition favors strong RPS galaxies. All the above makes the signal in gas-phase metallicity difference at different outflow velocities we are searching for very small (of an order of 10% or less).

If the mixing-driven stripping continues beyond the immediately stripped tails we model here, there must be a well-defined trend in the mixed gas fraction as a function of distance across very long tails of RPS galaxies (often dubbed jellyfish galaxies). Indeed, global RPS galaxy simulations forming such long tails (>100 kpc) show a correlation between clouds' distance and ICM mass fraction (Tonnesen & Bryan 2021). With an assumption that the mixing rate is constant over time, Tonnesen & Bryan (2021) laid out a simple model for the ICM mass fraction as a function of distance, which qualitatively agrees with the increasing ICM mass fraction in clouds farther away in their simulations. The recent analysis of MUSE observations of RPS galaxies shows that warm ionized gas metallicities decrease as a function of distance from stellar disks (Franchetto et al. 2021). The stripped gas, in reality, would experience much more complicated dynamical and thermal evolution, including deceleration by gravity, evaporation/fragmentation, and perhaps recondensation/growth by cooling in the mixing layers. The simple extrapolation of the clouds' velocity and ICM mass fraction at very far distances may not work well in predicting velocity and metallicity correlations quantitatively. Still, potentially illuminating results in Tonnesen & Bryan (2021, see their Figure 9) are that the slope in the cold clouds' velocity and ICM fraction correlation remains nearly linear over a large range of distances. Again, high-precision measurements of metallicity across velocity channels to measure the slope in the v and Z correlation will be the most direct ways to confirm whether the mixing-driven acceleration is the dominant mechanism for the RPS.

RPS galaxies often show star formation activity outside the main, old stellar disk (e.g., Kenney & Koopmann 1999; Sun et al. 2007; Poggianti et al. 2016). The mixing of the ICM also creates an imprint on the metallicity of stars formed in the extraplanar region. In the strong ICM models, the extraplanar star formation in the stripped tails occurs 2–3 kpc above the stellar disk, creating star clusters with lower metallicity (by 0.05–0.1 dex) than those formed in the disk (see Figure 14). Observationally, the relative difference of stellar metallicities between young stars from the unstripped inner part and stripped extraplanar region can be compared. If the intrinsic metallicity gradient of galaxies is subtracted, stacking analysis may enhance a potential signal.

6.1.2. Imprints in Gas Phases

RPS not only simply strips the cool ISM as is but also involves significant phase transition from cool to hotter phases (see Section 4). In the regions that experience strong RPS, the shredded cool gas escapes the simulation domain (stripped from galaxies) before it cools back (e.g., ICM-P14). This is manifested in H i deficiency in RPS galaxies (Chung et al. 2009; Ramatsoku et al. 2019, 2020). At the same time, the mass-loaded hot gas gets brighter in X-rays. The fate of such stripped gas is not traced in our simulations, but it is possible to cool back the stripped gas and form H i/Hα tails (e.g., Ramatsoku et al. 2019, 2020). In fact, large-scale simulations do show late-time cooling at more than tens of kiloparsecs away from the disk (Tonnesen & Bryan 2012, 2021; Lee et al. 2022), which might be responsible for the long, extended tails seen in Hα and CO (e.g., Vollmer et al. 2008; Lee et al. 2017; Jáchym et al. 2017, 2019).

Recently, sensitive, high-resolution observations of molecular gas tracers reveal the prevalence of extraplanar molecular gas in RPS galaxies (e.g., Vollmer et al. 2008; Jáchym et al. 2017, 2019; Moretti et al. 2018b). Regarding the origin of extraplanar molecular clouds, it is unclear whether they are remnants of directly stripped molecular gas from the ISM disk or destroyed and reformed in the extraplanar region. Lee & Chung (2018) used high-resolution ALMA observations of NGC 4522 and detected ¹³CO in the extraplanar molecular clumps at ≲a few kiloparsecs above the stellar disk near truncation radii. Given the relatively short formation time of molecular gas (∼10 Myr) compared to the stripping timescale (∼100 Myr), they concluded that both scenarios are feasible.

Our simulations lack the resolution and physical processes to follow the molecular species explicitly in the simulation. Instead, if we consider the dense gas (n_H > 10 cm⁻³; see Figure 13) as a proxy of the molecular gas, we find that in the ICM-P7h model, where the later-time extraplanar star formation occurs, (1) the dense gas fraction is comparable to the noICM model and (2) most of the cold gas is located at 1−3 kpc above the midplane within which stars form at late times. This cold, dense gas is not completely shredded and recondensed, but a significant fraction of the ICM has been mixed in, which is evidenced by the lower metallicity of star clusters formed in the extraplanar regions (Figure 14), up to ${f}_{\mathrm{ICM}}^{}\sim 0.15$. Our simulations suggest that the extraplanar molecular gas (not those in the stripped tails farther than a few tens of kiloparsecs from the disk) mostly originates from the gas directly stripped from the disk. Still, the ICM mixing is important to accelerate the molecular gas (Figure 10). Given that the ICM-P14 model quickly runs out of its dense gas, as the remaining enthalpy flux from the hot ICM after cooling is large enough to evaporate many small cold clouds, the marginally strong ICM condition can be optimal for pushing the molecular gas outward without destroying it. This translates into an expectation that the extraplanar molecular gas in active RPS galaxies can be most abundant near the truncation radii, which seems to be consistent with observations (Lee & Chung 2018; Moretti et al. 2018b).

Another potentially strong observable signature of RPS can be an enhanced diffuse X-ray brightness (e.g., Kouroumpatzakis et al. 2020; Campitiello et al. 2021). The diffuse thermal X-ray emission from the hot ISM is expected to correlate with SN rates and hence SFRs. From the observations of star-forming galaxies, Mineo et al. (2012) obtained a linear scaling relation between the diffuse (excluding resolved X-ray binaries) soft X-ray luminosity in 0.5–2 keV (L_{X,0.5–2 keV}) and SFR (${\dot{M}}_{* }$) as ${L}_{{\rm{X}},0.5\mbox{--}2\,\mathrm{keV}}/{\dot{M}}_{* }=8.3\times {10}^{38}\,\mathrm{erg}\,{{\rm{s}}}^{-1}{\,({M}_{\odot }\,{\mathrm{yr}}^{-1})}^{-1}$ (see also Ranalli et al. 2003; Mineo et al. 2014). Assuming a canonical SN energy of 10⁵¹ erg, the SN energy injection rate is ${\dot{E}}_{\mathrm{SN}}=3.3\,\times {10}^{41}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,({\dot{M}}_{* }/({M}_{\odot }\,{\mathrm{yr}}^{-1}))$ for the standard IMF (1 SN per 100 M_⊙; e.g., Kroupa 2001). Then, the observed relation means the SN energy to soft X-ray conversion efficiency (or "soft X-ray efficiency" in short) of ${\epsilon }_{\mathrm{SN}\to {\rm{X}}}\equiv {L}_{{\rm{X}},0.5-2\,\mathrm{keV}}/{\dot{E}}_{\mathrm{SN}} \tilde 0.25 \%$. ⁸ The analysis of the TIGRESS suite (Kim et al. 2020b) shows similar soft X-ray efficiencies of ${\epsilon }_{\mathrm{SN}\to {\rm{X}}} \tilde 0.1 \% -0.2 \%$, with higher efficiency for higher SFR surface density (C.-G. Kim et al., in preparation; see Figure 15).

Zoom In Zoom Out Reset image size

Here we calculate the X-ray surface brightness for the ICM models and compare them with the TIGRESS suite. We first obtain the X-ray emissivity in the soft X-ray band (0.3−2 keV) for each cell using the apec.v2 table (Foster et al. 2012) adopted in yt (Turk et al. 2011). The gas metallicity is fixed to solar metallicity. We then integrate total soft X-ray luminosity and divide it by the area to get the mean X-ray surface brightness. Figure 15 plots as colored points the soft X-ray surface brightness ${{ \mathcal S }}_{{\rm{X}},0.3-2\,\mathrm{keV}}$ as a function of Σ_SFR,40 (Equation (19) with t_bin = 40 Myr, the duration of SNe in each cluster). The weak ICM models are all consistent with the relation with soft X-ray efficiency of 0.1% (lower dashed line) as in the noICM model and the TIGRESS suite. However, the strong ICM models show the enhancement of X-ray surface brightness consistent with a lower X-ray efficiency for the ICM inflows as ${{ \mathcal S }}_{{\rm{X}},0.3-2\,\mathrm{keV}}/({\dot{E}}_{\mathrm{ICM}}/A)\sim 0.05 \%$ (horizontal dotted lines), where the ICM total energy flux is ${\dot{E}}_{\mathrm{ICM}}/A\,=0.5{\rho }_{\mathrm{ICM}}{v}_{\mathrm{ICM}}({v}_{\mathrm{ICM}}^{2}+5{c}_{s,\mathrm{ICM}}^{2})$. We note that the X-ray surface brightness of the pure ICM is much lower than that from the ICM−ISM interaction, mainly due to the ICM's low density. Somewhat lower X-ray efficiency of the ICM can be understood because of a larger mixing area involved in the ICM−ISM interaction than that of superbubbles driven by SNe.

In RPS galaxies, the diffuse X-ray brightness can then be enhanced by an order of magnitude compared to that expected purely from SNe before the majority of the shock-heated gas is stripped away (shown as decreasing X-ray at lower SFRs in the strong RPS models).

The detailed look of the RPS process as multiphase gasdynamical interaction is reminiscent of a collection of shock/wind−cloud interactions. The main question of shock/wind−cloud interaction studies is how cold gas can be accelerated before it is completely shredded. In adiabatic cases, the drag/acceleration timescale is always shorter than the cloud-crushing timescale (Klein et al. 1994). In other words, the energy transferred from the hot shock/wind to the clouds is fully retained to heat up the clouds, while surface instabilities shred them (e.g., Klein et al. 1994; Mac Low et al. 1994; Cooper et al. 2009; Scannapieco & Brüggen 2015; McCourt et al. 2015). In order to prolong the cloud lifetime, several mechanisms have been proposed, including radiative cooling (e.g., Cooper et al. 2009) and magnetic fields (in both wind and cloud; e.g., Dursi & Pfrommer 2008; McCourt et al. 2015; Sparre et al. 2020). Recently, it has been realized that when clouds are large enough and cooling is strong, all the enthalpy flux can be radiated away, while the significant mass and momentum of the hot phase are mixed and added into the cool phase without completely shredding the clouds. This allows the clouds to keep growing even while they are being accelerated by shock/wind−cloud interactions (e.g., Armillotta et al. 2016; Gronke & Oh 2018, 2020; Kanjilal et al. 2021; Li et al. 2020; Sparre et al. 2019, 2020). It is certainly true that in our simulations the chunk of the ISM that is facing the ICM inflows is large (a few hundreds of parsecs). The critical size above which the cool clouds can grow by cooling of the hot gas proposed by Gronke & Oh (2018) is

$$\begin{eqnarray}&&{R}_{\mathrm{crit}}\approx 2\mathrm{pc}\displaystyle \frac{{T}_{\mathrm{cl},4}^{5/2}{{ \mathcal M }}_{\mathrm{wind}}}{{P}_{3}{{\rm{\Lambda }}}_{\mathrm{mix},-21.4}}\displaystyle \frac{\chi }{100},\end{eqnarray} \tag{ 21 }$$

where T_cl,4 ≡ T_cl/10⁴ K is the cloud temperature, P₃ ≡ P/(10³ k_B cm⁻³ K) is the ambient pressure, ${{ \mathcal M }}_{\mathrm{wind}}$ is the hot wind Mach number, Λ_mix,−21.4 = Λ(T_mix)/(10^−21.4 erg s⁻¹ cm³) is the cooling coefficient at the temperature of the mixed gas, and χ is the density contrast between wind and clouds.

The critical size ⁹ is of order of a few to tens of parsecs at the typical conditions of our simulations with P/k_B ∼ 10⁴−10⁵ cm⁻³ K, χ = 10²⁻³, and trans-to-subsonic wind Mach number (note that the ICM Mach number at injection was supersonic, but it is quickly thermalized and becomes subsonic at the time of interaction near the midplane). The bulk ISM from the first interaction cannot be completely shredded, while there is continuous shredding at the interfaces as the ICM penetrates through low-density channels (see Figure 4). In the later-time evolution, the strong ICM models successfully stripped the majority of the cool ISM from the disk midplane. There are smaller, fragmented cold cloudlets embedded in the ICM inflows (see Figure 5), which are vulnerable to shredding/evaporation. This results in a broad cold-to-hot phase transition layer seen in Figures 6(d) and (e). However, these clouds' wakes meet other cool gas and add their mass back to the cool phase. This evolution is more equivalent to that seen in shock−multicloud interactions (Banda-Barragán et al. 2021, 2020) rather than the growth of cool clouds in idealized shock/wind−cloud interaction simulations where the mass is added from an infinite hot reservoir via cooling of the mixed gas.

The compression by the ICM can enhance SFRs by 30%–50% for a short period (a few tens of Myr), while the enhanced star formation is sustained in the weak ICM models (or inner part of an RPS galaxy) as the ISM remains compressed. In our strong models, representing the outer region of an RPS galaxy, star formation is quenched at timescales of ∼100 Myr.

Many previous simulations of RPS stripping galaxies including star formation recipes commonly show the enhanced star formation activity before quenching (e.g., Kronberger et al. 2008; Steinhauser et al. 2012; Ruggiero & Lima Neto 2017). Despite the qualitative agreements in the roles of RPS in star formation, global star formation enhancement found in many of these simulations is usually higher than ours by a factor of a few and persists longer (Kronberger et al. 2008; Steinhauser et al. 2012; Ruggiero & Lima Neto 2017). Keep in mind that earlier simulations in this category adopt a parameterized model for the ISM (cannot directly follow the gas phase cooler than 10⁴ K) and star formation (e.g., Cen & Ostriker 1992; Springel & Hernquist 2003) to model an entire galaxy in a wind tunnel (e.g., Kronberger et al. 2008; Steinhauser et al. 2012) or in a galaxy cluster (e.g., Ruggiero & Lima Neto 2017). Therefore, the SFRs obtained in previous global simulations can be sensitive to the adopted star formation recipes, although the global nature of such models (e.g., ICM wind inclination) can also be a reason for the difference (see Section 6.4).

Recently, Lee et al. (2020) presented simulations of an RPS galaxy with varying ICM inflow strengths and directions. Combined with adopted higher resolution (adaptive mesh refinement down to 20pc) and explicit ISM cooling and heating treatments, star formation in this work occurs in the cold, dense gas at number density above 100 cm⁻³, representing self-gravitating clouds, as in our simulations. In their moderate ICM inflow model, star formation at the outer region becomes suppressed for early 150 Myr, while the central region of the galaxy shows an enhancement of SFRs for several hundred Myr, compared to their NoWind case. The quantitative agreements with our models are encouraging and indicative of the importance of high-resolution modeling of star formation in the multiphase ISM.

Star formation enhancement prior to the quenching has been observed in RPS galaxies (e.g., Crowl & Kenney 2006; Kenney et al. 2014). Recently, Vulcani et al. (2018) reported a systematic enhancement of the SFR (0.2 dex) for 42 RPS galaxies compared to the counterpart galaxies. The spatially resolved SFRs have been estimated for some of those galaxies, showing signs of central SFR enhancement before quenching (Vulcani et al. 2020). In addition, Roberts & Parker (2020) identified 41 RPS candidate galaxies in the Coma Cluster and reported enhanced SFR (0.3 dex) of them. Meanwhile, Crowl & Kenney (2006, 2008) measured the age of the youngest stellar population at the H i truncation radii of RPS galaxies to estimate a quenching timescale—how long ago star formation has been quenched since the H i gas stripping. Crowl & Kenney (2008) derived the quenching timescale of a few hundred Myr for the Virgo RPS galaxies that appear to be currently undergoing active RPS. These results are broadly consistent with our results, while more spatially resolved analyses in observations are warranted for more quantitative comparisons.

In this work, we had to limit our simulation domain to a kiloparsec-size box to achieve high resolution with explicit treatments of ISM physics (Kim & Ostriker 2017, 2018). We thus cannot cover an entire galaxy or model a galaxy orbiting within a realistic ICM. Consequently, we missed a few important physical processes involved in RPS.

First of all, we have to fix the ICM inflow direction perpendicular to the disk, i.e., face-on interaction. The ICM inflow inclination can be arbitrary for galaxies infalling/orbiting in a cluster. If the interaction is more edge-on, the ICM may preferentially compress the inflow-side ISM, while stripping the extraplanar gas more easily (Roediger & Brüggen 2006; Lee et al. 2020).

Second, the local model cannot capture global geometrical effects, which may be especially important in the stripping process at the truncation radius. After rapid stripping of gas outside the truncation radius, continuous stripping occurs through global hydrodynamical instabilities (e.g., Roediger & Hensler 2005; Tonnesen & Stone 2014), in addition to local instabilities introduced by the penetrating ICM. Another interesting global effect is the inward radial migration of the stripped gas; the inner disk protects the tails from further interactions with the hot ICM. Such gas that is still bound to the galaxy will fall back (e.g., Schulz & Struck 2001; Tonnesen & Bryan 2009; Tonnesen & Stone 2014). In the future, more realistic, time-varying ICM inflows can be modeled, although global, cosmological models are needed for modeling of realistic variation of ICM ram pressure strengths and angles including the change of the orbits (e.g., Tonnesen 2019).

Third, although we vary the ICM ram pressure, we only consider a representative ISM disk condition similar to the solar neighborhood. It is generally expected that the ratio of the ICM ram pressure to the ISM anchoring pressure ${P}_{\mathrm{ICM}}/{{ \mathcal W }}_{\mathrm{GG}}$ is the main control parameter that determines the dynamical impact of RPS. However, the microphysics of the ISM (e.g., chemistry and hence cooling and heating processes) will be particularly important for RPS in the multiphase ISM as the volume filling factors of different phase ISM vary over different conditions. The cooling rate in the mixing layer is one of the main parameters that determines the properties of mixing (Tonnesen & Bryan 2021; Fielding & Bryan 2022). In this regard, even for the same relative ram pressure strength, metallicity can change the efficiency of cooling and hence overall evolution of the multiphase RPS process.

The major advantage of the local framework used in this study is detailed modeling of ISM physics, which can be improved further by the future extensions of the TIGRESS framework with radiation and chemistry (J.-G. Kim et al. 2022, in preparation). These capabilities are critically important to understand the evolution of the cold molecular gas (Girichidis et al. 2021). Although the current model follows gas at cold temperature T < 100 K, within which star formation is modeled, the questions of the molecular cloud stripping remain unanswered. Are they intact during the journey to the far extraplanar region? How long do they survive? Can molecular clouds form again in the stripped tails? Future work using the new TIGRESS framework will shed light on these questions.

Finally, we point out the importance of thermal conduction, which is currently missing. As RPS should be viewed as the hydrodynamical interaction between gas phases at large temperature differences, conductive heat flux can be the dominant energy flux from the hot ICM to the ISM. The thermal conductivity of the ionized plasma increases steeply with the temperature (Spitzer 1962), but at the same time the conductive heat flux can be limited by the magnetic fields in the ISM that may wrap around the cool clouds as the ICM inflow sweeps up (e.g., Orlando et al. 2008). Direct numerical simulations including anisotropic conduction within the TIGRESS framework, where self-consistent magnetic field structure in the turbulent ISM is modeled, are vital to understanding the role of thermal conduction in RPS galaxies.