Taming the TuRMoiL: The Temperature Dependence of Turbulence in Cloud-Wind Interactions

Turbulent radiative mixing layers (TRMLs) play an important role in many astrophysical contexts where cool ($\lesssim 10^4$ K) clouds interact with hot flows (e.g., galactic winds, high velocity clouds, infalling satellites in halos and clusters). The fate of these clouds (as well as many of their observable properties) is dictated by the competition between turbulence and radiative cooling; however, turbulence in these multiphase flows remains poorly understood. We have investigated the emergent turbulence arising in the interaction between clouds and supersonic winds in hydrodynamic ENZO-E simulations. In order to obtain robust results, we employed multiple metrics to characterize the turbulent velocity, $v_{\rm turb}$. We find four primary results, when cooling is sufficient for cloud survival. First, $v_{\rm turb}$ manifests clear temperature dependence. Initially, $v_{\rm turb}$ roughly matches the scaling of sound speed on temperature. In gas hotter than the temperature where cooling peaks, this dependence weakens with time until $v_{\rm turb}$ is constant. Second, the relative velocity between the cloud and wind initially drives rapid growth of $v_{\rm turb}$. As it drops (from entrainment), $v_{\rm turb}$ starts to decay before it stabilizes at roughly half its maximum. At late times cooling flows appear to support turbulence. Third, the magnitude of $v_{\rm turb}$ scales with the ratio between the hot phase sound crossing time and the minimum cooling time. Finally, we find tentative evidence for a length-scale associated with resolving turbulence. Under-resolving this scale may cause violent shattering and affect the cloud's large-scale morphological properties.


INTRODUCTION
While scales and relevant physics may vary, interactions between regions of cooler gas and coherent flows of hotter gas are prominent in many contexts.These interactions are prevalent in the circumgalactic medium (CGM), such as high velocity clouds (e.g., Wakker & van Woerden 1997;Putman et al. 2012), ram-pressure stripping of infalling satellites (e.g., Emerick et al. 2016;Simons et al. 2020) and the resulting streams (e.g., Bland-Hawthorn et al. 2007;Bustard & Gronke 2022), or cooling flows from cosmic accretion (e.g.Mandelker et al. 2020).There are also instances of these interactions within the interstellar medium (ISM), like the stellarwind driven bubbles within star-forming clouds (e.g.Lancaster et al. 2021).They are also relevant to the ram-pressure stripping of cluster galaxies and star formation in the tails of jellyfish galaxies (e.g.Tonnesen & Bryan 2021).We take a particular interest in their role within galactic winds (e.g.Fielding & Bryan 2022).
Galactic winds are ubiquitous throughout cosmic time, and play a pivotal role in galaxy evolution; they regulate star formation and transport metals out of the interstellar medium (ISM) (Somerville & Davé 2015).Observations indicate that stellar-feedback-driven winds are inherently multiphase; they are composed of comoving gas phases that vary in temperatures by orders of magnitude (see Veilleux et al. 2005 andRupke 2018 for reviews of observational evidence).
Observations favor a model in which supernovae drive hot ≳10 6 K winds that accelerate and entrain clouds of cool ∼10 4 K gas from the ISM (e.g.Chevalier & Clegg 1985).This model is complicated by hydrodynamical instabilities that drive mixing of gas between the cloud and wind.Because the timescale for mixing to destroy the cloud (by homogenizing the gas phases) is shorter than the ram-pressure acceleration timescale, it's remarkably difficult to accelerate clouds before they're destroyed (Zhang et al. 2017).
Various ideas have been proposed to address this difficulty.Some, like magnetic shielding (e.g.McCourt et al. 2015;Grønnow et al. 2018;Cottle et al. 2020), may extend the cold-phase lifetime by reducing mixing (see also Forbes & Lin 2019, for other mechanisms).Others are alternative acceleration mechanisms like radiationpressure (e.g.Zhang et al. 2018) or cosmic rays (e.g.Wiener et al. 2019;Brüggen & Scannapieco 2020).Another idea suggests the remnants of destroyed clouds seed the in situ formation of clouds in cooling outflows (e.g.Thompson et al. 2015;Schneider et al. 2018;Lochhaas et al. 2021).
Radiative cooling is also known to extend the coldphase lifetime (e.g.Mellema et al. 2002;Fragile et al. 2004;Melioli et al. 2005;Cooper et al. 2009).This work focuses on the regime in which rapid cooling acts as a mechanism that facilitates cloud survival (e.g.Marinacci et al. 2010;Armillotta et al. 2016).In this regime, cooling in a thin layer of gas at the interface between the phases is able to overcome the destructive effects of mixing (Gronke & Oh 2018).As turbulent mixing feeds hot phase material into this layer, isobaric cooling removes the temperature differential in the new material (Fielding et al. 2020).This process facilitates the transfer of mass and momentum to the cold phase providing a powerful additional acceleration source and allowing cloud growth.Hereafter, we refer to this mechanism as turbulent radiative mixing layer (TRML) entrainment.
The literature largely agrees that the occurrence and efficacy of TRML entrainment is controlled by three principal dimensionless numbers: (i) the density contrast χ = ρ cl /ρ w between the cloud and the wind, (ii) the Mach number of the wind M w = v w /c s,hot , and (iii) the cooling efficiency ξ = τ mix /τ cool .Here, τ mix and τ cool specify the characteristic timescales for mixing and for cooling of the mixing layer.As in Fielding et al. (2020) we primarily consider ξ sh = t shear /t cool,min , where t shear = R cl /v w is the shear-time and t cool,min is the minimum cooling time.In practice, our choice for τ cool is similar to the popular option of using t cool,mix , the cooling time of gas within the mixing layer at T mix ∼ √ T cl T w and n mix ∼ √ n cl n w . 1 It has been suggested that the relevant cooling timescale is instead set by cooling in the hot, volume filling, wind phase (Li et al. 2020;Sparre et al. 2020).We reconcile differences between these cooling timescales in follow-up work (Abruzzo et al., in prep.).Despite the obvious central importance of turbulence mechanisms underlying the operation of TRMLs, we do not yet have a clear understanding of how the turbulent velocity v turb changes as the 3 principal dimensionless numbers (χ, M w , and ξ) are varied.This is closely related to two fundamental unanswered questions.
(i) What is the role of cooling in driving turbulence?Shear layer studies find no or very weak cooling time dependence of v turb (Fielding et al. 2020;Tan et al. 2021).In contrast, some cloud crushing simulations find that cooling induced pulsations may be the dominant driver of turbulence (Gronke & Oh 2020a,b).Reconciling these pictures requires a careful investigation of how v turb scales with ξ.
(ii) What is the timescale for turbulent mixing?Shear layer studies associate τ mix with the eddy turnover time at the outer scale, or τ mix ∼ L outer /v turb (T cl , ℓ = L outer ), where v turb (T cl , ℓ = L outer ) is a fixed fraction of the relative velocity for χ ≳ 100 (Fielding et al. 2020;Tan et al. 2021).This scales similarly to t shear .Wind-tunnel studies instead link τ mix with the cloud-crushing time, t cc = √ χR cl /v w .The Kelvin-Helmholtz and Rayleigh-Taylor instabilities have growth times of order t cc and destroy clouds over a few t cc , in the absence of cooling (Klein et al. 1994).Gronke & Oh (2018) predicts cloud survival when ξ GO = t cc /t cool,mix exceeds unity.While both choices give ξ a v w R −1 cl scaling, the latter introduces an extra dependence on χ −1/2 .This discrepancy could have profound impacts on cloud survival criteria and requires a careful understanding of how v turb scales with χ.
To address these questions, we investigate the turbulent properties that emerge in wind-tunnel simulations of cloud-wind interactions.While turbulence in TRMLs has traditionally been treated as homogeneous (e.g Begelman & Fabian 1990;Gronke & Oh 2018;Fielding et al. 2020), we will show that it depends not just on scale but also on phase.This has important implications for mixing and hence cloud survival.Although most previous work on TRML entrainment has focused on cloud-wind density contrasts of χ = 100 − 300 (see however Sparre et al. 2020;Gronke & Oh 2018, 2020a), galactic winds are expected to have χ ≳ 10 4 (Fielding & Bryan 2022).Furthermore, we have preliminary evidence for important changes to the dynamics and clumping structure for χ ≫ 102 (Gronke & Oh 2020b).In this work we, therefore, place particular emphasis on higher χ results.
In § 2, we describe the suite of simulations used in this investigation.Videos of these simulations can be found at http://matthewabruzzo.com/visualizations/.In § 3, we describe and compare three approaches for characterizing multiphase turbulence, followed by a description of the results from applying these methods to our simulation suite in § 4. Subsequently, we describe implications of our results and detail our conclusions in § 5 and § 6.

SIMULATIONS
We ran a suite of 3D uniform grid hydrodynamical simulations using the enzo-e 2 code, which is a rewrite of enzo (Bryan et al. 2014) built on the adaptive mesh refinement framework cello (Bordner & Norman 2012, 2018).Our simulations employed the van Leer integrator (without constrained transport) (Stone & Gardiner 2009) with second order reconstruction and the HLLC Riemann solver.
Our simulations begin with a motionless spherical cloud embedded within a hot, uniform, laminar wind in the x direction.We imposed an inflow condition on the upstream boundary (positive x) and outflow conditions for the other boundaries.The cloud and wind material are initialized with p/k B = 103 K cm −3 .The cloud density ρ cl in all of our simulations is chosen such that T cl = 5010 K.This roughly corresponds to the temperature where heating starts to dominate over cooling (without self-shielding).The wind density is then determined by the desired value of χ.
To model radiative cooling, we use the grackle 3 library (Smith et al. 2017), assuming solar metallicity and no self-shielding.Specifically, we use the tabulated heating and cooling rates for optically thin gas in ionization equilibrium with the z = 0 Haardt & Madau (2012) UV background.We turn off cooling in gas with T > 0.6T w in our simulations with χ ≤ 10 3 .This helps to avoid complications from cooling in the hot wind in our χ = 100 simulations in which the ratio of the cooling time of the hot wind to the cooling time of the mixing layer is t cool,w /t cool,mix ∼ 40.In higher χ simulations cooling of the wind fluid is so slow that this ceiling has no discernible impact. 4For simplicity, we also turned off heating/cooling below T cl .
To break initial symmetries, we initialized the density of each cell, within the cloud, to the average of ρ(x), where ρ(x) ρ cl = 1 + 0.099 For each i, we drew a random unit vector êi and values for λ i and ϕ i from [R cl /8, R cl ] and [0, π).Cells on the cloud edges were initialized with subsampling; each subcell had a width of R cl /128.
Our simulations have resolutions of R cl /∆x = 4, 8, 16, 32, 64.Unless stated otherwise, results are presented for R cl /∆x = 16.By default, the wind-aligned dimension and transverse dimensions for most of our simulations' domains had sizes of 120R cl and 20R cl , respectively corresponding to a 1920 × 320 2 grid at our fiducial R cl /∆x = 16 resolution.The sizes were somewhat larger (360R cl and 30R cl ) for our χ = 10 4 simulations in order to minimize the impact of the bow shock reflections, prevent dense material from leaking out of the transverse boundaries in cases of shattering, and to give room for tail formation.While the default dimensions are adequate for determining whether our clouds survive in runs with M w ≥ 3 or χ = 10 3 , we find that boundary effects can impact later time measurements.Thus for such cases, with radiative cooling and R cl /∆x = 4, 8, 16, we present results from runs with a wind-aligned length of 240R cl .In all cases, the cloud was initialized at the center of the domain and we employed a frame-tracking scheme that updated the reference frame every t cc /16 such that the mass-weighted velocity for cells with ρ ≥ √ ρ cl ρ w was zero.
Table 1 presents a list of our simulations.As we will discuss in § 3, our measurements of v turb involve averages over velocity properties.Thus, leakage of material from the domain could plausibly bias our measurements.However, the generality of the scaling relations derived in this work, which apply to runs that do and do not leak material, suggests that overall effects on our v turb measurements are probably mini-4 We only explicitly checked the effects of a cooling wind in the χ = 100, Mw = 1.5 runs and χ = 10 3 , ξ sh = 27.8 run.We expect only minimal late time complications in our χ = 1000, ξ sh = 83.4run since it has t cool,w /tcc ∼ 51 (Abruzzo et al. 2022).While our χ = 300, ξ sh = 3.2 run has the same t cool,w /tcc, complications may be significant since that run is close to the survival threshold.
Complications are likely significant in larger χ = 300 runs.

CHARACTERIZING TURBULENCE
The primary goal of this paper is to characterize the turbulent properties of the turbulent radiative mixing layer that mediates mixing and cooling between the hot wind and cold cloud.Although much effort has been devoted to understanding turbulence in single phase media, there has been considerably less work for multiphase systems (e.g., Mohapatra et al. 2022;Gronke et al. 2022;Gronke & Oh 2022).The potential dependence of turbulent properties on both scale and the gas' local thermodynamic state complicates the interpretation of conventional methods for characterizing v turb .There are a number of possible ways to extend existing turbulence measures; however, their novel nature means that they can be difficult to interpret and their robustness is unclear.In order to get around this difficulty, in this paper we consider three distinct methods for characterizing our multiphase turbulent simulations.These are built around three different ideas based on (i) a filter-based Note-Unless otherwise noted, all runs were initialized with T cl = 5010 K, where t cool,cl /t cool,min ∼ 105.All simulations were initialized with an initial thermal pressure of p/k B = 10 3 cm −3 K.In each run, t cool is minimized at T = 1.83 × 10 4 K with a value of 75.5 kyr; the sound speed at this temperature is 18.6 km/s.The cooling length, ℓ cool = cst cool , is minimized at T = 1.70 × 10 4 K with a value of 1.43 pc.
a Denotes whether clouds survive (i.e. if the cold phase mass ever drops to 0). "Borderline" indicates cases where the line between survival vs. destruction and rapid subsequent precipitation is fuzzy b the cold phase mass, ρ > √ χρ cl , dropped to ∼0.01%, 0.07% of the initial value in the R cl /∆x = 8, 16 runs before growth.At these times, there is no mass denser than ρ cl /3.c the mass of gas denser than √ χρ cl (ρ cl /3) drops to 16% (6%) of it's original value and begins monotonic growth after 21.5tcc (24tcc).In an alternate version of the same run, where the domain dimensions are 120 × 10 2 , the mass instead drops to 32% (8%) of its initial value and starts growing after 12.5tcc (11tcc).
technique, (ii) a geometric approach, and (iii) classic structure function ideas.
We describe these approaches below, and to supplement our description of the methods, we apply each to Shows the phase dependence of v turb , measured via filtering, of the R cl = 64∆x run of our χ = 1000, ξ sh = 27.8,Mw = 1.5 simulation at 2.5tcc.The top panel includes contributions from all three velocity components.The bottom panel just includes contributions from the components transverse to vwind ; this is consistent with how v turb from filtering is measured throughout the remainder of this work.The solid orange line denotes the median while the dashed orange lines bound values between the 15th and 85th percentile.The dotted-dashed line shows v turb magnitudes that are equal to the sound-speed.The steep drop-off in v turb near Tw, is an artifact of the fact that wind is initially laminar.

Filtering
In our first approach, we attempt to explicitly remove the bulk flows by filtering out the large-scale bulk velocities.Specifically, we estimate v turb by applying a high-pass Gaussian filter, with density weighting, to each component of the velocity.The size of the filter is chosen to correspond to scales on which the bulk flow is varying, that is approximately the cloud radius.We use density weighting (which corresponds to smoothing the momentum) in order not to be dominated by the volume-filling hot gas component.
More precisely, in this approach, the i-th component of the turbulent velocity is given by where f σ (x) is the formula for a normalized, separable, three-dimensional Gaussian.In short, the convolution of f σ (x) with v i (r) (i.e. the fraction term) estimates the laminar part of v i , and subtracting if from v i gives the turbulent part.
Throughout this work, we use a Gaussian filter with a standard deviation of R cl /4; this was chosen after extensive experimentation to visually pick out turbulent regions with a minimal "bleed" into the laminar regions.Our results do not depend qualitatively on the exact choice of the filtering scale as long as it is on the order of the cloud size.
The rightmost two panels in Fig. 1, illustrate the hi-pass filtered transverse velocity components for the aforementioned simulation at 4.5t cc , and the center panel shows the combined magnitude of these components, v trans,hi .The left two panels show the density and specific thermal energy slices.Note that here, as elsewhere in this paper, we use e ≡ (p/ρ)/(γ − 1) to denote the specific thermal energy of the gas.This quantity is closely related to temperature, but is easier to compare among runs with different χ values since e wind = χe cl (due to variations in mean molecular mass T w < χT cl ).The inset panels make it readily apparent that the turbulent velocity has a clear phase dependence.
In Fig. 2 we show the phase dependence explicitly (at 2.5t cc ), plotting the 2D distribution of mass as a function of temperature and (top) hi-pass filtered velocity including all components, v tot,hi , and (bottom) just the transverse components, v trans,hi for this snapshot.Because of the spatial gradients that persist in the downstream velocity component (see Appendix A), we use v trans,hi to estimate v turb in the remainder of this work.
This approach is sensitive to turbulence on scales below the high-pass filtering limit; since this is approximately the driving scale (of order the cloud radius), we expect this to be a good measure of the turbulent properties, although it may also remove some of the contribution to the turbulence on scales just below the driving scale.One possible downside of this approach is the contribution of the bulk flow in scales at and below the filtering scales; we have explored alternate weighting schemes and find only minor differences.Although we do not have detailed scale information (except the re-moval of large scales), this approach does permit a very fine examination of the turbulent properties with phase (i.e., specific energy) as seen in Fig. 2.
Indeed, this figure clearly shows a different dependence on specific energy below and above log 10 (e/e cl ) ≈ 0.7, which corresponds to the peak of the cooling curve we adopt.We return to this point in § 4. Fig. 1 qualitatively shows that spatial variations in turbulence are largely explained by the spatial variations in gas phase.The main exception is the hottest phase, which is "contaminated" by unmixed, laminar gas (this is reflected in Fig. 2).
Because measuring phase information isn't as seamless for our other approaches, we define a set of nominal coarse phase bins to be used with them.We define the bin edges in terms of log(e/e cl )/ log χ to ease comparisons between runs with different χ values.The bin edges are −∞, 1/12, 3/12, 5/12, 7/12, 9/12, 11/12, which are illustrated by the vertical dotted lines in Fig. 2.

Geometric
Our second approach uses the geometry of isosurfaces in the flow to characterize the turbulence.To motivate this, consider a toy model in which a cloud's geometry is a sphere or a cylinder.The cloud is oriented such that the azimuthal angle, ϕ, measures the angle on the plane transverse to v wind .While cloud acceleration and accretion (e.g. by a TRML, Fielding et al. 2020) can drive steady coherent flows along the wind and radial directions, turbulence is the only source of motion along φ.In other words, we can characterize v turb with v ϕ .Tan et al. (2021) drew a similar conclusion in shearing box simulations about the utility of the dispersion of the velocity component perpendicular to shear and inflow directions.
Despite their more complex morphology, we can apply the same logic to real clouds.For a given snapshot, we employ the Lewiner et al. (2003) marching cubes algorithm to construct five topologically correct meshes of triangle facets that trace specific internal energy isosurfaces using values that coincide with the centers of the closed bins mentioned in § 3.1.We supplement these with additional isosurfaces at values near the peak of the cooling curve (we vary the precise locations based on the χ value of the simulation).Fig. 3a shows a cutaway visualization of several of these iso-surfaces at 2.5t cc for our χ = 1000, ξ sh = 27.8 simulation.
For each facet, we define v ϕ−like ≡ v • (v wind × n), where v is the linearly interpolated velocity and n is the outward normal vector.Finally, we estimate v turb for an iso-surface with the area-weighted standard deviation of v ϕ−like (excluding facets with vwind × n = 0).
Fig. 3b-c shows area-weighted distributions of v r−like ≡ v • n and v ϕ−like for the previously mentioned simulation.The distribution for v r−like shows a negative mean for each iso-surface, which is consistent with netinflow of gas.In contrast, the distributions for v ϕ−like is centered on 0, which is exactly what we expect.
While this approach does not provide any scale information about turbulence, it can be used to provide detailed phase information.For example, Fig. 3d illustrates qualitatively similar phase-dependence to the filtering measurements.However, in contrast to the filtering technique, this approach requires generation of a separate surface for each phase to be probed and so is much more computationally intensive.As is discussed in Appendix A, the main advantage of this approach is that it provides the most accurate early-time measurements.

Velocity Structure Function
Our final turbulence measure is the velocity structure function, which has the advantage of explicitly exploring the dependence on length scale ℓ, but comes with some uncertainty due to the potential influence of gradients in the large-scale bulk flows.
To compute this measure, we consider a velocity vector field that is sampled at a collection of points.Let |δv| denote the magnitude of the velocity difference between a pair of points.We define the first and second order velocity structure functions, ⟨|δv|⟩(ℓ) and ⟨(δv) 2 ⟩(ℓ), as the average values of |δv| and (δv) 2 for all pairs of points separated by a distance ℓ.Except where otherwise noted, the velocity differences only include components orthogonal to the wind direction (see Appendix A for further explanation).
Given the obvious phase dependence in our other turbulence metrics, we compute the structure function for individual phases of the gas in our simulations, using the same bins defined in § 3.1.All structure function calculations in this work are computed using all pairs of points from individual phase bins.We note that both points in each pair always comes from the same phase bin, and we leave consideration of cross-phase terms to future work.We omit the hottest phase-bin from our analysis because a large fraction is laminar (contaminating the signal) and it is computationally expensive to compute.
We also use discrete bins of ℓ, which depend on the cell width, ∆x, in our simulations.The ith ℓ bin is centered on ℓ = i∆x and has a width of ∆x.However, for i = 0 and i = 1, we have adjusted the bins such that they only Illustrates iso-temperature surfaces and derived v turb measurements for our the R cl /∆x = 64 run of our χ = 1000, ξ sh = 27.8,Mw = 1.5 simulation at 2.5tcc.Panel a shows a cut-away of five nested iso-surfaces measured at log χ e/e cl = 1/6, 1/3, 1/2, 2/3, 5/6 (for this system, T = 1.3 × 10 4 K, 3.3 × 10 4 K, 10 5 K, 3.3 × 10 5 K, 3.3 × 10 4 K, 10 6 K).The arrow illustrates the φ direction measured in the plane transverse vw.Panels b and c respectively show the normalized areaweighted distributions of the v normal and v ϕ−like velocity components measured on the iso-surfaces pictured in a. Panel d shows the standard deviation of the distributions from panel c (colored diamonds), as well as data derived from other isosurfaces (gray circles), plotted as a function of log χ e/e cl .
contain values for pairs of cells that share a face and an edge, respectively.In other words, the i = 0 bin (i = 1 bin) only contains values for cells exactly separated by ℓ = ∆x (ℓ = √ 2∆x).Throughout this work, we largely focus on ⟨(δv) 2 ⟩(ℓ) because it has a similar magnitude to our other v turb metrics.The top panel of Fig. 4 shows ⟨(δv) 2 ⟩(ℓ) measured for each phase bin of our χ = 1000, ξ sh = 27.8 simulation at 2.5t cc .The peak in ⟨(δv) 2 ⟩(ℓ) at ℓ ∼ R cl , present in all phases (in some cases it manifests as a change in slope), is expected since the outer scale should be of order the cloud size R cl ; although the complicated cloud structure at late times is unlikely to correspond to a narrow range for the injection of turbulence.We leave investigation of the behavior above R cl to future work.
We generally observe a weaker dependence on ℓ than the ∝ ℓ 1/3 scaling expected for idealized, Kolmogorov turbulence (for ⟨(δv) 2 ⟩(ℓ)), although this depends somewhat on phase.We caution that the precise ℓ scal-ing at intermediate (inertial) scales may not be a robust measurement due to the bottleneck effect, which arises for under-resolved turbulent cascade (e.g.Rennehan 2021; Mohapatra et al. 2022).
In agreement with the other measures, we also see a general decrease in the turbulent velocity with temperature.

Comparison
We have shown that much more information about the phase, scale, and spatial dependence of turbulence can be gleaned from these simulations when using metrics beyond the standard root-mean-square approach.We now compare these more refined turbulent metrics to each other.
The top row of Fig. 5 shows the v turb phase dependence, measured with each approach, in a R cl /∆x = 16 of the run of the previously mentioned simulation at t = 1.0, 5.5, 9.5, 17.5t cc (see § 4.3 for a discussion of how resolution affects our measurements).The differing approaches achieve remarkable qualitative agreement about the magnitude, phase dependence, and temporal dependence of v turb .In particular, all measures show that v turb increases rapidly with temperature at early times, before transitioning to a flatter slope at later times.In addition, all approaches show very similar amplitudes.However, the approaches are clearly not interchangeable.Indeed, this plot demonstrates the utility of computing all three turbulence measures, allowing us to ascertain the robust results without over-interpreting features that are not seen in at least two of the techniques.
When considering the volume-averaged properties of the entire system, our geometric approach offers the most robust measurements because it is most resilient to biases that may arise from the gradients in the laminar component of the flow at early times (see Appendix A for more details).
In the context of phase-dependence, the filtering approach clearly is the most convenient metric because it naturally provides turbulence as a function of phase.
However, unlike the other approaches, the filtering approach does not examine the turbulence of one phase in isolation to the others, which may introduce "artifacts" in this type of comparison.We will show in § 4.3, that the negative slope at large T , at late times may be a resolution effect.
The ⟨(δv) 2 ⟩(ℓ) approach captures much of the same phase dependence while also opening a window into the scale dependence of the turbulence.With that said, it is the most computationally expensive.

RESULTS
Having established the robustness and relative merits of our turbulence metrics, we now examine what they tell us about the cloud-wind interaction.We start (in § 4.1) by describing the phase dependence of v turb , its scaling with dimensionless parameters, and timedependence.Then, in § 4.2 we briefly discuss the driving scale before turning to an evaluation of numerical convergence in § 4.3.
For the purpose of this discussion and subsequent sections, we define the cold phase as all gas with densities of at least √ ρ cl ρ w (i.e. the density of the mixing layer).We also define the relative velocity, v rel , as the difference between v w (at the inflow boundary) and the mass-weighted velocity of the cold phase (initially v rel is v w but declines as the gas is entrained).

Turbulent Properties
Throughout this group of subsections, we will compare simulations with different parameters and ∆x = R cl /16.We first consider the phase dependence of v turb , then show how v turb scales between simulations, and finally describe the time-evolution of v turb .

Phase dependence
We begin by presenting the phase dependence in two limiting cases of the χ = 1000 and M w = 1.5 cloudwind interaction.These two cases are: (i) a run without cooling (ξ sh = 0) and (ii) a run where cooling is sufficient for entrainment (ξ sh = 27.8).
The bottom row of Fig. 5 shows the non-radiative run.In this case, the scaling of v turb with e (or T ) is consistent with a power-law where M turb is constant (i.e.v turb ∝ c s ∝ √ e) throughout the cold phase's lifetime.The amplitude of the turbulence decreases as v rel drops but there is no indication that the scaling with phase changes (note that at late times, the cold gas is entirely absent, due to mixing with the hot phase and so we cannot measure its turbulent properties).As expected, these trends are unaffected by our choice of turbulent metric.

No Cooling
Each row shows phase-dependence of v turb measured with different metrics (large panels), and bulk property evolution (smaller panels) for a χ = 1000, Mw = 1.5, ∆x = R cl /16 simulation.The top row shows our ξ sh = 27.8 simulation, (the cloud is entrained) and the bottom row shows the same simulation without cooling (the cloud is destroyed).The v turb panels respectively show data measured with high-pass filtering, from isosurfaces, and using ⟨(δv) 2 ⟩(ℓ).The data is colored by the time at which they are measured (the vertical lines in the small panels denote the times) and the dashed black line indicates values of cs as a function of phase (assuming constant pressure).Note that the lowest temperature point on the blue and green curves of the isosurface panel in the lower row, are likely outliers: the relevant isosurfaces probably bound a small amount of mass.The bulk evolution panels respectively show the mass in the cold phase, (i.e.gas denser than √ ρ cl ρw) and the average velocity difference between the cold phase and vw as functions of time.
The top row of Fig. 5 shows the run with cooling.Compared to the constant-slope power-law phase dependence of turbulence in our non-radiative run, the case with cooling clearly has more complex behavior.We parameterize the phase dependence of systems with sufficient cooling for entrainment, at a given time, using a broken power-law, Figure 6.The top row compares the dependence of v turb with gas phase at four representative times in the clouds' evolutions (full-size panels on the left) and bulk property evolution (small panels on the right) for a separate collection of simulations.The top row compares 9 simulations with χ = 100, Mw = 1.5, and the same initial cloud temperature, but varying R cl (and so varying ξ sh ).The v turb panels shows filtering measurements after 1.5tcc (leftmost panel), when the relative velocity between the cold and hot phase are various fractions of its initial value (middle 2 panels), and after 20tcc (right panel).Data is only shown for a given simulation for 0 < log χ e/e cl < 0.9.The black dashed line shows v turb = cs ∝ √ e and the vertical grey line extends between the temperatures where the cooling length is minimized and t cool is minimized.The bulk property panels show the evolution of (upper) the cold phase's mass and (lower) the relative velocity in each simulation.The curves in these panels are annotated with dots to specify the snapshots displayed in the v turb panels.It's a little ambiguous whether the cloud survives in the χ = 100, t shear ∼ 0.57 run, or if its destruction seeds the prompt precipitation of cold phase material (see Table 1 for more details).The bottom row shows the same information, but for a set of 4 simulations that all have χ = 1000.The other difference is that the rightmost v turb panel in the bottom row compares v turb measurements at v rel /vw ∼ 0.2.We note that c s,hot t cool,min is 6.56 pc (20.7 pc) for simulations in the top (bottom) row.
Below e break , the scaling of v turb on e is constant in time.Above e break , the slope of the power-law dependence, α, has clear time-dependence.At very early times (≲t cc /8), geometric measurements provide some evidence (not shown) that α = 1/2; in this case Eq. 3 is equivalent to the scaling of our non-radiative run.As the system evolves, α decreases (i.e. the slope flattens above e break ).When the cloud is mostly entrained, α stabilizes at ∼ 0.6 While we don't show it, we note that similar behavior occurs in our χ = 1000, ξ sh = 2.78 run, but the cloud is destroyed long before α drops to 0.
This demonstrates an essential feature of the turbulence in systems with sufficient cooling for cloud survival and entrainment: the cold phase has a larger turbulent  Like the top row of Fig. 6 except that the pictured simulations primarily vary the cloud temperature.Each simulation has χ = 100 and Mw = 1.5.We expect at higher resolution that the power-law slope below e break in the purple curve will be closer to 0.5 (i.e. the slope of the dashed black line).6 except that the pictured simulations primarily vary in χ.We have made two compromises in presenting this data.First, we fix β to 0.25 for all panels.This is done as a simplification because β changes on a timescale related to χ.Second, the rightmost v turb panel compares simulations at a fixed value of v rel /vw rather than at a fixed time.The last panel typically compares the simulations at a point in evolution when v turb stabilizes (see § 4.1.3).However, that time seems to come much later in our χ = 10 4 simulation, after the simulation terminates.While we include the χ = 10 4 run for completeness, strong resolution dependence (see Table 1) and the atypical shape of the cool-phase mass evolution may indicate that it is not well-converged.As noted in § 2, some material that started in the cloud leaks out of the domain at 6.5tcc, which coincides with the large drop-off in cool-phase mass.
Mach number and turbulent kinetic energy than the hot phase.
4.1.2.Scaling with Cloud Parameters (χ, Mw, and ξ sh ) Now that we've established the behavior in these limiting cases, we discuss how the principal dimensionless numbers (χ, M w , and ξ sh ) affect the magnitude of v turb in simulations with rapid enough cooling to ensure cloud survival.At a given stage of a cloud's evolution (i.e. for a given value of v rel /v w or fixed time), we find that v turb satisfies the scaling, where ξ sh and M w both refer to values used to initialize the problem.This is equivalent to saying that v turb scales with the ratio between the hot-phase soundcrossing time (R cl /c s,hot ) and t cool,min .The best fit values for β are 0.25 and ∼ 0.5 at early and late times, respectively.This change in β seems to coincide with a  6 except that the pictured simulations primarily vary Mw.The solid (dotted) lines show data from simulations with R cl /Mw = 37.6 pc (376 pc) and ξ sh = 5.73 (57.3).We note that the c s,hot t cool,min is 6.56 pc for all simulations in this plot.As we will show in panels e-h of Fig. 10, v turb evolves more slowly in higher Mw runs.Consequently, the "late times" panel shows data from Mw = 0.75, 1.5 runs at 20tcc, and data from Mw = 3, 6 runs at 30tcc (we did not run the Mw = 6 simulation to late enough times or with a long enough domain for an optimal late-time comparison).
transition between temporal evolutionary stages, which we will discuss further in the next subsection and link to a change in the primary source of turbulent kinetic energy.
Figures 6-9 compare v turb measurements, adjusted to remove differences captured by this scaling, for different sets of simulations.In other words, the agreement between the curves in a given panel in these figures indicates the accuracy of the adopted scaling.Because the principal dimensionless numbers clearly affect the slope of v turb above e break , the reader should primarily consider agreement at e break (denoted by a vertical line) and in colder gas.Note that unlike previous plots, the black dashed line shows v turb ∝ c s (e) rather than v turb = c s (e).
First, we consider the scaling for runs with M w = 1.5.Fig. 6 shows the scaling on R cl ; the top (bottom) row shows runs with χ = 100 (χ = 1000).The impressive overlap of the curves in each panel demonstrates that the adopted scaling works remarkably well -there are occasional differences at high e, but the turbulence in the gas closest to the wind phase is the most challenging to accurately measure.The figure also clearly shows that the shape of the turbulence dependence with e changes over time, a point we will return to later.
Figures 7 and 8 provides evidence that v turb depends not just on R cl , but on the ratio R cl /c s,hot by comparing runs with different R cl and c s,hot values.The variation in c s,hot come from adopting different values for T cl and χ.Fig. 8 provides further confirmation that c s,hot is the correct sound-speed to include in this scaling because c s,hot has different χ-dependence from the sound speed in the (cold) cloud phase.
Finally, Fig. 9 demonstrates the v turb scaling for runs that vary in M w .It largely confirms the lack of M w dependence.
We provide a rough normalization for Eq. 4 when v rel /v w ∼ 0.75.In this case, we find that v turb (e break ) ∼ 0.4c s,break (R cl /(t cool,min c s,hot )) 1/4 .The precise normalization will change when using other techniques to measure v turb .
All of these results are computed with the filtering metric for turbulence.We note that these scalings are somewhat less clear for geometric and ⟨(δv) 2 ⟩(ℓ) measurements of the χ = 100 simulations (the scaling between χ = 1000 runs is clear for all metrics).For example, the geometric measurements show slightly different trends among the runs that initially lose mass, and suggest that β never changes from 0.25 in the χ = 100 runs.Although the latter quirk is difficult to explain, we are encouraged by the fact that the geometric approach does show the change in β for the χ = 10 3 simulations, and the fact that the ⟨(δv) 2 ⟩(ℓ) relation definitely supports β = 0.5 at late-times in our χ = 100 simulations.With that said, ⟨(δv) 2 ⟩(ℓ = R cl ) measurements do show more scatter than is present in the top row of Fig. 6.While this could be physical, the coarser phase bins may also contribute to this scatter.We defer further investigation to future work.

Temporal evolution
So far, we have focused on how the phase-dependence of the turbulence changes (at a set of different times) Figure 10.Each panel shows the temporal evolution of v turb (blue solid curve), the average radial inflow (olive curve), the surface area (violet dashed-dotted curve), and v rel (dashed red curve).The v turb , average inflow, and surfaces areas were all computed from an isosurface constructed at the temperature at which t cool is minimized (this is e/e cl = 4.8).More specifically, the average inflow is computed from the area-weighted average of the normal velocity component on each facet on the isosurface.In contrast, v rel measures the bulk relative velocity of all gas with ρ > √ ρ cl ρw and is normalized such that it starts at unity and approaches zero as a cloud is entrained.To denote that a cloud becomes entrained (is destroyed) we include a "✓" ("×") in the panel label.The top (bottom) row show runs that have χ = 100 (χ = 1000) and Mw = 1.5 while varying ξ sh .The middle row shows runs that have χ = 100 and ξ sh = 5.73 while varying Mw.Panels c and f show the same run.As noted in the Table 1, the cloud's fate in panel b is somewhat debatable, given that the cold-phase mass drops to 0.07% of it's initial value before growing.For this case, we have elected not to show data after the cloud starts growing.
with cloud properties.We turn our attention to how v turb changes with time in a single simulation.Given how the slope of the v turb broken power-law phase dependence is largely time independent for the cold phase up to the break, we focus on v turb at e = e min,cool ∼ e break .Fig. 10 shows, for a broad range of simulations, the time evolution of v turb , measured geometrically (we use this metric to isolate a narrow phase bin), the average inflow velocity, and surface area on the same isosurface.The figure also shows the time evolution of v rel .We do not show other types of v turb measures because they are less accurate at early times (see Appendix A), and do not distinguish between turbulence and gradients in inflowing gas as well as the geometric measurements.We start by considering the turbulent evolution in a characteristic case with cloud entrainment: Fig. 10c shows a χ = 100 run with M w = 1.5 and ξ sh = 5.73.v turb has two primary evolutionary stages.Initially, in the 'pre-entrained stage', v turb rapidly grows until it reaches a peak value; v turb is sustained near this peak for a short time, and then it starts to drop off, as the cloud becomes partially entrained.During the subsequent 'partially entrained' phase, v turb stabilizes at a smaller value (within a factor of ∼2 of the peak) that is sustained for the remainder of the run.
The primary source of turbulent energy during the pre-entrained stage appears to be the relative velocity.This would explain why v turb peaks within a few t cc : we expect large v rel to drive the Kelvin-Helmholtz and Rayleigh-Taylor instabilities, which have growth rates proportional to t cc (Klein et al. 1994).This also explains similar rapid turbulent growth during the initial stage of the non-radiative and the slow cooling simulations in panels a and b of Fig. 10.Furthermore, it explains why the drop in v turb , which indicates the transition between stages (and is most prominent in the radiative runs), follows the drop in v rel -this is presumably because v rel no longer provides enough turbulent energy to sustain the peak v turb .
The two stages of v turb evolution roughly coincide with the stages of areal growth identified in Gronke & Oh (2020a).The 'pre-entrained' stage coincides with the rapid surface area growth dominated by the formation of the cloud's tail.Likewise, the 'partially entrained' stage roughly corresponds to the slower isotropic areal growth that occurs once the cloud is entrained.It's also noteworthy that the average inflow velocity plateaus before the slower isotropic areal growth, which is consistent with findings from Gronke & Oh (2020a).
At face value it might seem surprising that there is net inflow in Fig. 10b even though we know that this run is losing mass during the first ∼10t cc (see the mass evolution of the χ = 100, R cl = 5.64 pc run in Fig. 6).However, this just illustrates that net inflow doesn't necessarily equate with mass growth; the inflowing gas will raise the temperature of the gas enclosed by an isosurface in the absence of sufficient cooling.
We now consider how the principal dimensionless numbers (ξ sh , χ, M w ), affect the v turb evolution with time.In general, we find that these parameters only minimally affect the overall trend, so we focus on the relatively small differences that do emerge.
First, we examine variation in ξ sh .Compared to panel c of Fig. 10, panel d illustrates that more efficient cooling can increase the maximum v turb as well as the value of v turb at late times.This is consistent with the scalings from the last subsection.In this panel, v turb approaches v inflow 's magnitude at late times.It's plausible that all entrained runs in the figure would show the same behavior if we had measurements for late enough times; it may just be most prominent in panel d because v inflow is elevated and the cloud is accelerated more quickly.This feature may suggest that v turb is dominated by the radial flow at late times.We also find that higher ξ sh simulations have a somewhat smaller surface area.
The bottom row of Fig. 10 shows data for a set of runs with χ = 10 3 , and varying entries of ξ sh .In simulations in which the cloud survives, the transition between evolutionary stages of v turb happens at larger v rel /v w when χ is larger.This transition appears to roughly coincide with the time at which the value of β, from Eq. 4, increases from 0.25.Differences in v turb 's magnitude are qualitatively consistent with the scaling given in that equation.Finally, the middle row of Fig. 10 compares runs with varying M w .Increasing M w appears to increase v turb 's initial growth rate, v turb 's magnitude, and the duration over which v turb 's maximum magnitude is sustained.There is also some indication that higher M w simulations may also have larger inflow rates and larger surface areas, even at late times.
Independent of ξ sh , χ, and M w , Fig. 10 illustrates that the acceleration timescale is tightly correlated with the stages of areal growth (the surface area and v rel curves feature abrupt slope changes at similar times).In contrast, the transition between v turb stages appears less tightly coupled with the acceleration timescale as the principal dimensionless numbers are changed.We attribute this mostly to the fact that the cold phase is not a rigid body with a single bulk velocity, but instead has different velocities at different spatial locations.
This differential acceleration is responsible for the cloud's head-tail morphology: downstream material moves faster than upstream material.Regions with larger v rel (compared to v w ) should generally have larger v turb , albeit with some scatter related to the local history of turbulent driving.This is illustrated for a high resolution version (to improve sampling) of our χ = 1000, ξ sh = 27.8 run in Fig. 11.Here, we explore the relation between our measured turbulence metric (at the cooling peak) and the relative velocity of the gas as a function of both time (colors) and location along the length of the cloud (different points with the same color).This demonstrates that there is a correlation between these quantities not just at different times for the whole cloud (as shown in Fig. 10), but also along a cloud at a given time, strengthening the case for a causative relation.
How does this relate back to the loose coupling seen between the v turb evolutionary stages and the acceleration timescale, when we vary the principal dimensionless numbers?Because entrained clouds in our various simulations have different wind-aligned lengths, we know that changes in these numbers alter the cloud's differential acceleration.Consider the temporal evolution of the volume-averaged v turb measurements for a narrow phase bin of a very coherently accelerated cloud and a less coherently accelerated cloud.One would naturally expect that that v turb measurements might spend more time near its maximum value in one of these cases.It's not much of a stretch to assume that v rel might be fairly different when v turb starts to decrease (i.e.begins transitioning between stages).Thus, we would find different coupling between v rel 's evolution and v turb 's evolution in these cases.

Evolution of the driving scale
We now briefly revisit the velocity structure function in order to investigate how the turbulent driving scale varies with time.The bottom panel of Fig. 4 shows ⟨(δv) 2 ⟩(ℓ) for the R cl /∆x = 64 run of our χ = 1000, M w = 1.5, ξ sh = 27.8 simulation when the cloud is mostly entrained in the wind (v rel /v w ∼ 0.27).Comparisons with the top panel (v rel /v w ∼ 0.94) reveal that the outer scale of turbulent driving, which coincides with the peak ⟨(δv) 2 ⟩(ℓ), does not change substantially from early to late times.Although we don't show it, we confirmed similar behavior in the R cl /∆x = 32 run of our χ = 100, M w = 1.5, ξ sh = 5.73 simulations for similar values of v rel and at times when the cloud is more entrained.
We note that it's unclear why the 9/12 ≤ e/e cl < 11/12 phase bin's ⟨(δv) 2 ⟩(ℓ ∼ 0.3R cl ) measurement, from the lower panel, is smaller than comparable measurements for other phase bins.This feature also appears in the R cl /∆x = 32 version of this simulation.In contrast, this feature is absent from the aforementioned χ = 100 run; in that case ⟨(δv) 2 ⟩(ℓ) is always larger for a given ℓ in hotter gas.

Convergence
In this section, we discuss how numerical resolution impacts our various measurements.We primarily compare the measurements among different resolution runs of our χ = 1000, ξ sh = 27.8 simulation, varying R cl /∆x from 4 to 64.

Turbulence Metrics
The large panels in the top row of Fig. 12 compare the phase dependence of v turb using filtering measurements of v turb at various points in the cloud's lifetime.The figure shows that resolution appears to slightly affect the magnitude and the slope of the v turb phase dependence above e break .Importantly, the figure also suggests that the occurrence of a negative slope of the phase dependence is likely a resolution effect.The full phase dependence of v turb is well converged for R cl /∆x ≥ 32.These same conclusions apply to our other turbulence metrics (shown in the other rows).
Fig. 13 shows how resolution affects the temporal evolution of various quantities.The top panels show convergence in the total cold phase mass7 and v rel ; the only noteworthy feature is that rapid growth begins slightly sooner at higher resolutions.However, the surface area measurements are not converged at all; it increases more rapidly for higher resolution runs.These results are consistent with the findings of Gronke & Oh (2020a) for a χ = 10 simulation.
There are some differences in the v turb evolution.While low resolution runs have a strong, sharp peak followed by a flat region, higher resolution runs have a moderate peak with a gradual descent.With that said, there seems to be convergence for R cl /∆x ≳ 16, and all of the runs qualitatively agree with our picture that there are two stages of v turb evolution.The average inflow velocity measurements are similar overall but do show some significant differences -its somewhat unclear what the relevant trends are.We defer further investigation of inflow velocity convergence to future work.

Phase structure
Resolution strongly affects the 2D thermodynamic e− p phase-space distribution.Previous work (e.g.Fielding et al. 2020;Abruzzo et al. 2022)   The top row compares the v turb phase dependence of different resolution runs of a simulation at 4 selected stages of evolution (full-size) and bulk property evolution (small panels) of each simulation.The illustrated simulations all have χ = 10 3 , ξ sh = 27.8,Mw = 1.5.The v turb panels shows filtering measurements after 1.5tcc (left panel) and when the relative velocity between the cold and hot phase are various fractions of its initial value (other panels).Data is only shown for a given simulation for 0 < log χ e/e cl < 0.9.The black dashed line shows v turb = cs ∝ √ e and the vertical grey line extends between the temperatures where the cooling length is minimized and t cool is minimized.The bulk property panels respectively show the evolution of (top) the cold phase's mass and (bottom) the relative velocity in each simulation.The curves in these panels are annotated with dots to specify the values during the snapshots displayed in the v turb panels.The second and third rows of v turb panels are the same as the top row, except that they respectively display geometric and ⟨(δv) 2 ⟩(ℓ = R cl ) measurements.The R cl /∆x = 32, 64 runs make use of a smaller simulation box than the other displayed runs.While visual inspection of our simulations, lead us to expect boundary effects to bias measurements in those cases, when v rel /vw ∼ 0.25, this doesn't seem to be an issue for this exercise.
isobar that is bounded by the properties of the cloud and wind.However, a pressure decrement emerges in the phase distribution at points along this isobar where cooling is not resolved.
Each point in the internal energy-pressure (e-p) phase space has an associated cooling length-scale c s t cool .Fig. 14 shows that the size of the pressure decrement scales inversely with how well c s t cool is resolved.The figure also suggests that resolving the minimum cooling length scale (i.e., ∆x ≲ ℓ cool ∼ min(c s t cool )), which is equivalent to the "shattering" length scale (McCourt et al. 2018), is adequate to largely remove the pressure decrement for χ ∼ 100, which is consistent with results from prior works (e.g.Abruzzo et al. 2022).
Under-resolved cooling is not the sole reason for the gas distribution's deviations from the pressure isobar.Ji et al. (2019) previously argued that it is actually the sum of the turbulent pressure, ρv 2 turb , and thermal pres-sure that should match the external pressure.Fig. 14 illustrates the median turbulent pressure as a function of e with dashed lines.In both χ cases, the turbulent pressure shows clear convergence in our higher resolution runs.The turbulent pressure's lack of dependence on e below e break and inverse correlation with e above e break (at the pictured time) are consistent with the scaling described in Eq. 3. The factor of ∼3 difference in the maximum values (i.e. at e ∼ e break ) of the turbulent pressures between the two χ cases helps explain why the χ = 1000 case has larger deviations in the thermal pressure from the external pressure.This difference is consistent with the scaling from Eq. 4. For context, we expect v 2 turb at e break to be a factor of ∼4.85 2β larger in this χ = 1000 case, although the value of β is ambiguous; Fig. 6 suggests that these particular χ = 1000 and χ = 100 runs should have β = 0.5 and β = 0.25 at 11.5t cc .
Figure 14.Solid colored lines show the median thermal pressure as a function of temperature for multiple resolutions of our Mw = 1.5 simulations with χ = 100, ξ sh = 5.73 (left) and χ = 10 3 , ξ sh = 27.3 (right) at 11.5tcc.The background shaded region shows the size of the cooling length associated with a point in phase-space measured relative to the simulation's R cl .There is not an associated length scale below e cl or above ∼ew because we have disabled cooling and heating at these temperatures.At low pressures, just above e cl there isn't an associated length scale because heating dominates.For the sake of comparison, the dashed lines show the median turbulent pressure, ρv 2 turb (derived from the filtering approach).
At the finite resolutions of our simulations, there is a decrement in the total pressure in our simulations.However, at infinite resolution it is plausible that the total pressure of the gas is constant.In short, the minimum c s t cool along the segment of the pressure isobar, connecting the cloud and the wind phase properties, specifies the grid-scale requirement for fully resolving phasestructure.Remarkably, the degree to which we resolve ℓ cool appears to have minimal impact on the 1D e phase distribution.This is shown in Fig. 15b (we will discuss the rest of the figure in the next section).

Turbulent structure and Cloud Morphology
Our results hint that under-resolving turbulence might influence various properties of these interactions.To illustrate this, we turn to Fig. 15, which shows measurements taken from different resolution runs of our χ = 10 3 , ξ sh = 27.8 simulation at 7.5t cc .Each panel in the top rows displays lines of data that comes from each  resolution.Subsequent rows just show measurements taken from a single resolution.
Fig. 15a shows the first-order velocity structure function, ⟨|δv|⟩(ℓ), measured for gas in the 1/12 ≤ log χ e/e cl < 3/12 (8.1 × 10 3 K ≤ T ≤ 2.0 × 10 4 K) phase bin at various resolutions.Values are divided by the bin's maximum sound speed and the gray region denotes the width of the bin.⟨|δv|⟩(ℓ) specifies the average magnitude of the velocity differences8 for pairs of points separated by a length-scale ℓ.
As the separation ℓ decreases so does the velocity difference.On scales comparable to the cloud radius, ℓ ∼ R cl , the slope and normalization of ⟨|δv|⟩(ℓ) are re-markably well converged.9As the separation approaches the grid scale the velocity differences are damped by numerical dissipation.Where this numerical dissipation kicks in relative to the sound speed appears to have a major impact on the morphology of the system.In reality the true physical viscosity of these systems is uncertain, but is likely to be much less than the effective numerical viscosity even in our highest resolution simulation.
On large scales the velocity differences are greater than the sound speed, but at small enough separations the velocity differences become subsonic.We define the turbulent sonic length, ℓ turb,sonic , as the scale at which ⟨|δv|⟩(ℓ) passes through the point ⟨|δv|⟩(ℓ turb,sonic ) = c s .By extrapolating the slope from large separations we can estimate ℓ turb,sonic in the limit of infinite resolution (and very small viscosity), which, in this case, falls around 0.07R cl .This ℓ turb,sonic is not resolved by the simulations with R cl /∆x = 4 or 8, is marginally resolved by the R cl /∆x = 16 simulation, and is fairly well resolved by the R cl /∆x = 32 and 64 simulations.When ∆x > ℓ turb,sonic , the average velocity difference, in a given phase bin, can be supersonic at the viscous scale (i.e. between adjacent cells).
Panels d, f, h, j, and l show the distribution of velocity difference magnitudes, in the previously mentioned phase-bin, measured at ℓ = ∆x; the average values of these distributions give the leftmost points of the curves in panel a.These panels illustrate that as ℓ turb,sonic /∆x decreases, fewer pairs of cells have supersonic velocity differences.They also show that some grid-scale supersonic velocity differences persist when ℓ turb,sonic is barely resolved.
We now investigate the question: What are the consequences of under-resolving ℓ turb,sonic ?Panels e, g, i, k, and m of Fig. 15 show the projected density of these simulations.The dramatic differences in these maps suggest that the degree to which ℓ turb,sonic is resolved may be linked to morphological differences between simulations.We find that the cold phase in higher resolution simulations is composed of more large-scale structures and has a narrower transverse extent, whereas in lower resolution simulations the cold phase is clumpier and more dispersed.The cold phase in the low resolution simulations has effectively shattered while in the higher resolution simulations that have ℓ turb,sonic /∆x > 1 the cold phase remains more intact (McCourt et al. 2018).These effects support a picture in which under-resolving turbulence intensifies shattering by enabling the presence of supersonic velocity differences on the grid scale.This also naturally explain the wider dispersal of cold gas in low resolution simulations since the most intense shattering cause explosive breakup of clouds (Gronke & Oh 2020b).Physically, supersonic grid-scale velocity differences will lead to large pressure imbalances that will in turn promote the dispersal as opposed to coagulation of cold cloudlets (Gronke & Oh 2022).
Using Eq. 4, which captures how v turb scales with v rel /v w and t cool,min , we can write a rough scaling relation for ℓ turb,sonic .Assuming that ⟨|δv|⟩(ℓ) ∝ ℓ ζ , we find for cold gas with e cl ≤ e ≤ e break that For sake of convenience we take ζ = 1/3 (the scaling for Kolmogorov turbulence), which is close to what is found in the simulations (see Fig. 15a).At early times in the 'pre-entrained' stage (e.g., when v rel /v w ∼ 0.75) β = 1/4, which yields a precise prediction for the turbulent sonic length (6) The normalization is measured empirically in our χ = 1000, ξ sh = 27.8 simulation.Note that we focus on the ⟨|δv|⟩(ℓ) measurements from the same phase-bin that includes e break since, as we have shown above, this is the region of phase space where these scalings are robust, but the general trends will be the same for other bins below e break .Due to the fact that we have only measured the length where ⟨|δv|⟩(ℓ) is equal to the maximum sound-speed of the phase bin, this relation should be considered an upper-limit on ℓ turb,sonic .
It's more intuitive to compare this relation against other known length-scales, like the minimum radius for cloud survival, R cl,crit or the minimum cooling length.
For fixed cloud properties, we find that This demonstrates that the turbulent sonic length tends to be more difficult to resolve in runs with larger χ and 10 This assumes that R cl,crit has the scaling from the Li et al.
(2020)/ Sparre et al. (2020) survival criterion, since this does an accurate job predicting cloud survival (see § 5.5).Survival criteria have the generic form, τ cool /tcc < q and thus R cl,crit ∝ τ cool Mw/q.In this case, τ cool ∼ t cool,w and q ∝ R 0.3 cl M −2.5 w n 0.3 w v 0.6 w , or equivalently q ∝ M −1.9 w R 0.3 cl p 0.3 µ 0.3 w .For T ≳ 10 5 K, t cool roughly scales as e 2.7 p −1 and the mean molecular weight, µ, is constant.Putting this together yields R 1.3  cl,crit ∝ χ 2.7 e 2.7 cl M 2.9 w p −1.3 c s,cold when Tw ≳ 10 5 K.
higher M w .If we assume that ℓ cool = min(c s t cool ) ∼ c s,break t cool,min and c s,break ∼ c s,cold , we find that ℓ turb,sonic /ℓ cool ∼ ξ ulations, it should be clear that ℓ turb,sonic exceeds ℓ cool in all of our entrained runs.We find that this relation reproduces the value of ℓ turb,sonic measured from the R cl /∆x ∼ 32 run of our χ ∼ 100, M w = 1.5, ξ sh = 5.73 simulation, to within ∼50%.The lower resolution runs of that simulation all resolve ℓ turb,sonic , and we are encouraged that none of them shows signs of shattering (the transverse extent is fairly consistent among runs).In the R cl /∆x ∼ 16 run of our, χ ∼ 10 4 , ξ sh = 172.8simulation, we find that ℓ turb,sonic is smaller than the grid scale, when v rel /v w ∼ 0.75, which is consistent with the relation's prediction.We note that both resolution runs of this simulation clearly shatter.We performed a few spotchecks with a handful of our other runs and the relation seems accurate to within a factor of a few, but more careful modeling is required since ℓ turb,sonic is close to R cl /8 and R cl /16 in many of our runs.
While we primarily presented this analysis for the phase bin containing e break , we also find evidence (not shown) suggesting that the general results can be extrapolated to lower temperature bins.This is intuitive, given our earlier finding that v turb (e)/c s (e) is roughly constant for e cl ≤ e ≤ e break .
Although our association of these large-scale morphological changes with ℓ turb,sonic requires further investigation, it presents an attractive way to understand several outstanding related questions, namely, when do clouds shatter (Gronke & Oh 2020b), and why do higher M w simulations require so much higher resolution to achieve convergence (Gronke & Oh 2020a; Bustard & Gronke 2022, we elaborate further in section 5.7).Although this need not be true in the general case (e.g. if there are external drivers of turbulence), ℓ turb,sonic > ℓ cool in all of our runs.Thus, the resolution effects on large-scale morphology may be more closely related to under-resolved cooling.
We clarify that resolving small scale structure (e.g.surface area and number of clumps) has other conditions unrelated to resolving ℓ turb,sonic .Sparre et al. (2019) and Gronke & Oh (2020a) each show that convergence of such properties is very weak in high resolution simulations that resolve ℓ cool (in both studies, ℓ turb,sonic > ℓ cool ).

Phase Dependence of turbulence
We have demonstrated for the first time that the turbulent velocity, v turb , in a mixing layer follows a bro-ken power law dependence on temperature or internal energy.A major implication of this finding is that the turbulent kinetic energy density is not constant across gas phase.Consider the ratio of the turbulent kinetic energy densities in the hot and cold phases, or ϵ = ρ w v turb (T w ) 2 /(ρ cl v turb (T cl ) 2 ).Per Eq. 3, this evaluates to ϵ = (e break /e w ) 1−2α .We remind the reader that α, the power-law slope above e break , starts out near 1/2 at the earliest times and decreases to ∼0 at a rate that depends on the principal dimensionless numbers.Thus, during the bulk of the cloud-wind interaction, the cold phase has a larger turbulent kinetic energy density (i.e.ϵ < 1).This contradicts (explicit and implicit) assumptions that ϵ = 1 in multiple works on TRMLs.
For example, we consider the arguments that lead to the expression for the temperature of the mixing layer, T mix ∼ √ T cl T w (Begelman & Fabian 1990;Gronke & Oh 2018).This relation derives from the average of the cold and hot phase temperatures, weighted by the mass flux from each phase into the mixing layer.The derivation assumes that each phase's mass flux scales with the respective v turb values.Because the derivation involves arguments equivalent to assuming ϵ = 1, it overestimates the hot phase's v turb and consequently the mass flux when compared against the values for the cold phase.Thus, √ T cl T w overestimates T mix and the size of the discrepancy is inversely correlated with ϵ.Because the value of t cool is commonly monotonic between t cool,min and t cool ( √ T cl T w ) (e.g.see figure 14 of Abruzzo et al. 2022), typical calculations overestimate t cool,mix by an amount also negatively correlated with ϵ.
In another case, Fielding et al. (2020) explicitly assumes that ϵ = 1.The only practical implication is that their quoted measurement of f turb = v turb /v rel is too large by a factor of √ ϵ.Thus, f turb might have a weak dependence on the shape of the cooling curve.In their analysis of clouds in a turbulent medium, Gronke et al. (2022) also assumes ϵ = 1, but this may be valid since they consider externally driven turbulence.

Observable Predictions
It may be possible to observe v turb 's broken power-law phase dependence in real-world systems.For example, previous studies have already placed constraints on temperature and nonthermal motion in the circumgalactic medium of other galaxies by measuring the widths of absorption lines for elements with different atomic masses (e.g.Rudie et al. 2019;Qu et al. 2022).Similar measurements may also be possible for high velocity clouds, for which there an abundance of absorption (e.g.Fox et al. 2004) and emission line data (e.g.Tufte et al. 1998;Hill et al. 2009).One could also imagine using 21 cm emis-sion or Mgii absorption to extend such an analysis to probe the turbulent properties down to lower temperatures, where gas is atomic (e.g.Marchal et al. 2021).
It may also be possible to perform a similar exercise for gas in multiphase galactic outflows (e.g.Strickland & Heckman 2009;Reichardt Chu et al. 2022).
Additionally, one can perform more straight-forward comparisons against observational measurements of turbulent measurements in ∼10 4 K gas.However, given the simplifying assumptions in this work (described further in § 5.7) and the fact the drivers of turbulence may vary between different systems, such comparisons must be interpreted with great caution.Nevertheless, we find it encouraging that there is evidence that the Perseus molecular cloud has transonic turbulent Mach number (Burkhart et al. 2015), just like we see in a fair number of our simulations.We also find it encouraging that studies of CGM clouds (e.g Rudie et al. 2019;Qu et al. 2022) recover non-thermal broadening measurements within a factor of a few of 10 km s −1 , which nicely matches the turbulent velocities in our simulations.We leave further comparisons to future work.

What drives mixing?
We now return to one of the motivating questions, the origin of turbulence in the flow.From the results in this paper, the short answer appears to be that both shear and cooling drive the turbulence responsible for mixing.As we conclude in § 4.1.1,shear is the primary driver of turbulence at early times.After the cloud becomes partially entrained, v turb falls off before stabilizing at a lower value.The long-term support of a non-zero v turb value, as v rel goes to zero, suggests that some form of "cooling-induced mixing" mechanism takes over.To put this another way, the primary source of turbulent kinetic energy changes with time.At early times, turbulent kinetic energy primarily comes from the large relative shear velocities between fluid elements.At late times, it instead comes from the radial kinetic energy of inflowing material.
Possible origins for the late-time turbulence include rapid cooling driven pulsations in the cloud, (Gronke & Oh 2020a) 11 , or simply the net radial inflow driven by the initial shear-driven turbulence.This later explanation is supported by the correlation of v turb 's late-time magnitude with v inflow , which itself correlates with a run's cooling efficiency.We plan to provide a detailed analysis of the temporal evolution of v turb and its dependence on v rel in a follow up work.
A few other features are consistent with this conclusion.First, we find the rapid growth of surface area, when shear primarily drives mixing, and subsequent stabilization at a roughly constant value, when mixing is primarily driven by pulsations or radial inflow, to be consistent.Second, the minimal variance in the driving scale, as the cloud is elongated, is also consistent.At early times the driving scale is linked with the length of the wind-aligned axis of the cloud, of order R cl .Because the cloud's transverse extent doesn't change much with time, the typical radial separation between opposite inflow 'fronts' of the clouds should still be of order R cl at late times.Finally, the saturation of the inflow velocity after cooling-driven mixing has fully developed fits into this picture since the shear-driven contribution will have become subdominant.
Gronke & Oh (2020a) noted that the anti-correlation between the cold cloud mass growth rate and v rel might suggest that shear-driven turbulence from the KH instability might not fuel mass growth, and instead might be a competing destructive process.However, our most efficiently cooling M w = 1.5 runs with χ = 100, 300, 1000 have significant v rel when they start monotonically growing.In other words, mass growth at early times in these runs should primarily arise from shear-driven turbulence.With that said, mass growth is still negatively correlated with v rel since the surface area is still increasing.
The evolution in the v turb phase dependence is also consistent with this picture.When shear primarily drives turbulence at early times, turbulent kinetic energy is roughly constant with phase (as in non-radiative simulations where shear is the only turbulent driver).In contrast, when cooling drives turbulence, it does so primarily in regions with short cooling times, which explains why turbulence in the hot phase drops off.

What is the mixing timescale?
The canonical estimates for the characteristic mixing timescale are t cc and t shear .We find that the turbulent velocity scales with R β cl c −β s,hot t −β cool,min , where β is 0.25 at early times and 0.5 at late times.Notably, it has no dependence on M w for most of the cloud's evolution.Therefore, the characteristic mixing time has no v rel dependence.
With that said, the initial value of M w does affect the temporal evolution of v turb .Fig. 10 also provides some indications that the magnitude of v turb may have some dependence on M w at very early times.Comparing panels g to h (as well as f to g) reveal that the peak values of v turb , when v rel > 0.8, is larger in the higher M w run by more than the factor of √ 2 expected by Eq. 4 from differences in R cl .

Survival Criterion
There has been great interest in the literature about the minimum radius for cloud survival (e.g.Gronke & Oh 2018;Li et al. 2020;Sparre et al. 2020;Kanjilal et al. 2021;Abruzzo et al. 2022;Farber & Gronke 2022).We will provide more firm conclusions about this topic in an upcoming work (Abruzzo et al., in prep.).However, we do note that the our results are most consistent with the predictions of Li et al. (2020) with the corrections described by Sparre et al. (2020) for supersonic winds.

Convergence
What does it mean to resolve the cloud-wind interaction?The obvious ideal is to achieve point-wise convergence, but this is generally prohibitively computationally expensive except in rare cases (e.g., Lecoanet et al. 2016).Short of this ultimate goal there are lesser gradations of convergence that depend on the question at hand.The easiest quantity to achieve convergence in is the net mass growth of the cold phase.We show in Fig. 15c that the mass growth is fairly well converged for resolutions of R cl /∆x ≳ 8.This likely corresponds to some minimum threshold to resolve any turbulent mixing, and is consistent with previous findings (e.g.Gronke & Oh 2020a).The hardest quantity to achieve converge in is the 2d p − e phase distribution, which requires resolving the minimum cooling length (ℓ cool ; also known as the shattering length).Therefore, if one is interested in simply capturing the total amount of mass in the cold phase then the resolution requirements are much less onerous than if one is interested in capturing the detailed phase structure (or cloud morphology).The details of the phase structure can be extremely important for comparisons to observations since the pressure decrement that develops in under-resolved simulations occurs in precisely the region traced by commonly observed ions, such as Mgii (e.g., Nelson et al. 2021;Burchett et al. 2021).
Here we propose an intermediate convergence criterion for the large-scale morphology of cold structures which requires resolving the turbulent sonic length ℓ turb,sonic by several cells.This is in general less stringent than the requirement to resolve the minimum cooling length.At face value, the difficulty of resolving ℓ turb,sonic in galaxyscale simulations suggests that the detailed morphological properties of cool (∼10 4 K) gas, involved in TRML entrainment, within galactic outflows and the circumgalactic medium are unlikely to be correct.However, the implications of accurately capturing the morphol-ogy may be more complex in more realistic systems because of the way cloud shape and size couples to other physical process absent in our simulations.For example, in systems in which the hot phase is itself turbulent, such as in galactic wind simulations (e.g.Schneider et al. 2020), under-resolving ℓ turb,sonic may lead to artificially shattered clouds which will in turn be more likely to be destroyed than if they were able to remain coherent.Therefore, having ∆x < ℓ turb,sonic may prove to be essential for determining the overall phase structure and evolution of turbulent multiphase flows that are ubiquitous in and around galaxies.
This discussion about large-scale morphological convergence of cool gas in larger-scale models deserves elaboration on two finer points.First, it assumes applicability of our results about the emergent turbulent properties in the cloud-wind interactions; we discuss how the equilibrium T cl and shape of t cool (T ) affect this in the next subsection ( § 5.7).Second, we are extrapolating from simulations of isolated clouds, whereas larger-scale models often include multiple clouds in an outflow (e.g.Cooper et al. 2008;Kim & Ostriker 2018;Schneider et al. 2020).This is not an issue when the inter-cloud spacing is large enough for clouds to be treated individually, albeit with a hot phase that is already turbulent from upstream interactions.However, more work is required to make predictions when the inter-cloud separation is small (such work might use a multi-cloud setup akin to Alūzas et al. 2012;Banda-Barragán et al. 2020).

Comparison to prior work
At early times, when the KH instability is the primary driver of mixing, one might expect similarities between our runs and the TRML simulations of Fielding et al. (2020) and Tan et al. (2021).Unfortunately, it's difficult to draw direct comparisons since those works highlight properties after reaching a quasi-steady state.In contrast, our runs never reach such a state since v rel evolves with time.More meaningful comparisons could be made if the cloud was in a potential that was tuned to maintain v rel at late times.Additionally, Tan et al. (2021) point out that we would likely expect different v turb scaling to be dependent on geometry.Nevertheless, we find the presence of inflowing gas at early times to be encouraging (especially when juxtaposed with our adiabatic runs that don't have net inflow).The fact that v turb and the inflow velocity show signs of scaling with cooling efficiency is also encouraging.
Likewise, we expect similarities with Gronke & Oh (2020a) at late times when turbulence is driven by "cooling-induced mixing" .Although we broadly see similar qualitative evolution in the surface area, detailed comparisons of other properties are challenging.While both works measured v inflow , we expect differences in our methodologies will complicate comparisons of these quantities at late times.Gronke & Oh (2020a) used v inflow ∼ ṁcold /(Aρ w ) while we directly measure the velocity component normal to the e mix isosurface (the scaling doesn't change much if we use the e break isosurface).In other words, their measurements are weighted by mass flux and ours are weighted by surface area.We expect that this difference in methodology explains why our results indicate that inflow starts much earlier in our runs; early time inflow that doesn't correspond to mass growth won't be picked up by their measurements.Because our work focused on measuring v turb , rather than v inflow , we defer detailed scaling of v inflow to followup work.
Gronke & Oh (2020a) found that cold phase mass evolution's convergence in a M w = 6 simulation run at R cl /∆x = 8, 32 to be quite poor.In contrast we found that the cold phase mass evolution in our R cl /∆x = 8, 16 runs of our M w = 6 simulation to be fairly well converged.While it's possible that we could see differences at higher resolution, it's plausible this difference could arise from differences in the cloud temperature.The clouds in Gronke & Oh (2020a) had a temperature of T cl = 4 × 10 4 K.This translates to values of t cool,min and e cl that are factors of ∼ 5 and ∼ 11.7 larger.Consequently, we expect ℓ turb,sonic /R cl,crit to be 7.3 times smaller in their simulations, which means they could be under-resolving R cl according to our new resolution criterion.
More generally, one might ask "How does the choice of T cl affect our results?"given that the equilibrium T cl varies12 greatly among cloud-crushing and galactic outflow studies.For context, this work focuses on runs with T cl ∼ 5 × 10 3 K, while other works commonly include simulations with T cl ∼ 10 4 K (e.g.Li et al. 2020;Kanjilal et al. 2021;Abruzzo et al. 2022;Schneider et al. 2020) or T cl ∼ 4 × 10 4 K (e.g.Gronke & Oh 2018, 2020a;Abruzzo et al. 2022).We expect the applicability of our results is more-strongly tied to the shape of t cool (T ) over T cl ≲ T ≲ T w than the precise value of T cl .Fig. 7 suggests our results are minimally affected when t cool (T cl ) exceeds the minimum value of t cool computed over the temperature range.However, the applicability is less clear when t cool (T ) is minimized at T cl (i.e. if T cl ≳ 2 × 10 4 K for p/k B = 10 3 K cm −3 , Z ⊙ , z = 0) or at a value of T exceeding √ T cl T w .Finally, we note that some works also consider conditions with T cl < 500 K (e.g.Banda-Barragán et al. 2021;Farber & Gronke 2022).Further investigation is required to understand the applicability of our results in this context, but our above discussion about t cool (T )'s shape is relevant.
We next draw comparisons with works that studied multiphase gas in turbulent box simulations.For example, Gronke et al. (2022) initialized a pressure-confined cool (T cl = 4×10 4 K) cloud in a hot ambient background and studied how the system evolved while driving turbulence in the hot phase.Mohapatra et al. (2022) studied the turbulent properties of multiphase gas (comparable to ICM conditions) that emerged from driven turbulence and radiative cooling in a box of initially hot (T = 4 × 10 6 K) gas.These studies respectively observed that the amplitude of the first and second order velocity structure functions (⟨|δv|⟩(ℓ) and ⟨(δv) 2 ⟩(ℓ)) have lower amplitudes in the cold-phase gas than in the other phases, which is in good qualitative agreement with our results.We note that the sub-Kolmogorv scaling of our ⟨(δv) 2 ⟩(ℓ) measurements are more consistent with the hydrodynamic volume-weighted heating run from Mohapatra et al. (2022) than the mass-weighted run.However, as mentioned in § 3.3, the driving scale is not sufficiently resolved to remove the bottleneck effect's influence on the slope of ⟨(δv) 2 ⟩(ℓ).To be more concrete, we note that Mohapatra et al. (2022) illustrated that the driving scale must be resolved by more than 192 cells, in a non-radiative turbulence simulation, to remove the bottleneck effect's influence on the slope.For that reason, we refrain from making detailed comparisons.

Caveats
This work made a number of simplifying assumptions and omitted a variety of potentially relevant physical effects that could potentially modify our results.Future work should consider: Other sources of turbulence: We only analyzed the turbulence that emerged from two phases that initially had coherent velocities without turbulence.In reality, external processes, like supernovae, can drive turbulence in the wind; this likely alters the interaction's evolution and makes survival more difficult (e.g.Schneider et al. 2020).Additionally, differences in the initial cloud structure, due to turbulent driving before encountering a wind, can affect the rates at which mixing destroy clouds (e.g.Schneider & Robertson 2017;Banda-Barragán et al. 2019).
Thermal Conduction: The omission of thermal conduction from our simulations will certainly affect the morphology of the cold-phase (e.g.Brüggen & Scanna-pieco 2016;Li et al. 2020).However, we take solace in the fact that mass transfer through the TRML will be minimally affected in simulations where cooling is fast relative to the mixing time (Tan et al. 2021).
Magnetic fields: It is well known that magnetic fields can extend the lifetime of clouds (e.g.Dursi & Pfrommer 2008;McCourt et al. 2015).Banda-Barragán et al. (2018) showed that magnetic fields have a stabilizing effect on initially turbulent clouds embedded in a laminar wind.While realistic magnetic field strengths don't seem to strongly affect the criteria for survival through rapid cooling, they do have a number of other effects that will almost certainly affect the system's turbulent properties (Gronke & Oh 2020a).Among others, such effects include non-thermal support, which could alter cooling properties, suppression of the KH instability and alteration of cloud morphology, leading to higher surface areas (Gronke & Oh 2020a).
Cosmic Rays: Cosmic rays were also omitted from our simulations.They are a known sources of non-thermal pressure support, which may alter cooling properties (Butsky et al. 2020).They can also provide another mechanism for accelerating clouds (Wiener et al. 2019;Huang et al. 2022).
Gravity: Our simulations neglected the effects of gravity because we generally expect our χ ≤ 10 3 simulations to be Jeans stable.However, one could imagine that external gravitational fields could sustain an elevated shear velocity (Tan et al. 2023) and consequently influence the system's turbulent properties.
More realistic cooling: All of our simulations assume simplified equilibrium cooling and neglect self-shielding.However, given that all our simulations where the cloud survives have N Hi > 10 17.2 cm −2 , self-shielding may be relevant.Including more realistic cooling could modify our results (Farber & Gronke 2022), but we leave that for future work.
Viscosity: Our simulations do not have explicit viscosity (Li et al. 2020;Jennings & Li 2021).This may affect turbulent properties near the scale of turbulent dissipation.

CONCLUSION
We have investigated the multiphase turbulent properties that emerge from interactions between cool clouds and hot supersonic flows (or winds).The relative efficiency of turbulent mixing and radiative cooling in mixing layers govern the outcome of such interactions.To address difficulties associated with characterizing multiphase turbulence, our analysis employed three distinct methods to measure v turb .We found the following primary results for simulations, in which cooling is suffi-cient for the cloud to survive the interaction and become entrained: • Radiative cooling dramatically changes the v turb temperature13 scaling.In non-radiative simulations v turb has a scaling consistent with the sound speed's temperature scaling: v turb ∝ c s ∝ √ T .In runs with sufficient cooling for entrainment, this scaling only applies for gas colder than T break , the temperature where t cool is minimized.Above T break , the power-law slope starts near 0.5 and flattens to ∼0.Consequently, cold gas generally has larger turbulent Mach number and turbulent kinetic energy than hot gas.
• v turb has two stages of temporal evolution.
The shear velocity initially drives rapid growth of v turb at early times in the "pre-entrained" phase.
As the cloud becomes partially entrained, v turb drops off before stabilizing at a lower value, one that is of comparable magnitude to the average inflow velocity.
• The driving scale is of order the cloud radius throughout the cloud's entire evolution.
• The grid scale should exceed the minimum cooling length, ℓ cool ∼ min(c s t cool ) to resolve 2D phase structure.The 1D temperature phase structure is remarkably well-converged at lower resolutions.
• Our simulations suggest the existence of a minimum length scale for resolving turbulence, ℓ turb,sonic , for clouds with an equilibrium temperature of 5 × 10 3 ≲ (T cl /K) ≲ 2 × 10 4 .Under-resolving this scale seems to artificially amplify the violence of shattering.When this scale is resolved, the entrained cool phase is composed of larger clouds.
We thank M. Gronke for useful discussions and for sharing some sample code to compute the velocity structure function.We are grateful to James Bordner, Mike Norman, and the other enzo-e developers.GLB acknowledges support from the NSF (AST-2108470, XSEDE), a NASA TCAN award, and the Simons Foundation.DBF is supported by the Simons Foundation through the Flatiron Institute.Software: numpy (Harris et al. 2020), matplotlib (Hunter 2007), yt (Turk et al. 2011), scipy (Virtanen et al. 2020) Our approaches for characterizing v turb all build on the idea that a velocity field can be decomposed into a laminar part and a turbulent part.Consider an ideal turbulent flow in which the laminar part of the velocity field is uniform.In this scenario, the magnitude of the laminar part sets the average of the velocity field and the turbulent part sets the dispersion in the velocity values.For this reason, our methods for measuring a spatially averaged v turb (in a given gas phase) all measure this dispersion in one way or another.
Unfortunately, the flows considered in this work are more complex: the laminar portion of the flow has spatial gradients.Figure 16a illustrates these gradients for several velocity components measured on the log χ e/e cl = 1/6 isosurface of our χ = 1000, ξ sh = 27.8 simulation at 0.5t cc .In more detail, the panel shows the conditional distributions14 of multiple velocity components as a function of cos θ spherical , where θ spherical is the polar angle measured from the center of the inflow boundary.
Unless they are removed, such gradients can dominate or inflate the dispersion of the global velocity distribution, which can bias our measurements of v turb .Fig. 16b, suggests that this is less of an issue after early times (once v turb has had time to grow) because the dispersion from turbulence is larger relative to the laminar variations.However, it's clear that these gradients still remain problematic in the wind-aligned velocity component.Fig. 11b shows that large variations in the wind aligned velocity persist to later times, even as the cloud is accelerated.
We expect our v turb measurements from our geometric approach to be unaffected by this issue because it estimates v turb from the dispersion in v ϕ−like , which maintains a mean of zero at all times.However, the laminar variations will bias the measurements using our other approaches at early times.While one might expect our filtering measurements to be resilient to this effect, because it uses a local estimate of the laminar flow, at least some bias will remain given that these early-time gradients are most naturally described in spherical components.Throughout this work, we elect to just focus on turbulence in velocity components orthogonal to the wind direction, in our filtering and ⟨(δv) 2 ⟩(ℓ) measurements, in order to avoid biases from the wind-aligned velocity component.
As an aside, the resilience of our geometric approach to these biases are related to the definition of the velocity components.Consider ûr−like , which we define the unit vector parallel to the specific internal energy gradient (i.e.ûr−like = ∇e/||∇e||).Because this vector is always normal to the specific internal energy isosurfaces, we can define

Figure 1 .
Figure 1.Slice of a χ = 1000, ξ sh = 27.8,Mw = 1.5 simulation at 4.5tcc.This simulation has a resolution of 64 cells per cloud radius and the cloud is eventually entrained in the wind.The left two panels show the density and specific internal energy, which is T /µ scaled by physical constants.The right two panels show the high-pass filtered components of the velocity field transverse to vwind and the center panel shows the combined magnitude of these values.The insets highlight how the turbulent velocity has a clear temperature dependence.
Figure 2.Shows the phase dependence of v turb , measured via filtering, of the R cl = 64∆x run of our χ = 1000, ξ sh = 27.8,Mw = 1.5 simulation at 2.5tcc.The top panel includes contributions from all three velocity components.The bottom panel just includes contributions from the components transverse to vwind ; this is consistent with how v turb from filtering is measured throughout the remainder of this work.The solid orange line denotes the median while the dashed orange lines bound values between the 15th and 85th percentile.The dotted-dashed line shows v turb magnitudes that are equal to the sound-speed.The steep drop-off in v turb near Tw, is an artifact of the fact that wind is initially laminar.
Figure 3.Illustrates iso-temperature surfaces and derived v turb measurements for our the R cl /∆x = 64 run of our χ = 1000, ξ sh = 27.8,Mw = 1.5 simulation at 2.5tcc.Panel a shows a cut-away of five nested iso-surfaces measured at log χ e/e cl = 1/6, 1/3, 1/2, 2/3, 5/6 (for this system, T = 1.3 × 10 4 K, 3.3 × 10 4 K, 10 5 K, 3.3 × 10 5 K, 3.3 × 10 4 K, 10 6 K).The arrow illustrates the φ direction measured in the plane transverse vw.Panels b and c respectively show the normalized areaweighted distributions of the v normal and v ϕ−like velocity components measured on the iso-surfaces pictured in a. Panel d shows the standard deviation of the distributions from panel c (colored diamonds), as well as data derived from other isosurfaces (gray circles), plotted as a function of log χ e/e cl .
Figure 7.Like the top row of Fig.6except that the pictured simulations primarily vary the cloud temperature.Each simulation has χ = 100 and Mw = 1.5.We expect at higher resolution that the power-law slope below e break in the purple curve will be closer to 0.5 (i.e. the slope of the dashed black line).

Figure 8 .
Figure8.Like the top row of Fig.6except that the pictured simulations primarily vary in χ.We have made two compromises in presenting this data.First, we fix β to 0.25 for all panels.This is done as a simplification because β changes on a timescale related to χ.Second, the rightmost v turb panel compares simulations at a fixed value of v rel /vw rather than at a fixed time.The last panel typically compares the simulations at a point in evolution when v turb stabilizes (see § 4.1.3).However, that time seems to come much later in our χ = 10 4 simulation, after the simulation terminates.While we include the χ = 10 4 run for completeness, strong resolution dependence (see Table1) and the atypical shape of the cool-phase mass evolution may indicate that it is not well-converged.As noted in § 2, some material that started in the cloud leaks out of the domain at 6.5tcc, which coincides with the large drop-off in cool-phase mass.

Figure 9 .
Figure9.Like the top row of Fig.6except that the pictured simulations primarily vary Mw.The solid (dotted) lines show data from simulations with R cl /Mw = 37.6 pc (376 pc) and ξ sh = 5.73 (57.3).We note that the c s,hot t cool,min is 6.56 pc for all simulations in this plot.As we will show in panels e-h of Fig.10, v turb evolves more slowly in higher Mw runs.Consequently, the "late times" panel shows data from Mw = 0.75, 1.5 runs at 20tcc, and data from Mw = 3, 6 runs at 30tcc (we did not run the Mw = 6 simulation to late enough times or with a long enough domain for an optimal late-time comparison).
Figure 11.Points of a given color show v turb and (vw−v isosurface )/vw measurements for different sections of the e break isosurface from R cl /∆x = 64 run of our χ = 10 3 , ξ sh = 27.8 simulation.The points' colors indicate the simulation time that the measurement is associated with.The isosurfaces are split into bins based on each facet's position along the vwind .Each bin has a width of R cl ; there are more bins when the cloud is more elongated.The averages and standard deviations are all weighted by the area of each facet.While it is not shown, we have evidence indicating the data's slope may change when the principal dimensionless numbers are varied.
established that gas in χ ≤ 100 simulations is roughly distributed along the

Figure 13 .
Figure13.Illustrates how the evolution of various quantities in our χ = 1000, ξ sh = 27.8 simulation are affected by resolution.The top two rows show evolution of the coldphase mass and of the relative velocity between the cold and hot phases.Subsequent rows show evolution of quantities computed from the e break isosurface including surface area, average inflow velocity, and the turbulent velocity.
Figure 15.Illustrates resolution effects on the χ = 1000, ξ sh = 27.8 simulation at 7.5tcc.The top row shows ⟨|δv|⟩(ℓ) measurements for gas with 1/12 ≤ log χ e/e cl < 3/12 (panel a), the projected 1D phase distribution (panel b), and the bulk mass evolution for cold gas with ρ > ρmix (panel c).The dotted black line in panel a shows ⟨|δv|⟩(ℓ) ∝ ℓ 1/3 , the scaling expected for Kolmogorov turbulence.The brown shaded region in panel b denotes the gas phases considered in ⟨|δv|⟩(ℓ) while the vertical dotted line indicates the location of t cool,min .While the top row shows measurements from all resolutions, subsequent rows only show data for individual simulations.Panels d, f, h, j, and l shows the distribution of velocity difference magnitudes at the grid scale for gas with 1/12 ≤ log χ e/e cl < 3/12 (the average of this distribution is ⟨|δv|⟩(ℓ = ∆x)).The region enclosed by the grey dashed lines in these panels and panel a denote the range of cs values for the selected phase bin.Panels e, g, i, k, and m shows the density projection for each run.

Figure 16 .
Figure16.The probability density functions of several velocity components (in the cloud's rest-frame), as measured on the log χ e/e cl = 1/6 iso-surface for our χ = 1000, ξ sh = 27.8,R cl /∆x = 64 simulation at multiple times.The contours bound the region containing the most frequently occurring 68.4% of values at a given cos(θ spherical ).The fluctuations in a distribution's mode arises from the mostly-spherical laminar flow at early times.The dotted lines show the mean values of v r−like and v ϕ−like as functions of cos(θ spherical ).The vertical extent of a contour arises from turbulence (and are somewhat inflated by asymmetries in the flow).At early times, estimating v turb from the variance in any velocity component, other than v ϕ−like , without explicitly accounting for these laminar variations, will yield over-estimates.

Table 1 .
Table of simulations.χ Mw R cl (pc) tcc/t cool,mix t shear /t cool,min break = e min,cool , which coincides with the minimum of t cool . 5For now, we're just interested in α; the following subsections will discuss v turb,break .