A simple model for mixing and cooling in cloud-wind interactions
DDraft version January 27, 2021
Typeset using L A TEX twocolumn style in AASTeX63
A simple model for mixing and cooling in cloud-wind interactions
Matthew W. Abruzzo , Greg L. Bryan ,
1, 2 and Drummond B. Fielding Department of Astronomy, Columbia University, New York, NY 10027, USA Center for Computational Astrophysics, Flatiron Institute, 162 5th Ave, New York, NY 10003, USA (Received ?; Revised ?; Accepted ?)
Submitted to ApJABSTRACTWe introduce a simple entropy-based formalism to characterize the role of mixing inpressure-balanced multiphase clouds, and demonstrate example applications using enzo-e (mag-neto)hydrodynamic simulations. Under this formalism, the high-dimensional description of the sys-tem’s state at a given time is simplified to the joint distribution of mass over pressure ( P ) and entropy( K = P ρ − γ ). As a result, this approach provides a way for (empirically and analytically) quantifyingthe impact of different initial conditions and sets of physics on the system evolution. We find thatmixing predominantly alters the distribution along the K direction and illustrate how the formalismcan be used to model mixing and cooling for fluid elements originating in the cloud. We further confirmand generalize a previously suggested criterion for cloud growth in the presence of radiative cooling,and demonstrate that the shape of the cooling curve, particularly at the low temperature end, can playan important role in controlling condensation. Moreover, we discuss the capacity of our approach togeneralize such a criterion to apply to additional sets of physics, and to build intuition for the impactof subtle higher order effects not directly addressed by the criterion. Keywords:
Astrophysical fluid dynamics (101) — Galaxy evolution (594) — Interstellar medium (847)— Circumgalactic medium (1879) — Galaxy winds (626) INTRODUCTIONStellar feedback driven galactic outflows play a crit-ical role in galaxy formation and evolution. They areimportant for regulating star formation and transport-ing metals out of galaxies (White & Rees 1978; Dekel &Silk 1986; White & Frenk 1991). To reproduce observedgalaxy properties, large-scale cosmological simulationsoften assume that these galactic winds effectively accel-erate cool gas, with mass loading factors of ∼ −
100 (e.g.Pillepich et al. 2018; Dav´e et al. 2019). Moreover, obser-vations of outflows serve as direct evidence for the exis-tence of stellar feedback. A robust test of simulations istheir ability to reproduce these observations (Somerville& Dav´e 2015).Observations have revealed that these winds are in-herently multi-phase (for reviews of observations see
Corresponding author: Matthew W. [email protected]
Veilleux et al. 2005; Rupke 2018); there is cool gas co-moving with hot gas. This multi-phase nature presentsa challenge to the conventional model that the outflowsare driven by hot ( (cid:38) K) winds produced by super-novae. There has been a considerable effort to deter-mine whether ram pressure acceleration of cool ( ∼ K) ISM, by these winds, can produce a co-moving multi-phase flow (e.g. Klein et al. 1994; Cooper et al. 2009;Scannapieco & Br¨uggen 2015; Schneider & Robertson2017; Sparre et al. 2019). Difficulties arise because hy-drodynamical instabilities (e.g. Kelvin-Helmholtz andRayleigh-Taylor) grow from the initial velocity differencebetween the cloud and wind, and drives mixing betweenthe phases; the phases can be homogenized before thecloud is entrained.Klein et al. (1994) showed that this destruction of coolclouds is roughly characterized by the cloud-crushingtime scale . For a non-radiative cloud with density ρ cl and radius R cl , initially at rest with respect to a hotwind, with density ρ w = ρ cl /χ ( χ ∼ a r X i v : . [ a s t r o - ph . GA ] J a n Abruzzo et al. locity v w , this time-scale is given by t cc = χ / R cl v w . (1)Because the cloud is destroyed within a few t cc , and t cc is a factor of ∼ χ / smaller than the ram pressureacceleration time-scale, it’s challenging for hot winds toentrain the cool gas before it’s destroyed.Subsequent studies have modeled additional physicaleffects in attempts to delay cloud disruption long enoughfor them to be embedded within the wind. The mostcommon additional set of physics is radiative cooling(e.g. Cooper et al. 2009; Scannapieco & Br¨uggen 2015;Schneider & Robertson 2017). However, the generalconsensus was that the cloud’s lifetime is not prolongedenough for it to be fully entrained in the wind. More-over, Zhang et al. (2017) compellingly showed (withsemi-analytic methods) that real observations cannot bereproduced by ram-pressure accelerated cold clouds inthe interval of time before they are destroyed.Another approach for extending the cloud’s lifetimehas been the inclusion of magnetic fields. Certain con-figurations can inhibit mixing and provide an additionaltension force that resists destruction. McCourt et al.(2015) demonstrated that the presence of a tangled mag-netic field in the cloud gave promising results for χ = 50(with the inclusion of radiative cooling). Unless mag-netic pressure dominates ( β <
1) in the wind, however,magnetic fields alone don’t appear to inhibit the disrup-tion of higher density contrast clouds ( χ ∼ in situ cloud forma-tion within a cooling outflow (e.g. Thompson et al. 2015;Schneider et al. 2018; Lochhaas et al. 2020). These al-ternatives have achieved varying degrees of success, butno single model appears to apply in all cases.Recent work (Armillotta et al. 2016; Gronke & Oh2018) has shed new light on cloud acceleration by a hotwind in the limit of rapid cooling. Gronke & Oh (2018,2020a) showed that if mixed gas cools sufficiently fast,then it becomes a part of the colder cloud phase be-fore it is further homogenized with the wind. We here-after refer to this process as turbulent radiative mixinglayer (TRML, which we pronounce as “turmoil”) en-trainment. This process not only inhibits the depletionof the cloud mass, but also transfers mass and momen-tum to it from the wind (similar to an inelastic collision). Gronke & Oh (2018) argued that this process occurswhen the mixing time-scale, ∼ t cc , exceeds the coolingtimes-scale of the mixing layer, t cool , mix ; the mixinglayer has a temperature T mix ∼ √ T cl T w and numberdensity n mix ∼ √ n cl n w . They recast this criterion, t cc > t cool , mix , as a radius requirement for a sphericalcloud. Clouds should survive when their radius exceeds R cl , crit ∼ v w t cool , mix χ / ≈ T / , M w P Λ mix , − . χ , (2)where T cl , ≡ T cl / (10 K), M w is the mach numberof the wind, P ≡ nT / (10 cm − K), and Λ mix , − . ≡ Λ( T mix ) / (10 . erg cm s − ) .This picture has been bolstered by recent related stud-ies of individual shear layers (Ji et al. 2019; Fieldinget al. 2020; Tan et al. 2020). These simulations lack theoverall cloud geometry but are able to reach significantlyhigher resolution. The key finding from these studies isthat the crucial parameter that determines the rate ofcooling and the rate at which wind material is advectedinto the cloud is the ratio of the cooling time to the eddyturn over time. The eddy turn over time is comparableto the cloud crushing time.Recently, Li et al. 2020 found a different survival cri-terion, based on the cooling time of the wind, t cool , w .They argue that clouds survive, when t cool , w /t cc <
10 ¯ f ,while ¯ f is an order-unity term with weak power-law de-pendence (exponents range from ∼ . ∼ .
6) on R cl , n w , and v w . Sparre et al. (2020) reached a similar re-sult (with additional M w dependence), while Kanjilalet al. (2020) found results in support of the Gronke &Oh (2018) criterion. Therefore, this disagreement overthe survival criterion remains unresolved.Prior works have clearly assembled models for howvarious circumstances and physical effects modify thecloud-wind interaction. Unfortunately, simple charac-terizations of how different effects influence the interac-tion are often incompatible with one other. The complexmultidimensional nature of this process is a barrier tomaking comparable and composable characterizations.We, therefore, explore higher order characterization ofthe mixing and cooling processes in these systems in or-der to isolate the competing effects. Following the arguments from Begelman & Fabian (1990), (when µ has ρ and p dependence) we find that T mix ∼ √ χµ mix T cl /µ cl , n mix ∼ µ cl n cl / ( √ χµ mix ), and R cl , crit scales with µ − / . Using grackle (Smith et al. 2017), in tabulated mode, we find slightlydifferent characteristic values for solar metallicity, χ = 100, M w = 1, T cl , = 1, and P = 1 Under these conditions, µ cl = 0 .
827 and R cl , crit ∼ mix is unchanged). ixing & Cooling in Cloud-Wind Interactions p − K ) phase distribution broadlycapture the cloud-wind interaction’s evolution. Our ap-proach is particularly conducive for comparing the im-pact of radiative cooling against mixing. Furthermore, itnaturally complements the Gronke & Oh (2018, 2020a)physical model for turbulent radiative mixing layer en-trainment.Our paper is organized as follows: in § § § §
5, we describeour results from two example applications that demon-strate how the formalism both: (i) broadly captures thesystem’s evolution and (ii) can quantitatively character-ize how different conditions and physical effects mod-ify the system’s evolution. These applications include anon-radiative parameter study ( §
4) and a more detailedstudy involving radiative cooling ( § § § MODEL OVERVIEW2.1. p − K Representation
The cloud-wind interaction is fundamentally an inter-action between a finite pool of colder, dense gas (thecloud) and a large reservoir of hotter, more diffuse gas(the wind). We are interested in understanding howdifferent physical conditions affect the interaction’s out-come. The two outcomes are: (i) the homogenizationof the colder phase within the more abundant hotterphase (cloud destruction) or (ii) the long-term survivaland coexistence of both phases (entrainment). It’s in-structive to specify the system’s state purely in terms ofthis thermodynamic description.Figure 1 helps illustrate this premise for a non-radiative simulation (see § R cl / ∆ = 64 NR-X100 simulation). The left column il-lustrates snapshots of the system’s morphological evolu-tion while the center-left column shows the n − T phasespace evolution. The initial properties of each gas phaseare denoted by black circles in the phase diagrams. Atearly times, the initial shock introduces pressure pertur-bations and slightly elevates the entropy of some fluidelements originating in the colder phase. As the inter-action progresses, adiabatic mixing drives gas from thecolder dense phase towards the hotter, diffuse phase. At a given snapshot, the phase distribution clearly encodesinformation about the system’s state.Although number density, n , and temperature, T , arefamiliar thermodynamic quantities, we choose to baseour mixing model on p − K phase space (center-right col-umn of Figure 1). The quasi-isobaric nature of the prob-lem makes pressure, p , an intuitive choice for a phasespace axis. Although the supersonic wind seeds smalltransient pressure perturbations, the second phase spacedimension effectively indexes the continuum of proper-ties gas can have between the two initial states. Follow-ing convention, we pair p with an entropy-like quantity, K = pρ − γ (hereafter, we refer to K as entropy) andwe take γ = 5 / p ( K ) that increase by factor of 10.The choice of p − K space over p − T space is somewhatdiscretionary. Because T directly characterizes thermalenergy, it more directly governs heat flow (e.g. cooling).However, it’s easier to characterize scale-free, dynam-ical effects in terms of K ; relating T to hydrodynamicquantities requires scale-dependent knowledge about themean molecular mass, µ . Additionally, the entropy of afluid element is unchanged by compression or expansion,which explains why the spread in gas at a given pressurelies along K contours in the center column of Figure 1.Only irreversible processes change entropy: shocks andmixing increase it while radiative cooling decreases it.Moreover, fluid elements have continuous trajectoriesthrough p − K space in the absence of shocks becausemixing and cooling modify K smoothly.We now take a more careful look at the p − K repre-sentation and evolution for the non-radiative, hydrody-namic cloud-wind interaction. At initialization, all gashas a single pressure p , the colder phase lies at K cl = pρ − / , and the hotter phase lies at K w = χ / K cl . Asthe system evolves, hydrodynamical instabilities drivemixing of the two phases. Mixing increases the entropyof the colder phase gas and initially decreases the en-tropy of the hotter phase gas. By the time the cloudis destroyed, all gas has an entropy of K ∼ K w . Thegas pressure remains relatively constant throughout thisprocess.The precise p − K evolution depends on the initialconditions and the simulated physics. Changes to eithermay modify how the distribution evolves. Herein liesthe true value of this thermodynamic description of thesystem’s state: it provides a low-dimensional domain forcharacterizing mixing that is well-suited for comparingdifferent physical processes.2.2. Mixing Model
Abruzzo et al. [ cl R cl ] 10 M cl,0 d Md log nd log T [dex ] 10 M cl,0 d Md log Pd log K [dex ]5.02.50.02.55.0 y [ R c l ] t / t cc = 0.5wind5.02.50.02.55.0 y [ R c l ] t / t cc = 2.5wind5.02.50.02.55.0 y [ R c l ] t / t cc = 4.5wind5.02.50.02.55.0 y [ R c l ] t / t cc = 6.5wind15 10 5 0 5 10 15 x [ R cl ]5.02.50.02.55.0 y [ R c l ] t / t cc = 8.5wind 10 T [ K ] T [ K ] T [ K ] T [ K ] n [cm ]10 T [ K ] K [ c m g / s ] passivescalarall mass K [ c m g / s ] K [ c m g / s ] K [ c m g / s ] m i x i n g P / k B [Kcm ]10 K [ c m g / s ] M cl,0 dMdlogK [dex ] Figure 1.
Illustration of the cloud-wind interaction until 8 . t cc for an initially spherical cloud, a wind with Mach number M w = 1 . χ = 100 at a resolution of R cl / ∆ x = 64. The left column illustrates the cloud surface density. Wenote that simulated domain is far longer than these panels depict (the shock remains in the domain over the entire simulation)and the axes labels don’t describe the position in the domain. The center columns depict Mass-weighted n − T (center-left) and P − K (center-right) phase diagrams of all gas in the simulation. The n − T diagrams are computed assuming a fixed meanmolecular mass of µ = 0 .
6. The black circles near the top (bottom) of the panels denote the initial location of fluid elementsoriginating in the wind (cloud). The dashed (dotted) gray lines in the n − T diagram denote logarithmically spaced lines ofconstant P ( K ). The rightmost column depicts the one dimensional entropy distribution. The red histogram traces all gas inthe domain, while the blue traces just the fluid elements initialized in the cloud. Consider the motion of fluid elements originating inthe colder phase through p − K space; this is illustratedby the red arrow in Figure 1. We anticipate factors thatinhibit mixing, like strong magnetic fields, to deceleratethe rate that fluid elements increases in entropy. Con-versely, we expect factors that hasten mixing, such aslarger values of M w , to accelerate that rate. Thus, wecan empirically model mixing by characterizing the fluidelements’ motion through p − K space. For simplicity, we largely ignore pressure perturba-tions in the context of mixing. Motion along the pres-sure dimension may be particularly relevant in caseswith strong sources of non-thermal pressure support.The right column of Figure 1 illustrates just the entropyevolution of fluid elements.In this work, we focus on the motion of fluid elementsthat originate in the cloud. Because we trace these fluidelements with a passive scalar, we refer to their totalmass as M ps . We defer analysis of the p − K evolution ixing & Cooling in Cloud-Wind Interactions ( p / k B ) [K cm ]282930313233 l o g K [ c m g / s ] c m c m c m c m c m K K K l o g t c oo l = l o g ( K / K c oo l ) [ M y r ] Figure 2.
Illustration of contours of − K/ ˙ K cool (equal to t cool for a fixed composition, perfect gas) in pressure-entropyspace. The dashed (dotted) denote logarithmically spacednumber density (temperature) contours. In the lower-leftwhite region, heating dominates cooling. for fluid elements originating in the hot phase to a fu-ture work. These fluid elements encode information thatis most relevant at times and locations in phase spacewhere the motion of fluid element from the clouds aren’trepresentative of all fluid elements at that location. Webriefly revisit this point in § t can be quantitatively described by its initialpressure ( p ), the distribution of the initial cloud fluidelements with respect to entropy, ( dM ps /dK )( K, t ), and˙ K ( K, t ), the ensemble averaged Lagrangian derivative of K for all initial cloud fluid elements with a given valueof K . The form of ˙ K ( K, t ) dictates the outcome of theinteraction. It’s dependent on the initial conditions andmodeled physics (in full generality, the notation resem-bles ˙ K ( K, t ; p , γ, χ, R cl , M w , β, . . . )).For purely non-radiative (magneto)hydrodynamic in-teractions, the only sources of entropy are the initialshock and mixing. Because we expect mixing to dom-inate outside of highly supersonic flows, we refer to˙ K ( K, t ) as ˙ K mixing ( K, t ) for these simulations. Becauseof the scale-free nature of such an interaction, the valueof ˙ K mixing ( aK cl , bt cc ) / ( K cl t − ), where a and b are arbi-trary positive values, is constant for any choice of ρ cl , p ,and R cl (as long as γ , χ , M w , β , and the initial geometryremain unchanged). When ˙ K mixing ( K, t ) is known, it can be used to pre-dict the outcome of interactions involving radiativecooling through comparisons against expected contri-butions from cooling, ˙ K cool ( p, K ). For optically thingas, ˙ K cool = K ˙ e/e = − K/t cool and Figure 2 illustrates t cool ( p, K ). When ˙ K mixing ( K, t ) + ˙ K cool ( p , K ) (cid:29) K ∈ [ K cl , K w ], we generally expect the colder phase tobe destroyed. Conversely, the existence of a large sub-interval over [ K cl , K w ], where the sum is far less thanzero, suggests long term survival of a cool phase. Theoutcome is ambiguous when the sum is close to zero be-cause the sum does not directly give the total ˙ K ( K, t )for an interaction with cooling (hereafter ˙ K total ( K, t )).This reasoning is reminiscent of the comparisons be-tween t cc and t cool , mix that underlies the Gronke & Oh(2018) survival condition. In fact, we can apply anal-ogous arguments in the context of our mixing modelto derive a criterion comparable to t cc > t cool , mix . Forsimplicity, suppose ˙ K mixing , char ( K ) gives the character-istic mixing rate for a non-radiative interaction at en-tropy K at times when dM ps /dK ( K, t ) > § K mixing , char ( K ) is indeed a well-conceivedquantity). Then, following their logic, we expect turbu-lent radiative mixing layer entrainment to occur when˙ K mixing , char ( K mix ) < (cid:12)(cid:12)(cid:12) ˙ K cool ( p , K mix ) (cid:12)(cid:12)(cid:12) , where the en-tropy at the mixing layer is K mix ∼ √ K cl K w . METHODS3.1.
Simulations
To run magnetohydrodynamical (MHD) simulationsfor this work, we make use of the enzo-e code. Thiscode is a rewrite of enzo (Bryan et al. 2014) that tar-gets exascale computing and is built-on the distributed,scalable, adaptive mesh refinement (AMR) framework, cello (Bordner & Norman 2012, 2018). Although enzo-e is still under active development, it has maturedenough that it can be used for basic scientific studies.For this work we implemented the second-order accu-rate unsplit VL + CT (van Leer + Constrained Trans-port) algorithm presented by Stone & Gardiner (2009).This is a predictor-corrector scheme that employs theConstrained Transport (CT) method (Evans & Haw-ley 1988). Each simulation uses second order recon-struction and the HLLD Approximate Riemann Solver(Miyoshi & Kusano 2005). We provide a brief assess-ment on how the choice of integrator affects the evolu-tion of the cloud-wind interaction in Appendix A. http://cello-project.org The prediction step always uses first order reconstruction
Abruzzo et al.
Table 1.
Table of simulations (additional simulations are discussed in the appendices). All simulations were initialized with aninitial thermal pressure of p/k B = 10 cm − K.name χ M w a β b R cl (pc) T cl c (K) R cl / ∆ x Cooling d t cool , mix /t cc R cl /(cid:96) cool e t cool , mix /t cool , cl NR-X100
100 1.5 ∞ ×
8, 64 N
NR-X300
300 1.5 ∞ × NR-X1000 ∞ × NR-X100-M0.75
100 0.75 ∞ × NR-X100-M3
100 3 ∞ × NR-X100-M4.5
100 3 ∞ × NR-X100-B10
100 1.5 10 1 4 × NR-X100-B100
100 1.5 100 1 4 × NR-X100-B1000
100 1.5 1000 1 4 × SlowRC-T1e4
100 1.5 ∞
64 Y 8.982 0.668 1.044
FastRC-T4e4
100 1.5 ∞ ×
64 Y 0.2165 388 63.06
FastRC-T1e4
100 1.5 ∞
57 10
64 Y 0.1576 38.1 1.044
BPLawRC-1
100 1.5 ∞ BPLawRC-2
100 1.5 ∞ BPLawRC-6
100 1.5 ∞ BPLawRC-60
100 1.5 ∞ a sonic Mach number of the wind b plasma beta (thermal pressure divided by magnetic pressure) c For non-radiative simulations, T cl is computed assuming µ = 0 . d Indicates whether radiative cooling is included. Simulations with the value “BPL” used a custom broken power-law coolingcurve with shaped given by the specified η mix , cl = t cool , mix /t cool , cl and Equation 7 (and a constant µ of 0.6). e (cid:96) cool is the “cooling length,” (cid:96) cool = min( t cool c s ) (McCourt et al. 2018) In our cloud-wind simulations, we solve the ideal,adiabatic MHD equations on a fixed, uniform three-dimensional Cartesian grid. We initialize each simula-tion with a spherical cloud of radius R cl embedded in asteady wind with M w = 1 . χ lower than the cloud. The cloud initiallyhas no bulk velocity and is in pressure equilibrium withthe wind. We also initialize a passively advected scalarthat traces the gas initially confined to the cloud. In asubset of our non-radiative simulations we also initializemagnetic fields transverse to the wind that are constantthroughout the entire domain.Our runs including radiative cooling employ the grackle chemistry and cooling library (Smith et al.2017) and assume solar metallicity. Our main coolingruns use the Haardt & Madau (2012) UV backgroundmodel and don’t use the self-shielding approximation. In § https://grackle.readthedocs.io/ For simplicity, we restrict cooling to only occur be-tween ∼ T cl and ∼ . T w ; we modified the tables tohave Λ /m H = 10 − erg cm s − g − in the restricted re-gions . We assess the consequences of restricted coolingin Appendix A.The simulation domain extends 100 R cl downwind ofthe cloud’s initial location and 12 R cl along each trans-verse dimension. Gas with wind properties (includingmagnetic fields, if β is finite) flows into the domain 20 R cl upwind of the cloud’s initial center of mass, We enforceoutflow conditions for the other boundaries. Table 1 provides a summary of our simulation proper-ties. We primarily employ low-resolution ( R cl / ∆ x = 8)simulations for our non-radiative parameter study. Theprimary simulations we use to assess the cloud-wind in- Due to an oversight, CMB Compton cooling occurs (and actuallydominates) in the restricted regions for simulations without bro-ken power-laws. For these cases t cool , restrict (cid:38) . n H / ( n e µ ) × s (cid:38) t cc . Thus, radiative losses are minimal in these re-stricted regions (simulations are run for (cid:46) t cc ). Enforcement of outflow conditions, on transverse boundaries,maintains ∇ · B = 0 at roughly the same precision as periodicboundaries. ixing & Cooling in Cloud-Wind Interactions l o g ( K / K c l ) m ps loss = 100 0.50 t cc m ps loss t cc m ps loss t cc m ps loss t cc m ps loss t cc l o g ( K / K c l ) m ps loss = 300 m ps loss 0% m ps loss 14.9% m ps loss 46.0% m ps loss10 P / k B ( K cm )0.00.51.01.52.0 l o g ( K / K c l ) m ps loss = 1000 P / k B ( K cm )0% m ps loss 10 P / k B ( K cm )0.2% m ps loss 10 P / k B ( K cm )31.6% m ps loss 10 P / k B ( K cm )60.1% m ps loss 10 l o g M c l , d M p s d l o g P d l o g K [ d e x ] Figure 3.
Mass-weighted Pressure-entropy phase evolution for gas originating in the cloud (in contrast Figure 1 shows all gasin the simulation) for three non-radiative, hydrodynamic simulations. Each row depicts a simulation with a different χ (densitycontrast) and simulation time increases from left to right. The fraction of the gas originating in the cloud (i.e. the passive scalar)that has exited the simulation domain is recorded in each panel. We scale the entropy axes in terms of log χ ( K/K cl ) to facilitateeasier comparisons between simulations, since the colder phase gas is initialized with K = K cl and will have K ∼ K w = χ / K cl when it homogenizes with the hotter gas (these locations are indicated by black circles). Finally, the black square denotes K mix = √ K cl , K w and the dashed line line denotes ρ cl / teractions including radiative cooling have R cl / ∆ x = 64.We briefly investigate the impact of resolution in Appen-dices A and C.Finally, we note that our simulations employ a refer-ence frame tracking scheme. Unfortunately, an imple-mentation bug caused our R cl / ∆ x <
64 simulations toeffectively have no frame tracking (compared to v w , theframe velocity is near-zero). However, our R cl / ∆ x = 64simulations used an improved version of the code inwhich the frame velocity was properly updated every0 . t cc . In Appendix B we provide more details aboutthe scheme and show that the differences have a negli-gible impact on our result.3.2. ˙ K Calculation
As a post-processing step, we estimate ˙¯ K as a func-tion of K and t for each of our simulations. Recall that˙ K ( K, t ) is the Lagrangian entropy derivative averagedover all fluid elements originating in the cloud with en- tropy K . ˙¯ K is simply ˙ K averaged over some time inter-val ∆ t .Because mass (or in this case, passive scalar mass)is conserved, we can write an analog to the continuityequation: ∂∂t (cid:18) dM ps dK (cid:19) + ∂∂K (cid:18) ˙ K dM ps dK (cid:19) = 0 . (3)This describes the changes in passive scalar mass profileas a function of K for a pair of snapshots measured at t n and t n +1 .Consider a set of discrete K bins where the i th bin hascenter K i , width δK i , and encloses a passive scalar massof m ps ,i . Integrating equation 3 in time from t n to t n +1 and over the i th K -bin (from K i − / = K i − . δK i to K i +1 / = K i + 0 . δK i ) yields m n +1ps ,i − m n ps ,i t n +1 − t n = (cid:18) ˙ K dM ps dK (cid:19) n +1 / i − / − (cid:18) ˙ K dM ps dK (cid:19) n +1 / i +1 / , (4) Abruzzo et al. l o g ( K / K c l ) m ps loss = 10 0.50 t cc m ps loss t cc m ps loss t cc m ps loss t cc m ps loss t cc l o g ( K / K c l ) m ps loss = 100 m ps loss 0% m ps loss 0% m ps loss 33.6% m ps loss10 P / k B ( K cm )0.00.51.01.52.0 l o g ( K / K c l ) m ps loss = 1000 P / k B ( K cm )0% m ps loss 10 P / k B ( K cm )0% m ps loss 10 P / k B ( K cm )0.1% m ps loss 10 P / k B ( K cm )10.8% m ps loss 10 l o g M c l , d M p s d l o g P d l o g K [ d e x ] Figure 4.
Mass weighted pressure-entropy evolution (for gas originating in the cloud) as in Figure 3 except that χ is fixed at100 and the simulations include transverse magnetic fields. In each of the simulations, the initial magnetic fields is initializedwith constant β throughout the entire domain. where ( ˙ KdM ps /dK ) n +1 / is time-averaged between t n and t n +1 . By selecting a minimum bin, K , such that˙ K / = 0 at all t , we can compute ( ˙ KdM ps /dK ) n +1 / atall bin interfaces from changes in the measured profiles.Finally, if we assume that ( dM ps /dK ) i +1 / is near con-stant between t n and t n +1 then˙¯ K n +1 / i − / ≈ (cid:18) ˙ K dM ps dK (cid:19) n +1 / i − / (cid:44)(cid:18) dM ps dK (cid:19) n +1 / i +1 / . (5)We approximate ( dM ps /dK ) ni +1 / via linear interpola-tion of the average dM ps /dK from adjacent bins, andthen average the values from t n and t n +1 to estimate( dM ps /dK ) n +1 / i +1 / . To enforce our assumption, we focuson ˙¯ K measurements where f , given by f = | ( dM ps /dK ) n +1 i +1 / − ( dM ps /dK ) ni +1 / | min(( dM ps /dK ) n +1 i +1 / , ( dM ps /dK ) ni +1 / ) , (6)is less than 0.25 and omit measurements altogetherwhere f >
2. While our measurements of ˙¯ K may besomewhat biased, the overall dependence on K and t isstill useful, particularly when the dependence is stablein time. In practice, we estimate ˙¯ K from pairs of snapshotssatisfying t n +1 = t n + t cc /
2. If the passive scalar ad-vects into or out of the simulation domain, then Equa-tion 3 should include an additional source or sink term.Thus, we only consider snapshots at times before 1% ofthe initial passive scalar mass escapes the domain. Wenote that vorticity near the transverse outflow bound-aries can introduce artificial passive scalar inflow. Inpractice, this is only an issue for
NR-X100-M0.75 andwe conservatively discard all data from that run mea-sured after 7 t cc .Appendix C includes a brief study on how resolutionaffects measurements of the passive scalar mass profileand our measurements of ˙¯ K ( K, t ). Because the profile’sevolution is most sensitive to resolution for K -bins hold-ing under 1% of the total passive scalar mass M cl , , oursubsequent analysis primarily focuses on ˙¯ K n +1 / i − / mea-surements where m ni − , m ni , m n +1 i − , and m n +1 i are all ≥ . M cl , .For simulations with no cooling or inefficient cool-ing, the dependence of our ˙¯ K measurement on K and t are remarkably robust with respect to resolution. Atthe same time, measurements for simulations with rapid ixing & Cooling in Cloud-Wind Interactions NON-RADIATIVE PARAMETER STUDYRESULTS4.1. P − K Phase Space
In this section, we apply our formalism to a suite of(magneto)hydrodynamic non-radiative simulations thatprobe a wide variety of properties.Figure 3 illustrates the time evolution of the passive-scalar weighted P − K distribution for three hydrody-namic simulations with varying χ ( NR-X100 , NR-X300 , NR-X1000 ). Unlike the panels in the center-right col-umn of Figure 1, these only show phase distributionsfor fluid elements originating in the cloud, depict lowerresolution simulations, and have rescaled entropy axes.Throughout this work, we plot entropy as log χ ( K/K cl )to remove most of its χ dependence. With this defini-tion, it spans values from 0 through γ = 5 / K cl to K w (denoted byblack circles). Fluid elements that have already exitedthe domain should generally have comparable or larger K to the remaining ones since they mixed faster. Thehigher χ simulations lose fluid elements more quickly be-cause they have larger v w . The dotted black line denotes ρ cl /
3, a density threshold commonly used to identifycloud mass (e.g. Scannapieco & Br¨uggen 2015; Schnei-der & Robertson 2017; Gronke & Oh 2018). The figureshows that the vast majority of fluid elements cross thisthreshold by ∼ . t cc , (as discussed in Appendix C thereis some resolution dependence), which is consistent withprior work (e.g. Sparre et al. 2019). The log χ scalingclearly removes most of the χ dependence dependencein the entropy evolution.Each simulation’s P evolution follow a common evo-lution, largely independent of χ . In each case, the shockinitially produces large P perturbations. By 2 . t cc , themotions giving rise to the under-pressured gas have beenslightly damped. While not shown, the distribution’s ex-tent is stable between ∼ . t cc and ∼ . t cc , albeit withminor fluctuations in the minimum. The low pressuregas at these times is presumably supported by the vortic-ity produced by the initial shock, the post-shock flow in the shearing layer (at the cloud boundary), and the for-mation of the vortex rings (Klein et al. 1994). Between3 . t cc and 4 . t cc the minimum pressure drops once moreand subsequently all pressure perturbations damp away.The χ -dependence manifests in the P distribution intwo main ways. First, the minimum pressure at 4 . t cc is larger for χ = 100 than it is in the other cases. Whileit’s unclear how robust this difference is, we note thatKlein et al. (1994) reported that the post-shock flowwas the primary generator of vorticity in their χ = 10, M w = 0 .
9, 2D ellipsoidal cloud simulation and arguedthat it scales with ∼ χ / .The other difference, is that the mode of the χ = 1000pressure distribution has a positive offset at 8 . t cc (themode is roughly double the initial pressure at 9 . t cc ).This is an unexpected artifact caused by the reflectionof waves and discontinuities off of the transverse outflowboundaries. We reran these simulations with three timeslarger transverse widths and found that this artifact isfirst noticeable at ∼ t cc . Furthermore, while all threesimulations were affected, the magnitude strongly scaledwith χ . We don’t expect this effect to significantly influ-ence our results because subsequent analysis of NR-X180 , NR-X300 , and
NR-X1000 (our only χ >
100 simulations)ignores data at t/t cc > , . , and 4 . χ on the gas distri-bution, we now consider the effect of magnetic fields.Figure 4 illustrates the phase evolution of the fluid el-ements originating in the cloud for three simulationswith transverse magnetic fields of different strengths( NR-X100-B10 , NR-X100-B100 , NR-X100-B1000 ). Incontrast to the pure hydro cases, cloud destruction pro-ceeds far more slowly when β = 10. It takes longerfor fluid elements to cross the ρ cl / t = 8 . t cc , a much larger fraction of the fluid elementshave K < K w (and lie within the simulation domain).As β increases, the clouds are more readily destroyedand the phase distributions bear greater resemblanceto those of the pure hydro simulations. The minimumpressure appears to correlate with the initial β at earlytimes. This suggests that its supported by magneticstress. Because the magnetic fields are initially trans-verse, we expect the shock that propagates through thecloud to transfer energy to the magnetic field, therebyelevating the magnetic pressure. We expect the pressureto be less supported by vorticity (than in the purely hy-drodynamical case) because magnetic fields impede itsgrowth.Both figures clearly illustrate that much of the inter-esting evolution of the cloud-crushing problem occursover the K dimension. While there is variation in the0 Abruzzo et al. l o g M c l , d M p s d l o g K = 100 w = 1.5= = 300 w = 1.5= = 1000 w = 1.5=10 l o g M c l , d M p s d l o g K = 100 w = 0.75= = 100 w = 3.0= = 100 w = 4.5=0 1 2log ( K / K cl )10 l o g M c l , d M p s d l o g K = 100 w = 1.5= 10 0 1 2log ( K / K cl )= 100 w = 1.5= 100 0 1 2log ( K / K cl )= 100 w = 1.5= 1000 0.52.54.56.58.5 t / t cc Figure 5.
Temporal evolution of the mass profile as a func-tion of entropy for fluid elements originating in the cloud, dM ps /d log χ K . The profiles are normalized by the initialcloud mass. Each panel depicts a different simulation; χ , M w , and β vary across the top, middle, and bottom rows.Colors denote measurement times and dash lines denote pro-files measured after 1% of the passive scalar leaves the do-main. Gray vertical dashed lines denote K cl and K w (theinitial entropy of fluid elements in the cloud and wind). Thescaling of time (in terms of t cc ) and entropy largely removesthe effects of χ and M w (where M w ≥ .
5) for pure hydro-dynamic mixing. However, it does not account for differencesin subsonic winds. At late times, the distributions’ extentsconvey that increasing magnetic field strengths more effec-tively impede mixing. pressure, its both a transient effect that largely fadesaway at late times, and is smaller than the variation in K . 4.2. Passive Scalar Mass Distribution over K Having qualitatively established that the bulk motionof cloud fluid elements through P − K space both occursprimarily along K and reflects the initial conditions andmodelled physics, we now consider a more quantitativeparameterization of mixing.If we assume that pressure perturbations are broadlyunimportant for the system’s evolution, we can integrateover the phase distribution’s pressure dependence to get dM ps /dK . We effectively trade the information encodedin the pressure perturbations for a dimensionality reduc-tion. Recall that ( dM ps /dK ) dK specifies the mass ofall fluid elements with entropy between K and K + dK .Figure 5 illustrates the time evolution of ( dM ps /dK ) dK for a selection of times for each non-radiative simulationlisted in Table 1 with a resolution of R cl / ∆ x = 8. At t = 0 . t cc , each simulation’s profile has a peak at K cl and a long tail extending to K w . Over time, mixingincreases the entropy of the fluid elements near the loweredge of the histogram, K min . By t = 3 . t cc , K min startsto increase, indicating that all of the fluid elements fromthe cloud have started mixing. We largely ignore differ-ences in the profiles at intermediate and late times thatare depicted by dashed lines because different fractionsof passive scalar remain in the simulation domain whenthose are measured.For the hydrodynamical simulations (top two rows ofFigure 5), the narrow histograms near K w at t = 8 . t cc reflect how the initial cloud fluid elements have largelyhomogenized with the wind. The figure shows that thescaling of our K -bins and t in terms of log χ ( K/K cl )and t cc almost entirely captures the evolution’s χ de-pendence. This t scaling also largely removes the M w dependence for supersonic simulations; however, the his-tograms’ upper edges do scale weakly with M w .The subsonic run, NR-X100-M0.75 , has the mostunique evolution. In this case, mixing appears to morerapidly increase entropy below log χ ( K/K cl ) ∼ . t (cid:46) . t cc ) the profile evo-lution is largely the same as before, but by t = 3 . t cc thesuppression of mixing by the magnetic fields causes theevolutionary paths to diverge. Since stronger fields (in agiven configuration) more strongly suppress mixing, therate at which K min increases scales with increasing β .4.3. Mixing Rate Estimation
In this section, we use this distribution evolution toestimate the rate at which fluid elements from the cloudmix. Figure 6 illustrates our measurements of ˙¯ K ( K, t ),averaged over ∆ t = t cc /
2, as functions of K for eachof our simulations. See § § K mixing ( K, t ) because mixing is the dominant entropygeneration mechanism.In each hydro simulation, ˙¯ K mixing / ( K cl t − ) broadlyhas a power law relationship in terms of K/K cl , witha near-unity slope for 0 . (cid:46) log χ K/K cl (cid:46) .
4. Thesection below log χ K/K cl ∼ . ∼ . NR-X100 ) while the upper section’s slope maybe slightly shallower ( ∼ . NR-X100 ) and could betime dependent. Note that the top row of Figure 24from Appendix C suggests that these trends are robustto resolution effects. ixing & Cooling in Cloud-Wind Interactions l o g K m i x i n g [ K c l t cc ] = 100 w = 1.5= = 300 w = 1.5= = 1000 w = 1.5=0.00.51.01.52.0 l o g K m i x i n g [ K c l t cc ] = 100 w = 0.75= = 100 w = 3.0= = 100 w = 4.5=0 1 2log ( K / K cl )0.00.51.01.52.0 l o g K m i x i n g [ K c l t cc ] = 100 w = 1.5= 10 0 1 2log ( K / K cl )= 100 w = 1.5= 100 0 1 2log ( K / K cl )= 100 w = 1.5= 1000 0.51.52.53.54.55.56.57.5 t / t cc Figure 6.
Passive scalar mass weighted and time aver-aged values of ˙ K mixing ( K, t ) computed from the histogramsin Figure 5. Each ˙ K curve is averaged over 0 . t cc . Mea-surements have been omitted (made transparent) at any binedges where the interpolated dM ps /dK changes by more thana factor of 3 (1.25) over this interval. Measurements are alsomade transparent when they are adjacent to histogram binscontaining under 1% of the initial mass. The curves are onlycomputed at times when more that 99% of the fluid elementsoriginating in the cloud lies in the simulation domain. The higher χ simulations appear to have slightlysteeper slopes than NR-X100 , but the dearth of highquality measurements make this comparison tenuous,especially at high K . The higher M w simulations havegreater variance in their measurements than NR-X100 (possibly due to their stronger initial shocks), but areotherwise broadly consistent. Without better measure-ments, we’re unable to make any comparisons with
NR-X100-M0.75 .Next, we consider the MHD simulations. As in Fig-ure 5, the β ≤
100 ˙¯ K mixing ( K, t ) measurements onlystart diverging from the
NR-X100 measurements at t ∼ . t cc . The suppression of mixing gives ˙¯ K mixing ( K )a shallower slope. As the initial field strength de-creases, the measurements more closely resemble thosefrom NR-X100 .Figure 7 illustrates K/ ˙¯ K mixing ( K, t ), which is the timethat mixing takes to double K , in units of t cc . Itmakes the slope variations in ˙¯ K mixing , above and be-low log χ K/K cl ∼ .
8, more apparent. Additionally, K/ ˙¯ K mixing ( K, t ) is generally within a factor of ∼ . t cc in each hydro simulation. This implies that K / K m i x i n g [ t cc ] = 100 w = 1.5= = 300 w = 1.5= = 1000 w = 1.5=10 K / K m i x i n g [ t cc ] = 100 w = 0.75= = 100 w = 3.0= = 100 w = 4.5=0 1 2log ( K / K cl )10 K / K m i x i n g [ t cc ] = 100 w = 1.5= 10 0 1 2log ( K / K cl )= 100 w = 1.5= 100 0 1 2log ( K / K cl )= 100 w = 1.5= 1000 0.51.52.53.54.55.56.57.5 t / t cc Figure 7.
Time scales, for mixing to double the K of afluid element that originated in the cloud. For the hydrody-namical simulations, the time scales are relatively consistentacross χ times and they are of the same order as t cc . Al-though the timescales for the MHD simulations are similarat early times, they get longer at late times. The data selec-tion is the same as in Figure 6. ˙¯ K mixing ( K, t ) and t − share similar M w and χ depen-dence.The results in the section broadly indicate that˙¯ K mixing ( K, t ) robustly characterizes cloud destructionthrough mixing. For idealized conditions, our resultsfurther suggest that ˙¯ K mixing doesn’t have a strong t de-pendence and we can approximate ˙¯ K mixing ( K, t ) witha time-independent function, ˙ K mixing , char ( K ). In thepresence of additional physical effects (e.g. the presenceof magnetic fields), ˙¯ K mixing ( K, t ) shows stronger timedependence, and improves on the description of clouddestruction offered by t cc . This is conveyed in Figures 6and 7 for NR-X100-B10 ; ˙¯ K mixing ( K, t ) clearly capturesthe decreasing destruction rate, presumably caused bythe tangling of magnetic fields. RESULTS WITH COOLING5.1. P − K Phase Evolution
Next, we apply our mixing model to hydrodynamicsimulations with radiative cooling. We consider twomain regimes of cooling: slow, t cool , mix ∼ t cc , andfast, t cool , mix ∼ . t cc . Per Gronke & Oh (2018), thecloud should be destroyed in the former case and survivein the latter. For the fast cooling regime, we consider2 Abruzzo et al. [ cl R cl ] 10 M cl,0 d Md log nd log T [dex ] 10 M cl,0 d Md log Pd log K [dex ]5.02.50.02.55.0 y [ R c l ] t / t cc = 0.5wind5.02.50.02.55.0 y [ R c l ] t / t cc = 2.5wind5.02.50.02.55.0 y [ R c l ] t / t cc = 4.5wind5.02.50.02.55.0 y [ R c l ] t / t cc = 6.5wind15 10 5 0 5 10 15 x [ R cl ]5.02.50.02.55.0 y [ R c l ] t / t cc = 8.5wind 10 T [ K ] T [ K ] T [ K ] T [ K ] n [cm ]10 T [ K ] K [ c m g / s ] passivescalarall mass K [ c m g / s ] K [ c m g / s ] K [ c m g / s ] P / k B [Kcm ]10 K [ c m g / s ] M cl,0 dMdlogK [dex ] Figure 8.
The same as Figure 1 except this shows
FastRC-T4e4 , which shows significant cloud growth. In this simulation, µ is a function of n H and T . As in Figure 1, mixing initially drives material from the colder phase to the hotter phase, but rapidcooling slows the transfer rate. Shortly before 4 . t cc , cooling causes this transfer to reverse: the cool phase starts to accretemass. The stable phase diagrams are a manifestation of this growth because there’s a limitless supply of hot phase gas. Thegrowth is even more obvious in the right-column; by 8.5 t cc , the mass of the gas with K < × cm g − / s − has more thandoubled. two separate initial cloud temperatures: T cl = 10 K( FastRC-T1e4 ) and T cl = 4 × K (
FastRC-T4e4 ).Figure 8 depicts the latter case. However, we onlypresent one slow cooling simulation with T cl = 10 K( SlowRC-T1e4 ) because T cl has minimal impact in thisregime. We compare these simulations against the non-radiative simulation NR-X100 .Figure 9 depicts how radiative cooling modifies the P − K phase space distribution for fluid elements orig-inating in the cloud. The key takeaway is that coolingslows the spread of cloud material into the backgroundhigh entropy phase. This suppression is stronger forhigher cooling rates, which reflects the fact that inter- mediate entropy material cools to low entropy prior tomixing with high entropy material.Unsurprisingly, the evolution of our slow cooling caseis minimally changed from the non-radiative case; therate at which gas migrates to the high entropy phaseis slower. Rapid cooling more significantly modifies thedistribution. Consistent with our expectation of entrain-ment, a reservoir of gas is always present at ( p , K cl )throughout the system’s evolution. Interestingly, afteran early transient phase, which has a large scatter in p , the distribution approaches a near steady state. Inthis state, the conditional pressure distributions have re- ixing & Cooling in Cloud-Wind Interactions l o g ( K / K c l ) m ps loss t cool,mix t cc = 0.50 t cc m ps loss t cc m ps loss t cc m ps loss t cc m ps loss t cc m ps loss t cc m ps loss t cc l o g ( K / K c l ) m ps loss T cl = 10 K t cool,mix t cc = 9.0 m ps loss 0% m ps loss 0% m ps loss 0% m ps loss 0% m ps loss 0.7% m ps loss0.00.51.01.52.0 l o g ( K / K c l ) m ps loss T cl = 4 × 10 K t cool,mix t cc = 0.22 m ps loss 0% m ps loss 0% m ps loss 0% m ps loss 0% m ps loss 0% m ps loss10 P / k B ( K cm )0.00.51.01.52.0 l o g ( K / K c l ) m ps loss T cl = 10 K t cool,mix t cc = 0.16 P / k B ( K cm )0% m ps loss 10 P / k B ( K cm )0% m ps loss 10 P / k B ( K cm )0% m ps loss 10 P / k B ( K cm )0% m ps loss 10 P / k B ( K cm )0% m ps loss 10 P / k B ( K cm )0.4% m ps loss 10 l o g M c l , d M p s d l o g P d l o g K [ d e x ] Figure 9.
Mass weighted pressure-entropy evolution for fluid elements originating in the cloud. This is similar to Figure 3except that each simulation has χ = 100 and R cl / ∆ x = 64. Instead, the radiative cooling effectiveness differs between rows.The top row has no cooling, the second row has slow cooling, and the bottom rows have fast cooling. At t ≥ . t cc , thesedistributions are qualitatively similar to the phase distributions that include all gas in the domain at low and intermediate K . duced scatter and a mode that that lies mostly along the p isobar, but has a decrement near log χ ( K/K cl ) ∼ . FastRC-T1e4 , where (cid:96) cool (see Appendix C)is actually resolved.The over-dense diagonal line, in Figure 9, intersecting( p , K cl ) lies along the isotherm corresponding to thecooling curve’s temperature floor. In the absence of thisfloor, the gas would cool to lower K . This is shown inAppendix A.Figure 10 illustrates the bulk motion of the fluid ele-ments originating in the cloud along K . For both fastcooling simulations, it shows a bi-stable medium withlong-lived cold and hot gas. During the early stages ofthe interaction, the colder phase loses mass to the hotterphases, but after some time this reverses. While the ex- act timescale depends on T cl , cooling gradually becomesmore effective at opposing cloud destruction. Eventu-ally, it is effective enough that it not only prevents loss ofadditional mass but also recaptures the lost mass. Thisbehavior manifests over a notably shorter timescale for FastRC-T4e4 .5.2.
Turbulent Radiative Mixing Layer Entrainment
Figure 11 shows how the different cooling regimes af-fect the bulk property evolution of the colder, denserphase (gas with ρ > ρ cl / § Abruzzo et al. t / t cc log ( K / K cl )10 l o g M c l , d M p s d l o g K T cl = 10 K t cool,mix t cc = 9.00 1 2log ( K / K cl )10 l o g M c l , d M p s d l o g K t cool,mix t cc = log ( K / K cl ) T cl = 4 × 10 K t cool,mix t cc = 0.220 1 2log ( K / K cl )10 l o g M c l , d M p s d l o g K T cl = 10 K t cool,mix t cc = 0.16 Figure 10.
The one-dimensional entropy distribution asin Figure 5 except that some of the simulations include ra-diative cooling, and they all have χ = 100, M w = 1 .
5, and R cl / ∆ x = 64. The simulations are in the same order as forFigure 9. The bottom and middle panels depict the evolutionof the velocity and purity fraction (i.e. the cold phasemass fraction of fluid elements initialized in the cloud).The correlation in the evolution of velocity and purityfraction reflects an inelastic collision in the fast cool-ing limit; this is expected for turbulent radiative mixinglayer entrainment (Gronke & Oh 2018; Schneider et al.2020, Tonnesen & Bryan, in prep.). The sustained highpurity fraction signals that a different process, proba-bly ram pressure, dominates acceleration in the non-radiative and weak cooling regime. Note that the minoroffset in the velocity and purity fraction evolution sug-gests that ram pressure could play a subdominant rolein the fast cooling limit.Interestingly, Figure 11 also indicates that the t cool , mix /t cc criterion alone doesn’t fully specify thecloud evolution. In the fast cooling limit, the rate ofcloud growth depends on T cl ; FastRC-T1e4 takes atleast twice as long as
FastRC-T4e4 to show growth de-spite having nearly identical t cool , mix /t cc ratios. Thisdepressed cloud growth in FastRC-T1e4 is accompaniedby a delay in the time at which the cloud is entrained.Figure 12 underscores the significance of this differ-ence in growth. The green and brown curves show themass growth of simulations that are respectively identi- M ( > c l / ) / M c l , t cool, mix / t cc = t cool, mix / t cc = 9.0, T cl = 10 K t cool, mix / t cc = 0.22, T cl = 4 × 10 K t cool, mix / t cc = 0.16, T cl = 10 K P u r i t y F r a c t i o n t / t cc v / / v w i n d Figure 11.
Evolution of total mass (top), purity fraction(middle; fraction of mass originating in the cloud), and av-erage velocity (bottom; ∆ v = v w − v cl ) of cells satisfying ρ > ρ cl / v/v w indicate thatmixing drives acceleration. cal to FastRC-T4e4 and
FastRC-T1e4 , except that theyhave t cool , mix /t cc ∼ The green curve shows nearlyidentical growth to
FastRC-T1e4 (shown in cyan), de-spite the difference in t cool , mix /t cc . As we’ll conclude This difference in t cc is achieved by reducing R cl by a factor offive. ixing & Cooling in Cloud-Wind Interactions t / t cc M ( > c l / ) / M c l , t cool, mix / t cc = 0.22, T cl = 4 × 10 K t cool, mix / t cc = 1.1, T cl = 4 × 10 K t cool, mix / t cc = 0.16, T cl = 10 K t cool, mix / t cc = 1.0, T cl = 10 K Figure 12.
Dependence of cold phase mass growth on T cl . To make the comparison as fair as possible, all datawas measured from simulations with χ = 100, M w = 1 . p/k B = 10 cm − K, and R cl / ∆ x = 16. The magenta andcyan curves are measured from lower resolution versions of FastRC-T4e4 and
FastRC-T1e4 , while the green and browncurves have radii of 1000 pc and 9 pc. Note, we have veri-fied that the brown curve goes to zero (the y-axis starts at0.1). This demonstrates that the t cool , mix /t cc criterion alonedoesn’t fully specify the cloud evolution. below, this difference in growth arises from differencesin t cool /t cc between T cl and T mix (for reference the min-imum t cool /t cc for the green curve is half of that for FastRC-T1e4 ). Moreover, the fact that the brown curvegoes to zero, despite having a comparable t cool , mix /t cc tothe green curve, illustrates that this difference can evenmodify the survival cloud survival criterion.To interpret these results, we consider them in termsof our mixing model. For each cooling case, Fig-ure 13 compares standalone non-radiative ˙¯ K mixing ( K, t )measurements and the ˙ K cool ( K ) prediction against the˙¯ K total ( K, t ) measurements from the simulations with ra-diative cooling.Given our result from § K/ ˙ K mixing ( K mix , t ) ∼ t cc for most times when ( dM ps /dK )( K mix ) >
0, the sur-vival criterion t cc > t cool , mix can be directly visualizedin terms of this plot. Turbulent radiative mixing layerentrainment is expected when − ˙ K cool exceeds the char-acteristic value of ˙ K mixing at K mix = χ / K cl . Given˙ K mixing , char ( K ) ∼ K/t cc and | ˙ K cool | ’s mostly inversedependence on K for realistic ISM conditions, satisfac-tion of the survival criterion basically guarantees that˙ K mixing , char ( K ) < ˙ K cool ( p , K ) for K ∈ [ K cl , K mix ].Thus, we predict a negative ˙ K total ( K, t ) over that in-terval, and, by extension, entrainment.The outcome of the slow cooling case ( t cool , mix /t cc ∼
10) is clear-cut. Because | ˙ K mixing , char ( K ) | (cid:29)| ˙ K cool ( p , K ) | , ˙ K total ( K, t ) resembles ˙ K mixing ( K, t ) andthe cloud is destroyed. Likewise, the ultimate fates in K m i x i n g [ K c l t cc ] t cool,mix t cc = T cl = 10 K t cool,mix t cc = 9.010 K m i x i n g [ K c l t cc ] t cool,mix t cc = T cl = 4 × 10 K t cool,mix t cc = 0.220 1 2log ( K / K cl )10 K m i x i n g [ K c l t cc ] t cool,mix t cc = 0 1 2log ( K / K cl ) T cl = 10 K t cool,mix t cc = 0.16 0.52.03.55.06.58.09.511.012.514.0 t / t cc Figure 13.
Comparison of passive scalar mass weightedand time averaged ˙ K ( K ) measured with and without radia-tive cooling for three sets of physical conditions with χ = 100, M w = 1 . R cl / ∆ x = 64. The left panels compare mea-surements for just adiabatic mixing ( ˙ K mixing ( K )) against thepredicted contributions from cooling ( ˙ K cool ( K )), while the right panels depict measurements involving both mixing andcooling. The red dashed line shows − ˙ K cool ( K ) for the ini-tial pressure and the shaded regions show variations from 0.5dex pressure perturbations. The other curves are computedfrom the histograms in Figure 10 (only one set of ˙ K mixing ( K )curves are shown). Curves are omitted when at least 1% ofthe fluid elements originating in the cloud have left the do-main. Curve segments adjacent to bins with under 1% of theinitial cloud mass are transparent. the fast cooling cases are also predictable. However, thedetailed shape and temporal evolution of ˙ K total ( K ) isless straightforward.The ˙ K total ( K, t ) evolution for the cool cases directlyreflects the discussion from the previous section ( § K total ( K, t ) ispositive because radiative cooling is unable to preventinitial mixing of the phases. As the process continues,˙ K total ( K, t ) gradually decreases with time, which indi-cates that cooling becomes more effective at combatingmixing. Eventually, the opposition of cooling to mix-ing becomes so effective that it reverses the transfer ofgas between phases, which causes ˙ K total ( K, t ) to becomenegative.6
Abruzzo et al. T / T mix l o g t c oo l / t c oo l , m i x T cl = 10 K T cl = 4 × 10 KBPL mix, cl = 1BPL mix, cl = 2BPL mix, cl = 6BPL mix, cl = 60 1 0 1210 l o g ( T ) / ( T m i x ) K / K cl (for broken power law curves) Figure 14.
Comparison of the truncated cooling curves,measured at p/k B = 10 K cm − and normalized by theirproperties at T mix , used in our ( χ = 100) radiative coolingsimulations. The upper axis denotes the entropy for just thebroken power-law cooling curves, which assume a fixed µ of0.6. In all displayed curves, K cl (= 0 . γ K mix ), K mix , and K w (= 10 γ K mix ) always coincide with T cl , T mix , and T w , bydefinition. The gray dashed (dotted) lines denote the relativelocations of T cl and T w when T cl = 10 K (4 × K). Foreach broken power law cooling function, T cl = 0 . T mix and T w = 10 T mix . It is around this time that our measurements of˙ K total ( K, t ) lose meaning, since an increasing fractionof the colder phase is composed of gas originating in thehot phase. This is discussed in further in Appendix C.Nevertheless, we have included the measurements be-cause they illustrate, if imprecisely, the expected behav-ior.We now consider why
FastRC-T4e4 begins rapidgrowth in roughly half the time as
FastRC-T1e4 . Fig-ure 13 suggests that this difference derives from the localshape of the cooling curve. The major difference is that | ˙ K mixing ( K ) / ˙ K cool ( K ) | is roughly an order of magnitudesmaller in 0 . (cid:46) log χ ( K/K cl ) (cid:46) . FastRC-T4e4 .Phrased another way, the value of t cool near the entropy(or temperature) of the colder phase appears to deter-mine how rapidly the cloud grows.5.3. Cooling Curve Variations
Motivated by the impact of cooling curve shape forfixed t cool , mix /t cc ratios, in this subsection, we more sys-tematically investigate how simple variations in the local K / K cl )10 M c l , d M p s d l o g K T cl = 4 × 10 K T cl = 10 K mix, cl = 1 mix, cl = 2 mix, cl = 6 mix, cl = 60 Figure 15.
Late time entropy distributions of fluid ele-ments originating in the cloud for simulations run with thecustom broken power-law cooling functions (see Figure 14).The shaded region depicts entire temporal variation during t/t cc ∈ [10 , .
5] and the colored lines denote the median.For comparison, distributions are also shown for R cl / ∆ x = 8versions of FastRC-T4e4 and
FastRC-T1e4 . At least 99% ofthe passive scalar remain in the domain through 13 . t cc innearly all cases. However, this is only true for BPLawRC-1 ,through 11 t cc . shape of the cooling curves near T cl (or K cl ) modify theonset of growth.To approximately match the realistic cooling curve,but allow us to the change the lower section in a system-atic way, we model the cooling-curve with a three-piecebroken power law. The shape at low T is controlledby a single parameter: η mix , cl = t cool , mix /t cool , cl . For T = T /T mix , the cooling time is given by t cool ( T ) t cool , mix = η − , cl T ≤ η . , cl T η . , cl < T ≤ α − / α T . α − / < T . (7)The value of α is 0.138259 and is set by the intersectionof the middle and upper segments . For simplicity, we The upper power law segment mimics properties of the cool-ing curve used for
FastRC-T4e4 . They share the same( t cool ( T ) /t cool , mix ) at T = 10 . and have comparable power-law slopes from there to T = 10 . . ixing & Cooling in Cloud-Wind Interactions M ( > c l / ) / M c l , mix, cl = 1 mix, cl = 2 mix, cl = 6 mix, cl = 60 t / t cc v / / v w i n d Figure 16.
Like Figure 11 except that the data is forsimulations using custom power law cooling function (seeFigure 14) with χ = 100, M w = 1 .
5, and R cl / ∆ x = 8. Allsimulations have the same t cool , mix /t cc , but show significantlydifferent evolution, demonstrating the importance of coolingbelow T mix . define Λ( T ) such that the above equation is satisfied atall pressures and for constant µ .We run four simulations which resemble the fast cool-ing simulations, but use constant µ and employ brokenpower law cooling functions with η mix , cl ∈ [1 , , , BPLawRC-1 , BPLawRC-2 , BPLawRC-6 , BPLawRC-60 ). Forconsistency with earlier simulations, the cooling func-tions are truncated below T cl and above 0 . T w . Fig-ure 14 illustrates the shapes of the cooling curves nor-malized by the properties at the mixing layer.Figure 15 shows the dependence of the late-timepassive scalar mass-entropy profiles on η mix , cl . Thehigh η mix , cl simulations resemble FastRC-T4e4 , and as η mix , cl decreases, they begin to more closely resemble FastRC-T1e4 . The temporal stability of these profilesindicates that each simulation has a bistable medium.However,
BPLawRC-1 ’s cold phase may not yet be stable;its peak near K cl ( K w ) decreases (increases) by a factor Ordinarily, µ drops by ∼
21% between 10 K and 4 × K.Above that, it only drops by ∼ FastRC-T1e4 , we confirmed that these µ variations have negligi-ble impact on mass evolution. This simulation used R cl / ∆ x = 8,a fixed µ , a re-scaled Λ( ρ, T ) such that t cool ( ρ, T )’s shape is un-changed, and a R cl that maintained t cool , mix /t cc = 0 . K m i x i n g [ K c l t cc ] t cool = mix,cl = 110 K m i x i n g [ K c l t cc ] t cool = mix,cl = 210 K m i x i n g [ K c l t cc ] t cool = mix,cl = 60 1 2log ( K / K cl )10 K m i x i n g [ K c l t cc ] t cool = 0 1 2log ( K / K cl ) mix,cl = 60 0.52.03.55.06.58.09.511.012.514.0 t / t cc Figure 17.
Like Figure 13 except that the measured ˙ K curves were computed from the histograms shown in Fig-ure 15. The histograms were measured from simulations us-ing custom power law cooling function (see Figure 14) with χ = 100, M w = 1 .
5, and R cl / ∆ x = 8. of ∼ t cc and 14 t cc . Note that FastRC-T4e4 ’sprofiles in Figure 10 suggest that for high η mix , cl cases,there may be more variability at intermediate K athigher resolutions.Figure 16 establishes a clear trend: as η mix , cl increasesand the break in the cooling curve moves to lower T ,rapid cloud growth sets in more quickly. This figure alsosupports our expectation that BPLawRC-1 is just startingto grow in mass between 11 t cc and 14 . t cc . The over-taking of the mass and velocity growth in BPLawRC-60 by BPLawRC-6 may be a resolution effect. Figure 22 ofAppendix C shows that velocity evolution is particularlysensitive to resolution.Our results clearly demonstrate that while t cool , mix isimportant for identifying the conditions under whichcooling occurs, the shape of the cooling curve below T mix is important for determining when rapid growthcommences.8 Abruzzo et al.
To gain some intuition for why these different cool-ing times matter, we examine the ˙¯ K total ( K, t ) measure-ments in Figure 17. Unsurprisingly, these measurementsresemble the fast cooling simulations; it’s most obvi-ous when considering measurements from simulations ofcomparable resolution (see Figure 24 from Appendix C).As expected, the η mix , cl ≤ T ) resemble FastRC-T1e4 , while the η mix , cl ≥ FastRC-T4e4 . Comparing therelative magnitudes of ˙¯ K mixing ( K, t ) and ˙ K cool ( p , K ),in the middle left two panels provide some intuition forwhy there is such a big difference between η mix , cl = 2and η mix , cl = 6.However, the precise property of the cooling curve onthe interval T cl (cid:46) T (cid:46) T mix that controls the growthrate remains somewhat ambiguous. We speculate thatthe crucial quantity is some kind of (weighted) averageover the interval, possibly related to (although not ex-actly equal to) t cool , cl or min t cool , and hereafter refer toit as the characteristic cooling time of the cold phase˜ t cool , cl . Regardless of ˜ t cool , cl ’s true nature, Figure 2clearly illustrates that it must be considerably smallerfor FastRC-T4e4 than it is for
FastRC-T1e4 . Therefore,rapid growth commences more quickly in
FastRC-T4e4 .The onset of rapid growth may coincide with the tran-sition between the “tail growth” and “entrained phases”of the cloud’s areal growth (Gronke & Oh 2020a). Be-cause areal growth is more rapid in the earlier phase,this transition likely corresponds to the point when thesystem is able to reach equilibrium and cooling is ableto balance the destructive mixing effects (Fielding et al.2020). If true, then the delayed transition may implythat a larger area is required when η mix , cl is smaller.Such an interpretation would be consistent with the mix-ing layer being linked to both ˜ t cool , cl and t cool , mix , ratherthan just the latter. DISCUSSION6.1.
What does our model offer?
Our entropy evolution mixing model provides threemain benefits. First, it offers a condition-agnosticmethod for the characterization and quantificationof the cloud-wind interaction’s evolution. Our non-radiative parameter study showed that the model mean-ingfully captures the processes of cloud destruction. Foridealized, hydrodynamic interactions, it reproduces thewell-known destruction time-scale t cc . However, the in-dependence of these characterizations with respect tothe interaction’s physical conditions warrants emphasis.Our model offers a robust approach for describingcloud destruction in circumstances where the t cc descrip-tion breaks down. We have already demonstrated that it quantitatively captures the well-documented effectsthat magnetic draping have on extending the cloud’slifetime (e.g Dursi & Pfrommer 2008; McCourt et al.2015; Banda-Barrag´an et al. 2018; Gronke & Oh 2020a).Another interesting application might be characterizingthe evolution of networks of small clouds where the ideaof having a monolithic cloud with a well-defined R cl doesnot really apply.Second, our model facilitates comparisons betweenthe effects of radiative cooling and empirical charac-terizations of other effects that affect the system’s evo-lution. We describe a complementary relationship tothe t cc < t cool , mix survival criterion from Gronke & Oh(2018) at length in § § T cl = 10 K case. However, our mixing model re-vealed that these factors are sensitive to the character-istic cooling time of the colder phase ˜ t cool , cl . We deferdiscussions of this finding’s significance to § p at eachvalue of K , the model is conducive to layering addi-tional quantities atop p − K space; one could imagineconstructing manifolds in higher dimensional space. Forexample, one could supplement p − K space with thewind-aligned velocity (similar to Schneider & Robert-son 2017; Kanjilal et al. 2020). This would also connectour model to the established relation between a fluid ele-ment’s wind-aligned velocity and the fraction of its massthat originated in the hot phase (e.g. Melso et al. 2019;Schneider et al. 2020, Tonnesen & Bryan, in prep).As mentioned earlier, it would also be useful to con-sider the entropy flow for fluid elements originating inthe hotter phase in addition to the fluid elements fromthe colder phase.6.2. Limitations and Missing physics
The omission of geometric information may be a limi-tation of our model. Consider a non-radiative hydrody- ixing & Cooling in Cloud-Wind Interactions K mixing , char ( K ) and thus pre-dict stricter conditions for turbulent radiative mixinglayer entrainment. In reality, Gronke & Oh (2020a)showed that turbulent clouds not only survive under thesame conditions as spherical clouds, but initially growfaster because they have larger surface areas. Addition-ally, it’s unclear how well the model captures the evolu-tion of a system in which each phase has different levelsof non-thermal pressure support.We note that these proposed limitations are entirelyhypothetical. Simulations are needed to assess whetherthere are actually issues in these scenarios. Regardless,we are unaware of any alternative models with similarpredictive power that are devoid of these issues.This work entirely neglected relevant physical effectslike viscosity, conduction, and cosmic rays. It also didn’tconsider magnetic fields at the same time as radiativecooling. Additionally, we artificially prevented coolingbelow T cl , which appears to have a large impact onentrainment (c.f. Appendix A; Gronke & Oh 2018).Moreover, we only considered idealized initial condi-tions. The influence of metallicity variations, differentmagnetic field configurations, and the presence of tur-bulence warrant attention in future work. However, weemphasize that our model is well-equipped for character-izing how each of these conditions modify the conditionsfor turbulent radiative mixing layer entrainment.6.3. Turbulent Radiative Mixing Layer Entrainment
Our mixing model is conducive to applications re-lated to cloud survival through turbulent radiative mix-ing layer entrainment. It naturally provide a generalcondition under which this entrainment mechanism isexpected (i.e. ˙ K mixing ( K, t ) + ˙ K cool ( p , K ) (cid:28) K ∈ [ K cl , K w ], see § K mixing ( K, t )from non-radiative simulations to predict the interac-tion’s fate. Therefore, our mixing model complementssuch criteria, and can be used to help improve them.6.3.1.
Relevant timescale
Our results in § T cl through T mix . While the Gronke & Oh (2018) sur-vival criterion, which is based on t cool , mix , appears toaccurately predict cloud survival, it doesn’t fully spec-ify the interaction’s evolution. Specifically, the charac-teristic cooling time of the colder phase, ˜ t cool , cl , affectshow quickly rapid cloud growth commences (the delayfrom the start of the simulation appears correlated with˜ t cool , cl /t cool , mix ).The delay is significant because it provides additionalopportunity for other processes (e.g. externally-driventurbulence in the wind) to destroy the cloud. This raisesa broader point. Although the distinction between cloudsurvival and destruction is of primary interest, knowinghow close a surviving cloud comes to being destroyed(or how quickly rapid growth commences) would be in-sightful. We discuss how the delay in rapid growth mayaffect the prevalence of turbulent radiative mixing layerentrainment in § t cool , w , which un-derlies the Li et al. (2020) and Sparre et al. (2020)criteria, is not the dominant cooling time-scale. Fig-ure 3 of Gronke & Oh (2018) and Figure 20 in Ap-pendix A show that switching cooling on and off above ∼ . T w has minimal effect on the mass evolution for FastRC-T4e4 and
FastRC-T1e4 , respectively. The mainconsequence of the wind cooling is that the system’sequilibrium pressure drops (by (cid:46)
30% for
FastRC-T4e4 and
FastRC-T1e4 ), which does not appear to be signif-icant. While this does affect the cooling function, wedon’t expect it to be significant in most cases (see Ap-pendix A for further discussion).The fact that the difference in ˜ t cool , cl between FastRC-T1e4 and
FastRC-T4e4 so efficiently accountsfor the variations in mass growth rates reinforces ourconclusion that t cool , w is not the dominant time-scale.If t cool , w were dominant, we would expect it to explainthe difference.Finally, we address Sparre et al. (2020)’s proposed ex-planation for why t cool , w could be important to turbulentradiative mixing layer entrainment: they suggest that afluid element’s temperature evolution from T w to T cl israte-limited by an initial cooling phase near T w , set by t cool , w . While we acknowledge that initial cooling of thewind could possibly make entrainment easier, we expectthis to be high-order effect. In fact, the small impactthat switching cooling on and off above 0 . T w has onthe cold phase mass-growth suggests that the tempera-ture change at high T is dominated by mixing.6.3.2. Comparison with prior work
Our results are largely consistent with Gronke & Oh(2018, 2020a) and Kanjilal et al. (2020), but they need0
Abruzzo et al. to be reconciled with those of Li et al. (2020) and Sparreet al. (2020). Most works primarily considered a set of“typical” conditions with χ ≥ T cl = 10 K, andan p/k B = 10 K cm − . Li et al. (2020) and Sparreet al. (2020)’s results both support survival criteria thatrequire a larger minimum survival radius under theseconditions than the Gronke & Oh (2018) survival crite-rion. We largely attribute this difference to a combina-tion of choices, which include truncation of the coolingfunction, the simulation box size, and the standards foridentifying destroyed clouds.Under these “typical” conditions, every survival cri-terion implicitly requires that t cc (cid:38) t cool , cl . Becausegrowth takes a few t cc to develop, unrestricted coolingwill cause clouds to initially contract. Like Gronke &Oh (2018, 2020a) and Kanjilal et al. (2020), we explic-itly prevent gas below T cl from cooling. In contrast,Sparre et al. (2020) allows cooling down to T cl / × K. Because this contraction impedes growth(see Appendix A) and might make Sparre et al. (2020)’s M w ∼ . thisdifference may help to reconcile our results. However,Sparre et al. (2020)’s inclusion of magnetic fields couldplausibly inhibit these effects. While Li et al. (2020)didn’t truncate their cooling curve, their results prob-ably aren’t strongly affected because they allow theirinitial conditions to equilibrate before introducing thevelocity difference.Kanjilal et al. (2020) highlight two choices made byLi et al. (2020) and Sparre et al. (2020) that may fur-ther help to reconcile our results. First, Kanjilal et al.(2020) argue that Li et al. (2020)’s small box-size maycause misclassification of growing clouds. While plausi-ble, we note that Li et al. (2020) claimed to have verifiedtheir conclusions with longer boxes. Relatedly, the re-flection of shocks off of the transverse boundaries, likewe encountered for our non-radiative simulations, mayintroduce some artificial shock heating in Sparre et al.(2020)’s χ = 1000 simulations. However, it remains un-clear how significant this artifact is in simulations withcooling.Second, both studies choose standards for cloud sur-vival that implicitly place requirements on when growthcommences. Li et al. (2020) and Sparre et al. (2020) There are a few notable exceptions. Gronke & Oh (2018, 2020a)considered clouds with T cl = 4 × K. Li et al. (2020) considereda larger range of pressures and χ values. For χ (cid:38)
100 clouds with R cl > (cid:96) cool , Gronke & Oh (2020b)show that tripling the ρ cl ’s initial density from cooling-drivencontraction cause shattering. They further argue that the shockfrom oncoming winds with M w (cid:38) . T cl ). identify surviving clouds in cases where cloud growthcauses the total cold phase mass to never fall below10% of the initial mass and to have a positive derivativeat t = 12 . t cc , respectively. Thus, they classify cloudsdifferently that survive, but come closer to being de-stroyed (we expect the differences to be minimal for Liet al. 2020). Furthermore, many of Sparre et al. (2020)’s M w = 0 . M w = 1 . T mix . While we generally believe that t cool , mix and ˜ t cool , cl are the most important time scales for theprocess of entrainment, it’s plausible that process ofthermal conduction could increase the importance of t cool , w . 6.3.3. Prevalence
Our result that rapid cloud growth takes longer tocommence at larger ˜ t cool , cl /t cool , mix has important im-plications for the prevalence of turbulent radiative mix-ing layer entrainment. This effect is most pertinent forclouds at or near thermal equilibrium, where ˜ t cool , cl islargest. To give a concrete example with a realistic cool-ing curve, consider a system at p = 10 k B Kcm − inwhich M w = 1 . T cl ∼ × K and χ (cid:38) For such clouds, the delay in growth provides addi-tional opportunity for destructive processes (like mix-ing) to destroy the cloud. Therefore, these clouds re-quire a minimum survival radius that is somewhat largerthan the Gronke & Oh (2018) criterion predicts. Fig-ure 12 illustrates this effect for an idealized, initiallylaminar wind. However, this delay may be even morerelevant for clouds embedded in winds with turbulencedriven by external processes (e.g. supernovae) becausemixing may more efficiently destroy clouds In this sce- They each identify the cold phase with the density threshold √ ρ cl ρ w rather than ρ cl /
3. We don’t expect small differences inthresholds to make a significant difference. We conservatively chose this lower bound on χ to ensure thatthe equilibrium pressure drop from cooling above 0 . T w has nomore influence on cloud growth than it does for FastRC-T1e4 (seeAppendix A). Note, magnetic fields might directly mitigate this. Banda-Barrag´an et al. (2018) showed that magnetic field have a stabi-lizing effect on turbulent clouds in a laminar wind. They mightplausibly have a similar impact in a wind with externally-driventurbulence. ixing & Cooling in Cloud-Wind Interactions K floor (Schneider & Robertson 2018)make the delayed growth and entrainment, from large˜ t cool , cl /t cool , mix , relevant (as in FastRC-T1e4 ). We ex-pect that the combination of external turbulence anddelay potentially impedes growth near the galaxy. Atlarger radii, the hot phase’s ∼
10 times larger pressurethan the cold phase (in the simulation) may further ex-acerbate the effect. If the intermediate phase also has asomewhat elevated pressure, then ˜ t cool , cl /t cool , mix shouldbe larger because t cool has an inverse dependence onpressure for photo-ionized gas. Although Gronke & Oh(2020a) showed the rapid growth in expanding winds,their model explicitly assumes that the cold and hotphases are in sonic contact. CONCLUSIONWe have presented an entropy-based formalism for in-terpreting the cloud-wind interaction’s evolution. Thebasic premise of the approach is that information aboutthe system’s state is encoded in the evolution of its ther-modynamic phase space. We consider p − K phase spaceto take advantage of the system’s quasi-isobaric nature,and the conservation of a fluid element’s specific entropyin the absence of irreversible processes (like shocks, mix-ing, and heating/cooling). Thus, in the adiabatic limit,the gas distribution along K is primarily governed bythe history of mixing, the dominant cloud destructionprocess.We leverage the fact that mixing occurs in K -space tointroduce an empirical mixing model. We characterizemixing with the average rate of change in the entropyof fluid element’s originating in the cloud, ˙ K mixing ( K, t ).From this knowledge, we can define a mixing timescale, K/ ˙ K mixing ( K, t ) as a function of K and t . Additionally,the model provides the capability for making predic-tions about how radiative cooling will modify adiabaticmixing by facilitating comparisons of ˙ K cool ( p , K ) = K/t cool ( p , K ), with measurements of ˙ K mixing ( K, t ).We have considered two example applications, using enzo-e simulations, to demonstrate that this mixingmodel works as expected and provides useful insight.We enumerate our four main results below: 1. The timescale of cloud destruction from adiabaticmixing is well characterized by K/ ˙ K mixing ( K, t )for most entropy values ranging from the initialvalue in the cloud, K cl , to the initial value inthe wind. In fact, K/ ˙ K mixing ( K, t ) is comparableto t cc for idealized, non-radiative, hydrodynamicalinteractions.2. In addition, the model can characterize the changein destruction rate due to other physical pro-cesses. For example, we have demonstrated thatthe model reflects the reduction in the cloud de-struction rate from the tangling of initially trans-verse magnetic fields.3. These characterizations are well suited for compar-isons against the effects of radiative cooling. Ananalogous form of the Gronke & Oh (2018) sur-vival criterion, t cc > t cool , mix , can be formulatedin terms of this model.4. We used our model to show that the local shape ofthe cooling curve can influence the process of cloudentrainment via the rapid cooling of gas that hasmixed with the wind. Independent of the coolingtime at the mixing layer, variations in the char-acteristic cooling time of the cold phase can morethan double the elapsed time required for cloudsto commence rapid growth and become entrained.ACKNOWLEDGMENTSThe authors are grateful to James Bordner, MikeNorman, and the other enzo-e developers. Weprimarily performed simulations and analysis usingthe NSF XSEDE facility. GLB acknowledges fi-nancial support from the NSF (grant AST-1615955,OAC-1835509, AST-2006176) and computing supportfrom NSF XSEDE and the Texas Advanced Com-puting Center (TACC) at the University of Texas atAustin. We also acknowledge computing resources fromColumbia University’s Shared Research Computing Fa-cility project, which is supported by NIH Research Fa-cility Improvement Grant 1G20RR030893-01, and as-sociated funds from the New York State Empire StateDevelopment, Division of Science Technology and In-novation (NYSTAR) Contract C090171, both awardedApril 15, 2010. Software: yt (Turk et al. 2011), Launcher Utility(Wilson & Fonner 2014),2
Abruzzo et al. t / t cc M ( > c l / ) / M c l , t cool,mix / t cc = 0.22 T cl = 4 × 10 K t / t cc t cool,mix / t cc = 0.16 T cl = 10 K fiducial: no cool T T cl & T T w unrestricted coolingno cool T T cl no cool T T w Figure 18.
Comparison of how cooling curve restrictions affect the mass evolution of the cold phase (gas with ρ > ρ cl /
3) forlow resolution versions ( R cl / ∆ x = 8) of both FastRC-T4e4 (left) and
FastRC-T1e4 (right).
APPENDIX A. SIMULATION ROBUSTNESSA.1.
Cooling Curve Restrictions
Figure 18 illustrates how switching cooling on and off at temperatures above ∼ . T w and below T cl , affects the coldphase mass evolution for FastRC-T4e4 and
FastRC-T1e4 . Switching cooling on and off above ∼ . T w is relativelyinsignificant. However, when cooling is allowed below T cl , the initial period of mass loss is reduced and is followed bya period of slower growth. These results are consistent with Gronke & Oh (2018), who argue that the cooling below T cl affects mass evolution because it causes the cloud to contract.Figure 19 shows how the difference between the fiducial restricted cooling curves (used in the bulk of this work) andthe unrestricted cooling curves impact the p − K evolution for all gas in FastRC-T4e4 and
FastRC-T1e4 . There aretwo main differences. First, cooling below T cl , drives mass to much lower K (or equivalently, T ). Second, cooling ofthe wind slightly decreases (by (cid:46) FastRC-T4e4 and
FastRC-T1e4 , when the cool-ing curves are unrestricted. Because t cool , w /t cc for FastRC-T1e4 is ∼
46 and is ∼ . FastRC-T4e4 ,allowing cooling above ∼ . T w causes a slightly larger pressure drop and has slightly more affect on mass growth for FastRC-T1e4 . At the same time,
FastRC-T4e4 has a significantly larger density increase than
FastRC-T1e4 ; by t/t cc = 2 . ∼ ρ cl and ∼ ρ cl , respectively. This makes sense given thatthe former’s t cc is both ∼
39 times larger than the latter’s t cc and (cid:38)
100 times larger than the time needed for gas tocool (whether the cooling proceeds isobarically or isochorically) between the respective T cl .Interestingly, when cooling is unrestricted, FastRC-T4e4 has faster cold phase growth than
FastRC-T1e4 , eventhough its gas contracts more. Although increased resolution may affect growth rate, this suggests that the earlytime value of ˜ t cool , cl (see §
5) may be more important for setting the properties of cloud growth. Future work mustinvestigate cloud evolution using the full cooling curve.A.2.
Influence of Hydrodynamical Integrator
Figure 20 illustrates how differences in the hydrodynamical integrator modify the cold phase mass evolution. Specif-ically we compare the VL+CT integrator (which is used in the rest of this paper) with the ppm integrator, which waspreviously ported from enzo (Bryan et al. 2014). The defining differences are that the ppm solver is dimensionallysplit and employs third order spatial reconstruction. Additionally, while the VL+CT integrator employs a predictor-corrector scheme, the ppm solver updates the grid in a single pass. More minor differences include implementationchoices for the dual energy formalism and the choice of the HLLD (Two-Shock) Riemann solver for our simulationswith VL+CT (ppm). When magnetic fields are zero the HLLD solver reduces to anHLLC solver. ixing & Cooling in Cloud-Wind Interactions l o g ( K / K c l ) T cl = 4 × 10 Kno cool
T T cl & T T w t cc t cc t cc t cc t cc l o g ( K / K c l ) T cl = 4 × 10 Kunrestrictedcooling l o g ( K / K c l ) T cl = 10 Kno cool
T T cl & T T w P / k B ( K cm )101 l o g ( K / K c l ) T cl = 10 Kunrestrictedcooling P / k B ( K cm ) 10 P / k B ( K cm ) 10 P / k B ( K cm ) 10 P / k B ( K cm ) 10 l o g M c l , d M d l o g P d l o g K [ d e x ] Figure 19.
Comparison of how restricting cooling affects the mass weighted pressure-entropy evolution. The top row andsecond row from the bottom are similar to the bottom two rows of Figure 9. The two differences are that this figure showsthe distribution of all fluid elements in the simulation domain (not just the fluid elements originating in the cloud) and thesimulation resolution is lower ( R cl / ∆ x = 8). The other rows show simulations with same initial conditions, but have unrestrictedcooling. The gray region indicates wherever t cool > t cc or heating dominates. We primarily consider how the different integrators affect the mass growth in
FastRC-T4e4 and
FastRC-T1e4 whenusing our standard restricted cooling curves (bottom row of Figure 20). The simulations using the PPM integratoreach show elevated mass growth rates, compared to the VL+CT simulations. However, we find solace in the way thatthe simulations using the PPM integrator appear to trend towards the converged curves (discussed in Appendix C) aswe increase resolution.We also compare how the difference in integrators affect the mass growth when cooling is unrestricted (top row ofFigure 20). Interestingly, the
FastRC-T1e4 mass evolution has slower growth when using the PPM curve. Nevertheless,simulations with both integrators indicate that
FastRC-T4e4 shows faster growth than
FastRC-T1e4 . Our resultssuggest that while the precise values measured in simulations with different integrators may differ, the trends betweenthe values measured in different simulations (with a single integrator) are fairly robust. B. FRAME TRACKING SCHEME COMPARISONTwo different developmental versions of enzo-e were employed in this work. The main difference between them isin the reference frame tracking scheme. The earlier version (used for simulations with R cl / ∆ x <
64) has a bug which4
Abruzzo et al. C oo li n g C u r v e R e s t r i c t i o n M ( > c l / ) / M c l , t cool,mix / t cc = 0.22 T cl = 4 × 10 K t cool,mix / t cc = 0.16 T cl = 10 K fiducial: no cool T T cl & T T w no cool T T cl & T T w , PPM integratorunrestricted coolingunrestricted cooling, PPM integrator t / t cc V a r y i n g R e s o l u t i o n M ( > c l / ) / M c l , R cl / x =8 R cl / x =16 t / t cc fiducial: no cool T T cl & T T w no cool T T cl & T T w , PPM integrator Figure 20.
Comparison of how different integrators affect the mass evolution for the cold phase (gas with ρ > ρ cl /
3) for lowerresolution versions of both
FastRC-T4e4 (left column) and
FastRC-T1e4 (right column). The top (bottom) row further showshow the the mass evolution in each integrator is affected by cooling curve restrictions (differences in resolution). All simulationsin the top row have a fixed resolution of R cl / ∆ x = 8. We note that all that the solid orange and blue curves are identical tothe curves in Figure 18. l o g ( K / K c l ) Broken FrameTracking 0.50 t cc t cc t cc t cc t cc t cc t cc P / k B ( K cm )0.00.51.01.52.0 l o g ( K / K c l ) Fixed FrameTracking P / k B ( K cm ) 10 P / k B ( K cm ) 10 P / k B ( K cm ) 10 P / k B ( K cm ) 10 P / k B ( K cm ) 10 P / k B ( K cm ) 10 l o g M c l , d M d l o g P d l o g K [ d e x ] Figure 21.
Comparison of a cloud-wind interaction’s mass weighted pressure-entropy evolution simulated with two versionsof enzo-e . These simulations’ initial conditions are the same as those for bottom row of Figure 9, except that R cl / ∆ x = 16.Unlike Figure 9, this figure shows the distribution of all fluid elements in the simulation domain. The simulation in the toprow (like all other simulations with R cl / ∆ x <
64) effectively has no frame tracking while the bottom row has frame tracking (aslightly refactored Riemann Solver). The differences between these simulations are minimal. only allows the frame velocity to be updated once, immediately after the very first update cycle. Thus, the framevelocity remains near zero.The bug is fixed in the later version (used for simulations with R cl / ∆ x = 64). Every 0 . t cc , the frame velocity(measured in the cloud’s initial reference frame) is updated to match the minimum velocity of cells with a passivescalar density of at least ρ cl / ixing & Cooling in Cloud-Wind Interactions t / t cc M ( > c l / ) / M c l , t cool, mix / t cc = t cool, mix / t cc = 9.0, T cl = 10 K t cool, mix / t cc = 0.22, T cl = 4 × 10 K t cool, mix / t cc = 0.16, T cl = 10 K t / t cc P u r i t y F r a c t i o n R cl / x = 8 R cl / x = 16 R cl / x = 32 R cl / x = 64 t / t cc v / / v w i n d Figure 22.
Comparison of how resolution affects the evolution of total mass, purity fraction, and average velocity. Theevolution clearly converges at high resolution.
This strategy was selected to ensure that the bow shock remained in the simulation domain . The later version ofthe code also features a more efficient Riemann Solver implementation.To assess how the code differences affect our results, we compare the results of two simulations using a single set ofinitial conditions but with the different code versions. We used initial conditions matching FastRC-T1e4 , but with aresolution of R cl / ∆ x = 16. Figure 21 provides a comparison of the phase space evolution for each simulation. It isclear that both versions of the code produce consistent results. C. RESOLUTION STUDYIn this appendix we briefly assess how resolution affects our measurements of the cloud wind-interaction. To dothis we consider the primary four initial conditions discussed in section § NR-X100 , SlowRC-T1e4 , FastRC-T4e4 , FastRC-T1e4 ) at the resolutions R cl / ∆ x = { , , , } .C.1. Relevance of Shattering
A relevant length scale for our convergence study is the so-called “cooling length”, (cid:96) cool ∼ min( c s t cool ) (McCourtet al. 2018). McCourt et al. (2018) first showed in 2D simulations that large clouds with sizes exceeding (cid:96) cool are More aggressive strategies exist that both satisfy this criterionand increase the time that the cold phase remains in the domain Abruzzo et al. l o g M c l , d M p s d l o g K t cool,mix t cc = R cl / x = 8 R cl / x = 16 R cl / x = 32 R cl / x = 64 l o g M c l , d M p s d l o g K T cl = 10 K t cool,mix t cc = 9.0 l o g M c l , d M p s d l o g K T cl = 4 × 10 K t cool,mix t cc = 0.22 K / K cl )10 l o g M c l , d M p s d l o g K T cl = 10 K t cool,mix t cc = 0.16 K / K cl ) 0 1 2log ( K / K cl ) 0 1 2log ( K / K cl ) 0.52.54.56.58.510.512.514.5 t / t cc Figure 23.
Comparison of dM ps /d log χ K evolution with respect to resolution. Different rows illustrate different simulationsand resolution varies between columns. prone to fragmenting into a swarm of cloudlets of size ∼ (cid:96) cool . Thus, simulations where clouds “shatter” may not havewell-converged properties when (cid:96) cool is not converged. For reference, Table 1 lists each simulation’s (cid:96) cool . The lengthscale is well resolved at all resolutions of SlowRC-T1e4 and barely resolved ( (cid:96) cool = 1 . x ) for FastRC-T1e4 when R cl / ∆ x = 64. However, it’s not resolved in any other simulations with cooling.More recently, Gronke & Oh (2020b) considered cloud shattering in 3D simulations and linked the shattering ofclouds to cloud growth through cooling. They demonstrated that all clouds with R cl > (cid:96) cool that are over pressurizedcompared to the ambient medium and have density inhomogeneities undergo some degree of shattering. However, thecloud’s fate depends on how the density contrast (when the cloud is over-pressurized), compared to χ crit ∼ χ crit , the cloud breaks apart. Otherwise, the cloud re-coagulates, and has the opportunity toacrete material from the cooling ambient medium.Although Gronke & Oh (2020b) only studied simulations in which the thermal instability made clouds over-pressurized (the contraction leads to overshooting pressure equilibrium), they argued that similar conditions arisefrom the shock that supersonic winds drive through clouds. Because they predict that clouds should only shatterwhen M w (cid:38) . χ = 100 and gas can’t cool below T cl ), we don’t expect the clouds in our simulations to shatter.Nevertheless, resolution of (cid:96) cool could be important for convergence of simulation properties because of its link togrowth. C.2. Measurement sensitivity to Resolution
Figure 22 illustrates how the evolution of the cloud’s bulk properties (survival fraction, purity fraction, and bulkvelocity) vary with resolution. The figure illustrates a remarkable level of convergence which seems to imply that thenet effects of mixing generally have only a weak dependence on resolution.Although the net effect of mixing doesn’t change significantly, the microscopic details can and do change withresolution. While we might not expect the average time derivative of K for all fluid elements to vary much with ixing & Cooling in Cloud-Wind Interactions K m i x i n g [ K c l t cc ] t cool,mix t cc = R cl / x = 8 R cl / x = 16 R cl / x = 32 R cl / x = 64 K m i x i n g [ K c l t cc ] T cl = 10 K t cool,mix t cc = 9.0 K m i x i n g [ K c l t cc ] T cl = 4 × 10 K t cool,mix t cc = 0.22 K / K cl )10 K m i x i n g [ K c l t cc ] T cl = 10 K t cool,mix t cc = 0.16 K / K cl ) 0 1 2log ( K / K cl ) 0 1 2log ( K / K cl ) 0.52.03.55.06.58.09.511.012.514.0 t / t cc Figure 24.
Comparison of passive scalar mass weighted and time averaged ˙ K ( K ) at different resolutions. As in Figure 23,different rows illustrate different simulations and resolution varies between columns. resolution, the derivative for individual fluid elements can vary wildly. For this work, we’ve measured the time-averaged ˙¯ K ( K, t ) just for fluid elements originating in the cloud. Therefore, we expect the ˙¯ K ( K, t ) measurements tobe fairly robust for K -bins in which the majority of the fluid elements originated in the cloud. However, when a largefraction of the fluid elements in a bin originated in the wind (and the fluid elements originating from the cloud are nolonger representative of all fluid elements in the bin), ˙¯ K ( K, t ) should be treated with care.Thus, we expect our ˙¯ K ( K, t ) measurements for low to intermediate K to be fairly robust for NR-X100 (and non-radiative simulations in general) and
SlowRC-T1e4 since the purity fraction remains high. However, for the fast coolingcases where the purity fraction drops, we expect more variation in ˙¯ K ( K, t ) at increasing K and over a larger range in K at later times. These are generally reflected in the convergence properties of dM ps /dK and ˙¯ K , which are illustratein Figures 23 and 24.We note that dM ps /dK ( K, t ) appears to have the most variability in bins with under 1% of the initial mass. Thus, wefocus our assessment of ˙¯ K ( K, t ) throughout this work on values computed from K bins that include at least 1% of thecloud’s initial mass. Because the measurements that don’t satisfy this condition can still be instructive (particularlyfor simulations with rapid cooling), we still show the other measurements in our figure, as translucent lines.We further note that calculation of the average ˙ K of all fluid elements in the system would improve substantiallyupon the reliability of our measurements. While doing this is possible, we consider our current measurements to beadequate for the purposes of conveying the premise of our mixing model.REFERENCES Armillotta, L., Fraternali, F., & Marinacci, F. 2016,MNRAS, 462, 4157, doi: 10.1093/mnras/stw1930 Banda-Barrag´an, W. E., Federrath, C., Crocker, R. M., &Bicknell, G. V. 2018, MNRAS, 473, 3454,doi: 10.1093/mnras/stx2541 Abruzzo et al.
Begelman, M. C., & Fabian, A. C. 1990, MNRAS, 244, 26PBordner, J., & Norman, M. L. 2012, in Proceedings of theExtreme Scaling Workshop, BW-XSEDE ’12(Champaign, IL, USA: University of Illinois atUrbana-Champaign), 4:1–4:11.http://dl.acm.org/citation.cfm?id=2462077.2462081Bordner, J., & Norman, M. L. 2018, arXiv e-prints,arXiv:1810.01319. https://arxiv.org/abs/1810.01319Bryan, G. L., Norman, M. L., O’Shea, B. W., et al. 2014,ApJS, 211, 19, doi: 10.1088/0067-0049/211/2/19Cooper, J. L., Bicknell, G. V., Sutherland, R. S., &Bland-Hawthorn, J. 2009, ApJ, 703, 330,doi: 10.1088/0004-637X/703/1/330Dav´e, R., Angl´es-Alc´azar, D., Narayanan, D., et al. 2019,MNRAS, 486, 2827, doi: 10.1093/mnras/stz937Dekel, A., & Silk, J. 1986, ApJ, 303, 39,doi: 10.1086/164050Dursi, L. J., & Pfrommer, C. 2008, ApJ, 677, 993,doi: 10.1086/529371Evans, C. R., & Hawley, J. F. 1988, ApJ, 332, 659,doi: 10.1086/166684Fielding, D. B., Ostriker, E. C., Bryan, G. L., & Jermyn,A. S. 2020, ApJL, 894, L24,doi: 10.3847/2041-8213/ab8d2cGronke, M., & Oh, S. P. 2018, MNRAS, 480, L111,doi: 10.1093/mnrasl/sly131—. 2020a, MNRAS, 492, 1970, doi: 10.1093/mnras/stz3332—. 2020b, MNRAS, 494, L27, doi: 10.1093/mnrasl/slaa033Haardt, F., & Madau, P. 2012, ApJ, 746, 125,doi: 10.1088/0004-637X/746/2/125Ji, S., Oh, S. P., & Masterson, P. 2019, MNRAS, 487, 737,doi: 10.1093/mnras/stz1248Kanjilal, V., Dutta, A., & Sharma, P. 2020, arXiv e-prints,arXiv:2009.00525. https://arxiv.org/abs/2009.00525Klein, R. I., McKee, C. F., & Colella, P. 1994, TheAstrophysical Journal, 420, 213, doi: 10.1086/173554Li, Z., Hopkins, P. F., Squire, J., & Hummels, C. 2020,MNRAS, 492, 1841, doi: 10.1093/mnras/stz3567Lochhaas, C., Thompson, T. A., & Schneider, E. E. 2020,arXiv e-prints, arXiv:2011.06004.https://arxiv.org/abs/2011.06004McCourt, M., Oh, S. P., O’Leary, R., & Madigan, A.-M.2018, MNRAS, 473, 5407, doi: 10.1093/mnras/stx2687McCourt, M., O’Leary, R. M., Madigan, A.-M., &Quataert, E. 2015, MNRAS, 449, 2,doi: 10.1093/mnras/stv355Melso, N., Bryan, G. L., & Li, M. 2019, ApJ, 872, 47,doi: 10.3847/1538-4357/aafaf5Miyoshi, T., & Kusano, K. 2005, Journal of ComputationalPhysics, 208, 315, doi: 10.1016/j.jcp.2005.02.017 Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS,473, 4077, doi: 10.1093/mnras/stx2656Rupke, D. 2018, Galaxies, 6, 138,doi: 10.3390/galaxies6040138Scannapieco, E., & Br¨uggen, M. 2015, ApJ, 805, 158,doi: 10.1088/0004-637X/805/2/158Schneider, E. E., Ostriker, E. C., Robertson, B. E., &Thompson, T. A. 2020, ApJ, 895, 43,doi: 10.3847/1538-4357/ab8ae8Schneider, E. E., & Robertson, B. E. 2017, TheAstrophysical Journal, 834, 144,doi: 10.3847/1538-4357/834/2/144—. 2018, ApJ, 860, 135, doi: 10.3847/1538-4357/aac329Schneider, E. E., Robertson, B. E., & Thompson, T. A.2018, ApJ, 862, 56, doi: 10.3847/1538-4357/aacce1Smith, B. D., Bryan, G. L., Glover, S. C. O., et al. 2017,MNRAS, 466, 2217, doi: 10.1093/mnras/stw3291Somerville, R. S., & Dav´e, R. 2015, ARA&A, 53, 51,doi: 10.1146/annurev-astro-082812-140951Sparre, M., Pfrommer, C., & Ehlert, K. 2020, arXive-prints, arXiv:2008.09118.https://arxiv.org/abs/2008.09118Sparre, M., Pfrommer, C., & Vogelsberger, M. 2019,MNRAS, 482, 5401, doi: 10.1093/mnras/sty3063Stone, J. M., & Gardiner, T. 2009, NewA, 14, 139,doi: 10.1016/j.newast.2008.06.003Tan, B., Oh, S. P., & Gronke, M. 2020, arXiv e-prints,arXiv:2008.12302. https://arxiv.org/abs/2008.12302Thompson, T. A., Fabian, A. C., Quataert, E., & Murray,N. 2015, MNRAS, 449, 147, doi: 10.1093/mnras/stv246Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, TheAstrophysical Journal Supplement Series, 192, 9,doi: 10.1088/0067-0049/192/1/9Veilleux, S., Cecil, G., & Bland-Hawthorn, J. 2005,ARA&A, 43, 769,doi: 10.1146/annurev.astro.43.072103.150610White, S. D. M., & Frenk, C. S. 1991, ApJ, 379, 52,doi: 10.1086/170483White, S. D. M., & Rees, M. J. 1978, MNRAS, 183, 341,doi: 10.1093/mnras/183.3.341Wiener, J., Zweibel, E. G., & Ruszkowski, M. 2019,MNRAS, 489, 205, doi: 10.1093/mnras/stz2007Wilson, L. A., & Fonner, J. M. 2014, in Proceedings of the2014 Annual Conference on Extreme Science andEngineering Discovery Environment, XSEDE ’14 (NewYork, NY, USA: ACM), 40:1–40:8,doi: 10.1145/2616498.2616534Zhang, D., Davis, S. W., Jiang, Y.-F., & Stone, J. M. 2018,ApJ, 854, 110, doi: 10.3847/1538-4357/aaa8e4 ixing & Cooling in Cloud-Wind Interactions29