The Launching of Cold Clouds by Galaxy Outflows. V. The Role of Anisotropic Thermal Conduction

The full system of MHD equations with radiative cooling and conduction solved in our simulations is

$$\begin{eqnarray}&&{\partial }_{t}\rho +{\rm{\nabla }}\cdot (\rho {\boldsymbol{u}})=0,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\rho [{\partial }_{t}{\boldsymbol{u}}+({\boldsymbol{u}}\cdot {\rm{\nabla }}){\boldsymbol{u}}]+{\rm{\nabla }}\left(p+\displaystyle \frac{1}{2}| {\boldsymbol{B}}{| }^{2}\right)-({\boldsymbol{B}}\cdot {\rm{\nabla }}){\boldsymbol{B}}=0,\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{{\rm{\partial }}}_{t}E+{\rm{\nabla }}\cdot \left[\left(E+p+\displaystyle \frac{1}{2}|{\boldsymbol{B}}{|}^{2}\right){\boldsymbol{u}}-\left({\boldsymbol{B}}\cdot {\boldsymbol{u}}\right){\boldsymbol{B}}\right]\,+\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{n}^{2}{\rm{\Lambda }}(T)-{\rm{\nabla }}\cdot {\boldsymbol{q}}=0,\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{\partial }_{t}{\boldsymbol{B}}-{\rm{\nabla }}\times ({\boldsymbol{B}}\times {\boldsymbol{u}})=0,\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\boldsymbol{B}}=0,\end{eqnarray} \tag{ 14 }$$

where ρ is the density, u is the velocity, p = k_B T ρ/(μ m_p ) is the pressure, $E=p/(\gamma -1)+\tfrac{1}{2}\rho | {\boldsymbol{u}}{| }^{2}+\tfrac{1}{2}| {\boldsymbol{B}}{| }^{2}$ is the total energy density, and q is the electron thermal conduction.

Unlike in our previous papers, the simulations here were conducted using the performance-portable, open-source AthenaPK code. ³ AthenaPK implements various finite volume hydro- and magnetohydrodynamics algorithms in addition to Parthenon (Grete et al. 2022), which itself is a performance-portable AMR framework based on Athena++ (Stone et al. 2020), K-Athena (Grete et al. 2021), and Kokkos (Edwards et al. 2014; Trott et al. 2021). AthenaPK allows for efficient simulations with AMR on various devices including GPUs from different vendors and has demonstrated scalability with >90% parallel efficiency on up to 73,728 GPUs.

Here, we use an unsplit second-order, two-stage predictor-corrector time integrator, piecewise parabolic reconstruction of fluxes at cell boundaries (PPM), GLM divergence cleaning (Dedner et al. 2002), ⁴ and an HLL Riemann solver. The HLL family of Riemann solvers dates back to the work by Harten et al. (1983) and has the distinct advantage of being computationally simple but it only captures the fast magnetosonic waves. For convergence studies, we have also experimented with a piecewise linear (PLM) reconstruction scheme, which turned out to show a lower spectral bandwidth and smoother features (see Table 1 for an overview of runs).

In order to cope with the extreme conditions in the simulations, we implemented a first-order flux correction scheme. For cells in which the expected multidimensional update from the scheme above results in a negative temperature or density, all fluxes are recalculated using a forward Euler step with first-order reconstruction and LLF Riemann solver. This guarantees a positive, physically consistent update without the need to add artificial floors in the simulation.

We use the adaptive-mesh refinement capabilities of AthenaPK to refine the massless scalar where its concentration is ≥0.01 and de-refine if the concentration in the entire mesh block falls ≤0.001. The size of the computational box is 3000 × 1500 × 1500 pc³ resolved with 512 × 256 × 256 cells split into a mesh block of 32³ cells. It is essential that the box is large enough that it does not affect the flow around the cloud and that material does not escape laterally out of the box. The choice of our box ensures that this does not happen. In production runs, we used two levels of refinement, yielding an effective resolution of 1500 pc/256/4 = 1.46 pc, which corresponds to about 68 cells per cloud radius. In order to reduce computational costs, the adaptive refinement is disabled after 1 t_cc. This only pertains to the adaptive machinery so that the then static refinement, particularly around the cloud, remains in place. Finally, we use an adiabatic equation of state with a ratio of specific heats γ = 5/3 and tabulated cooling using the cooling table of Schure et al. (2009) for solar metallicity.

An important length scale to keep in mind is the Field length, ${\lambda }_{{\rm{F}}}=\sqrt{(\kappa (T)T/({n}^{2}{\rm{\Lambda }})}$ (Field 1965; Begelman & McKee 1990; McKee & Begelman 1990), the maximum length scale across which thermal conduction can dominate over radiative cooling in the unsaturated case. For 0.1 Spitzer conduction ${\lambda }_{{\rm{F}}}\approx 13\,\mathrm{pc}{T}_{7}^{7/4}{n}^{-1}{{\rm{\Lambda }}}_{-22}^{-1/2}$ cm. A second characteristic length scale is the cooling length λ_cool = c_s t_cool ≈ 2 × 10¹³ T₇ n⁻¹ Λ₋₂₂ cm. As described below in all our runs, R_c = 3.08 × 10²⁰ cm n = 1 cm⁻³ within the cloud and n = 0.001 cm⁻³ in the exterior medium. This means that in the ambient medium Δx ≪ λ_F, λ_cool, and inside the cloud our setup does not allow any cooling because at temperatures T < 10⁴ K, cooling is switched off. Finally, as soon as cloud material evaporates, the resulting high pressure will not allow for fragmentation due to cooling, and thus convergence is not likely to be an issue in our simulations.

New to our simulations is the ability to include anisotropic thermal conduction. Explicit integration of diffusive physics is usually hampered by a very restrictive time step limit that is inversely proportional to the square of the spatial resolution. When the diffusive terms are relatively large, or at very high resolution, this time step limit can prohibit simulations for sufficiently long times. Therefore, a Runge–Kutta–Legendre (RKL) super-time-stepping (STS) module Meyer et al. (2014) has been implemented. When STS is enabled, the diffusive equations are integrated forward in time by a separate super-time step in an operator-split update. Each super-time step is comprised of s stages and is equivalent to O(s²) times the explicit diffusive time step. Two operator-split super-time steps are required in a single (M)HD update for the second-order accurate RKL2 scheme.

Finally, to account for the different nature of classic and saturated fluxes (parabolic and hyperbolic, respectively), we follow Mignone et al. (2011) and use a smooth transition

$$\begin{eqnarray}&&{\boldsymbol{q}}=\displaystyle \frac{{q}_{\mathrm{sat}}}{{q}_{\mathrm{sat}}+{q}_{\mathrm{cl}}}{{\boldsymbol{q}}}_{\mathrm{cl}},\end{eqnarray} \tag{ 15 }$$

where q_sat is given by Equation (3) and Equation (5). We also employ limiters for calculating the temperature gradients following Sharma & Hammett (2007). This prevents unphysical conduction against the gradient, which may be introduced because the off-axis gradients are not centered on the interfaces.

As in previous papers in this series, we have carried out a suite of simulations that focuses on the conditions expected in galaxy outflows. The combination of supersonic flows, steep gradients in temperature and density, and the inclusion of radiative cooling, magnetic fields, and thermal conduction poses a formidable challenge to computational codes. Given the computational costs of carrying out such calculations, we have chosen only a limited set of parameters that are listed in Table 1. These runs were chosen to focus on a Mach number of 3.54, which is consistent with a cloud being impacted by a hot flow just outside of the driving radius of the galaxy, and a Mach 1 case to contrast these runs with a transonic case. Both parameter choices are also directly comparable to our results in previous papers of this series.

As before, we simulate a spherical cloud of radius 100 pc and temperature T_c = 10⁴ K in pressure equilibrium with the ambient gas. The cloud has a fiducial mean density of ρ_c = 10⁻²⁴ g cm⁻³, such that R_c n_i,c = 1.5 × 10²⁰ cm⁻². Similar to Grønnow et al. (2018) and Kooij et al. (2021), and unlike Papers I–III, we use a density profile that is smooth at the cloud edges and prescribed by

$$\begin{eqnarray}&&\rho (r)={\rho }_{{\rm{h}}}+\displaystyle \frac{1}{2}\left({\rho }_{{\rm{c}}}-{\rho }_{{\rm{h}}}\right)\times \left\{1-\tanh \left[s\left(\displaystyle \frac{r}{{R}_{{\rm{c}}}}-1\right)\right]\right\},\end{eqnarray} \tag{ 16 }$$

with steepness parameter s = 10. The Jeans length in the cloud is c_s(G ρ)^−1/2 ≈ 2 kpc, which means that the cloud is pressure confined rather than gravitationally bound, and justifies the neglect of gravity in our simulations. In order to mark the cloud material, we also include a massless scalar that is set to 1 in the initial cloud volume and set to 0 elsewhere.

The initial gas velocity is set to 0 within 1.3R_c and to a fixed speed, v_h, elsewhere. The gas flows in from the lower y-boundary with fixed v_h, density, and pressure, and leaves the computational domain on the upper y-boundary. This means that we have an inflow boundary condition on the lower y-boundary and outflow boundary conditions everywhere else. The magnetic fields are initialized either parallel to the flow (called parallel fields) or perpendicular to the flow (called perpendicular fields). We chose only to simulate weak magnetic fields as we are interested in the effect of the field direction on the anisotropy of thermal conduction and not on the dynamical effects of magnetic fields, although we also conduct two runs with a significant field strength to compare against our previous results, as discussed in more detail below. Hence, we set the magnetic field strength such that the ratio of the thermal pressure to the magnetic pressure, β, is 100 for the majority of our runs and β = 10 for our comparison runs. The full set of run parameters spanned by our simulations is described in Table 1.

Figure 1 shows the evolution of the dense cloud mass in our Mach 3.5 simulations. Here, as we have done in previous papers, we define cloud material as anything with a density greater than 1/3 of the initial density of the cloud. In the left plot of this figure, we contrast the evolution of the hydrodynamic case with two MHD cases with weak (β = 100) magnetic fields: one in which the field is initially parallel to the direction of the winds, and one in which the field is initially aligned in the perpendicular direction. Images of density, thermal pressure, velocity, and plasma-β from these simulations at t = 4t_cc are given in Figure 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In hydrodynamic interactions, the disruption of the cloud is caused by the impact of the exterior wind along both sides of the cloud and in the direction of the flow itself (e.g., Klein et al. 1994; Xu & Stone 1995; Fragile et al. 2004; Melioli 2005; Cooper et al. 2009; Marinacci et al. 2010, 2011; Banda-Barragán et al. 2019; Paper I). Along the sides of the cloud, the velocity shear leads to a mixing layer that grows at a rate of Δv_KH ∝ Δv χ^−1/2, where Δv is the velocity difference between the cloud and the hot medium (Chandrasekhar 1961), meaning that the cloud will be disrupted on a timescale of R/Δv_KH ∝ t_cc. This rate can be understood in terms of entrainment into a spatially growing shear layer made up of large-scale vortical structures convecting at a velocity v_c (Coles 1981; Dimotakis 1986), which exists even in the presence of cooling. However, the growth of the KH instability is slowed significantly if the exterior flow is supersonic (e.g., Chinzei et al. 1986; Papamoschou & Roshko 1988; Barre et al. 1994; Slessor et al. 2000) leading to a large increase in the cloud lifetimes in high-Mach number interactions (e.g., Paper I; Schneider & Robertson 2017).

In the streamwise direction, the most important effects are the buoyancy-driven Raleigh–Taylor (RT) instability (e.g., Chandrasekhar 1961), which operates effectively in the subsonic case, and the stretching of the cloud by pressure gradients, which operates effectively in the supersonic case (Paper I). If the exterior flow is subsonic, the acceleration of the dense cloud by the lower-density wind leads to the growth of the RT instability at a rate of h = α_b Agt², where g is the acceleration, A = (ρ_c − ρ_h)/(ρ_c − ρ_h) ≈ 1 is the Atwood number of the interaction, and α_b ≈ 0.03 − 0.07 is a constant that depends weakly on the initial geometry (Read 1984; Alon et al. 1995; Dimonte et al. 2004; Dimonte & Tipton 2006; Scannapieco & Brüggen 2008). Since $g\approx {\rho }_{{\rm{h}}}{v}_{{\rm{h}}}^{2}\pi {R}^{2}/({\rho }_{{\rm{c}}}\tfrac{4\pi }{3}{R}^{3})\approx {\chi }^{-1}{v}_{{\rm{h}}}^{2}/R$, setting h equal the size of the cloud gives a disruption time $\approx {\alpha }_{{\rm{b}}}^{-1/2}{\chi }^{1/2}R/{v}_{{\rm{h}}}$ or a few times t_cc.

In supersonic cases, such as the ones modeled here, the evolution of the cloud is dominated by the strong difference in pressure between the front and the back of the cloud (Paper I). In this case, the front of the cloud experiences an increase in pressure approaching the 1 + M² enhancement expected for a normal shock. On the other hand, the pressure increase downstream from the front of the clouds is more modest, scaling approximately as 1 + M due to the fact that this material passes through an oblique shock. As described in Paper I, this results in the stretching of the cloud on a timescale proportional to ${t}_{\mathrm{cc}}\sqrt{1+M},$ which serves as the dominant mechanism leading to cloud disruption in the supersonic case.

The presence of magnetic fields has a noticeable impact on these disruption processes, even at the moderate β = 100 values taken for these simulations. In the M3.5-par case in which the field lines are initially in the direction of the flow, the primary impact of the field is to stabilize the KH instability, such that the clumping of the material in the direction perpendicular to the flow is weaker than in the hydrodynamic case. Note, however, that as the KH instability is already weak in these supersonic runs, the impact of magnetic fields on the mass loss in the cloud is relatively minor.

In the case in which the field is oriented perpendicular to the flow, however, the impact of the magnetic field on the mass loss is much more significant. In the M3.5-trans simulation, the field has two major effects: the first stabilizes the RT and KH instabilities through magnetic draping (e.g., Dursi & Pfrommer 2008; Banda-Barragán et al. 2015; Kooij et al. 2021), leading to less clumping than seen in the other two cases. Second, and more importantly, it acts to counteract the stretching of the cloud by the pressure gradients. This is due to the fact that stretching in the streamwise direction leads to a separation between the magnetic field lines threading the cloud, and a corresponding drop in field strength due to flux freezing. As the magnetic pressure is proportional to B² this leads to a significant drop in the pressure within the cloud, which opposes the expansion decreasing the mass loss in the cloud. Finally, the compression of magnetic fields due to the bending of field lines leads to a magnetically dominated sheath behind the cloud, which is visible as the low-β regions surrounding the highly collimated tail visible in the last row of Figure 2.

The central panel of Figure 1 compares the mass evolution in three cases with perpendicular magnetic fields and different levels of conduction: no conduction (repeating the case in the left panel), isotropic thermal conduction, and anisotropic thermal conduction. Corresponding images of the density, thermal pressure, field lines, velocity, and plasma-β from the conducting simulations at t = 4t_cc are given in the top two rows of Figure 3.

Zoom In Zoom Out Reset image size

A direct comparison of the mass evolution and morphology between the M3.5-trans case without conduction to the M3.5-trans-iso case with isotropic thermal conduction reveals strong differences. As was observed in Paper II in the case with isotropic conduction, a strong evaporative wind develops, which compresses the cloud due to ram pressure (see also Mellema et al. 2002; Fragile et al. 2004; Orlando et al. 2005; Johansson & Ziegler 2013). The result is a much denser, compact system, which remains largely decoupled from the outside flow, losing mass at a rate that can be understood from a balance of impinging thermal energy and evaporation. The evaporative wind is almost completely dominated by ram pressure and thermal pressure, with magnetic fields playing only a minor role, as evidenced by the large plasma-β observed in this region.

In strong contrast to the isotropic conduction case, the mass-loss rates and morphology in the M3.5-trans-aniso run are almost identical to the run without thermal conduction. This is because the magnetic fields passing through the cloud are bent back into a magnetically dominated sheath behind the cloud, where the thermal pressure is small. This means that there is very little thermal energy available to be carried in the direction of the field lines and into the cloud, causing conduction to be strongly reduced.

Finally, the rightmost panel of Figure 1 compares the mass evolution in three cases with initially parallel magnetic fields and different levels of conduction. Corresponding images from the conducting simulations at t = 4t_cc are given in the bottom two rows of Figure 3. As in the M3.5-trans-iso case, the inclusion of isotropic thermal conduction in the M3.5-par-iso run results in a strong evaporative wind that compresses the cloud and largely decouples in from the outside flow. Furthermore, because the strength of this wind is such that the thermal and ram pressure dominate over magnetic forces, the overall evolution of the M3.5-trans-iso and M3.5-par-iso runs are very similar. Both clouds lose mass at a roughly constant rate, set by the energy input evaporating the cloud. Here, we should also comment on the bump in mass evolution at early times. This was not observed in our Papers I–IV. The difference stems from a difference in our initial conditions. In Papers I–IV the cloud was modeled as a hard sphere whereas here we taper the edges of the cloud according to Equation (16) in order to avoid unphysically steep gradients. As a result, the impinging wind piles up in front of the cloud leading to an initial rise of the cloud mass. However, this is only a transient effect that does not affect our results at later times.

Finally, in the M3.5-par-aniso case, the cloud mass-loss rate is even slower than that in the M3.5-par case without conduction. This is due to the fact that the radiative wind is largely suppressed except in small regions near the front and back of the cloud. Near the front of the cloud, electrons flow along fields lines from the upstream material powering an evaporative wind that insulates it from the development of the RT instability. At the same time, material at the back of the cloud is able to conduct heat from the wake region, which remains hot and thermally dominated in this case. This balances the cloud from expansion in the streamwise direction, also significantly extending its lifetime.

Next, we turn to the acceleration of the cloud by the supersonic wind. In Figure 4, we show the evolution of the velocity of the material with ρ > ρ_c/3, with lines and panels matching those used in Figure 1. As shown in previous work, in the hydrodynamic case, the hot wind is inefficient at accelerating the cold cloud. This also proves to be true for all the cases studied here, such that even after five cloud-crushing times, the cold material attains a speed between 70 and 120 km s⁻¹, which is less than 10% of that of the ambient medium.

Zoom In Zoom Out Reset image size

Comparing the various runs shows that in many ways the cloud velocities follow the same trends as the mass evolution. In cases without conduction, we find that the streamwise velocity of the cloud is highest in the purely hydrodynamic and M3.5-par runs, whereas the cloud is about 25% slower for the M3.5-trans runs (after 5 t_cc). Switching on anisotropic conduction has very little effect on these results, such that the velocity evolution in the M3.5-trans- and M3.5-par-aniso cases closely follow the cases without thermal conduction.

However, as also found in Paper II, in the runs with isotropic conduction, the clouds are only accelerated to half the speed as in the anisotropic case. This is because the evaporative pressure compresses the cloud such that it presents a smaller cross section to the impinging wind. This can be seen, for instance, in the streamwise projection of the density, which is shown in Figure 5.

Zoom In Zoom Out Reset image size

Note, however, that the velocity evolution for the runs without conduction is different from the case of strong magnetic fields. In Paper III, it was shown that the clouds in β = 1 perpendicular fields attained a higher velocity than the clouds in parallel fields (see Figure 12 in Paper III). Apparently, perpendicular fields lead to additional acceleration if the fields are strong. This higher acceleration was attributed to a hundred-fold increase in magnetic pressure at the front of the cloud as the perpendicular field lines are compressed by the flow.

We confirmed this by comparing fixed grid simulations of our setup at β = 100, 10, and 1 (not shown). These confirmed that the β = 1 case is markedly different from the lower field cases, in that it leads to a significantly stronger acceleration in the perpendicular field case.

In all cases, the magnetic tension exceeds the effect of the magnetic pressure in front of the cloud. When β becomes as low as 1, the tension becomes so strong as to change the flow around the cloud. The opening angle of the shock around the cloud increases and the stand-off distance of the shock increases. The wind thus needs to push against a magnetic wall that couples to the cloud and leads to a higher acceleration than in cases with higher β.

To compare the evolution of a cloud surrounded by a supersonic medium to one impacted by a transonic wind, we repeated our simulations reducing the ambient wind speed to 480 km s⁻¹, but leaving all other parameters untouched, which results in a Mach number of M ≈ 1. Here, we only looked at a single, purely hydrodynamic run (M1-hydro) and two runs with anisotropic conduction and different initial magnetic field orientations (M1-par-aniso and M1-trans-aniso). The mass and velocity evolution of these runs are shown in Figure 7, and the corresponding slices of key parameters at t = 3t_cc are shown in Figure 6.

Zoom In Zoom Out Reset image size

As seen in previous simulations (e.g., Paper I) at lower wind speeds, the pressure differences across the cloud are reduced, but the RT and KH instabilities become more efficient, leading to a more expedient mass loss. Thus, in the hydrodynamic case, we can see that by 3t_cc the cloud is highly disrupted, with much of the mass fragmented into a tail made up of a series of small clumps. Note that unlike the M3.5-hydro run, the strong KH that exists in transonic interactions means that the density sheared material drops well below the ρ_c/3 threshold that we use to define the boundary of the cloud. As a result, the effective cross section of the cloud drops well below π times the initial cloud radius squared.

Zoom In Zoom Out Reset image size

In the case M1-trans-aniso with perpendicular fields and anisotropic conduction, on the other hand, draping of the magnetic fields acts to preserve the cloud from instabilities, leading to a more ordered flow and lower overall mass-loss rate. Interestingly, because of the strong KH seen in the M1-hydro run, this run's cross section of dense material is greater in the case with perpendicular fields than in the hydro case. Finally, as in the higher Mach number M3.5-trans-aniso simulation, the field lines are bent behind the cloud, forming a low-beta region in which the thermal energy density is low. Also, this greatly reduces the ability of electrons to heat the cloud, suppressing the role of evaporation in the interaction (as in the M3.5-trans-aniso run).

On the other hand, in the M1-par-aniso run with anisotropic conduction and fields in the direction of the flow, electrons can flow freely from the hot upstream material. The resulting evaporative flow then plays two key roles in the overall cloud evolution. First, evaporation in the flow direction leads to the presence of a shock that precedes the cloud, reducing the velocity of the flow as it passes over the cloud itself. This reduces the growth of instabilities, preserving the cloud for longer times and reducing mass loss into the tail. However, a second effect of conduction also is to cause evaporation between the dense regions that arise from the suppressed instabilities. This causes clumps to move away from each other even in the direction perpendicular to the field lines. Interestingly, this causes the transverse size of the cloud M1-par-aniso to be even larger than in the hydro run. However, the overall mass loss is similar between the M1-par-aniso and M1-hydro runs.

A striking feature in the left panel of Figure 7 is the behavior of the purely hydrodynamical run. After 2.5 t_cc, the total mass with density above ρ_c/3 increases again. This is caused by condensation in the tail of the cloud as discussed in Gronke & Oh (2018). They find that in clouds with radii larger than a critical radius, radiative cooling will lead to eventual mass growth. For a cloud that is disrupted in $\tilde{t}$ cloud-crushing times, the temperature of the mixed gas is

$$\begin{eqnarray}&&{T}_{\mathrm{mix}}=\displaystyle \frac{{T}_{{\rm{c}}}{\dot{m}}_{{\rm{c}}}+{T}_{\mathrm{ps}}{\dot{m}}_{{\rm{w}}}}{{\dot{m}}_{{\rm{c}}}+{\dot{m}}_{{\rm{w}}}}\approx \tilde{A}\tilde{t}{\chi }^{-1/2}{T}_{\mathrm{ps}},\end{eqnarray} \tag{ 17 }$$

where ${\dot{m}}_{{\rm{c}}}=(4\pi /3){R}_{c}^{3}{\rho }_{h}\chi /(\tilde{t}\,{t}_{\mathrm{cc}})$, ${\dot{m}}_{{\rm{w}}}=\pi \tilde{A}{R}_{c}^{2}{\rho }_{h}{v}_{h}$, T_ps is the post-shock temperature, and $\tilde{A}$ is a factor that describes the area of the impinging gas that is mixed into the layer in units of $\pi {R}_{c}^{2}.$ With $\tilde{A}\tilde{t}\approx 1$, in the transonic case ${\chi }^{-1/2}{T}_{\mathrm{ps}}\approx \sqrt{{T}_{c}{T}_{h}}$ and this reduces to Equation (1) in Gronke & Oh (2018).

This critical radius is then found by setting the cooling time in this mixing gas to η times the cloud-crushing time. Again assuming $\tilde{A}\tilde{t}\approx 1$ this gives

$$\begin{eqnarray}&&{R}_{\mathrm{crit}}=\displaystyle \frac{{v}_{h}{t}_{\mathrm{cool},\mathrm{mix}}}{{\chi }^{1/2}}{\eta }^{-1}\approx 2\,\mathrm{pc}\,\displaystyle \frac{{T}_{{\rm{c}},4}^{5/2}M}{{P}_{3}{{\rm{\Lambda }}}_{-21.4}({T}_{\mathrm{mix}})}\displaystyle \frac{{\chi }_{0}}{100}{\eta }^{-1}\displaystyle \frac{{T}_{\mathrm{ps}}^{2}}{{T}_{{\rm{w}}}^{2}},\end{eqnarray} \tag{ 18 }$$

where T_c,4 is the cloud temperature in units of 10⁴ K, P₃ = n_h T/(cm³ K), Λ_−21.4 is the cooling rate in units of 10^−21.4 erg cm³ s⁻¹ and T_mix is the temperature of the mixed gas. This reduces to Equation (2) in Gronke & Oh (2018) in the transonic case.

Plugging in the parameters for our M = 1 simulation, we find that our cloud radius of 100 pc is safely above R_crit ≈ 2 pc. Hence, the growth of the cloud mass due to the condensation is not unexpected. However, we do not find this behavior in our M = 3.5 runs in which the post-shock temperature and mixing temperature are both increased by a factor of 1 + M², leading to a critical radius that is ≈1000 times greater.

We also do not see much condensation in the M = 1 runs with magnetic fields since the fields suppress the mixing that leads to radiative condensation. However, there is some cold gas in a thin filament in the wake of the runs with perpendicular magnetic fields (see Figure 3) that can be attributed to cooling. The exact critical sizes of clouds are also subject to debate (see, e.g., Li et al. 2020) and we are aware of work that looks at this problem for magnetized winds (Hidalgo-Pineda et al. 2023).

Finally, the right panel of Figure 7 shows the velocity evolution of the dense ρ = ρ_c/3 material in the three cases. As was true at higher Mach numbers, the evolution of the run with parallel fields and anisotropic conduction shows a velocity evolution similar to the hydro case. This is because while the M1-par-aniso run contains clumps that have moved away from each other due to evaporation, the overall cross section of the clumps remains similar to the M1-hydro case. Because the cross section of dense material is greater in the M1-trans-aniso case than in the M1-hydro run, the momentum transferred to the cloud exceeds that in the other two cases. This leads to a similar velocity evolution to the hydro case, even though there is much less mass loss in the M1-trans-aniso run.

From a technical point of view, there are multiple differences between this investigation and previous work. Most importantly, in Paper II where we studied the impact of thermal conduction in the purely hydrodynamic case, the saturation of the conductive flux was implemented by taking the minimum value of the classic and saturated flux. Here, we use a smooth transition to account for the parabolic and hyperbolic nature of the fluxes, cf., Equation eq:flux-transition. Moreover, thermal conduction was treated using the general implicit diffusion solver in Flash whereas we employ an explicit, second-order accurate RKL2 super-time-stepping scheme. Finally, we used different cooling tables (though both assume solar metallicity).

The difference in cooling tables also applies to the simulations in Paper III in which we studied the impact of strong magnetic fields (Paper III). While we do not expect significant differences between the divergence cleaning approach (a purely parabolic cleaning scheme in Flash versus the mixed hyperbolic-parabolic scheme in this paper), the spatial reconstruction methods matter. In Paper II we employed a second-order accurate scheme whereas we use a third-order scheme here resulting in an extended spectral bandwidth and more small-scale structure in the tail of the cloud. We reproduced the results of Paper III by reverting to a PLM reconstruction scheme. It is interesting to note that the third-order scheme employed here leads to a higher effective resolution in the simulation and also shows much better convergence properties than the runs with PLM reconstruction. Some results of runs with differing spatial resolutions are shown in the Appendix.

In addition to the reconstruction scheme, we also had to be careful about the choice of the Riemann solver. In previous work, we had used the Harten–Lax–van Leer discontinuities (HLLD) approximate Riemann solver developed by Miyoshi & Kusano (2005). In our setup, with supersonic flows with strong density contrasts, this shock-capturing scheme led to artificial heating in isolated cells and numerical instability that is related to the so-called carbuncle instability. This numerical instability affects the numerical capturing of shock waves and was first noticed by Peery & Imlay (1988) while simulating a high-speed flow around a blunt body. The nature of this instability is still not fully understood and is usually mitigated by adding numerical dissipation (see, e.g., Fleischmann et al. 2020; Minoshima & Miyoshi 2021). We tried the LHLLD solver by Minoshima & Miyoshi (2021) but still encountered numerical instabilities in a subset of the simulations. We found that the combination of PPM reconstruction schemes with HLL Riemann solvers yielded the best results in terms of effective resolution and numerical stability though at the cost of the Riemann solver not resolving the contact discontinuity.

Our work is complementary to recent work by Jung et al. (2023) who studied clumpy and uniform clouds in transverse magnetized, subsonic flows. They found that clumpiness has a large effect on the evolution of the cloud, but they did not include conduction. It will be interesting to study the effect of clumpiness also in our setup. It is also complementary to recent work by Jennings et al. (2023) who looked at the effect of anisotropic thermal conduction on clouds in the (hotter) intracluster medium.