© The Authors 2024

1. Introduction

Galaxy clusters represent ideal laboratories for studying galaxy evolution. Within a cluster environment, galaxies can undergo different interactions, namely with the gravitational potential of the cluster, with other cluster galaxies, or with the intracluster medium (ICM). Whereas gravitational interactions affect the stellar and gaseous content of a galaxy, interactions with the hot cluster atmosphere (ram-pressure stripping) only affect the gas. If a galaxy is on a rather radial orbit within the cluster, its velocity increases when approaching the cluster center (Boselli & Gavazzi 2006). At the same time, the ambient gas density of the ICM increases toward the cluster center. Thus, ram pressure, which is calculated as ICM density multiplied by the square of the galaxy velocity, can increase dramatically when the galaxy reaches its shortest distance from the cluster center.

Spiral galaxies that have undergone or are undergoing ram-pressure stripping show a truncated gas disk together with a symmetric stellar disk (e.g., VLA Imaging of Virgo in Atomic gas, VIVA, Chung et al. 2009). If the interaction is ongoing, a gas tail mainly detected in H I is present (e.g., Chung et al. 2007). The gas truncation radius is set by the galaxy’s closest passage to the cluster center via the criterion introduced by Gunn & Gott (1972):

$$\begin{array}{c}\hfill {\rho }_{\mathrm{ICM}}{v}_{\mathrm{gal}}^2=\pi \mathrm{\Sigma }{v}_{\mathrm{rot}}^2/R,\end{array}$$(1)

where ρ_ICM is the ICM density, v_gal the galaxy velocity, Σ the surface density of the interstellar medium (ISM), v_rot the rotation velocity of the galaxy, and R the stripping radius. If peak ram pressure occurred more than 400 − 500 Myr ago, the gas tails have disappeared and the truncated gas disk has become symmetric again (Vollmer 2009).

Under extreme conditions, ram-pressure-stripped galaxies often show extraplanar, one-sided optical and UV emission as well as important tails of ionized gas (e.g., Yagi et al. 2010; Poggianti et al. 2017). Because of their optical appearance, these objects are called jellyfish galaxies. A significant fraction of the Hα emission in these tails can be due to photoionization by massive stars born in situ in the tails (Poggianti et al. 2019). Extraplanar molecular gas is often found in jellyfish galaxies (Jáchym et al. 2014, 2019; Moretti et al. 2018), from which the ionizing stars are formed. The Virgo cluster spiral dwarf galaxy IC 3418 also shows a one-sided filamentary optical and UV tail (Hester et al. 2010). Kenney et al. (2014) coined the term “fireballs” to refer to the extraplanar H II regions with a tail of young stars detected in UV.

ESO 137−001, which is located in the Norma cluster, is one of the nearest jellyfish galaxies (D = 70 Mpc). Sun et al. (2006) discovered a bright 80 kpc-long X-ray tail to the northwest of the galaxy, pointing away from the cluster center. With a deeper Chandra exposure, Sun et al. (2010) discovered a fainter secondary X-ray tail to the south of the bright tail. The X-ray tail was also detected in Hα emission by Sun et al. (2007), Fumagalli et al. (2014), and Fossati et al. (2016). New deep MUSE observations lead to the detection of a faint Hα tail to the north of the bright X-ray tail (Sun et al. 2022; Luo et al. 2023). Jáchym et al. (2019) found a rich distribution of mostly compact CO regions extending to nearly 60 kpc in length and 25 kpc in width. In total, about 10⁹ M_⊙ of molecular gas was detected with ALMA in the tail, assuming the standard Galactic CO-to-H2 conversion factor. The 80 kpc gas tail of ESO 137−001 thus has a multi-phase nature, where the dense molecular gas is colocated with the warm and hot ionized gas (Fig. 1).

Fig. 1. ESO 137−001. Blue: X-rays (Sun et al. 2010). Red: Hα (Sun et al. 2022). Green: CO (Jáchym et al. 2019). The main features of the system are labeled. It appears that the X-ray stripping front is ahead of the Hα stripping front because of the large spatial smoothing scale of the X-ray image, but in reality the X-ray and Hα stripping fronts are at the same position, as shown in Fig. 2 of Luo et al. (2023). For a more detailed image, see Fig. 4.

ESO 137−001 is located at a projected distance of about 200 kpc from the cluster center. Its line-of-sight (LOS) velocity (4680 ± 71 km s⁻¹; Woudt et al. 2004) is close to the average cluster velocity (4871 ± 54 km s⁻¹; Woudt et al. 2008). This suggests that the galaxy orbit lies nearly within the plane of the sky. This is consistent with the large size of the tail. The Norma cluster is located close to the center of the Great Attractor, at the crossing of a web of filaments of galaxies (called the Norma wall; Woudt et al. 2008). The cluster is strongly elongated along the Norma wall, indicating an ongoing merger.

Boehringer et al. (1996) presented ROSAT data of the Norma cluster. The distribution of the X-ray emission of the ICM is elongated along the southeast–northwest direction. When a spherical X-ray emission distribution was subtracted from the image, an elongated southeastern structure and a northwestern ridge reminiscent of a large-scale shock became apparent. Moreover, the high-surface-brightness X-ray emission of the ICM detected by ROSAT and XMM-Newton is displaced toward the northwest of the central cD galaxy ESO 137−006 (see Fig. 2 of Sun et al. 2010). It thus appears that ESO 137−001 evolves in an asymmetric and highly dynamical ICM.

We present dynamical simulations of the ram-pressure stripping event to investigate the physics of the stripped gas and its ability to from stars. We also investigate whether or not the simulations can be used to predict H I maps and to constrain the orbit of ESO 137−001 within the Norma cluster. This article is structured in the following way: the dynamical model is presented in Sect. 2 with an emphasis on the stripping and mixing of the diffuse gas (Sect. 2.3) and the tested temporal ram-pressure profiles (Sect. 2.4). The already existing CO, Hα, and X-ray observations are briefly described in Sect. 3 followed by the results of our simulations (Sect. 4). The influence of the ICM–ISM mixing rate and the hot-gas-stripping efficiency on the simulation results are shown in Sects. 4.1 and 4.2. The existence of a third gas tail in the simulations is discussed in Sect. 5.1. The multi-phase gas masses, the dense molecular gas and CO emission, and the UV emission and star formation are inspected in Sects. 5.2, 5.3, and 5.4. Predicted H I maps and the Hα–X-ray correlation of the stripped gas are presented in Sects. 5.5 and 5.6. A possible orbit of ESO 137−001 within the Norma cluster is proposed in Sect. 5.8. Finally, we outline our conclusions in Sect. 6.

2. The dynamical model

For a proper treatment of the different gas phases during a ram pressure stripping event three-dimensional high-resolution hydrodynamic simulations are the first choice (e.g., Tonnesen & Bryan 2021; Lee et al. 2022). However, these simulations are complex and very much time-consuming. We think that, as a first attempt, simplified simulations with well-controlled sub-grid physics can help to understand the physics of the ram-pressure-stripped ISM and its ability to form stars.

We used the N-body code described in Vollmer et al. (2001), which consists of two components: (i) a non-collisional component that simulates the stellar bulge, stellar disk, and the dark halo and (ii) a collisional component that simulates the ISM as an ensemble of gas clouds (Sect. 2.1). Each particle of the collisional component or sticky particle corresponds to a gas cloud. A scheme for star formation was implemented, where stars are formed during cloud collisions and then evolve as non-collisional particles (Sect. 2.2). The effect of ram pressure is simulated as an additional force acting on gas clouds, which are not protected from the ram pressure wind by other clouds. Simulations with 19 different ram pressure profiles were calculated (Sect. 2.4).

Since our code is not able to treat diffuse gas of low density and high volume filling factor in a consistent way, we assume that warm gas clouds become diffuse if their densities fall below a critical density and if they are stripped out of the galactic plane, that is, if the stellar density drops below a given limit. When the clouds become diffuse their surface density Σ decreases and the acceleration by ram pressure p_ram increases (p_ram ∝ Σ⁻¹; Sect. 2.3).

Furthermore, the stripped warm gas gradually mixes with the ambient ICM. The mixing rate is given by an analytical model of a radiative turbulent mixing layer (Eq. 6; Fielding et al. 2020). Once the mixed mass is equal to the total cloud mass, the cloud temperature is set to 10⁷ K. At the same time the stripping efficiency is increased by a factor of 10 (or three in a second set of simulations) because of a decreased surface density of the heated clouds with respect to the warm clouds. Only for the calculation of the model X-ray and Hα maps the hot and warm gas mass fractions are taken into account (Sect. 2.5). Finally, the equation of motion of the hot gas clouds was modified by adding the acceleration caused by the gas pressure gradient a = −∇p_ISM/ρ_ISM using a Smoothed-particle hydrodynamics (SPH) formalism.

2.1. Halo, stars, and gas

The non-collisional component consists of 81 920 particles, which simulate the galactic halo, bulge, and disk. The characteristics of the different galactic components are shown in Table 1.

Our model was not initially intended to reproduce ESO 137−001 in detail. Its stellar content is twice as massive as that of ESO 137−001. The resulting stellar surface density profile and the model rotation curve are presented in Fig. 2. For comparison, we extracted a radial surface brightness profile from the HST H band image of ESO 137−001 and derived an exponential scale length of 2 kpc. The profile was scaled such that the disk is maximum at a radius of two times the scale length (gray line in the upper panel of Fig. 2). The model disk has a 1.5 times larger scale length and thus an about two times higher mass than observed. The model rotation velocity is ∼180 km s⁻¹ (lower panel of Fig. 2), and the rotation curve becomes flat at a radius of about 7 kpc. We also show the observed rotation curve (Luo et al. 2023) corrected for asymmetric drift using Eq. (9) of Iorio et al. (2017). For the calculation of the stellar velocity dispersion we used Eq. (B3) of Leroy et al. (2008) where we conservatively set the stellar scale length to 1 kpc. With L_K = 2.4 × 10¹⁰ L_⊙ this is at the lower end of the L_K–R_e, K relation of disk galaxies (Courteau et al. 2007). As noted by Jáchym et al. (2014), the corresponding L_K–v_rot relation yields v_rot = 110 − 120 km s⁻¹ for ESO 137−001. This is consistent with the asymmetric-drift corrected rotation curve. The model stellar rotation curve is about 50% higher than observed.

Fig. 2. Initial model conditions. Upper panel: Model stellar- (black solid line) and gas-mass (black dashed line) surface density profiles. Solid red line: Exponential stellar-mass surface density profile with a scale length of 2 kpc. Dashed red line: Analytical approximation to the gas surface density profile (see text). Lower panel: Model (black solid line) and observed asymmetric-drift-corrected (red line; Luo et al. 2023) rotation curve. Red dashed line: Observed rotation curve multiplied by a factor of 1.5.

We adopted a model where the ISM is simulated as a collisional component, that is, as discrete particles that each possess a mass and a radius and can have inelastic collisions (sticky particles). The advantage of our approach is that ram pressure can be easily included as an additional acceleration on particles that are not protected by other particles (see Vollmer et al. 2001).

The 20 000 particles of the collisional component represent gas cloud complexes that evolve in the gravitational potential of the galaxy. The gas surface density profile is presented in the upper panel of Fig. 2. It can be approximated up to a radius of ∼8 kpc by Σ = (1 + 4.5exp( − r/2.2 kpc))×10 M_⊙ pc⁻². The total assumed gas mass is $M_{\mathrm{gas}}^{\mathrm{tot}}=5.2\times {10}^9$ M_⊙, which corresponds to the total neutral gas mass before stripping. A radius is attributed to each particle, which depends on its mass assuming a constant surface density. The normalization of the mass-size relation was taken from Vollmer et al. (2012a). During each cloud-cloud collision, the overlapping parts of the clouds are calculated. Let b be the impact parameter and r₁ and r₂ the radii of the larger and smaller clouds. If r₁ + r₂ > b > r₁ − r₂, the collision can result in fragmentation (high-speed encounter) or mass exchange. If b < r₁ − r₂, mass exchange or coalescence (low-speed encounter) can occur. The outcome of the collision is simplified following Wiegel (1994). If the maximum number of gas particles or clouds (20 000) is reached, only coalescent or mass-exchanging collisions are allowed. In this way, a cloud mass distribution is naturally produced. The energy loss by partially inelastic cloud-cloud collisions results in an effective gas viscosity in the disk.

As the galaxy moves through the ICM, its clouds are accelerated by ram pressure a = ρ_ICMv_gal²/Σ. In addition, the gas clouds are accelerated by the gradients of the gravitational potential a = −∇ϕ. Within the galaxy’s inertial system, the galaxy’s clouds are exposed to a wind coming from a direction opposite to that of the galaxy’s motion through the ICM. The temporal ram pressure profile has the form of a Lorentzian, which is realistic for galaxies on highly eccentric orbits within the Virgo cluster (Vollmer et al. 2001). The effect of ram pressure on the clouds is simulated by an additional force on the clouds in the wind direction. Only clouds that are not protected from the wind by other clouds are affected. This results in a finite penetration length of the ICM into the ISM. Since the gas cannot develop instabilities, the influence of turbulence on the stripped gas is not included in the model. The mixing of the ICM into the ISM is very crudely approximated by the finite penetration length of the ICM into the ISM; in other words, up to this penetration length, the clouds undergo an additional acceleration caused by ram pressure.

2.2. Star formation

We assume that the star formation rate (SFR) is proportional to the cloud collision rate. At the end of each collision, a collisionless particle is created, which is added to the ensemble of particles. Since the mass of the newly formed stars is small compared to the total stellar mass, the newly created collisionless particles have zero mass (they are test particles) and the positions and velocities of the colliding clouds after the collision. These particles are then evolved passively with the whole system. The information about the time of creation is attached to each newly created star particle. We verified that the SFR of an unperturbed galaxy is constant within 1 Gyr. The simulations do not include stellar feedback. Clouds can lose kinetic energy via partially inelastic collisions. This energy loss does not lead to a significant decrease of the velocity dispersion of the gas clouds within an unperturbed galactic disk during 1 Gyr.

In the following, we link the star formation recipe based on cloud–cloud collisions to the recipes based on the local and global gas densities. Vollmer & Beckert (2003) expanded on an analytical model for a galactic gas disk that considers the warm, cold, and molecular phases of the ISM as a single, turbulent gas. This gas is assumed to be in vertical hydrostatic equilibrium, with the midplane pressure balancing the weight of the gas and the stellar disk. The gas is taken to be clumpy, so that the local density is enhanced relative to the average density of the disk. Using this local density, the free-fall time of an individual gas cloud (i.e., the fastest timescale for star formation) can be determined. The SFR is used to calculate the rate of energy injection by supernovae. This rate is related to the turbulent velocity dispersion and the driving scale of turbulence. These quantities in turn provide estimates of the clumpiness of gas in the disk (i.e., the contrast between the local and average density) and the rate at which viscosity moves matter inward. Vollmer & Leroy (2011) applied the model successfully to a sample of local spiral galaxies. Within the framework of the model, the SFR per unit volume is given by

$$\begin{array}{c}\hfill {\dot{\rho }}_{\ast }={\mathrm{\Phi }}_{\mathrm{V}}\rho t_{\mathrm{ff},\mathrm{cl}}^{-1},\end{array}$$(2)

where Φ_V and $t_{\mathrm{ff},\mathrm{cl}}^{-1}$ are the volume filling factor and free-fall time of self-gravitating gas clouds, respectively, and ρ is the global gas density. The volume filling factor is Φ_V = n_cl 4π/3 r_cl³, where n_cl is the number density and r_cl the radius of the self-gravitating gas clouds. For a self-gravitating cloud, the free-fall time equals the turbulent crossing time t_ff, cl = 2 r_cl/v_turb, cl. With the collision time given by t_coll = (n_cl π r_cl²v_turb)⁻¹, this leads to ${\dot{\rho }}_{\ast }=2/3\rho t_{\mathrm{coll}}^{-1}$. Thus, the SFRs of the two recipes are formally equivalent.

In numerical simulations, the star formation recipe usually involves the gas density ρ and the free-fall time $t_{\mathrm{ff}}=\sqrt{3\pi /(32G\rho )}$: ${\dot{\rho }}_{\ast }\propto \rho t_{\mathrm{ff}}^{-1}\propto {\rho }^{1.5}$. We verified that our star formation recipe based on cloud-cloud collisions leads to the same exponent of the gas density in a simulation of an isolated spiral galaxy. As a consequence, our code reproduces the observed SFR-total gas surface density, SFR-molecular gas surface density, and SFR-stellar surface density relations (Vollmer et al. 2012b).

2.3. Diffuse gas stripping and mixing

Our numerical code is not able to treat diffuse gas of low density and high volume filling factor in a consistent way. For a realistic treatment, 3D hydrodynamical simulations should be adopted. Nevertheless, we can mimic the action of ram pressure on diffuse gas by applying very simple recipes based on the fact that the acceleration by ram pressure is inversely proportional to the gas surface density, which in turn depends on the gas density for a gas cloud of constant mass. For the sake of simplicity, we do not modify the radii of the diffuse gas clouds for the calculation of the cloud-cloud collisions. The diffuse clouds are thus mostly ballistic particles under the influence of a ram-pressure induced acceleration. We divide the gas in our simulations into a dense and diffuse phase according to the local density. The gas and stellar densities are calculated via the 50 nearest neighboring particles.

Following Vollmer et al. (2021) we assumed that the warm (∼10⁴ − 10⁵ K) gas clouds become diffuse (i.e., their sizes and volume filling factor increase and their densities and column densities decrease) if they are stripped out of the galactic plane and if their densities fall below the critical density of $n_{\mathrm{crit}}^{\mathrm{warm}}=5\times {10}^{-3}$ cm⁻³. We further assume that the first condition is fulfilled if the stellar density at the location of the gas particle is lower than ρ_*^crit = 2.5 × 10⁻⁴ M_⊙ pc⁻³. For our model galaxy, this density is reached at a disk height of ∼3.5 − 4 kpc.

Once the gas has a high volume filling factor (i.e., it becomes diffuse), its volume increases and the surface densities decreases. The decrease in the surface density is taken into account by considering a gas cloud of a constant mass: for ρ_* < ρ_*^crit, the cloud size is proportional to n^−1/3, the cloud surface density is proportional to n^2/3, and the acceleration caused by ram pressure a = ρ_ICMv_gal²/Σ is increased by a factor of (0.044 cm⁻³/n_ISM)^2/3. Since it is assumed that the gas becomes diffuse, we set n_ISM to the global gas density. This last normalization and the critical densities were chosen such that they led to acceptable results for several observed galaxies undergoing ram pressure stripping. The critical stellar density is motivated by the fact that a significant change in the ISM properties only occurs once the ISM has entirely left the galactic disk. It turns out that this condition is necessary to avoid excessive gas stripping.

The initial temperature of the ISM is T_ISM = 10⁴ K. The hot (> 10⁶ K) diffuse gas is taken into account in the following way: we assume that once the stripped gas has left the galactic disk, it mixes with the ambient ICM. If (i) the gas density falls below the critical value of $n_{\mathrm{crit}}^{\mathrm{hot}}=5\times {10}^{-4}$ cm⁻³, (ii) the stellar density is below ρ_*^crit, and (iii) the ISM temperature is below 9 × 10⁶ K, ICM-ISM mixing begins to take place. On the other hand, if the stellar density exceeds its critical value (i.e., the gas is located within or close to the galactic disk), then the gas is assumed to cool rapidly and its temperature is set to T_ISM = 10⁴ K.

In Vollmer et al. (2021) we assumed that mixing occurs instantaneously and raises the temperature of the mixed gas clouds to

$$\begin{array}{c}\hfill T=\frac{n_{\mathrm{ICM}}6\times {10}^7\mathrm{K}+\sqrt{n_{\mathrm{ICM}}{n}_{\mathrm{ISM}}}{T}_{\mathrm{ISM}}}{n_{\mathrm{ICM}}+\sqrt{n_{\mathrm{ICM}}{n}_{\mathrm{ISM}}}},\end{array}$$(3)

where n_ICM is the density of the ICM and n_ISM is the density of the mixed ISM, which was continuously calculated for each particle. With n_ICM = 2 × 10⁻³ cm⁻³, n_ISM = 0.1 cm⁻³, and T_ISM = 10⁴ K the temperature of the mixed gas is T = 0.9 × 10⁷ K, consistent with the observed X-ray tail temperature of 0.8 keV (Sun et al. 2010).

In a new approach, we estimate the ICM-ISM mixing analytically based on the recipe of Fielding et al. (2020). These authors considered a radiative turbulent mixing layer in which cold and hot gas in pressure and thermal equilibrium move relative to each other. The Kelvin Helmholtz instability quickly develops turbulence that promotes mixing and populates a rapidly cooling intermediate-temperature phase. In quasi-steady state in the frame of the interface, radiative cooling losses are balanced by the advection of hot gas. Hot gas flows into the cooling layer at a speed v_in. The mass transfer rate from hot to cold gas is

$$\begin{array}{c}\hfill \dot{M}\sim {\rho }_{\mathrm{ICM}}{L}^2{v}_{\mathrm{in}},\end{array}$$(4)

where ρ_ICM is the ICM density and L the characteristic length of the mixing layer. According to Eq. (6a) of Fielding et al. (2020) the ratio between the inflow and the relative velocity between the hot and cold phases is

$$\begin{array}{c}\hfill \frac{v_{\mathrm{in}}}{v_{\mathrm{rel}}}={\left(\frac{{\rho }_{\mathrm{ISM}}}{{\rho }_{\mathrm{ICM}}}\right)}^{\frac{3}{8}}{\left(\frac{L}{v_{\mathrm{rel}}{t}_{\mathrm{cool}}}\right)}^{\frac{1}{4}}{\left(\frac{v_{\mathrm{turb}}}{v_{\mathrm{rel}}}\right)}^{\frac{3}{4}},\end{array}$$(5)

where the cooling time is t_cool = α/ρ_ISM, ρ_ISM is the ISM density, and v_turb is the turbulent velocity. Fielding et al. (2020) estimated $(\frac{v_{\mathrm{turb}}}{v_{\mathrm{rel}}})\sim 0.1-0.2$. We assume $(\frac{v_{\mathrm{turb}}}{v_{\mathrm{rel}}})=0.2$.

Combining Eqs. (4) and (5) yields

$$\begin{array}{c}\hfill \dot{M}={\rho }_{\mathrm{ICM}}^{\frac{5}{8}}{\left(\frac{M_{\mathrm{ISM}}}{4/3\pi }\right)}^{\frac{3}{4}}{\rho }_{\mathrm{ISM}}^{-\frac{1}{8}}{\alpha }^{-\frac{1}{4}}{v}_{\mathrm{rel}}^{\frac{3}{4}}{\left(\frac{v_{\mathrm{turb}}}{v_{\mathrm{rel}}}\right)}^{\frac{3}{4}}.\end{array}$$(6)

We call Ṁ the ICM-ISM mixing rate. With n_ICM = 10⁻³ cm⁻³, $\alpha =1.5\times {10}^7{n}_{\mathrm{ISM}}^{-1}$ yr, a molecular weight of μ_ICM = 0.6, and v_rel = 3000 km s⁻¹ we obtain

$$\begin{array}{c}\hfill \dot{M}={10}^{-3}{\left(\frac{M_{\mathrm{ISM}}}{{10}^5{\mathrm{M}}_{\odot }}\right)}^{\frac{3}{4}}{\left(\frac{n_{\mathrm{ISM}}}{{\mathrm{cm}}^{-3}}\right)}^{-\frac{1}{8}}{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1},\end{array}$$(7)

where M_ISM is the mass of an ISM cloud. The cloud mass M_ISM (in M_⊙) and the ISM density are calculated from the dynamical model at each timestep Δt. The total mixed mass is then updated by $M_{\mathrm{mix}}=M_{\mathrm{mix}}+\dot{M}\times \mathrm{\Delta }t$. Once the mixed mass is equal to the total cloud mass we set the cloud temperature to 10⁷ K which corresponds to the temperature given by Eq. (3). For a cloud mass of M_ISM = 10⁵ M_⊙ and a density of n_ISM = 1 cm⁻³ the ICM-ISM mixing time is about 100 Myr.

For the stripped hot gas clouds (> 10⁶ K), we modified the equation of motion by adding the acceleration caused by the gas pressure gradient a = −∇p_ISM/ρ_ISM using a Smoothed-particle hydrodynamics (SPH) formalism. For the sake of simplicity, an isothermal stripped ISM with a sound speed of c_s = 235 km s⁻¹, which corresponds to a temperature of 2.4 × 10⁶ K, is assumed for this purpose. We thus do not solve an explicit energy equation and the calculation of hydrodynamic effects is approximate. Its main purpose is to keep the hot diffuse ISM from clumping.

For the calculation of the acceleration caused by ram pressure a = ρ_ICMv_gal²/Σ, the acceleration is further increased by a heuristic factor of ten if the ISM temperature exceeds 10⁴ K and the stellar density is below ρ_* = 2.5 × 10⁻⁴ M_⊙ pc⁻³. Pressure equilibrium and a temperature ratio between the hot (∼10⁷ K) and the cool (∼10⁴ K) gas imply a density ratio of n_hot/n_cool = 10⁻³. The ratio of the surface densities $\mathrm{\Sigma }\propto {\rho }^{\frac{2}{3}}$ is then Σ_hot/Σ_cool = 10⁻². The heuristic factor is expected to be smaller than the inverse of this ratio. To investigate the effect of this heuristic factor, we recalculated all simulations with a factor of three instead of ten (Sect. 4.2). It turned out that the initial factor of ten used by Vollmer et al. (2021) gave the best results for ESO 137−001.

We note that without the inclusion of hydrodynamical effects, thin gas filaments cannot be produced by our simple model. It is assumed that the diffuse warm gas is heated and ionized by thermal conduction (see Sect. 2.5). Thus, the diffuse gas can significantly emit in the Hα line if its temperature is below 10⁵ K.

2.4. Parameters of the ram-pressure stripping event

Following Vollmer et al. (2001), we used a Lorentzian profile for the time evolution of ram pressure stripping:

$$\begin{array}{c}\hfill p_{\mathrm{rp}}=p_{\max }\frac{t_{\mathrm{HW}}^2}{{(t-t_{\mathrm{peak}})}^2+t_{\mathrm{HW}}^2},\end{array}$$(8)

where p_max is the maximum ram pressure occurring at the galaxy’s closest passage to the cluster center, t_HW is the width of the profile, and t_peak = 500 − 590 Myr is the time of peak ram pressure. The simulations were calculated from t = 0 to t = 800 Myr. We set p_max = 20 000, 36 000, 50 000, 80 000 cm⁻³(km s⁻¹)² and t_HW = 50, 750, 100 Myr. Following Nehlig et al. (2016) and Vollmer et al. (2018), we investigated the influence of galactic structure (i.e., the position of spiral arms) on the results of ram pressure stripping by varying the time of peak ram pressure between t_peak = 530 Myr and 590 Myr. For a different peak stripping the galactic spiral arms are not in the same place with respect to the leading edge of the interaction. We also used two different Gaussian profiles. The ram pressure profiles are specified in Table 2 and shown in Fig. 3. Each simulation took about four weeks on a single CPU.

Fig. 3. Model ram-pressure profiles (Table 2). The thick red line shows the highest-ranked model, and the vertical red lines show the range of time steps of the highest-ranked model.

The inclination angle of ESO 137−001 is i ∼ 65° (Sun et al. 2007). If the galaxy is moving in the plane of the sky in the opposite direction of its gas tail (close to the minor axis of the galaxy), the angle between the disk plane and the ram pressure wind is also about 65°. To investigate the influence of the angle between the disk plane and the ram pressure wind on the resulting gas distribution and velocity field, we set this angle to (50°, 60°, 75°). For the given parameter set we calculated 33 models.

In addition, we recalculated the model set (i) with a heuristic increase in the acceleration caused by ram pressure by a factor of three instead of ten and with a (ii) three-times-higher and (iii) three-times-lower ICM–ISM mixing rate; that is, a mass-inflow rate within the ICM–ISM mixing layer (Eq. 7). In total, we calculated 99 different models. We find that the models with the inflow rate of Eq. (7) give the best results.

2.5. Extraction of the observables

From the model calculations the following observables were extracted: H I, CO, Hα, NUV, and X-ray emission. Following Vollmer et al. (2008), the separation between the atomic and molecular gas phases is based on the gas density:

$$\begin{array}{c}\hfill \frac{M_{\mathrm{mol}}}{M_{\mathrm{tot}}}=\min (4.2\times {\rho }_{\mathrm{tot}},1),\end{array}$$(9)

where M_mol and M_tot are the molecular and total gas mass of the cloud and ρ_tot is the total gas density in M_⊙ pc⁻³. Elmegreen (1993) found the following relation between the molecular to atomic gas fraction and the gas pressure: M_mol/M_HI ∝ P^1.2. For a constant gas velocity dispersion the gas pressure is proportional to the gas density and M_mol/M_tot is about proportional to the gas density. The atomic gas mass of the cloud then is M_HI = 1 − M_mol/M_tot.

The calculation of the model X-ray emission map is based on the emission measure of the hot stripped ISM. For a temperature of T ∼ 10⁷ K the X-ray emissivity approximately yields

$$\begin{array}{c}\hfill {ϵ}_{\nu }\sim {10}^{-23}{n}_{\mathrm{e}}^2\mathrm{erg}{\mathrm{cm}}^{-3}{\mathrm{s}}^{-1},\end{array}$$(10)

(e.g., Fig. 34.2 of Draine 2011) where n_e is the electron density in cm⁻³. It is assumed that the gas is fully ionized. The X-ray luminosity of a hot gas particle or cloud is given by

$$\begin{array}{c}\hfill L_{\mathrm{X}}=4\pi R^3{ϵ}_{\nu }=5.3\times {10}^{35}\left(\frac{M_{\mathrm{cl}}}{1{\mathrm{M}}_{\odot }}\right)x_{\mathrm{hot}}\left(\frac{{\rho }_{\mathrm{hot}}}{1{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-3}}\right){\mathrm{erg}\mathrm{s}}^{-1},\end{array}$$(11)

where ρ_hot is the average density of the hot gas and x_hot its overdensity. To reproduce the observed X-ray luminosity of the stripped gas tail (L_{0.5 − 2 keV} ∼ 10⁴¹ erg s⁻¹; Sun et al. 2010) by the best-fit model we set x_hot = 70. For the calculation of the hot gas mass the hot portion of a gas cloud is taken into account. Hot gas within the galaxy, which is heated by supernova explosions, is not taken into account. Therefore, the observed strong X-ray emission emitted by the galactic disk is not reproduced by the model.

The model Hα images consist of two components: (i) the H II regions ionized by young massive stars and (ii) diffuse ionized gas that is ionized by the stellar UV radiation, heat conduction (Cowie & McKee 1977), or possibly by strong shocks induced by ram pressure stripping (as, e.g., in the diffuse Hα tail of NGC 4569; Boselli et al. 2016). The first component is modeled by the distribution of stellar particles with ages less than 20 Myr. This is the only component used for the images of the models without diffuse gas stripping. For the diffuse component we assumed that the gas with temperatures lower than 10⁶ K is ionized by thermal conduction: the mixing between the hot ICM and the warm stripped gas leads to a relatively dense (∼10⁻² cm⁻³; Sun et al. 2010) hot medium at a temperature of ∼10⁷ K (Sect. 2.3). Thermal electrons from this mixed ISM-ICM gas penetrate into the neutral warm stripped ISM clouds ionizing and heating them. Ultimately, this leads to the evaporation of the ISM clouds (Cowie & McKee 1977). The typical evaporation timescale for a cloud of 10²¹ cm⁻² is ∼100 Myr (Vollmer et al. 2001). In the case of a magnetic field configuration that inhibits heat flux (e.g., a tangled magnetic field), this evaporation time can increase significantly (Cowie et al. 1981) and might attain several 100 Myr. Stripped clouds of warm neutral gas can thus survive within the stripped gas tail. The influence of magnetic fields on thermal conduction in the turbulent stripped gas is discussed in Sect. 5.7.

For the determination of the emission measure of an evaporating gas cloud we followed Eq. (25) of McKee & Cowie (1977):

$$\begin{array}{c}\hfill \frac{\mathrm{d}(\mathrm{EM})}{\mathrm{d}lnT}={\left(\frac{n_{\mathrm{ICM}}}{n}\right)}^2\frac{\mathrm{d}N}{\mathrm{d}lnT}=3{\sigma }_0^{\frac{7}{8}}{n}_{\mathrm{ICM}}^2{\left(\frac{T}{T_{\mathrm{ICM}}}\right)}^{\frac{3}{2}}R,\end{array}$$(12)

where EM is the emission measure of the conduction front, T and T_ICM the temperature of the warm ISM and hot ICM, N and n are the column density and density of the conduction front, and ${\sigma }_0\simeq (\frac{T_{\mathrm{hot}}}{1.54\times {10}^7\mathrm{K}}{)}^2{n}_{\mathrm{ICM}}^{-1}{R}_{\mathrm{pc}}^{-1}$ is the saturation parameter. To also include the case of a classical evaporating cloud (Eq. 22 of McKee & Cowie 1977), we set σ₀ = 1 if σ₀ < 1. The gas cloud radius is given by R = (3 M_cl/(4 πρ))^1/3. Since the highest column densities occur at the transition between the inner classical and the saturated zone of the conduction front, we only used $\frac{\mathrm{d}N}{\mathrm{d}lnT}$ of the inner classical zone (Eq. 23 of McKee & Cowie 1977).

For all gas particles with temperatures lower than 10⁵ K the Hα luminosity was calculated with

$$\begin{array}{c}\hfill L_{\mathrm{H}\alpha }=\mathrm{EM}{\alpha }_{{\mathrm{H}}_2}^{\mathrm{eff}}\pi R^{2}h{\nu }_{\mathrm{H}\alpha }=2.7\times {10}^{36}{\sigma }_0^{\frac{7}{8}}{n}_{\mathrm{e}}^2{(xn)}^{-1}{\mathrm{erg}\mathrm{s}}^{-1},\end{array}$$(13)

where α_H2^eff = 1.17 × 10⁻¹³ cm³ s⁻¹ is the Hα effective recombination coefficient, h the Planck constant, ν_Hα the frequency of the Hα line, and n and x = 10³ are the average density and the overdensity of the warm stripped ISM. As for the X-ray emission, we set the local electron density of the hot gas to 70 times the mean ICM density. The overdensity of the warm gas was chosen such that the observed Hα luminosity of the gas tail (L_Hα ∼ 3 × 10⁴⁰ erg s⁻¹) is approximately reproduced by our best-fit model. For a justification of the chosen overdensities we refer to Sect. 5.7.

As explained in Sect. 2.2, the information about the time of creation is attached to each newly created star particle. In this way, the Hα emission distribution caused by H II regions can be modeled by the distribution of star particles with ages younger than 20 Myr. The UV emission of a star particle in the two GALEX bands is modeled by the UV flux of single stellar population models produced from the STARBURST99 software (Leitherer et al. 1999). The age of the stellar population equals the time since the creation of the star particle. The total UV distribution is then the extinction-free distribution of the UV emission of the newly created star particles.

The model H I, CO data cubes and the model NUV map were convolved with Gaussian kernels to the spatial resolutions of the actual observations. The X-ray map of Sun et al. (2010) was adaptively smoothed. We decided to convolve our model X-ray map with a Gaussian kernel with a half power width of 2.2 kpc, which lead to results well comparable to the adaptively smoothed X-ray map. The model distribution of the diffuse warm ionized gas was convolved with a Gaussian kernel with a half power width of 1 kpc, that of the H II regions with a kernel with a half power width of 0.6 kpc. The total Hα image was obtained by adding the images of the diffuse warm ionized gas and the H II regions. We then added the H II region component with a heuristic normalization to the model image of the diffuse warm ionized gas, which worked best to insure that extraplanar H II regions are well visible within the emission of the diffuse ionized component. In these images, the surface brightness of the disk emission with respect to the extraplanar emission is realistic in a qualitative, but not quantitative, sense. We used the following projection angles of ESO 137−001: inclination i = 66°, position angle PA = 10°.

3. Observations and model fitting

ESO 137−001 is one of the rare cases where CO, Hα, and X-ray data at decent spatial resolutions (≲1 kpc) are available. The X-ray observations were made by Chandra (Sun et al. 2010), the Hα observations by MUSE (Sun et al. 2022), and the CO observations by ALMA (Jáchym et al. 2019). All these observation are displayed in Fig. 4.

Fig. 4. Multiwavelength observations of ESO 137−001. Upper panel: MUSE Hα is shown in color (Sun et al. 2022; Luo et al. 2023); Chandra X-ray in black contours (Sun et al. 2010); and ALMA CO in white contours (Jáchym et al. 2019). It appears that the X-ray stripping front is ahead of the Hα stripping front because of the large spatial smoothing scale of the X-ray image, but in reality, the X-ray and Hα stripping fronts are at the same position, as shown in Fig. 2 of Luo et al. (2023). Lower panel: MUSE Hα velocity field (Sun et al. 2022). The color scale corresponds to the LOS velocities with respect to the systemic velocity in km s−1. At a distance of 70 Mpc 30″ correspond to about 10 kpc. The size of the images is 80 kpc × 60 kpc.

There are several ways to quantitatively compare our simulations to the observations: global properties as integrated fluxes of the disk and tail regions at a given wavelength or length of the tail, brightness profiles or mean LOS velocity along the axis of the system (along the tail), or a direct comparison of the resolved emission distributions. Global properties do not depend on the morphology of the emission in a given region (e.g., the tail region). Profiles along given axes are averaged over the perpendicular direction and thus depend on the global morphology along the axis. These quantities have the advantage that they mostly depend on the large-scale morphology and not much on particular realizations of a model, which might have different small-scale morphologies. In our case, different initial conditions of the galactic disk in terms of spiral arms and surface density profile will lead to a different morphology of the stripped gas tail (see, e.g., Vollmer et al. 2021). In the case of a limited number of simulations with only one initial condition, the comparison of global properties and profiles along the tail axis are the best choice. Since we made simulation with different ram pressure wind profiles and different initial conditions (Sect. 2.4), a direct comparison of the resolved model and observed emission distributions is appropriate. Because the emission distributions depend on the chosen projection (the inclination and position angles are given by the observations, the azimuthal angle has to be chosen), we produced model images with three different azimuthal angles. In addition, we allowed for a small possible rotation and a shrinking or expansion between the model and observed images. With about a hundred different simulations and three different projections, we are confident that it is possible to constrain the ram pressure profile and the time to peak ram pressure by a direct comparison of the resolved emission distributions.

For the search of the highest-ranked models we calculated the goodness of the fit for all timesteps of all models. This was also done for the velocity field. We define the goodness as the sum of the absolute differences between the model and observed image pixels. Before the calculation of the goodness of the model X-ray, Hα, and CO maps, the given model map was divided by the normalization factor Q described in Eq. (16) of Vollmer et al. (2020). The factor Q minimizes the χ² between the model and observed maps.

We clipped the model CO data at a surface density of 1 M_⊙ pc⁻² in a 10 km s⁻¹ channel. If a detection in three adjacent channels is required this yields a detection limit of 3 M_⊙ pc⁻², a value that is broadly consistent with the lowest contour of 6 M_⊙ pc⁻² in Fig. 4 of Jáchym et al. (2019).

The clipping values of the model X-ray and Hα surface brightness distributions were 10⁻⁶ erg cm⁻² s⁻¹. The model Hα clipping value corresponds to the MUSE 3σ Hα surface brightness of 1.6 × 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻² (Luo et al. 2023). The model X-ray clipping value corresponds to about one fourth of the Chandra 3σ X-ray surface brightness of ∼6 × 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻² for ∼16 kpc scales (Sun et al. 2010).

For the comparison of model and observed the Hα velocity fields we calculated the goodnesses for the clipped and unclipped model Hα data. The visual inspection of the best-fit models showed that the comparison based on the unclipped model Hα data lead to the best result, that is model Hα velocity fields, which reproduce the main characteristics of the observed Hα velocity field. The Hα and CO maps were aligned to the X-ray image (same center and same pixel size) before the calculation of the goodnesses.

Because of the unknown CO-H₂ conversion factor especially in the tail of ESO 137−001 we decided not to use the CO surface brightness maps but to produce binary maps of the model and observed CO maps. A pixel of the observed and model map was set to unity if the flux exceeds three times the rms and to zero otherwise. In this way the fitting procedure favored models whose qualitative morphologies resemble that of the observed CO distribution. In parallel, we also calculated the goodnesses using the CO surface brightness maps for the models with the nominal ICM-ISM mixing rate and stripping efficiency. It appeared that in this case model timesteps smaller than −50 Myr from peak ram pressure are preferred when the goodness is based on the CO, X-ray, Hα emission, and the Hα velocity field. However, the observed bifurcated structure of the tail was not well recovered by these models. Reassuringly, our highest ranked model (Fig. 5) was consistently among the preferred models for both types of goodnesses. The goodness distributions are smooth (Fig. A.1) with a change of slope after the first ∼100 highest-ranked models. For the sake of efficiency, we decided to inspect and present the 50 highest-ranked models. The visual inspection showed that the goodnesses based on the CO binary maps led to the best results in terms of model Hα and X-ray maps as well as the Hα velocity field. In particular, the emission of the model tails was recovered in the model CO images in a satisfactory way.

Fig. 5. Highest-ranked model A of ESO 137−001. Upper panel: Hα is shown in color; X-ray in dark gray contours; CO in white contours; and stellar content in black contours. The relative contours are the same as in Fig. 4, i.e., Hα: logarithmic transfer function with 15 contour levels; X-ray and CO: square root transfer function with 15 contour levels. Lower panel: Hα velocity field. The colors are the same as in Fig. 4. The corresponding time evolution and that of the mass surface densities of the different gas phases can be found be found in online movies.

Since we only calculated a restricted number of simulations, we took into account (i) a small possible rotation between the model and observed images from −15° to 15° in steps of 5° and (ii) a shrinking or expansion by a factor of 0.8 to 1.2 in steps of 0.1. These modifications were applied to the observed images. In addition, we added and subtracted 30° to and from the azimuthal projection angle.

Whereas the highest-ranked model at a single wavelength corresponds to the minimum goodness, the comparison of goodness at different wavelengths is not straight forward. We decided to rank the models at the different wavelengths to calculate the sum of the ranks at the different wavelengths. We define the highest-ranked model as the model with the smallest value of the sum of the ranks. The highest-ranked models were determined for (i) the combination of the different wavelengths, (ii) this combination plus the Hα velocity field, and (iii) only the Hα velocity field. Of course, a valuable model should reproduce the spatial distribution of the gas phases and their velocity fields if available. The 50 highest-ranked models of cases (i) to (iii) (Tables C.1–F.3) were inspected by eye, with only a small portion presented in the form of images in this article.

4. Results

The first 50 models with the smallest total ranks are presented in Tables C.1–C.3 for the comparison based on (i) the CO binary, Hα, X-ray images and the Hα velocity field and (ii) the CO binary image and the Hα velocity field, and (iii) the CO binary, Hα and X-ray images. After the visual inspection of all highest-ranked models, we chose three pre-peak models, for which there is a good resemblance between the model and observed multiwavelength images. The selection process is described in the supplementary material (Data availability).

Gaussian models (models 6 and 7) are excluded by our selection process. When the model selection is only based on the CO binary, Hα and X-ray maps (models with Δt < 0 in Table D.3), models 3 (high ram pressure and intermediate width, galaxy orbit of intermediate eccentricity) and 5 (highest ram pressure and smallest width or highly eccentric galaxy orbit) are preferred. When the model selection is based on the CO binary map and the Hα velocity field (models with Δt < 0 in Table C.2), the models with the lowest peak pressure and the largest width are preferred (model 1; less eccentric galaxy orbit). In addition, three models with the highest peak ram pressure and smallest width are present among the 50 highest-ranked models. When the model selection is based on the CO binary, Hα and X-ray maps and the Hα velocity field (models with Δt < 0 in Table C.1), the models with the highest peak pressure and the smallest width are preferred (model 5). Model 5 with a highly eccentric galaxy orbit is the overall preferred model because it is among the 50 highest-ranked models based on all three selection criteria. Shifting the time of peak ram pressure (models 5a and 5b) leads to a different morphology of the tail while the tail length and width are generally conserved.

It turned out that the observed extraplanar CO emission within 10 kpc of the galactic disk with the upturning extraplanar CO filament (Fig. 1 and upper panel of Fig. 4 of Jáchym et al. 2019) was very hard to reproduce. These features are only present in model 5 at Δt ∼ −20 to −40 Myr, which are amongst the highest-ranked pre-peak models of the comparison based on the CO/Hα/X-ray emission distributions and the Hα velocity field (Table C.1). We call this model the highest-ranked model A (Figs. 5 and 6).

Fig. 6. CO emission distribution (contours) on CO velocity field. Upper left: ALMA observations. The color levels are the same in all panels.

The parameters of the highest-ranked model A are presented in Table C.1. The CO, Hα, and X-ray emission distribution together with the Hα velocity field are presented in Fig. 5.

The peak-ram pressure model 1b is the highest-ranked model only based on the velocity field with Δt ≤ 0 Myr. We call this model the preferred model B (Fig. 7; Table C.2). Moreover, we visually identified model 1a with Δt = −10 Myr as the preferred model C (Fig. 8; Table C.2).

Fig. 7. Preferred model B of ESO 137−001. Upper panel: Hα is shown in color; X-ray in dark gray contours; CO in white contours; stellar content in black contours. Lower panel: Hα velocity field. See Fig. 5 for the description of the contours and colors. The corresponding time evolution and that of the mass surface densities of the different gas phases can be found be found in online movies.

Fig. 8. Preferred model C of ESO 137−001. Upper panel: Hα is shown in color; X-ray in dark gray contours; CO in white contours; and stellar content in black contours. Lower panel: Hα velocity field. See Fig. 5 for the description of the contours and colors. The corresponding time evolution can be found online (observations_modelC.gif). The time evolution of the mass surface densities of the different gas phases can be found online (gasphases_modelC.gif).

In models A, B, and C, the main tail and the southern secondary tail are present. As observed, the main tails of models A and C are significantly brighter than the secondary tails. In model B, the two tails have about the same surface brightness. The extents of the X-ray and Hα tails are reasonably reproduced by the models. Whereas the lengths of the secondary tails are well reproduced by the models, the sizes of model main tails are about 30% larger than the observed sizes. The observed broadening to the north of the main Hα tail at distances > 40 kpc is only present in model A. Overall, the morphology of the observed X-ray emission distribution is better reproduced by models A and C than that of model B. Especially, the observed local X-ray maximum ∼37 kpc from the galaxy center is only reproduced by model C. Furthermore, the Hα velocity field at distances > 30 kpc is better reproduced by models B and C. This is expected because they were selected based on the comparison of the velocity fields. The observed northern faint Hα tail is absent in all models, except for the slight broadening to the north of the main Hα tail at distances > 30 kpc in the highest-ranked model A.

The observed CO distribution is better reproduced by the highest-ranked model A. We observe two filamentary structures in the model CO emission distribution of the main tail: a northern convex and southern concave filament. The latter corresponds to the upbending extraplanar CO filament (Fig. 1). Whereas the southern filament is stronger than northern one in the ALMA observations, the model shows two filaments of similar column densities (see Sect. 5.3).

We conclude that the X-ray and Hα observations are best reproduced by model C, whereas the CO observations are better reproduced by model A. The existence of a two-tail structures is a common feature in our models. It is due to the combined action of ram pressure and rotation together with the projection of the galaxy on the sky. Magnetic fields might enhance the appearance of a tail that is bifurcated in the plane of the sky (Ruszkowski et al. 2014) but they are not included in our model. Too strong magnetic field might suppress thermal conduction, which is needed to explain the observed Hα tail emission (see Sect. 5.7).

4.1. The influence of the ICM–ISM mixing rate

We recalculated all models of Tables C.1–C.3 with a three times higher and lower ICM-ISM mixing rate. The ranking results for the comparisons (i) to (iii) are presented in Tables D.1–E.3.

As before, we focus on near ram-pressure peak models. For a higher mixing rate model 1c is preferred at Δt = −10 to 10 Myr for comparisons (i) and (iii) (Tables D.1 and D.3). Based on the comparison between the CO binary maps and the Hα velocity field models 1, 1a, and 1b at Δt = −50 to −10 Myr are preferred (Table D.2). All preferred models show strong X-ray main tails, but faint or absent outer Hα tails. The Hα extents of these model tails with a higher mixing rate are significantly smaller than those of the models with the nominal mixing rate. The reason for this behavior is the higher stripping efficiency (Sect. 2.3) of mixed hot gas together with the faster mixing. The outer parts of the tails in the model with the nominal mixing rate are already pushed to larger distances out of the field of view in the model with the three times higher mixing rate. In all preferred models the southern gas tail is barely visible in Hα emission.

The 50 highest-ranked models with a three times lower mixing rate are presented in Tables E.1–E.3. For a lower mixing rate the Hα and X-ray morphologies of the model gas tails are very different from the observed morphologies (Figs. E.1–E.3). Model emission Hα and X-ray is mostly seen up to distances of ∼20 kpc from the galaxy center. In addition, rare patchy emission regions with sizes of ∼10 kpc are present at larger distances. The gas distribution in the tail is much smoother with less overdensities than in the models with the nominal or a three times higher mixing rate. Because of the applied sensitivity limits only a small amount of emission is present in the model X-ray and Hα maps of the models with a lower mixing rate. None of the preferred model based on the comparisons (i) to (iii) reproduce the available observations.

We conclude that the highest-ranked models with the nominal ICM-ISM mixing rate reproduce observations significantly better than the models with a three times higher or lower ICM-ISM mixing rate.

4.2. The influence of the hot-gas-stripping efficiency

As stated in Sect. 2.3 the calculation of the acceleration caused by ram pressure is a = ρ_ICMv_gal²/Σ. This acceleration was further increased by a heuristic factor of ten if the ISM temperature exceeds 10⁴ K. To investigate the influence of this enhanced stripping efficiency of the hot gas on the model results, we recalculated all models of Tables C.1–C.3 with a three times lower stripping efficiency of the hot gas. The ranking results for the comparisons (i) to (iii) are presented in Tables F.1–F.3. The highest-ranked models with Δt ≤ 0 Myr are model 2 for comparison (i), models 1, 1a, and 1b for comparison (ii), and models 2 and 3 for comparison (iii).

The models with a lower stripping efficiency typically show stronger outer X-ray tails. This is expected because the gas in the outer tail is less accelerated than in the model with the nominal stripping efficiency and thus stays denser. Only the third and fourth highest-ranked model based on comparison (ii) (Fig. F.2) show two separate tails in X-ray and Hα emission with sizes larger than 30 kpc.

For a direct comparison with the models of the nominal stripping efficiency, we show the preferred model B with a three times lower stripping efficiency in Fig. 9. The gas truncation radius within the disk is larger than that of the preferred model B. Moreover, the X-ray and Hα tails are significantly stronger than those of the preferred model B. The emission of the northern and southern model tails are enhanced by about the same factor. Especially the model Hα emission of the southern tail is much stronger than it is observed. As expected, the velocity field is closer to the observed Hα velocity field than that of the highest-ranked model A.

Fig. 9. ESO 137−001 preferred models. Left panel: Model corresponding to the preferred model B of Fig. 7 with a three times lower stripping efficiency with respect to the nominal stripping efficiency. Right panel: Model close to the preferred model C (Δt = −20 Myr instead of Δt = −10 Myr) of Fig. 8. Upper panel: Hα is shown in color; X-ray in dark gray contours; CO in white contours; and stellar content in black contours. Lower panel: Hα velocity field. See Fig. 5 for the description of the contours and colors.

We conclude that overall the CO/Hα/X-ray highest-ranked models with a three times lower stripping efficiency reproduce the available observations less well than the highest-ranked models with the nominal stripping efficiency.

5. Discussion

After the detailed comparison between the models and observations based on different data (i) CO/X-ray/Hα emission and Hα velocity field, (ii) Hα velocity field, (iii) CO/X-ray/Hα emission we are left with three models, which reproduce the available observations in a satisfactory way: the highest-ranked model A (Fig. 5), the preferred model B (Fig. 7) and the preferred model C (Fig. 8). All these models have an ICM-ISM mixing rate, which is consistent with theoretical expectations (Sect. 2.3).

Model A has an angle between the galactic disk and the ram pressure wind of Θ = 60°, whereas models B to E have Θ = 75°. The LOS components of the model galaxy velocity unit vector of model A is tiny (0.07 − 0.09) and small of models B and C (0.33). Whereas the time to peak ram pressure of model A is Δt = −40 to −20 Myr (before peak ram pressure), it is Δt = 0 for model B and Δt = −30 to −10 Myr for model C. If we require that the galaxy is observed before peak ram pressure as deduced from its location within the Norma cluster and the direction of its gas tail (Sun et al. 2010), we prefer models A and C. If we also allow for peak ram-pressure, model B is also a preferred model.

In the following, we investigate the existence of a third gas tail in the simulation (Sect. 5.1) and inspect the stripped multi-phase model gas masses (Sect. 5.2), the model CO emission distribution and velocity field (Sect. 5.3), the model UV emission and star formation distributions (Sect. 5.4), and the Hα–X-ray correlation of the stripped gas tail (Sect. 5.6). Moreover, we give a possible orbit of ESO 137−001 within the Norma cluster (Sect. 5.8).

5.1. The faint northern gas tail

As discussed in Sect. 4 the faint northern gas tail (Fig. 1) is absent in our preferred models. However, such faint tail structures are present in our simulations but they are rare for the highest-ranked models. We found such structures in model 1b at timesteps 0 Myr ≤ Δt ≤ 40 Myr within the highest-ranked models when only the Hα velocity field is taken into account (Table C.2). The timestep Δt = 20 Myr is presented in Fig. 10. This model corresponds to model B with a 20 Myr later time and a somewhat different projection.

Fig. 10. ESO 137−001 model showing detached elongated Hα emission regions in the north and south of the two main tails.

Two detached elongated Hα emission regions in the north and south of the two main tails are visible at (40,  25) kpc and (25,   − 6) kpc (upper panel of Fig. 10). They are made of gas that resists the ram pressure wind because of its high density. This gas is constantly ablated by the ram pressure wind. The existence of such structures depends on the density structure of the stripped gas, which itself depends on the initial gas distribution and the temporal ram pressure profile. We conclude that our model is in principle able to produce structures as the observed faint northern gas tail, which is mainly detected in Hα emission.

5.2. Multi-phase gas masses

To provide a more quantitative comparison with observations, we extracted the tail gas masses within the regions of the model maps (90 kpc × 60 kpc; Table 3). We did this for the cold neutral medium (CNM), the warm neutral medium (WNM), the warm ionized medium (WIM), and the hot ionized medium (HIM). For the mass of the warm ionized medium of a gas particle we set

$$\begin{array}{c}\hfill M_{\mathrm{WIM}}=3{\sigma }_{0}70{\rho }_{\mathrm{hot}}{V}_{\mathrm{warm}}=3{\sigma }_{0}70{\rho }_{\mathrm{hot}}{M}_{\mathrm{cl}}{({10}^3{\rho }_{\mathrm{warm}})}^{-1},\end{array}$$(14)

where ρ_hot, warm are the average density of the hot and warm stripped gas.

The fractions of stripped (2 kpc ≤ x ≤ 80 kpc) gas in the different phases (CNM, WNM, WIM, and HIM) as a function of time are presented in Fig. 11.

Fig. 11. Fractions of stripped (2 kpc ≤ x ≤ 80 kpc) gas in the different phases. The time t = 0 Myr corresponds to peak ram pressure.

The WNM observed in the H I line represents ∼80% of the stripped gas at the beginning of the all three simulations. The WNM mass fraction then decreases to ∼20% after 200 − 300 Myr. Whereas it increases slightly towards peak ram pressure in model A, it decreases to ∼10% at peak ram pressure in models B and C. The CNM fraction varies between 10% and 20% in all three models. Only in model A it increases to ∼30% at peak ram pressure. The WIM fraction monotonically increases with time in all three simulations. It reaches a maximum of ∼5% in model A and ∼10% in models B and C. The HIM fraction increases to 60% after 200 − 300 Myr. It then decreases to ∼30% at peak ram pressure in model A and stays approximately constant in models B and C. Near peak ram pressure all three gas phases are equally present in the stripped gas in model A, whereas the HIM dominates the mass budget in models B and C.

For the CNM Jáchym et al. (2019) derived a mass of M_CNM ∼ 9 × 10⁸ M_⊙ within the disk and M_CNM ∼ 10⁹ M_⊙ in the tail region assuming a Galactic CO-H₂ conversion factor. Sun et al. (2007) estimated the mass of the WIM. Assuming a Galactic volume filling factor of 0.2 (Boulares & Cox 1990) yields M_WIM ∼ 2 × 10⁸ M_⊙. Sun et al. (2006, 2010) derived a HIM mass of M_HIM ∼ 10⁹ M_⊙ in the tail of ESO 137−001.

The model HIM masses vary between 2.5 an 5.9 × 10⁸ M_⊙. This is about half of the HIM mass derived from X-ray observations. The model WIM masses, which critically depend on the assumed volume filling factor of the WIM, range between 0.8 an 2.2 × 10⁸ M_⊙. This is quite close to the value derived from Hα observations. The WNM masses vary by about a factor of two between the different models (2.7 − 5.0 × 10⁸ M_⊙). The mass distribution of the different models is broadest for the CNM in the tail region¹ (0.4 − 4.3 × 10⁸ M_⊙). The CNM disk masses of all three models are comparable to the observed disk mass based on a Galactic CO-H₂ conversion factor. The CNM mass of the tail of model A is about half of the observed mass. The tail CNM masses of models B and C are a factor of 11 and 4 times smaller than that of model A, respectively.

Waldron et al. (2023) give a star formation rate of 1.2 M_⊙ yr⁻¹ in the disk of ESO 137−001. With a molecular gas mass within galactic disk of 7 × 10⁸ M_⊙ (Jáchym et al. 2019) the star formation efficiency (SFE) is 17 × 10⁻¹⁰ yr⁻¹. This is significantly higher than the average SFE in disks of local spiral galaxies (5.3 ± 2.5 × 10⁻¹⁰ yr⁻¹) and as high as the SFE in the centers of NGC 4736 and NGC 3351 (Leroy et al. 2008).

Leroy et al. (2013) stated that molecular gas in the central regions of spiral galaxies leads to more CO emission, appears more excited, and forms stars more rapidly than molecular gas further out in the disks. Sandstrom et al. (2013) found that most galaxies exhibit a lower conversion factor in the central kpc by factor of about two below the galaxy mean, on average. Moreover, the CO-H₂ conversion factor in the central part of the Galaxy is about four times lower than that in the disk (Bolatto et al. 2013). According to Sandstrom et al. (2013), NGC 4736 and NGC 3351 show a CO-H₂ conversion factor in the central kpc, which is more than five times lower than the galactic value. With a two times lower conversion factor than Galactic the CNM mass within the galactic disk of ESO 137−001 would be ∼4.5 × 10⁸ M_⊙ consistent with the molecular gas mass of model A. The SFE then is 8 × 10⁻⁹ yr⁻¹.

The molecular gas mass in the tail region is 9 × 10⁸ M_⊙ based on a Galactic CO-H₂ conversion factor. Comparison of fluxes of the ALMA+ACA observations with previous single-dish (APEX) observations (Jáchym et al. 2019) indicated that, in addition to the compact CO features, there is a substantial component (up to a factor of 3 to 4) of extended (scales > 6 kpc) molecular gas in the tail. Since a CO-H₂ conversion factor significantly lower than Galactic is expected in the tail, the missing ALMA flux might be compensated by the assumed Galactic CO-H₂ conversion factor. The star formation rate associated with this gas is about 0.5 M_⊙ yr⁻¹. Thus, whereas the CO flux in the disk region is as high as that of the tail, the associated star formation is about twice that of the tail. Since we integrated the entire star formation activity of our models without a sensitivity cutoff, our model SFR in the tail region represents an upper limit. The ratio SFR_disk/SFR_tail is thus a lower limit. When we take this limitation into account, model A best reproduces the SFR found by Waldron et al. (2023).

We conclude that models A best reproduces the observed CO emission and SFR fractions between the disk and tail regions.

5.3. The dense molecular gas and CO emission

The observed (Jáchym et al. 2019) and model CO velocity fields together with the CO emission distribution are presented in Fig. 6. The CO disk emission shows a north-south velocity gradient due to galactic rotation. The observed velocity amplitude is about 50 km s⁻¹ (see also Fig. 4 of Jáchym et al. 2019). LOS velocities > 50 km s⁻¹ are observed in the Hα velocity field of the southern tail (Fig. 3 of Luo et al. 2023). The upturning extraplanar CO filament mainly has velocities close to the systemic velocity. Only the outer northern part shows negative velocities <  − 20 km s⁻¹.

As already mentioned in Sect. 4 the morphology of the model dense gas tail is different from the observed one: model A shows a northern dense model gas filament with an extent of about 15 kpc. This filament is almost absent in the ALMA observations. The observed upturning southern filament is best reproduced by models A. Model C also shows a southern upturning filament but its extent is smaller than 5 kpc.

The observed velocity field of the model A southern filament is consistent with observations for distances > 5 kpc from the galactic disk. At smaller distances the model velocities are significantly higher than the observed ones. Higher model velocities compared to observations are expected as the model rotation velocities are a factor of 1.5 higher than observed (Fig. 2). In addition, the model gas extent within the disk plane is larger than observed leading to higher LOS velocities due to the radially increasing rotation curve. The observed negative velocities at the outer northern tip of the upturning filament is reproduced by model A.

We conclude that model A best reproduces the CO emission distribution and velocity for distances ≲15 kpc from the galactic disk.

5.4. UV emission and star formation

The observed HST F275W NUV emission distribution (see also Fig. 2 of Waldron et al. 2023) together with the CO emission distribution are presented in Fig. 12. The HST image clearly shows extraplanar NUV emission to the west of the galactic disk in the direction of the gas tail. This diffuse emission extends about 40″ ∼ 13.6 kpc to the west, except in the north where the extent is ∼25 kpc. The model with the strongest extraplanar NUV emission is model A, which qualitatively reproduces the observations. The observed most northern NUV filament is not present in model A. The observed southern extension of the NUV disk is present in model A albeit less extended. The extraplanar NUV emission distributions of models B and C are much less prominent than their observed counterpart.

Fig. 12. CO emission (red contours) on NUV emission distribution. Upper left: ALMA and HST observations.

We conclude that model A best reproduces the available NUV observations.

5.5. HI emission

The atomic hydrogen mass distribution can be determined by taking into account the total gas mass, the molecular fraction (Eq. 9), and the amount of mixed gas (Sect. 2.3). We convolved the model H I maps to a spatial resolution of 19″. The resulting model H I moment 0 maps are presented in Fig. 13. All model maps show an extended off-center maximum, which is elongated toward the northwest. The extent of the high-surface-density gas significantly varies between the models. The smallest extent (< 20 kpc) is observed in model C, the largest extent (∼30 kpc) in model A. Model B shows and intermediate extent. The extent of the low-surface-density tail is about proportional to that of the high-surface-density gas. The largest extent is present in model A (≳50 kpc) and the smallest extent in model C.

Fig. 13. Stellar disk (contours) on H I emission distribution. The contour levels are (1.5, 3, 6, 12, 24)×1019 cm−2.

5.6. The Hα–X-ray correlation of the stripped gas

Recently, Sun et al. (2022) found a strikingly linear correlation between the Hα and X-ray surface brightnesses of the diffuse stripped gas of cluster spiral galaxies at ∼10 − 40 kpc scales. ESO 137−001 is one of the prime examples in this work. Lee et al. (2022) found in their 3D hydrodynamical simulations that the ICM-dominant tail gas shows flux ratio F_X/F_Hα ∼ 1 − 20, while the gas in the disk vicinity (3 kpc < z < 10 kpc) exhibits a lower F_X/F_Hα of ∼1 (lower panel of their Fig. 15).

Triggered by these results, we established the Hα–X-ray correlation for the observations and our models. Sun et al. (2022) used uneven boxes to extract the mean surface brightnesses. For an objective comparison between observations and models we extracted the mean surface brightnesses on an equidistant grid with a grid size of 5 kpc. To do so we first convolved the Hα map with a Gaussian with a FWHM of 40 pixels = 8″ = 2.7 kpc. The resulting map was regridded to match the X-ray data in terms of pixel size (0.49″ = 0.17 kpc) and image size. The X-ray and Hα maps were then regridded to a common pixel size of 14.9″ = 5 kpc. This choice of a 5 kpc grid led to a linear slope for the observations (Figs. 14 and B.1), as it was derived by Sun et al. (2022). The Spearman rank coefficient is 0.77, meaning that the correlation is strong.

Fig. 14. Hα–X-ray correlation (right panels). The solid lines correspond to an outlier-resistant linear regression. Left panels: grayscale: X-ray emission distribution; contour: Hα emission distribution. The pixel size is ∼5 kpc. The observed Hα–X-ray correlation and the correlations of models B and C are presented in Fig. B.1.

For the models the Hα and X-ray maps were convolved with a Gaussian with a FWHM of 2.2 kpc. The Hα map was regridded to match the X-ray data in terms of pixel size (0.49″ = 0.17 kpc) and the X-ray and Hα maps were regridded to a common pixel size of 14.9″ = 5 kpc. Finally, the X-ray and Hα maps were clipped to yield the observed extents of the stripped gas tails. We verified that the resulting correlations are not sensitive to the FWHM of the convolutions.

The ranges of the Hα and X-ray surface brightnesses of all models are well comparable to the observed ranges. All three models show a clear correlation between the Hα and X-ray surface brightnesses (Spearman rank coefficients of 0.53, 0.69, and 0.62 for models A, B, and C, respectively). The correlations are weaker than the observed correlation. The slopes of the log–log correlations are 0.92 ± 0.11, 0.95 ± 0.10, and 1.08 ± 0.12 for models A, B, and C, respectively. The model slopes are thus consistent with the slope of the observed log–log correlation.

We conclude that our modeling of the Hα emission caused by ionization through thermal conduction is consistent with the results of Sun et al. (2022). We also tested a constant ionization fraction and an ionization fraction caused by the equilibrium between X-ray ionization and recombination. Both recipes led to a much shallower slope of the logarithmic Hα–X-ray correlation than observed. Thus, the ionization through thermal conduction sets the slope of the correlation. The model correlation scatter is mainly set by the variation of the model gas density. The main differences between the model and observed correlations are the higher Hα and X-ray brightnesses of the gas tail and the observed small scatter at high surface brightnesses. We can only speculate that the difference is caused by our approximate description of ram pressure stripping of diffuse gas (Sect. 2.3).

5.7. Heat transfer via turbulence and thermal conduction in the presence of a magnetic field

For the calculation of the emission measure of the stripped ISM due to thermal conduction (Eq. 12) the influence of the magnetic field was neglected. We assume that at the interface between the hot ICM and the warm stripped ISM Kelvin-Helmholtz instabilities develop, which lead to turbulent mixing of the hot ICM with the warm ISM (see Sect. 2.3). In this way turbulence is induced into the warm partially ionized ISM. Li et al. (2023) found a turbulent energy injection scale of L ∼ 1 kpc with an associated velocity of v_L ∼ 50 km s⁻¹. Following Eq. (5.32) of Sarazin (1988), the mean free path of electrons for T = 10⁷ K and n = 0.3 cm⁻³ is about λ ∼ 1 pc. If we further assume that the turbulence of the warm ISM is about sonic (Sonic Mach number M_S ∼ 1) with a temperature of 10⁴ K dictated by the cooling curve (v_turb ∼ 10 km s⁻¹), the Alfvénic Mach number is about 5.

Following Lazarian (2006) the scale at which the magnetic field gets dynamically important is $l_{\mathrm{A}}=L{M}_{\mathrm{A}}^{-3}\sim 8$ pc and thus l_A > λ. In this case Lazarian (2006) stated that heat conduction is decreased by a factor of three with respect to the classical Spitzer value (Sect. 2.2 of Lazarian 2006). According to Eqs. (7), (21), and (25) of McKee & Cowie (1977) the emission measure is proportional to the square root of the conductivity. A reduction of the emission measure by a factor of $\sqrt{(}3)$ due to the presence of magnetic fields is within the uncertainties of our calculations. We note that the global heat transfer is dominated by turbulent advective heat transfer (Fig. 1 of Lazarian 2006 with M_S > λ/(α β L) with β = 4 and $\alpha ={(m_{\mathrm{e}}/m_{\mathrm{p}})}^{\frac{1}{2}}$) and ($M_{\mathrm{A}}=5<{(L/\lambda )}^{\frac{1}{3}}=10$) as assumed in Sect. 2.3.

As a consistency check we can compare the density based on the observed mean emission measure to that based on the model overdensities of the warm and hot gas. The extent of the ionized region l of a stripped gas cloud can be estimated by setting the heat conduction timescale t_heat = l²/(λ v_e) to the turbulent timescale $t_{\mathrm{turb}}=L/(\sqrt{(}3)v_{\mathrm{turb}})\sim 11$ Myr, where v_e = 500 km s⁻¹ is the velocity of the thermal electrons. This yields l = 76 pc and 44 pc in the absence and presence of magnetic fields, respectively. An Hα surface brightness of the tail of 2 × 10⁻¹⁷ erg cm⁻² s⁻¹ arcsec⁻² (Sun et al. 2022) corresponds to an emission measure of EM ∼ 0.6 cm⁻⁶ pc, which leads to ionized gas densities of 0.09 and 0.12 cm⁻³, respectively. With an overdensity of 70 for the hot gas and 1000 for the warm gas (Sect. 2.5) the warm gas is ∼14 times denser than the hot gas. With an observed mean hot gas density of 10⁻² cm⁻³ (Sun et al. 2010) the warm gas density is thus 0.14 cm⁻³, well comparable to the value based on the observed emission measure.

5.8. ESO 137–001 within the Norma cluster

We adopt the systemic velocity of ESO 137−001 of 4647 km s⁻¹ from Luo et al. (2023). The cluster mean velocity is 4871 km s⁻¹ (Woudt et al. 2008). ESO 137−001 thus has a systemic velocity of 224 km s⁻¹ with respect to the cluster mean velocity. Assuming a LOS component of the unit galaxy velocity vector of −0.07 (model 5 of Table C.1) leads to a total galaxy velocity of ∼3200 km s⁻¹. Its projected distance to the cluster center is 180 kpc. Thus, the galaxy’s smallest distance to the cluster center and thus the ram pressure peak is expected to occur in about 50 Myr. This is somewhat longer but comparable to the time to peak ram pressure of the highest-ranked model A (Δt = −40 to −20 Myr; Table C.1). The associated ram pressure is 49 000 and 59 000 cm⁻³(km s⁻¹)². With a total velocity of 3200 km s⁻¹ the ICM density at the location of ESO 137−001 is n_ICM = 4.8 − 5.8 × 10⁻³ cm⁻³. This is about four times higher than the ICM density derived from X-ray observations (Sun et al. 2010).

As the next step we calculated possible orbits of ESO 137−001 within the Norma cluster. For the gravitational potential of the Norma cluster follow Jáchym et al. (2014) by assuming an NFW halo density profile with a Virial mass of M_vir = 10¹⁵ M_⊙ and scaling radius r_s = 346 kpc. For the ICM density distribution we used a β profile

$$\begin{array}{c}\hfill n_{\mathrm{e}}=n_{\mathrm{e},0}(1+{(r/r_{\mathrm{c}})}^2{)}^{-3/2\beta }\end{array}$$(15)

with n_e, 0 = 2.4 × 10⁻³ cm⁻³, r_c = 9.95′ = 200 kpc, and β = 0.555 (Boehringer et al. 1996). The observed X-ray emission of the Norma cluster ICM is displaced from the cluster center (ESO 137−006) to the northwest by about 150 kpc (Fig. 1 of Jáchym et al. 2014). We thus displaced the center of the distribution of Eq. (15) to (100 kpc, 100 kpc, 0). For the initial conditions we used the projected distance and the 3D velocity of ESO 137−001, which we derived from the highest-ranked model A. To simplify the model, we kept the orbit in the plane of the sky (x − y plane). x₀ = 300 kpc; y₀ = 100 kpc; z₀ = 0; v_x, 0 = −3180 km s⁻¹; v_y, 0 = −280 km s⁻¹; v_z, 0 = 0. The resulting galaxy orbit is presented in the upper panel, the resulting ram pressure profile in the middle panel of Fig. 15.

Fig. 15. A plausible orbit of ESO 137−001 within the Norma cluster. Upper panel: grayscale: ICM distribution. The timestep t = 0 Myr corresponds to the current location of ESO 137−001. The orbit is extrapolated into the future. The galaxy approaches the cluster center for the first time. Middle panel: ram pressure profile. Lower panel: solid line: last ram pressure stripping event of middle panel; dotted line: solid line multiplied by a factor of 2.5; dashed line: ram pressure profile of the highest-ranked model A (red line in Fig. 3).

This galaxy orbit is close to unbound and consistent with the orbital solutions found by Jáchym et al. (2014). If the velocity of the galaxy is increased by a factor of 1.2 the orbit becomes unbound. The initial total galaxy velocity for the unbound orbit is 3840 km s⁻¹ instead of 3200 km s⁻¹ derived from the highest-ranked model A. The last ram pressure maximum of the upper panel is enlarged in the lower panel of Fig. 15 (solid line) together with the ram pressure profile of the highest-ranked model A (dashed line). The maximum ram pressure of the model orbit is significantly smaller than that of model A. Modifications of y₀ did not lead to a significantly higher peak ram pressure. On the other hand, when the last ram pressure profile is multiplied by a factor of 2.5 (dotted line) it is consistent with the ram pressure profile of the highest-ranked model A.

This increase of the ram pressure profile can be obtained in two ways: (i) an increase of the ICM density or (ii) a moving ICM opposite to the motion of the galaxy within the cluster as it occurs for NGC 4522 in the Virgo cluster (Kenney et al. 2004; Vollmer et al. 2004). In the latter case the ICM velocity adds to the galaxy velocity and ram pressure increases quadratically with velocity. An ICM velocity of 3000 km s⁻¹ with respect to the cluster mean velocity in the opposite direction to the motion of ESO 137−001 leads to an increase of ram pressure by a factor of two. With a sound speed of $1263\sqrt{kT/6\mathrm{keV}}$ km s⁻¹ this corresponds to a Mach number of 2.4, comparable to that of the Bullet cluster (Markevitch et al. 2002). In a strongly perturbed galaxy cluster as the Norma cluster with an off-center ICM distribution the two possibilities (i) and (ii) are probable and plausible. If the moving ICM leads to a significant modification of the ram pressure a Lorentzian profile as a function of time might not be expected. In this context, Tonnesen (2019) studied the influence of galaxy orbits within a cluster, Roediger et al. (2014) the influence of instantaneous stripping by means of eight “wind-tunnel” hydrodynamical simulations. It is beyond the scope of this work to test different shapes of the ram pressure profile.

6. Summary and conclusions

The Norma cluster galaxy ESO 137−001 is one of the rare ram-pressure-stripped galaxies for which deep, kiloparsec-resolution observations are available in the CO line, Hα line, and X-ray emission. Its stripped X-ray (Sun et al. 2010) and Hα (Sun et al. 2022) tails are spectacular, extending ∼80 kpc behind the galactic disk. ESO 137−001 evolves in the highly dynamic ICM of the Norma cluster. The cold molecular, warm ionized, and hot ionized gas coexist within the gas tail as traced by the CO (Jáchym et al. 2014, 2019), Hα, and X-ray emission.

We made dynamical simulations of ESO 137−001 to constrain its 3D orbit within the Norma cluster and to investigate the physics of the ISM–ICM turbulent mixing process. Our numerical code is not able to treat diffuse gas in a consistent way. For a realistic treatment, 3D hydrodynamical simulations should be adopted. Nevertheless, we can mimic the action of ram pressure on diffuse gas by applying very simple recipes based on the fact that the acceleration by ram pressure is inversely proportional to the gas surface density, which in turn depends on the gas density for a gas cloud of constant mass (Sect. 2.3). The hot (> 10⁶ K) diffuse gas is taken into account by assuming that once the stripped gas has left the galactic disk, it mixes with the ambient ICM. In a new approach, we estimate the ICM–ISM mixing analytically (Eq. 7) based on the recipe of Fielding et al. (2020).

Following Vollmer et al. (2001), we used Lorentzian and Gaussian profiles for the time evolution of ram-pressure stripping (Sect. 2.4). Following Nehlig et al. (2016) and Vollmer et al. (2018), we investigated (i) the influence of galactic structure (i.e., the position of the spiral arms) on the results of ram-pressure stripping by varying the time of peak ram pressure, and (ii) the influence of the angle between the disk plane and the ram-pressure wind on the resulting gas distribution and velocity field. Moreover, we recalculated the initial model set with different ICM–ISM mixing rates (Sect. 4.1) and stripping efficiencies of the hot gas (Sect. 4.2).

For the modeling of the warm diffuse ISM, we assume that the gas with temperatures lower than 10⁶ K is ionized by thermal conduction (Sect. 2.5): thermal electrons from the hot ICM penetrate into the neutral, warm, stripped ISM clouds, ionizing and heating them. Ultimately, this leads to the evaporation of the ISM clouds (McKee & Cowie 1977).

To search for the highest-ranked models, we calculated the goodness of the fit for all time steps of all models (Sect. 3). This was also done for the velocity field. Whereas the highest-ranked model at a single wavelength corresponds to the minimum goodness, the comparison of goodness at different wavelengths is not straight forward. We decided to rank the models at the different wavelengths to calculate the sum of the ranks at the different wavelengths. We define the highest-ranked model as the model with the smallest value of the sum of the ranks.

Based on a detailed comparison between our dynamical models and the multiwavelength observation, we come to the following conclusions:

• The highest-ranked models with the nominal ICM–ISM mixing rate (Eq. 7) reproduce observations significantly better than the models with a three-times higher or lower ICM–ISM mixing rate (Sect. 4).

• We are not able to reproduce all observational characteristic with a single model. The X-ray and Hα observations are better reproduced by the preferred model C, whereas the CO observations are better reproduced by the highest-ranked model A (Sect. 4).

• The angle between the direction of the galaxy’s motion and the plane of the galactic disk is between 60° and 75°. Ram-pressure stripping thus occurs more face-on (Sect. 4).

• The two-tailed structure is a common feature in our models, and is due to the combined action of ram pressure and rotation together with the projection of the galaxy on the sky (Sect. 4).

• Structures like the observed northern faint Hα tail can in principle be reproduced by the model (Sect. 5.1). Their existence depends on the density structure of the stripped gas, which itself depends on the initial gas distribution and the temporal ram-pressure profile.

• Overall, the CO/Hα/X-ray highest-ranked models, with a three-times lower stripping efficiency reproduce the available observations less well than the highest-ranked models with the nominal stripping efficiency (Sect. 4.2).

• Model A best reproduces the observed CO emission and SFR fractions between the disk and tail regions (Sect. 5.2).

• Model A also best reproduces the CO emission distribution, velocity for distances of ≲20 kpc from the galactic disk, and the available NUV observations (Sects. 5.3 and 5.4).

The recently established linear correlation between the Hα and X-ray surface brightnesses (Sun et al. 2022) is reproduced by all three of the models studied here, albeit with a weaker correlation strength (Sect. 5.6). Our modeling of the Hα emission caused by ionization through thermal conduction is thus consistent with observations: the thermal electrons of the mixed ICM–ISM at T ∼ 10⁷ K penetrate into and ionize the warm stripped ISM.

Turbulent mixing is essential for modeling X-ray emission, and heat conduction is essential for modeling Hα emission. For the efficiency of heat conduction, knowledge of the magnetic field strength, which can be estimated via radio continuum observations, is important. We think that future 3D magnetohydrodynamic ram-pressure-stripping simulations should be able to resolve the turbulent mixing of diffuse gas in the tail (Tonnesen & Stone 2014; Ruszkowski et al. 2014) and account for thermal conduction.

We also predict the H I emission distributions for the different models (Fig. 13). The observed total H I mass will be a critical test of our models and an important constraint to add to them. Based on the 3D velocity vector derived from our dynamical model, we derive a galaxy orbit, which is close to unbound (Sect. 5.8). We argue that ram pressure is enhanced by a factor of ∼2.5 compared to that predicted for an orbit in an unperturbed spherical ICM. This increase can be obtained in two ways: (i) an increase in the ICM density or (ii) a moving ICM opposite to the motion of the galaxy within the cluster. In a strongly perturbed galaxy cluster, such as the Norma cluster, with an off-center ICM distribution, the two possibilities are probable and plausible.

Data availability

Appendices C–F can be found on Zenodo (https://doi.org/10.5281/zenodo.13837789). Movies associated to Figs. 5, 7, and 8 are available at https://www.aanda.org

Acknowledgments

The authors would like to thank J.D.P. Kenney for useful comments on the article and the anonymous referee for their constructive questions and suggestions.

References

Appendix A: Goodness distribution

Fig. A.1. Goodness distribution for the comparison based on the X-ray, Hα, and CO maps and the Hα velocity field. The straight solid line represents a fit to the distribution of the first 100 highest-ranked models. The dotted vertical lines correspond to the first 50 and 100 highest-ranked models, respectively.

Appendix B: Hα–X-ray correlation

Fig. B.1. Hα–X-ray correlation (right panels). The solid lines correspond to an outlier-resistant linear regression. Left panels: grayscale: X-ray emission distribution; contour: Hα emission distribution. The pixel size is ∼5 kpc.

All Tables

All Figures

Fig. 1. ESO 137−001. Blue: X-rays (Sun et al. 2010). Red: Hα (Sun et al. 2022). Green: CO (Jáchym et al. 2019). The main features of the system are labeled. It appears that the X-ray stripping front is ahead of the Hα stripping front because of the large spatial smoothing scale of the X-ray image, but in reality the X-ray and Hα stripping fronts are at the same position, as shown in Fig. 2 of Luo et al. (2023). For a more detailed image, see Fig. 4.

In the text

Fig. 2. Initial model conditions. Upper panel: Model stellar- (black solid line) and gas-mass (black dashed line) surface density profiles. Solid red line: Exponential stellar-mass surface density profile with a scale length of 2 kpc. Dashed red line: Analytical approximation to the gas surface density profile (see text). Lower panel: Model (black solid line) and observed asymmetric-drift-corrected (red line; Luo et al. 2023) rotation curve. Red dashed line: Observed rotation curve multiplied by a factor of 1.5.

In the text

	Fig. 3. Model ram-pressure profiles (Table 2). The thick red line shows the highest-ranked model, and the vertical red lines show the range of time steps of the highest-ranked model.
In the text

Fig. 4. Multiwavelength observations of ESO 137−001. Upper panel: MUSE Hα is shown in color (Sun et al. 2022; Luo et al. 2023); Chandra X-ray in black contours (Sun et al. 2010); and ALMA CO in white contours (Jáchym et al. 2019). It appears that the X-ray stripping front is ahead of the Hα stripping front because of the large spatial smoothing scale of the X-ray image, but in reality, the X-ray and Hα stripping fronts are at the same position, as shown in Fig. 2 of Luo et al. (2023). Lower panel: MUSE Hα velocity field (Sun et al. 2022). The color scale corresponds to the LOS velocities with respect to the systemic velocity in km s−1. At a distance of 70 Mpc 30″ correspond to about 10 kpc. The size of the images is 80 kpc × 60 kpc.

In the text

Fig. 5. Highest-ranked model A of ESO 137−001. Upper panel: Hα is shown in color; X-ray in dark gray contours; CO in white contours; and stellar content in black contours. The relative contours are the same as in Fig. 4, i.e., Hα: logarithmic transfer function with 15 contour levels; X-ray and CO: square root transfer function with 15 contour levels. Lower panel: Hα velocity field. The colors are the same as in Fig. 4. The corresponding time evolution and that of the mass surface densities of the different gas phases can be found be found in online movies.

In the text

	Fig. 6. CO emission distribution (contours) on CO velocity field. Upper left: ALMA observations. The color levels are the same in all panels.
In the text

	Fig. 7. Preferred model B of ESO 137−001. Upper panel: Hα is shown in color; X-ray in dark gray contours; CO in white contours; stellar content in black contours. Lower panel: Hα velocity field. See Fig. 5 for the description of the contours and colors. The corresponding time evolution and that of the mass surface densities of the different gas phases can be found be found in online movies.
In the text

Fig. 8. Preferred model C of ESO 137−001. Upper panel: Hα is shown in color; X-ray in dark gray contours; CO in white contours; and stellar content in black contours. Lower panel: Hα velocity field. See Fig. 5 for the description of the contours and colors. The corresponding time evolution can be found online (observations_modelC.gif). The time evolution of the mass surface densities of the different gas phases can be found online (gasphases_modelC.gif).

In the text

Fig. 9. ESO 137−001 preferred models. Left panel: Model corresponding to the preferred model B of Fig. 7 with a three times lower stripping efficiency with respect to the nominal stripping efficiency. Right panel: Model close to the preferred model C (Δt = −20 Myr instead of Δt = −10 Myr) of Fig. 8. Upper panel: Hα is shown in color; X-ray in dark gray contours; CO in white contours; and stellar content in black contours. Lower panel: Hα velocity field. See Fig. 5 for the description of the contours and colors.

In the text

	Fig. 10. ESO 137−001 model showing detached elongated Hα emission regions in the north and south of the two main tails.
In the text

	Fig. 11. Fractions of stripped (2 kpc ≤ x ≤ 80 kpc) gas in the different phases. The time t = 0 Myr corresponds to peak ram pressure.
In the text

	Fig. 12. CO emission (red contours) on NUV emission distribution. Upper left: ALMA and HST observations.
In the text

	Fig. 13. Stellar disk (contours) on H I emission distribution. The contour levels are (1.5, 3, 6, 12, 24)×1019 cm−2.
In the text

	Fig. 14. Hα–X-ray correlation (right panels). The solid lines correspond to an outlier-resistant linear regression. Left panels: grayscale: X-ray emission distribution; contour: Hα emission distribution. The pixel size is ∼5 kpc. The observed Hα–X-ray correlation and the correlations of models B and C are presented in Fig. B.1.
In the text

Fig. 15. A plausible orbit of ESO 137−001 within the Norma cluster. Upper panel: grayscale: ICM distribution. The timestep t = 0 Myr corresponds to the current location of ESO 137−001. The orbit is extrapolated into the future. The galaxy approaches the cluster center for the first time. Middle panel: ram pressure profile. Lower panel: solid line: last ram pressure stripping event of middle panel; dotted line: solid line multiplied by a factor of 2.5; dashed line: ram pressure profile of the highest-ranked model A (red line in Fig. 3).

In the text

	Fig. A.1. Goodness distribution for the comparison based on the X-ray, Hα, and CO maps and the Hα velocity field. The straight solid line represents a fit to the distribution of the first 100 highest-ranked models. The dotted vertical lines correspond to the first 50 and 100 highest-ranked models, respectively.
In the text

	Fig. B.1. Hα–X-ray correlation (right panels). The solid lines correspond to an outlier-resistant linear regression. Left panels: grayscale: X-ray emission distribution; contour: Hα emission distribution. The pixel size is ∼5 kpc.
In the text