Ram pressure stripping of the multiphase ISM: a detailed view from TIGRESS simulations

Ram pressure stripping (RPS) is a process that removes the interstellar medium (ISM) quickly, playing a vital role in galaxy evolution. Previous RPS studies have treated the ISM as single-phase or lack the resolution and physical processes to properly capture the full multiphase ISM. To improve this simplification, we introduce an inflowing, hot intracluster medium (ICM) into a self-consistently modeled ISM in a local patch of star-forming galactic disks using the TIGRESS framework. Our simulations reveal that the workings of RPS are not only direct acceleration of the ISM by ICM ram pressure but also mixing-driven momentum transfer involving significant phase transition and radiative cooling. The hot ICM passes through the low-density channels of the porous, multiphase ISM, shreds the cool ISM, and creates mixing layers. The ICM momentum is transferred through the mixing layers while populating the intermediate temperature gas and radiating thermal energy away. The mixed gas extends beyond galactic disks and forms stripped tails that cool back unless the ICM fluxes are large enough to prevent cooling until they escape the simulation domain. The mixing-driven momentum transfer predicts that the more ICM mixes in, the faster the ISM moves, resulting in the anti-correlation of outflow velocity and gas metallicity of the stripped ISM. The compression of the ISM disks due to the ICM ram pressure enhances star formation rates up to 50% compared to the model without ICM. With the ICM ram pressure higher than the disk anchoring pressure, star formation is quenched within ~100 Myr.


INTRODUCTION
The current concordance cosmology predicts that smaller-scale structures form first and large-scale structures form through hierarchical merging (e.g., Davis et al. 1985). In this scenario, the continuous merging of dark matter halos (and galaxies within them) is the chief channel of galaxy formation and evolution, which involves a variety of interactions of galaxies with each other and with their environments (Voit 2005). Among many types of interactions, the interstellar medium (ISM) in galaxies moving through the intrwoorak.c@gmail.com cgkim@astro.princeton.edu achung@yonsei.ac.kr acluster medium (ICM) can be removed by the ICM's ram pressure (Gunn & Gott 1972). This process, socalled ram pressure stripping (RPS), has been extensively studied due in part to its unique observational signatures identified by disturbed gaseous medium with undisturbed stellar components of galaxies (see Boselli & Gavazzi 2014;Boselli et al. 2021, for reviews). Also, RPS is known to effectively affect the ISM content of galaxies on a relatively short time scale (Abadi et al. 1999;Boselli et al. 2009), dramatically changing galaxy evolution in high-density environments.
The observational signatures of RPS are traced in multi-wavelengths. The atomic hydrogen (H I) 21 cm line has been a pioneering tool to identify RPS galaxies since H I gas is generally diffuse and extended well beyond the stellar disk, and hence vulnerable to the interaction with the surroundings. The early single-dish observations such as Haynes & Giovanelli (1984) found the cluster population to be overall deficient in H I compared to the field counterpart. More direct signatures of RPS such as gas truncation into the stellar disk and/or gas tails have been reported by many H I imaging studies (e.g., Cayatte et al. 1990;Kenney et al. 2004;Chung et al. 2009;Scott et al. 2010;Ramatsoku et al. 2020). RPS galaxies also show common signs including locally enhanced synchrotron radiation or ionized gas tails which can be observed through radio continuum, optical lines, or X-ray emission (e.g., Gavazzi et al. 2001;Vollmer et al. 2010;Sheen et al. 2017;Sun et al. 2010;Poggianti et al. 2017).
On the other hand, the case of molecular gas content, which is generally present in the inner part of galaxies with a higher density compared to the other ISM phases, does not show clear evidence of RPS (e.g., Kenney & Young 1989). Depending on the sample selection and the observational strategy, only a handful number of studies find molecular gas deficiencies (e.g., Fumagalli et al. 2009;Corbelli et al. 2012;Chung et al. 2017), and the impact of ram pressure on the molecular gas disk has long been under debate. However, more recent high-resolution radio observations begin to show that the morphological characteristics of RPS seen in H I disks such as asymmetry and compression are shared by CO disks (e.g., Lee et al. 2017;Lee & Chung 2018;Brown et al. 2021) and have revealed the enhancement of molecular gas (e.g., Moretti et al. 2018aMoretti et al. , 2020aCramer et al. 2021). These data imply that the molecular ISM is essentially affected by ram pressure in similar ways as more diffuse components although whether it is stripped along with atomic gas or not is still unclear. All these observations indicate that RPS is the process in which the multiphase ISM from cold molecular gas to hot ionized gas is involved.
This naturally raises a question on the effect of RPS in star formation, whether RPS is only effective in stripping the diffuse gas limiting the supply of gas for future star formation (e.g., Koopmann & Kenney 2004;Crowl & Kenney 2008) or even enhance star formation by compressing the gas further (e.g., Merluzzi et al. 2013;Kenney et al. 2014;Vulcani et al. 2018). There have been extensive observational studies showing that both scenarios are possible. Therefore more complete understandings of the impact of ram pressure and its consequences on star formation and galaxy evolution require the studies of the multiphase ISM responding to the ICM ram pressure.
A large number of theoretical studies of RPS have been conducted mainly using numerical simulations in two contexts. One is self-consistent cosmological simulations within which galaxies experience RPS as they are moving in clusters (e.g., Steinhauser et al. 2016;Ruggiero & Lima Neto 2017;Jung et al. 2018;Yun et al. 2019), and the other is more controlled simulations of a single galaxy interacting with an inflowing ICM (so called wind-tunnel simulations, e.g., Vollmer et al. 2006;Kronberger et al. 2008;Jáchym et al. 2009;Tonnesen & Bryan 2009Steinhauser et al. 2012;Tonnesen & Stone 2014;Ruszkowski et al. 2014;Ramos-Martínez et al. 2018;Lee et al. 2020;Tonnesen & Bryan 2021). Such simulations reproduce long tails seen in observations and show overall agreement with the prediction for the effectiveness of RPS when the ICM ram pressure is stronger than the ISM anchoring pressure (e.g., Roediger & Brüggen 2007;Tonnesen & Bryan 2009). In both cases, the outer dimensions covered had to be larger than a few tens of kpc, banning the use of pc-scale, high-resolution required for explicit modeling of the ISM physics. Instead, subgrid models of multiphase ISM, star formation, and feedback (e.g., Springel & Hernquist 2003) are often adopted.
There exists a few RPS galaxy simulations that include gas cooling down to ∼ 100 K, but most of them has not particularly focused on modeling of the full multiphase (cold, warm, and hot) ISM properly (e.g., Tonnesen & Bryan 2010Bekki 2014). The radiative heating by photoelectric effect of FUV on small grains is ignored, which is the major heating source of the warm and cold ISM (Wolfire et al. 1995(Wolfire et al. , 2003. Simulations with not enough resolution (Tonnesen & Bryan 2012;Bekki 2014) cannot resolve the Sedov-Taylor stage of SNe that are critical in driving turbulence and creating the hot gas (Kim & Ostriker 2015;Kim et al. 2017a;Steinwandel et al. 2020). Such simulations tend to overcooling the gas and confine the ISM in very thin, unresolved disks. The mass and volume distributions of the multiphase ISM that the ICM is interacting are severely compromised. To our best knowledge, only Lee et al. (2020) have marginally high resolution and set of physics to treat the full range of the multiphase ISM and star formation and feedback explicitly.
Nevertheless, Tonnesen & Bryan (2021) conducted the RPS simulation with a single galaxy with a full range of cooling function and claimed that RPS occurs via mixing between the ICM and ISM. This interesting result qualitatively agrees with the recent in-depth studies of radiative mixing layers (Fielding et al. 2020;Tan et al. 2021) in context of shock/wind-cloud interaction simulations (Gronke & Oh 2018Gronke et al. 2022;Li et al. 2020;Kanjilal et al. 2021;Sparre et al. 2020;Abruzzo et al. 2022) and starburst-driven galactic winds (Schneider et al. 2020) that emphasize the mixing-driven momentum transfer as the major acceleration mechanism for cooler, denser gas.
To enable global galaxy simulations with a large dynamic range, the usual practice is to adaptively refine the resolution elements to achieve constant mass resolution (Lee et al. 2020;Tonnesen & Bryan 2021). Given typical ∼ 2 decades temperature contrast between cold, warm, and hot phases in pressure equilibrium, the spatial resolution of adjacent thermal phases differs by a factor of 5. The interaction between hot and cooler phases and the mixing layers produced by such interaction can be severely altered by large differences in spatial resolution of interacting phases. Simulations with uniformly high resolution are thus necessary to model different phases and their interactions more robustly (e.g., Kim & Ostriker 2018;Schneider et al. 2020). Since the mixing-driven momentum transfer is key physical process in general multiphase hydrodynamical interactions, the multiphase RPS deserves more careful studies using self-consistent, multiphase ISM models that interacts with the ICM.
To this end, we conduct a new suite of numerical simulations focusing on a smaller section of galactic disks with an inflowing ICM. Our numerical models build on the TIGRESS framework developed to model the starforming ISM self-consistently (Kim & Ostriker 2017). TIGRESS solves ideal magnetohydrodynamics (MHD) in a shearing-box with Athena (Stone et al. 2008) and includes additional ISM physics including optically-thin cooling at full temperature range, self-gravity of gas and newly formed stars, star cluster formation in gravitationally bound objects using sink particles, and massive star feedback in the forms of supernovae (SNe) and farultraviolet (FUV) radiative heating. The original closed box model has been used to study the internal regulation of star formation rates (SFRs) and driving of multiphase outflows (Kim et al. 2020a;Kim & Ostriker 2018 among other implications. In this paper, we choose a particular ISM model representing the solar neighborhood condition. We take a snapshot of self-consistently modeled multiphase ISM in a quasi-steady state and conduct controlled numerical experiments with different ICM ram pressures covering relatively weak and strong pressure regimes compared to the ISM anchoring pressure by stellar disks. Our chosen parameters approximately represent the conditions of NGC 4522, a prototypical RPS galaxy in Virgo, at different radii .
As the first of a kind study using local models, this paper focuses on fostering an in-depth understanding of the inner workings of multiphase RPS. In addition, with the help of self-consistent star formation and feedback models of TIGRESS, we investigate how RPS affects star formation in and out of the galactic disks. In the future, we further study the role of magnetic fields in RPS, especially on the dense, molecular gas.
The structure of this paper is as follows. In Section 2, we summarize the TIGRESS framework and introduce the ISM and ICM models. In Section 3, we first overview our simulations using the time evolution of horizontallyaveraged and globally-integrated quantities. Section 3.3 then delineate a variety of physical properties in two representative models. In Section 4.1, we analyze the mass, momentum, and energy transfers between thermal phases. We then check the prediction of the mixingdriven momentum transfer in Section 4.2. Section 5 presents the impact of RPS on SFRs and the extraplanar star formation. Section 6 discusses the main observational imprints from RPS by the mixing-driven momentum transfers. Also, we discuss our results in context and caveats. Finally, the main conclusions are summarized in Section 7.

NUMERICAL METHODS AND MODELS
In this section, we begin by summarizing the numerical methods employed in the TIGRESS framework to simulate the multiphase ISM with star formation and feedback in Section 2.1 for completeness. We then explain the evolution of the ISM without ICM inflows in Section 2.2. The readers who are familiar with the local ISM simulations and only interested in the ICM-ISM interaction can skip the first two subsections. Section 2.3 explains the ICM inflow setup. The tracer fields and gas phases used throughout the paper are defined in Section 2.4.

TIGRESS framework
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 ∼kpc 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 ( , while z is the vertical coordinate. The simulation domain is corotating with galactic rotation speed at the center of the simulation domain, Ω 0 ≡ Ω(R 0 ), arising inertial forces including the Coriolis force and the tidal potential in the momentum equation. The flat rotation curve is assumed, d ln Ω/d 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 FFT 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 selfgravitating collapse as indicated by the Larson-Penston density threshold at ρ LP ≡ (8.86/π)c 2 s /G∆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 −3 for 8 pc resolution and ∼ 300 cm −3 for 4 pc resolution. Note that the typical mass of sink particles is in a range between a few 10 3 M to 10 5 M , representing star clusters rather than individual stars. We treat each particle as a population of stars with a fully-sampled initial mass function of Kroupa (2001).
We use the STARBURST99 stellar population synthesis model to obtain 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). If M SNR /M sf < 1, where M sf = 1540 M (n amb / 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 51 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 SNR p SNR = 2.8 × 10 5 M km s −1 (n amb / cm −3 ) −0.17 as calibrated in Kim & Ostriker (2015). The total and metal mass of SN ejecta, M ej = 10 M and Z SN M ej = 2 M with Z SN = 10Z , 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. (2020a). 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 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 in prep.), 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.

ISM disk model
In this work, we make use of the solar neighborhood model of the TIGRESS simulation suite, setting the parameters for gravitational potential of stars and dark matter (see below; Equation 4). We adopt the angular velocity of galactic rotational Ω 0 = 28 km s −1 kpc −1 , giving rise to the orbit time t orb = 2π/Ω 0 = 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 2017Kim et al. 2020a). 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, Kim & Ostriker (2018); Vijayan et al. (2020) for galactic fountains and winds, and 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 −2 . We introduce initial velocity perturbation and set thermal pressure to ensure that the disk is in Figure 1. Space-time diagrams of horizontally averaged (a) hydrogen number density nH , (b) outgoing mass flux ρvzsgn(z), (c) turbulent pressure ρv 2 z , (d) thermal pressure P , and (e) magnetic pressure PB ≡ B 2 /(8π) for the noICM model. We only show the evolution during a self-regulated state over t ∼ 250 − 500 Myr as a reference that can be directly compared with the models with the ICM. The horizontal dotted line demarks the midplane (z = 0). 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 time scales of 40 − 50 Myr (see Kim et al. 2020a). 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 horizontallyaveraged physical quantities in the space (vertical coordinate z) and time (t) plane as defined by 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) ρv 2 z , (d) thermal pressure P , and (e) magnetic pressure P B ≡ B 2 /(8π). Here, the outward vertical velocity is defined by v out ≡ v z 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.
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 5−6 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 4 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 ten micro Gauss 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. (2020a).

ICM models
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 −2 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.1Z (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 −1 . The total pressure of the ICM at injection is which is dominated by the ram pressure for our chosen v ICM ≥ 10 3 km s −1 (or Mach number of the ICM M ICM > 3.3). 1 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 1 Note that the adopted ICM sound speed (or T ICM ∼ 4 × 10 6 K) is about a factor of two smaller than that of the ICM in the Virgo cluster (T ICM ∼ 2 × 10 7 K; Shibata et al. 2001), which is still smaller than the inflow velocity v ICM so that the results are expected to be qualitatively unchanged.
the external gravity 2 is where the functional form of the external gravitational potential is We adopt the parameters representing solar neighborhood conditions: galactocentric distance R 0 = 8 kpc, stellar surface density Σ * = 42 M pc −2 , stellar scale height z * = 245 pc, and midplane dark matter density ρ dm = 6.4 × 10 −3 M pc −3 (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 W ext ≈ W GG , where This "restoring" force per area (often called as "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;Köppen et al. 2018;Jaffé 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, W GG slightly underestimates the maximum W 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, 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). We consider four different ICM conditions (with additional two higher resolution runs), covering P ICM /k B ∼ 1 − 14 × 10 4 cm −3 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) have P ICM /W GG < 1, and ICM-P7(h) and ICM-P14 have P ICM /W GG > 1. Throughout the paper, the former two models are referred to the weak ICM models, and the latter two models are referred to the strong ICM models, respectively. Our parameter choice brackets the relative strength of the ISM-ICM interaction seen in NGC 4522, a prototypical galaxy undergoing ram pressure stripping in the Virgo cluster 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 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.

Tracer Fields and Gas Phases
In the TIGRESS framework, the gas is divided into five thermal phases based on its temperature, corre- sponding typical discriminators of the three-phase ISM (but including thermally unstable phases; McKee & Ostriker 1977). Each cell is exclusively assigned as cold, unstable, warm, intermediate (warm-hot ionized medium), and 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 . 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 SN , and ICM mass fraction f 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, ρf SN , and ρf ICM are solved with the velocity field of gas. For each SN event, we add the total and metal density of SN ejecta, ρ ej ≡ M ej /V SNR and Z SN ρ 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 ICM = 1 is then added and followed by another passive scalar. Also, the ICM metallicity is set to Z ICM = 0.1Z 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: • genuine ISM -initial gas from the beginning of the noICM model (Z ISM,0 = Z ), • SN ejecta -gas added in the feedback region by SNe (Z SN = 10Z ), and • ICM -gas added from the bottom boundaries as the ICM inflow (Z ICM = 0.1Z ).
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 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 v cool 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).

OVERALL EVOLUTION
In this section, we provide an overview of the RPS process in our simulations and visual impressions using a variety of quantities at different times in two representative models.

Overview of RPS in Simulations
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 ICM of the hot gas, f hot ICM ≡ ρf ICM hot / ρ 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 3rd and 5th rows. We plot the reference line of f hot ICM = 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.
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 pc and z = 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 lowdensity 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 ICM < 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 a quite 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 with 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 gas disk. This is in part because we only model the face-on, plane-parallel interaction, ideal for the simple criteria to work best.

Time Evolution of Masses
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. 2020a). The mean SFR and mass outflow rate over the simulation duration of 250 − 500 Myr are Σ SFR = 3.1 × 10 −3 M kpc −2 yr −1 andΣ gas,out = 7.7 × 10 −4 M kpc −2 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 the 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 Σ gas (t) = Σ gas,max and Σ gas,max /2 is ∼ 130 Myr and 60 Myr for the ICM-P7 and ICM-P14 models, respectively. At around similar time scale, 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.

Morphological 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. Also, significant penetration of the ICM (see (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 ρv 2 z (panel (g)) show a sharp change between the left and right sides. The enhanced ICM mass fraction f ICM (panel (f)) and the relatively large outward velocity (panel (d)) in 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 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 as bright orange color 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 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)).
At this time, the ICM-P7h model already loses about ∼ 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. 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 in to transfer momentum and maintain 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.

STRIPPING OF THE MULTIPHASE ISM
The acceleration and stripping of the ISM is a generic feature of the disk interacting with the ICM. To provide a more quantitative view, we first investigate when and where the hot ICM exchanges its mass, momentum, and energy with the cool ISM. We then seek evidence of mixing-driven momentum transfer.

Mass, Momentum, and Energy Transfers between Thermal Phases
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 ∈ (z min , z max )) to obtain a set of conservation equationṡ where q = M , p, and E for mass, vertical momentum, and total energy, which are respectively defined as M ≡ zmax zmin ρdV, and p ≡ zmax zmin ρv z dV, and The horizontally-averaged fluxes (the square bracket term in the left lend side of Equation 9) are respectively defined by for mass, vertical momentum, and total energy as and . Note that we do not include the gravitational potential term in the energy flux such that the work done by gravity is appeared as a sink term in the right hand side of Equation 9, similarly to the momentum sink term by weight.
In the mass conservation equation, we have mass sink due to star formation (Ṁ * ) and mass source due to SN ejecta (Ṁ SN =Ṅ SN M ej ) if stars are born and SNe are exploded in the particular volume of interest. For our chosen stellar population synthesis model (STAR-BURST99 with the Kroupa IMF; Leitherer et al. 1999), 1 SN ejects 10 M for ∼ 100 M of new star formed, implyingṀ SN ∼ 0.1Ṁ * on average. We note that we did not subtract the mass of stars exploding as SNe in the simulations. The gravitational weight acts as sink of vertical momentum,ṗ sink = WA and where Φ = Φ ext + Φ sg + Φ tidal includes the external potential (Equation 4), the self-gravity potential returned from the solution of the Poisson's equation for both gas and star cluster particles, and the tidal potential arising from the local rotating frame Φ tidal = −qΩ 2 x 2 . Finally, there are energy sources from SN energy injection and shearing-box stressĖ source =Ṅ SN E SN + w 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Ė cool = (n 2 H Λ(T ) − n H Γ)dV and the work done by gravitẏ 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 the cool phase. We include SN mass and energy source to the hot phase; only a small fraction ( 10%) of unresolved SNe injects mass and energy into the form of the cooler phases. By integrating Equation 9 over a time interval ∆t, we obtain phq ph net = 0, wherė 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 By fully accounting for the temporal changes of each quantity, fluxes through vertical faces, and source/sink in each phase,q ph net represents any loss/gain of mass, momentum, and energy through phase transitions within space-time 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 17) 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. 3 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 non-conservation errors.

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 space-time 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 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 space-time 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 space-time bins.
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 creates 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 ∼ kpc. It can be both 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. 2017b;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 with 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 higherz 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-tohot 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 Figure 6. Net mass change rates per unit area for each phase due to phase transition. The space-time bin is separated by ∆t = 10 Myr and ∆z = 400 pc. Mass sink (source) by star formation (SN ejecta) is calculated in each space-time bin and added (subtracted) to the cool (hot) phase to isolate the gain and loss solely due to phase transition. The dashed lines in each panel show the ICM-ISM interface as defined in Figure 2. The phase transition driven by the ICM shows a general trend (clearer with stronger ICM pressure) described by the gain (loss) of the hot (cool) phase near the ICM-ISM interface followed by the loss (gain) of the cool (hot) phase in the extraplanar region.
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 2018Sparre et al. 2020;Kanjilal et al. 2021;Li et al. 2020;Abruzzo et al. 2022).
As RPS is much more efficient in ICM-P14, the hotto-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. 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 as well as the loss of the downward momentum flux. The noICM model shows the change of sign in each phase 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).

Momentum Transfer
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 further understanding of vertical momentum transfer in detail, Figure 8 shows a decomposition of each term in Equation 18 in 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) due 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.  There is no explicit source term, while we show (b) the time dependent term and (f) sink term due to gravity. The flux term is further decomposed into the (c) kinetic, ρv 2 z , (d) thermal, P , and (e) magnetic, Pmag − B 2 z /4π, flux terms. The momentum transferred from hot to cool phase is mostly used to provide appropriate support against the increased weight. At later times (active stripping stage), the cool phase gains kinetic flux.
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 (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, where 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.

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. But, 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. 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 (e) is much larger than the residual energy flux, which is then shared by the actual kinetic flux gain in (c) and work done against gravity in (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 (panel (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).

Mixing Driven Acceleration and Stripping
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 still in play.
As we show in Section 4.1.1, there is a significant phase transition 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 8; 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 cool z plane (first and third columns) and the ICM mass fraction f cool ICM and vertical velocity v cool z plane (second and fourth columns). The solid lines in the second and fourth columns are corre-   Generally, the outflow velocity is above the simple linear prediction, indicating 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 react to the ICM as a whole. At late times, cool clouds are more fragmented so that the crossection to the ICM inflows gets smaller while the surface area of the mixing layers gets larger.
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 cool ICM is generally enhanced. The correlation between f cool ICM and v cool z becomes less clear in the ICM-P3h model due 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 due 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, the 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, we in Figure 11 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,SN ). 4 Throughout the simulation, the contribution 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 is a 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 % of the injected fluxes (e.g., Kim et al. 2020a) 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).
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.
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 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 is not a major driver of outflows, except very early 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 cool ICM shown in Figure 10 shows f cool ICM < 0.1 − 0.2. Having demonstrated that the mixing is the main mechanism to drive outflows (or stripping), the limited range of the f cool ICM 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.

IMPACT OF THE ICM ON STAR FORMATION
We take advantage of self-consistent modeling of star formation and feedback implemented in the TIGRESS framework to study the impact of the ICM ram pressure on star formation in and out of the ISM disks. We first present the changes in overall SFRs and their links to dense gas in the simulations in the presence of the ICM inflows. We then take a detailed look at extraplanar star formation.

Enhancement and quenching of star formation
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: 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% to 50% in the weak ICM models for more than 200 Myr explored in this paper. Also, the temporal modulation of Σ SFR in these models gets stronger with higher peaks. 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. Figure 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 −3 . 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 free-fall time of gas at n H = 10 cm −3 (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 due 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).

Extraplanar Star Formation
One of the intriguing properties commonly found among the RPS galaxies is star-forming patches outside ) and (c) dense gas surface density Σ dense ≡ Σgas(nH > 10 cm −3 ). The colored solid lines correspond to the models with the different ICM pressure, while the black dashed line is the noICM model. Right: box and whisker plots of (b) ΣSF R and (d) Σ dense for early (t < 330 Myr) and late (t > 330 Myr) periods. Boxes include 25th to 75th percentile with median (white dashed horizontal line) and mean (red circle). Whiskers represent 5th to 95th percentile with outliers shown as diamonds.
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.
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 highmetallicity 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 is 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 results 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 gener-ally true for the extraplanar star formation within a few kpc 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. 2017Jáchym et al. , 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.

Imprints in Metallicity of Stripped Tails
The main observational imprint of the mixing-driven acceleration model is on the anti-correlation 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. 5 Assuming f SN = 0, Equation 7 and Equation 8 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 −1 , this implies the metallicity reduction in the mixed cool gas is 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 −1 . 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 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 −1 even in ICM-P14 with v ICM = 2000 km s −1 . 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 10s to 100s kpc scales. Although the stripped tails traveled very far from the disk can have much larger f ICM of clouds up to 0.8, the range of f ICM at a particular position is limited to ∼ 0.1 − 0.2, translating into the velocity difference < ∼ 200 − 400 km s −1 . 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. 2020b), 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 condition favors strong ram pressure stripping 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 ram pressure stripping.
RPS galaxies often show star formation activity outside the main, old stellar disk (e.g., Kenney & Koop-mann 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 dexes) than those formed in the disk (see Figure 14). Observationally, the relative difference of stellar metallicities between young stars from 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.

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. 2019Ramatsoku et al. , 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. 2019Ramatsoku et al. , 2020. In fact, large scale simulations do show late time cooling at more than tens of kpc away from the disk (Tonnesen & Bryan 2012Lee 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. 2017Jáchym et al. , 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;Moretti et al. 2018b;Jáchym et al. 2017Jáchym et al. , 2019. The origin of extraplanar molecular clouds 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 13 CO in the extraplanar molecular clumps at a few kpc above the stellar disk near truncation radii. Given the relatively short formation time of molecular gas (∼ 10 Myr) compared to the stripping time scale (∼ 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 −3 ; 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 frac-tion 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 ICM ∼ 0.15. Our simulations suggest that the extraplanar molecular gas (not those in the stripped tails farther than a few tens kpc from the disk) is mostly originated by 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 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).
Here, we calculate the X-ray surface brightness for the ICM models and compare them with the TIGRESS Figure 15. Soft X-ray surface brightness, S X,0.3−2 keV , as a function of SFR surface density in the past 40 Myr, ΣSFR,40. The ICM models presented in this paper are shown as colored points, while the TIGRESS suite results (C.-G. Kim et al. in prep) are shown as gray points with their mean values as black circle for references. The two diagonal dashed lines denote the SN to soft X-ray efficiencies of SN→X = 1 and 0.1%. The horizontal dotted lines are for the ICM to soft X-ray conversion efficiency ICM→X of 0.05% for given ICM total energy flux.
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 (http://www.atomdb.org; 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 with the area to get the mean X-ray surface brightness. Figure 15 plots as colored points the soft X-ray surface brightness S X,0.3−2 keV as a function of Σ SFR,40 (Equation 20 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 the noICM model in 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 S X,0.3−2 keV /(Ė ICM /A) ∼ 0.05% (horizontal dotted lines), where the ICM total energy flux isĖ ICM /A = 0.5ρ ICM v ICM (v 2 ICM +5c 2 s,ICM ). 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).

Ram pressure stripping and shock/wind-cloud interactions
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 time scale is always shorter than the cloud crushing time scale . 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;Mc-Court et al. 2015;Sparre et al. 2020). Recently, it is 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 2018Kanjilal et al. 2021;Li et al. 2020;Sparre et al. 2019Sparre et al. , 2020. It is certainly true that in our simulations, the chunk of the ISM that is facing the ICM inflows is large (a few 100 pc). The critical size above which the cool clouds can grow by cooling of the hot gas proposed by Gronke & Oh (2018) is where T cl,4 ≡ T cl /10 4 K is the cloud temperature, P 3 ≡ P/(10 3 k B cm −3 K) is the ambient pressure, M wind is the hot wind Mach number, Λ mix,−21.4 = Λ(T mix )/(10 −21.4 erg s −1 cm 3 ) is the cooling coefficient at the temperature of the mixed gas, and χ is the density contrast between wind and clouds. The critical size 7 is of order of a few to tens of parsec at the typical conditions of our simulations with P/k B ∼ 10 4−5 cm −3 K, χ = 10 2−3 , 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 are 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 Figure 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 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.

Star formation in RPS galaxies
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 time scales 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 4 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 star formation rates 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) presents simulations of an RPS galaxy with varying ICM inflow strengths and directions. Combined with adopted higher resolution (adaptive mesh refinement down to 20 pc) and explicit ISM cooling and heating treatments, star formation in this work occurs in the cold, dense gas at number density above 100 cm −3 , representing self-gravitating clouds, as in our simulations. In their moderate ICM inflow model, star formation at the outer region becomes suppressed during 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 measured the age of the youngest stellar population at the H I truncation radii of RPS galaxies to estimate a quenching time scalehow long ago star formation has been quenched since the H I gas stripping. Crowl & Kenney (2008) derived the quenching time scale of a few hundred Myr for the Virgo RPS galaxies which 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.

Caveats and Future Perspectives
In this work, we had to limit our simulation domain to a kpc-size box to achieve high-resolution with explicit treatments of ISM physics (Kim & Ostriker 2017. We thus cannot cover an entire galaxy nor 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 strips 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 solar neighborhood. It is generally expected that the ratio of the ICM ram pressure to the ISM anchoring pressure P ICM /W 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 parameter that determines the properties of mixing (Tonnesen & Bryan 2021;. 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. in prep). 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 TI-GRESS 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.

CONCLUSIONS
We conduct high-resolution MHD simulations of the ICM-ISM interaction to understand how the ICM strips the multiphase ISM from the galactic disk and how star formation changes in and out of the disk. We model the star-forming ISM using the TIGRESS numerical framework. We solve the ideal MHD equations in a local shearing-box with gas and (fixed) stellar gravity, optically thin cooling, star formation, and massive star feedback in the form of SNe and FUV radiative heating (Kim & Ostriker 2017). We take a snapshot of a fully developed ISM from the solar-neighborhood model and simulate it with hot ICM inflows from the bottom boundaries to model face-on interactions of the disk ISM moving in a cluster. We adopt four different strengths of the ICM ram pressure P ICM = ρ ICM v 2 ICM while the ISM condition is fixed. The relative strength of the ICM ram pressure to the ISM anchoring pressure W GG = 2πGΣ gas Σ * covers a range of conditions representing the inner and outer radii of the truncation radius of a galaxy experiencing ram pressure stripping.
Our main findings are as follows: 1. We find that the simple RPS condition comparing P ICM and W GG (Gunn & Gott 1972) works well to predict overall stripping of the ISM disk even in our simulations with the multiphase ISM. Although the porous multiphase ISM structure allows the ICM to penetrate the disk through lowdensity channels and pollute the upper region of the disk regardless of the ICM strength, the effect of the ICM in accelerating the bulk ISM remains insignificant in the weak ICM models with P ICM /W GG < 1. In this case, the majority of the ICM stays below the disk midplane (Figure 2), and the gas remains within the ISM disk over the entire simulation duration (∼ 250 Myr). However, the ICM-ISM interface marches toward the other side of the ISM disk in the strong ICM models with P ICM /W GG > 1. The ICM quickly strips the ISM in a half-mass stripping time scale of 60-130 Myr (Figure 3).
2. In the strong ICM models, the mixing-driven momentum transfer from the ICM to the ISM plays an essential role in RPS (Section 4). At the ICM-ISM interface, the hot ICM inflow first shreds the cool ISM, adding mass into the hotter phases while all phases gain kinetic energy by ram pressure. In the stripped tails, the hot and intermediate phases (genuine ICM and shredded ISM) mix into the cool gas continuously. Most of the hot gas energy is radiated away, while mass and momentum are transferred to the cool phase. These hydrodynamical interactions between hot ICM (energy reservoir) and cool ISM (mass reservoir) result in accelerated cool gas after significant mixing of the hot ICM.
3. The same momentum transfer process also occurs in the weak ICM models. But, the amount transferred to the ISM, together with the SN injected momentum fluxes, is simply used to support the deformed, one-sided disk with increased weight. There is not enough excess momentum and energy to drive strong, continuous outflows (or RPS) as in the strong ICM models.
4. RPS via the mixing-driven momentum transfer imprints on the metallicity of the stripped tails. We find that star clusters formed in the stripped gas (z > 1 kpc from the midplane) show the metallicity lower than the new stars in the disk by ∼ 0.1 dex ( Figure 14). Furthermore, we find a linear relationship between velocity and ICM mass fraction in the stripped cool gas as expected in the mixingdriven momentum transfer, giving rise to an equivalent anti-correlation between velocity and metallicity.
5. Star formation is enhanced (30-50%) in all ICM models at the early epoch of the simulation compared to the noICM model. This enhancement persists in the weak ICM models for the entire simulation time (∼ 250 Myr), while the SFR is greatly reduced after ∼ 100 Myr in the strong ICM models.
As the first results from novel RPS simulations using the local TIGRESS framework, we focus on general responses of the ISM disk with varying ICM inflow strengths. In a forthcoming paper, we will delve deep into the role of magnetic fields in the marginally strong model and the differential stripping of the cold and warm ISM.