aa r X i v : . [ a s t r o - ph . GA ] A ug Draft version August 28, 2018
Preprint typeset using L A TEX style emulateapj v. 11/10/09
MIXING IN SUPERSONIC TURBULENCE
Liubin Pan & Evan Scannapieco
School of Earth and Space Exploration, Arizona State University, P.O. Box 871404, Tempe, AZ, 85287-1404.
Draft version August 28, 2018
ABSTRACTIn many astrophysical environments, mixing of heavy elements occurs in the presence of a supersonicturbulent velocity field. Here we carry out the first systematic numerical study of such passive scalarmixing in isothermal supersonic turbulence. Our simulations show that the ratio of the scalar mixingtimescale, τ c , to the flow dynamical time, τ dyn (defined as the flow driving scale divided by therms velocity), increases with the Mach number, M , for M ∼ <
3, and becomes essentially constantfor M ∼ > . This trend suggests that compressible modes are less efficient in enhancing mixing thansolenoidal modes. However, since the majority of kinetic energy is contained in solenoidal modes at allMach numbers, the overall change in τ c /τ dyn is less than 20% over the range 1 ∼ < M ∼ <
6. At all Machnumbers, if pollutants are injected at around the flow driving scale, τ c is close to τ dyn . This suggeststhat scalar mixing is driven by a cascade process similar to that of the velocity field. The dependenceof τ c on the length scale at which pollutants are injected into flow is also consistent with this cascadepicture. Similar behavior is found for the variance decay timescales for scalars without continuingsources. Extension of the scalar cascade picture to the supersonic regime predicts a relation betweenthe scaling exponents of the velocity and the scalar structure functions, with the scalar structurefunction becoming flatter as the velocity scaling steepens with Mach number. Our measurements ofthe volume-weighted velocity and scalar structure functions confirm this relation for M ∼ < , but showdiscrepancies at M ∼ >
3, which arise probably because strong expansions and compressions tend tomake scalar structure functions steeper. INTRODUCTION
Understanding mixing in turbulent flows is essentialto interpreting a wide range of observations including:the metallicity dispersion in open clusters (e.g., Friel &Boesgaard 1992; Quillen 2002; DeSilva et al. 2006), thecluster to cluster metallicity scatter (Twarog et al. 1997),the abundance scatter along different lines of sight in theinterstellar medium (ISM) (Meyer et al. 1998, Cartledgeet al. 2006), the scatter in metallicities and abundanceratios in coeval field stars (e.g., Edvardsson et al. 1993;Nordstrom et al. 2004; Reddy et al. 2003, 2006). Thedegree of chemical inhomogeneity in these objects is con-trolled by the efficiency at which interstellar turbulencetransports and mixes metals from supernova explosionsand stellar winds (Scalo & Elmgreen 2004). This mix-ing process is also essential to the pollution of primordialgas by the first generation of stars in early galaxies (e.g.,Pan & Scalo 2007), and the transition from PopulationIII to Population II star formation (e.g., Scannapieco etal. 2003).Mixing in incompressible turbulence has been exten-sively studied in the fluid dynamics literature, where itis usually referred to as “passive scalar turbulence” (e.g.,Shraiman & Siggia 2000). The classic picture, knownas the Obukohov-Corrsin phenomenology, is a cascadeof scalar fluctuations, which is caused by the stretchingof the concentration field by the velocity field. This pro-duces structures of progressively smaller size, down to thediffusion scale at which molecular diffusion homogenizesand erases the fluctuations. In fact, molecular diffusionis the only process responsible for true mixing, but it op-erates very slowly at the large scales at which the scalarsources are injected. By providing a path from the scalarinjection scale to the diffusion scale, the cascade greatly accelerates the mixing process.Turbulent stretching is faster in smaller eddies, andthe mixing timescale is essentially given by the timescaleof eddies at the scalar injection scale, which is indepen-dent of the molecular diffusion coefficient. Scalar struc-tures at small scales are determined by the velocity struc-tures that control the scalar cascade. For an incompress-ible flow with a Kolmogorov spectrum, the Obukohov-Corrsin phenomenology predicts a − / § §
2, wereview the physics of turbulent mixing, and point outthe important issues that need to be addressed for a fullunderstanding of mixing in supersonic turbulence. In § § §
5, and our major results are summarized in theconclusion section ( § TURBULENT MIXING AND KINETIC ENERGYDISSIPATION
Physics of Mixing
We study mixing of pollutants in a turbulent flow v ( x , t ) with a density field ρ ( x , t ). The mixing processcorresponds to the homogenization of the concentrationfield C ( x , t ), defined as the ratio of the pollutant density to the local flow density. The evolution of the concen-tration field is described by the advection-diffusion equa-tion, ∂ t C + v i ∂ i C = 1 ρ ∂ i ( ρκ∂ i C ) + S ( x , t ) , (1)where κ denotes the molecular diffusivity, and the term S ( x , t ) represents the scalar sources. To understand thephysics of turbulent mixing, it is crucial to recognize therole of the turbulent velocity. Intuitively, a velocity fieldredistributes the concentration field and changes its spa-tial configuration, but does not affect the mass fractionat a given concentration level, since the tracer particlessimply follow the flow elements. This is true for any ve-locity field no matter how complex. Therefore, the veloc-ity field itself does not mix, because true mixing meansspatial smearing of pollutants, which can be caused onlyby molecular diffusion.The observation that the velocity field is not responsi-ble for true mixing can be formally confirmed from theequation for the variance of the concentration fluctua-tions, which is usually referred to as the scalar energyin the turbulence literature. To correctly reflect thisintuition in compressible flows, one needs to considerthe density-weighted variance, h ˜ ρC i , where the density-weighting factor ˜ ρ is the ratio of ρ to the average flowdensity, ¯ ρ , in the flow, and h· · ·i denotes an ensemble av-erage, which is equal to the average over the flow domainfor statistically homogeneous flows, as in our simulations.The equation for the density-weighted variance can bederived from eq. (1) with the help of the continuity equa-tion, yielding ∂ t h ˜ ρC i + ∂ i h ˜ ρC v i i = 2 h ∂ i (˜ ρκC∂ i C ) i− h ˜ ρκ ( ∂ i C ) i + 2 h ˜ ρSC i . (2)The second term on the left hand side of this equationis an advection term, which corresponds to the spatialtransport of concentration fluctuations between differentregions. In the statistically homogeneous case, this termvanishes, and without it, eq. (2) does not have an ex-plicit dependence on the velocity field. This proves theintuition that the velocity field does not truly mix. Infact, this is also true in an inhomogeneous flow. Whilethe advection term is not zero in this case, it is a surfaceterm, and thus its integral over the entire flow domain iszero. This means that the advection term conserves theglobal scalar variance .On the right hand side of eq. (2), the last term isthe source term, which corresponds to the variance in-crease due to the injection of new pollutants. The othertwo terms on the right hand side of this equation arisefrom the diffusion term in eq. (1). The first of theseis again a surface term, which vanishes in the homo-geneous case and conserves the global scale variance inthe inhomogeneous case. The term − h ˜ ρκ ( ∂ i C ) i rep-resents scalar dissipation. This term is negative defi- Without density-weighting, there would be a term −h C ∂ i v i i in the equation for the (volume-weighted) scalar variance. Thisterm is not zero in general, and represents the change in the volumefraction of the flow elements with a given concentration due tocompressions and expansions. This change in the volume fractionhas nothing to do with real mixing, and it is more appropriateto use the density-weighted variance, which is conserved by theadvecting velocity field. an & Scannapieco 3nite and thus monotonically reduces the scalar variance,corresponding to homogenization of the scalar fluctua-tions by molecular diffusion. We denote the dissipationrate as ¯ ǫ c ≡ h ˜ ρκ ( ∂ i C ) i and define a scalar dissipationtimescale τ c ≡ h ˜ ρC i / ¯ ǫ c to characterize it. The lineardependence on κ in the definition of ¯ ǫ c may leave an im-pression that the scalar dissipation rate strongly dependson κ . However, that is not true because the scalar gra-dient in the definition also depends on κ . With decreas-ing κ , larger scalar gradients can develop, which maycompensate the decrease in κ . As seen from the phys-ical picture below and in § .The molecular diffusivity κ is tiny in most practical en-vironments, thus the dissipation rate is significant onlyif the scalar gradient is large, meaning that true mixingoccurs only at very small scales. As mentioned earlier,at the injection scale of scalar sources, the moleculardiffusion is usually very slow, and a significant mixingrate needs a velocity field to produce small-scale scalarstructures, or equivalently, large concentration gradients.In an incompressible turbulent flow, the velocity fieldstretches, contracts and folds the scalar field, leading toscalar structures at smaller and smaller scales. The sizeof the scalar structures keeps decreasing toward the dif-fusion scale where the molecular diffusion efficiently ho-mogenizes the structures. Since κ is small, a wide scaleseparation usually exists between the scalar source sizeand the diffusion scale.The diffusion scale is essentially the scale at which thediffusion term in eq. (1) becomes larger than the ad-vection term, and thus can derived by comparing thesetwo terms accounting for their scale dependences. Theestimate of the diffusion scale depends on the Schmidtnumber, Sc , defined as the ratio of viscosity to diffusiv-ity, which determines the diffusion scale relative to theviscous dissipation scale of the velocity field (the scaleat which kinetic energy is removed by viscosity). Thederivation of the diffusion scale at different Sc valuescan be found in Monin and Yaglom (1975) and Lesieur(1997). For heavy elements in a neutral interstellarmedium, the diffusivity is expected to be smaller thanthe viscosity due to the larger cross section of heavy el-ements. However, for relatively light atoms like oxygen,the diffusivity is probably of the same order as the vis-cosity, and the Schmidt number is close to unity.A second important dimensionless number for passivescalars is the Peclet number, P e . It is defined as
U L/κ where U and L are characteristic velocity and lengthscale of the turbulent flow. The Peclet number is relatedto the Reynolds number, Re , by P e = ScRe . Therefore,for diffusing species with Sc ≈
1, the Peclet number isclose to the Reynolds number, ∼ > , in typical inter-stellar clouds.In summary, a turbulent velocity field acts as a catalystthat enhances the mixing efficiency by feeding molecu-lar diffusion with large gradient structures. The mixingtimescale, τ c , is thus determined by the rate at which theturbulent velocity field brings the scalar structures downto the diffusion scale. It is essentially independent of κ The mixing timescale may have a weak dependence on κ if theSchmidt number (see definition below) is much larger than 1. for the case with Sc ∼ < Scalar Cascade
The classic picture for the generation of scalar struc-tures at small scales in incompressible turbulence is acascade similar to the energy cascade of the velocity field(e.g., Lesieur 1997). A constant flux of the scalar en-ergy occurs along the cascade over the convective range,defined as the range of scales between the characteris-tic length scale of the scalar sources and the diffusionscale. In this range, the scalar fluctuations are regulatedprimarily by the advecting velocity. The flux feeds thedissipation at the diffusion scale and is thus expected tobe equal to the average scalar dissipation rate. Assum-ing the timescale for a cascade step is of the order of theeddy turnover time at a given scale leads to the follow-ing scaling for the concentration difference, δC ( l ), over aseparation l , δC ( l ) ≃ ¯ ǫ c lδv ( l ) , (3)where δv ( l ) is the velocity difference over l , and ¯ ǫ c is thescalar dissipation rate.In incompressible turbulent flows, we have the Kol-mogorov scaling for the velocity difference, δv ( l ), in theinertial range, i.e., δv ( l ) ≃ ¯ ǫ v1 / l / , where ¯ ǫ v is the av-erage dissipation rate of kinetic energy. With this ve-locity scaling, eq. (3) gives δC ( l ) ≃ ¯ ǫ − / ¯ ǫ c l / in theinertial-convective range (i.e., the intersection of the iner-tial range for the velocity field and the convective rangefor the scalar field). This scaling for the scalar differ-ence, usually referred to as the Obukhov-Corrsin scalinglaw, corresponds to a k − / scalar spectrum. This spec-trum has been confirmed for passive scalars in isotropicincompressible turbulence by both experiments and nu-merical simulations (Monin and Yaglom 1975; Sreeni-vasan 1996; Mydlarskiand and Warhaft 1998; Yeung etal. 2002; Watanabe & Gotoh 2004).On the other hand, the scalar structure function andspectrum have not been studied in supersonic turbulence.Therefore it is unknown whether eq. (3) originally pro-posed for passive scalar mixing in incompressible turbu-lence is valid in supersonic turbulence. The scaling ofthe velocity difference in supersonic turbulence has beenshown to be significantly steeper than the Kolmogorovscaling (e.g., Kritsuk et al. 2007). Thus, if eq. (3) holdsfor turbulent flows at any Mach number, one would ex-pect the scalar structure functions to become flatter withincreasing Mach number. We will measure the scalingexponents, ζ v and ζ c , of the 2nd order velocity structurefunctions ( h δv ( l ) i ∝ l ζ v ) and the scalar structure func-tions ( h δC ( l ) i ∝ l ζ c ) from our simulations at differentMach numbers. Eq. (3) predicts the following relationbetween the two exponents, ζ c ≃ − ζ v / , (4)which will be tested against our simulation data.Eq. (3) also suggests that, if the scalar field andthe flow are driven at similar scales, the dissipationtimescales for scalar fluctuations and for kinetic en-ergy will be similar. They are both determined by theturnover time of largest eddies at the driving/sourcescale, which is of the order of the dynamical timescale, Mixing in Supersonic Turbulencebecause the cascade becomes progressively faster atsmaller scales, and the most time is spent at largescales. Numerical simulations confirmed that these twotimescales are similar in incompressible turbulence (see,e.g. Donzis, Sreenivasan, & Yeung et al. 2005 and refer-ences therein) and that the scalar dissipation timescaleis about half the energy dissipation timescale. Further-more, in the supersonic case, the energy dissipation ratehas been found to be similar to that in incompressibleturbulence, i.e., close to the dynamical timescale (Stoneet al. 1998; Mac Low 1998; Padoan and Norlund 1999),implying an energy cascade similar to that in incompress-ible turbulence. However, in the supersonic case, thedissipation timescale for forced scalars has never beencomputed. Kinetic Energy Dissipation
Our simulations also give us an opportunity to explorethe dissipation of kinetic energy in supersonic turbulence.The momentum equation in a driven turbulent flow isgiven by, ∂ t v i + v j ∂ j v i = − ∂ i pρ + 1 ρ ∂ j ( σ ij ) + f i ( x , t ) , (5)where p ( x , t ) is the pressure, f ( x , t ) is the driving force,and σ ij is the viscous stress tensor. For an ideal gas,the bulk viscosity is negligible and the viscous tensor σ ij = ρν ( ∂ i v j + ∂ j v i − ∂ k v k δ ij ) where ν is the kinematicviscosity.From the momentum equation, eq. (5), and the conti-nuity equation, we can derive an equation for the averagekinetic energy per unit mass, h ˜ ρv i , ∂ t h
12 ˜ ρv i = 1¯ ρ h p∂ i v i i + h ˜ ρf i v i i− * ˜ ρν (cid:18) ∂ i v j + ∂ j v i − ∂ k v k δ ij (cid:19) + (6)The three terms on the right hand side of eq. (6) repre-sents the pdV work,the energy injection from the driv-ing force, and the viscous dissipation of kinetic energy,respectively.The pdV term, which vanishes in incompressible flows,is a two-way energy exchange between kinetic energy andthermal energy. Although the conversion by pdV workbetween the two energy forms is reversible, the exchangerates in the two directions can be asymmetric. In fact,the term preferentially converts kinetic energy to heatdue to an anticorrelation between the density (or pres-sure) and the velocity divergence. Consider, for exam-ple, an expanding region and a converging region withthe same |∇ · v | . In the converging region, the density istypically larger, and thus the rate at which the pdV workconverts kinetic energy to heat is faster than the conver-sion from thermal to kinetic energy in the correspondingexpanding region. Therefore, in supersonic turbulence,the pdV work provides another mechanism for kineticenergy loss in addition to viscous dissipation. We de-note the rate of kinetic energy loss by this mechanismas ¯ ǫ p ≡ −h p ¯ ρ ∂ i v i i , and point out that this effect has notbeen explicitly considered in previous studies.Viscous dissipation, corresponding to the second termon the right hand side of eq. (6), is an irreversible con- version of kinetic energy to heat. If ρν is constant, whichis true if the flow temperature is roughly constant, theoverall average kinetic energy dissipation rate in the flowcan be written as ¯ ǫ ν = ˜ ρν [ h ( ∇ × v ) i + h ( ∇ · v ) i ] us-ing integration by parts. We define a viscous dissipationtimescale, τ diss , as h ˜ ρv i / ¯ ǫ ν , and point out that thistimescale for kinetic energy loss purely by viscous dissi-pation, to our knowledge, has also never been evaluatedin supersonic turbulence.In fact, the only timescale for kinetic energy lossestimated in earlier studies is the overall dissipationtimescale τ v that characterizes the total rate for the ki-netic energy loss by both pdV work and viscous dissipa-tion, τ v = h ˜ ρv i / (¯ ǫ ν + ¯ ǫ p ) (Stone et al. 1998; Mac Low1998; Padoan & Norlund 1999; Lemaster & Stone 2009).For example, the timescale for kinetic energy loss, t diss defined in Lemaster & Stone (2009) corresponds to thetimescale τ v defined here. Below we will evaluate both τ diss and τ v , separating out the relative importance ofkinetic energy loss by pdV work and viscous dissipationas a function of Mach number. METHODS
To study mixing in supersonic turbulence, we sim-ulated hydrodynamic turbulent flows at various Machnumbers using the FLASH code (version 3.2), a mul-tidimensional hydrodynamic code (Fryxell et al. 2000)that solves the Riemann problem on a Cartesian gridusing a directionally-split Piecewise-Parabolic Method(PPM) solver (Colella & Woodward 1984; Colella &Glaz 1985; Fryxell, M¨uller, & Arnett 1989). The hydro-dynamic equations and the advection-diffusion equationwere solved on a fixed 512 grid for a domain of unitsize with periodic boundary conditions. We adopted anisothermal equation of state in all our simulated flows,which was obtained by setting the ratio of specific heatsto be 1.0001. We also performed simulations at the res-olution of 256 for the purposes of checking convergencyand examining the effect of different driving schemes.Unless explicitly stated otherwise, the following descrip-tion of the methods and results is based on the 512 simulations. Velocity Forcing
To drive turbulence in our simulation, we made use ofthe FLASH3 “Stir” package (Benzi et al. 2008), which weheavily modified for our purposes. The flow velocity wasdriven by the forcing term, f ( x , t ) in eq. (5), where theamplitude was adjusted to achieve different Mach num-bers. In our 512 simulations, we drove the turbulencewith solenoidal modes only, i.e., ∇ · f = 0 and took f tobe a Gaussian random vector with an exponential tempo-ral correlation. A finite correlation timescale is adoptedin the forcing scheme (e.g., Padoan and Nordlund 1999).This is different from studies using an infinitesimal cor-relation timescale with independent forcing at each timestep (e.g., Lemaster & Stone 2009) or an infinite corre-lation timescale with a fixed driving force (e.g., Kritsuket al. 2007). The driving scheme can be summarizedas h f i ( k , t ) f j ( k , t ′ ) i = P f ( k ) (cid:16) δ ij − k i k j k (cid:17) exp (cid:16) − t − t ′ t f (cid:17) , where t f is the forcing correlation timescale. The forcingwave numbers are chosen to be in the range 1 ≤ k/ π ≤ simulations with a different driving scheme where we set1/3 of the forcing energy to be in compressible modes.The driving force for these flows was otherwise the sameas for the 512 simulations.To characterize the large scale properties of the flow,we define a forcing length scale, L f , from the forcingspectrum, L f = R πk P f ( k ) d k / P f ( k ) d k . For the forc-ing spectrum chosen here, we find L f = 0 .
46 in units inwhich the domain size is 1. Lemaster & Stone (2009)adopted a forcing spectrum, P f ( k ) ∝ k exp( − k/k pk ),that strongly peaks at a wave number k pk . Insertingthis spectrum into our definition for L f , we obtain that L f is equal to 2 π/k pk , which is exactly the same as thecharacteristic flow length scale used in that study. Foreach simulation in our study, we define a dynamical time, τ dyn ≡ L f /v rms , where v rms is the 3D density-weightedrms velocity.We simulated turbulent flows at six different Machnumbers ranging from transonic to highly supersonic. Ateach Mach number, the simulation covered more than 10dynamical times. After an initial phase of steady in-crease, the rms velocity saturated in a few dynamicaltimes, and a statistically steady state was reached. Inthe steady state, the rms velocity was not exactly con-stant, but had small variations with time, with an (rms)amplitude of about 2% to 5%, depending on the Machnumber.Following Lemaster & Stone (2009), we used thedensity-weighted (3D) rms velocity in our definition forthe Mach number, i.e., M = h ˜ ρv i / /c s where c s is thesound speed. For each flow, we averaged the density-weighted rms velocity over about 10 dynamical times,and the average Mach numbers in the six flows were 0.9,1.4, 2.1, 3.0, 4.6 and 6.1. If the volume-weighted rmsvelocity was used, the Mach number was only slightlylarger (by about 7% except for the M = 0 . L f was thesame at all Mach numbers, and the dynamical timescale τ dyn = L f /v rms was proportional to 1 /M . Forced and Decaying Scalars
We studied both forced scalars and decaying scalars.For the forced case, we continuously added pollutantsin the flow by the source term S ( x , t ) in the advection-diffusion equation (eq. 1). Similar to the velocity driving, S was taken to be a Gaussian variable . with a spectrumgiven by h S ( k , t ) S ( k , t ′ ) i = P s ( k ) exp (cid:16) − t − t ′ t s (cid:17) . We con-sidered several different forcing schemes for P s ( k ). In thefirst scheme, P s ( k ) was chosen to have the same shapeas the velocity forcing spectrum, i.e., P s ( k ) ∝ P f ( k ).With exactly the same forcing pattern as the velocity, We point out that scalar sources are probably non-Gaussian inpractical situations. However, it is a common practice to set thescalar source to be Gaussian in studies of passive scalar physics,which is useful in revealing non-Gaussian scalar features that arisepurely from mixing. this scheme allows for a direct comparison of the dissipa-tion timescale of the scalar variance to that of the kineticenergy. As discussed below, this reveals important dif-ferences in the physics of scalar and kinetic energy dis-sipation, especially concerning the role of compressiblemodes and shocks.Like the rms velocity, the scalar variance exhibitedtemporal variations in the steady state. In order to re-duce their amplitude, we used three independent scalarswith the same source spectrum with identical scalarand flow forcing. Averaging over the three indepen-dent scalars resulted in a much smaller amplitude (inthe range of 2%-5%) in the temporal variations of therms concentration than that of each individual scalar.To test the dependence of the mixing timescales onthe characteristic length scale of the scalar sources, weconsidered three other schemes in which the scalars wereforced at larger wave numbers. In these cases, we chosethe source wave numbers to be 3 . ≤ k/ π ≤ . , . ≤ k/ π ≤ . , and 16 . ≤ k/ π ≤ . , respec-tively. Again, each independent mode in these rangeswas given the same forcing power. Similar to the lengthscale, L f , of the flow driving force, we defined a sourcescale, L s = R πk P s ( k ) d k / P s ( k ) d k . The source lengthscales are 0 .
25, 0 .
124 and 0 .
061 the size of the computa-tion box for the three cases, respectively.In our simulations, we take the scalar source termto have zero mean and do not consider the evolutionof the mean concentration, which is simply determinedby the rate at which pollutants are injected. By defini-tion, the concentration is positive and bounded in therange 0 ≤ C ≤
1. The negative concentration valuesin our simulations should be understood as relative tothe mean concentration. The amplitude of the scalarfluctuations in our simulations is arbitrary (due to thearbitrarily chosen amplitude for the source term). Thisdoes not affect the study of mixing physics, because theadvection-diffusion equation is linear with the concentra-tion.In all the simulations, we took both t f and t s to be1/4 of the sound crossing timescale, defined as the boxsize divided by the sound speed. We also tried differentvalues of t f and t s (but with a resolution of 256 ), andfound that the results were not affected by the specificchoice of the forcing/source correlation timescales.For the decaying case, we set an initial configurationfor the scalar field when the turbulent flows were welldeveloped, and let it evolve with no continuing sources,i.e., subsequently setting S ( x , t ) to zero. For each Machnumber, we studied four decaying scalars with differentinitial spectra. The initial scalar spectra were taken tobe the same as the source spectra chosen for forced cases,i.e., the wave number ranges were set to be 1 ≤ k/ π ≤ . ≤ k/ π ≤ . , . ≤ k/ π ≤ . , and 16 . ≤ k/ π ≤ .
5, respectively.Although the main motivation for considering scalarsforced at different scales is theoretical for a full under-standing of mixing, they are also of practical interest.For example, if the main energy source for the ISMturbulence is supernova explosions, then, for mixing ofnew metals from SNe, the source length scale would besame as the flow driving scale. On the other hand, thescalar injection scale would be different from the flowforcing scale if the pollutant sources and the flow en- Mixing in Supersonic Turbulenceergy sources are different. An example is self-enrichmentin star-forming clouds, where the scalar source lengthscale may be smaller than the flow driving scale if theturbulent energy is mainly from the cascade of the in-terstellar turbulence at larger scales. The decaying casemay be applicable to mixing in galaxies with an episodicstar-formation history: the mixing process for elementsproduced by core-collapse SNe in between star bursts inthese galaxies would be a decaying case.Finally, we point out that neither the viscous term, ρ ∂ j σ ij , nor the diffusion term, ρ ∂ i ( ρκ∂ i C ), was explic-itly included in our simulations. Rather, kinetic energydissipation and scalar dissipation were both controlledby numerical diffusion. Therefore, both the viscous dis-sipation scale and the diffusion scale are expected to bearound the resolution scale. The effective viscosity anddiffusivity are thus similar, with the effective Schimdtnumber, Sc ≈
1, and the Peclet number about equalto the Reynolds number. For the case with identicalscalar and flow forcing, the convective scale range coin-cides with the inertial range, while the scalars forced atlarger wave numbers has a smaller convective range thanthe inertial range.
Calculation of the Timescales
To calculate the scalar dissipation timescale, τ c , weneed the scalar dissipation rate, ǫ c , which could not bedirectly computed in our simulations because the dissi-pation was controlled by numerical diffusion. However,in the case of forced scalars, after the scalar fluctuationsreach a statistically stationary state, we have a balancebetween the dissipation term and the source term in eq.(2), i.e., ǫ c = ǫ s ≡ h ˜ ρSC i . The computation of the sourceterm is straightforward, and we calculated the scalar dis-sipation timescale as τ c = h ˜ ρC i / h ˜ ρSC i . For the schemewith identical scalar and flow forcing, we used the averageof the source term over the three independent scalars.Similarly, the viscous dissipation of kinetic energy wasobtained using energy balance, ¯ ǫ f − ¯ ǫ p − ¯ ǫ ν = 0, in thestatistically stationary state. The energy injection by thedriving force per unit time, ¯ ǫ f , and the contribution frompdV work, ¯ ǫ p were calculated by h ˜ ρf i v i i and ¯ ρ − h p∂ i v i i in eq. (6), respectively. In the stationary state, this al-lowed us to evaluate the viscous dissipation timescale as τ diss = h ˜ ρv i / (¯ ǫ f − ¯ ǫ p ), and the timescale for the over-all kinetic energy loss including the contribution from thepdV work as τ v = h ˜ ρv i / ¯ ǫ f .The kinetic energy loss by both the pdV work and theviscous dissipation was stored as thermal energy. Wepoint out that the isothermal equation of state adoptedin our simulations was used to imitate a roughly constanttemperature due to the effect of efficient radiative coolingin interstellar clouds. In nature, the heat converted fromkinetic energy would be radiated away.We find that, if the forcing correlation timescale, t f ,was much smaller than the typical timescale for the ki-netic energy loss to thermal energy, τ v , the energy injec-tion rate goes like ¯ ǫ f ≃ t f R P f ( k ) d k . On the other hand,if t f is larger than τ v , as is the case for our choice of τ f ,¯ ǫ f ∝ τ v R P f ( k ) d k because only the forcing energy withina timescale τ v remains as kinetic energy in the flow. Thesame applies to the variance input from scalar sources,i.e., ¯ ǫ s ≃ t s R P s ( k ) d k if t s ≪ τ c or ≃ τ c R P s ( k ) d k for t s ≫ τ c .Similar to the velocity and scalar variances, the otherquantities needed for the calculation of timescales alsoshow temporal variations to different degrees. Thesequantities include the kinetic energy injection by thedriving force, the pdV work, and the scalar variance in-put by the sources. In our calculations, we used the tem-poral average of all these quantities over about 10 dynam-ical times. Comparing with results from 256 simulationsfor the corresponding flows, we find that the timescaleshave converged at our resolution of 512 zones (see alsoLemaster & Stone 2009). RESULTS
The Turbulent Flows
In the left panels of Fig. 1, we show the X-component ofthe velocity field on a slice (X-Y plane) of the simulationgrid at one snapshot for two Mach numbers, M = 0 . M = 6 .
1. The snapshots are taken at t = 8 τ dyn and5 . τ dyn for M = 0 . M = 6 .
1, respectively. The con-centration field shown here is from a scalar forced at theflow driving scale. The velocity field exhibits extremelydifferent features in these two limits. In the M = 0 . M = 6 . M = 6 .
1. For example, strong shears andvortices are usually found behind narrow postshock re-gions. As we will show later, the solenoidal modes playimportant roles in both energy dissipation and mixing atall Mach numbers.In Table 1, we list the overall properties of the flowas a function of Mach number. In the 2nd column, wegive the ratio of the density-weighted divergence varianceand vorticity variance, h ˜ ρ ( ∇ · v ) i / h ˜ ρ ( ∇ × v ) i . This ra-tio characterizes the contribution of compressible modesto viscous dissipation relative to solenoidal modes, andit is smaller than 1 even for M = 6 .
1, suggesting thatthe majority of kinetic energy is dissipated by solenoidalmodes even in highly supersonic flows (Kritsuk et al.2007). However, we point out that this ratio is onlya rough estimate for the relative contribution from thetwo modes, as it is not clear how the kinetic energy isdissipated by numerical diffusion by the high-order PPMsolver. In the 3rd column, we also give the same ratios,but with no density weighting. For both cases, the ra-tio increases with the Mach number at small M . Thedensity-weighted ratio saturates at M ≈
3, while thevolume weighted one saturates at a little larger M, andthe density-weighted ratios are all smaller because of ananti-correlation between the density and divergence.The 4th column lists the fraction of kinetic energy incompressible modes. We computed this fraction by de-composing the energy spectrum into a solenoidal partand a potential part and comparing the kinetic en-ergy contained in each. Although the flow is driven bysolenoidal modes only, a finite fraction of kinetic energyappears in the compressible modes, which are convertedfrom the solenoidal modes. The fraction increases froman & Scannapieco 7 Fig. 1.—
The flow velocity in X direction (left) and the scalar field (right) on a slice of the simulation grid, with a linear color scale goingfrom the minimum value (blue) to the maximum value (red) on the slice. The top two panels are for the M = 0 . t = 8 . τ dyn , andthe bottom two panels correspond to the M = 6 . . τ dyn .
5% at M = 0 . M = 6 . k/ π ≈
4. Thus we cal-culated the kinetic energy contained in the compressiblemodes and the solenoidal modes in the inertial range byintegrating the potential and solenoidal spectra over allwave numbers k/ π ≥ , which also includes the negli-gible contribution from the dissipation range. The 5thcolumn in Table 1 gives the fraction of kinetic energyin compressible modes at k/ π ≥
4, which increases atsmall M and saturates at about 1/3 for M ∼ >
3. Thefraction of 1/3 implies an energy equipartition betweenthe solenoidal modes and the potential modes in the in- ertial range, which is reached at M ≈
3. We find theslope of the compressible spectrum is close to that of thesolenoidal spectrum in the inertial range, meaning thata equipartition is established at each wave number in theinertial range for M ≈
3. We find this equipartition isuniversal regardless the compressible fraction in the flowdriving (see below).Note that the fractions given in columns 4 and 5 do notaccount for density-weighting, which cannot be perfectlyapplied to the fraction of kinetic energy in compressiblemodes at large M . Although the velocity field can be de-composed into a solenoidal part, v s , and a potential part, v p , the density-weighed ratio h ˜ ρv i / ( h ˜ ρv i + h ˜ ρv i ) can-not be interpreted as the energy fraction at high Machnumbers. This is because the denominator is not equalto the total energy, which is given by h ˜ ρv i + 2 h ˜ ρ v p · v s i + h ˜ ρv i . In fact, we find that the density-weighted corre-lation term h ˜ ρ v p · v s i is not zero. Rather, it is negative, Mixing in Supersonic Turbulence TABLE 1Flow properties at different Mach numbers M h ˜ ρ ( ∇ · v ) i / h ˜ ρ ( ∇ × v ) i h ( ∇ · v ) i / h ( ∇ × v ) i f comp f comp ( k/ π ≥ f pdV decreases with increasing M, and becomes significant at M ∼ > −
3. Therefore, density-weighting for the fractionof energy in the compressible modes cannot be well de-fined at large Mach numbers. However, there is no suchproblem for the volumed-weighted fraction because it canbe shown that h v p · v s i is exactly zero in the homogeneouscase, and we thus only show the volume-weighted frac-tion in Table 1.The last column in Table 1 lists the ratio of the kineticenergy loss by pdV work to the total kinetic energy loss.This will be discussed in detail in § ) simulations covering a similar range of Machnumbers but with 1/3 forcing energy contained in thecompressible modes. We found that, with this differ-ent driving scheme, most quantities listed in Table 1 re-main essentially unchanged, except the overall fractionof energy in the compressible modes (column 4). Thisis expected since the overall compressible fraction de-pends on the driving scales, which contain most of thekinetic energy, while the rest of the quantities depend onthe statistics at smaller scales, which should be driving-independent.The change in the overall compressible fraction is slightfor M ∼ <
2, only a few percent larger than the case withsolenoidal driving, but there is a significant increase inthe fraction at M ∼ >
3. With the compressible driv-ing scheme, the fraction saturates at ≈
24% for M ∼ > ≈
18% given in Table 1. Wenote that the overall compressible fraction in the flowturns out to be smaller than that (1/3) in the drivingforce at all Mach numbers, indicating a preferential con-version from compressible modes to solenoidal modes atthe driving scales. This is different from the case withsolenoidal driving where a fraction of kinetic energy insolenoidal modes is converted to compressible modes.Of particular interest is that, in the inertial range, theequipartition between solenoidal and compressible modesis also found in the M ∼ > The Concentration Field
In the right panels of Fig. 1, we show the concentra-tion on slices taken from simulations with Mach num-bers M = 0 . M = 6 .
1. As was the case for thevelocity field, there are large differences between the tworuns. At M = 0 . , the scalar field shows numerous small-scale structures, which are clearly produced by stretch- ing and shear in vortices. There are many structureswith sharp concentration contrasts, which are usually re-ferred to as “cliffs and ramps” in the literature for pas-sive scalar in incompressible turbulence (e.g., Shraiman& Siggia 2000). These cliffs and ramps are produced byturbulent stretching, which rearranges the scalar gradi-ents that initially exist only at large scales and bringsfluid elements with very different concentration levelsnext to each other. The “cliffs and ramps” have largescalar gradients, and are thus strongly dissipative.On the other hand, the scalar field at M = 6 . M = 0 . §
5, where we also compare the efficiency of the compress-ible modes relative to the solenoidal modes in amplifyingthe scalar gradients.Note that there is a fundamental difference betweenscalar edges and velocity shocks. Although the scalaredges in Fig. 1 appear to be discontinuous, they areactually continuous. This can be shown by consider-ing the concentration change across a shock front in thelimit of infinitesimal shock thickness. From the forma-tion mechanism of the scalar edges, the concentrationchange over a shock front can be roughly estimated bythe product of three factors: the preshock gradient, thedensity jump factor (i.e., the squeezing factor) and theshock thickness. Therefore, when the shock thicknessapproaches zero (corresponding to decreasing viscosity),the concentration change across the shock front wouldalso approach zero, since the density jump factor wouldremain constant and finite as the shock width decreasesto infinitesimal. This means that the concentration isessentially continuous, unlike the velocity shocks, whichare exact discontinuities in the limit of infinitesimal vis-cosity. This argument is confirmed by simulations for 1Dsupersonic flows at very high resolutions.The continuity of the scalar field across a shock of in-finitesimal width is also expected from the conservationof the tracer mass. Mass conservation implies that thetracer density has the same jump as the flow densityacross a shock, and thus the concentration, which is theratio of the two densities, would be continuous across thefront.The reason the edges appear to be discontinuous inour simulation is due to the limited resolution. A finiteand significant shock width in our simulations means thatthe materials across the shock front could come from wellseparated regions (especially for strong shocks), makinga considerable scalar change across the shock front pos-sible. With increasing resolution, we would expect thescalar structures around the shocks appear to be morecontinuous, while the shocks would clearly remain dis-continuous. If a shock had an infinitesimal width, itseffect on the scalar field is just to amplify the scalar gra-dient in the postshock region.The difference between the scalar edges and the veloc-ity shocks has an interesting implication. Every shockis intrinsically a strong dissipative structure for kineticenergy, but it is not true that every shock produces astrong dissipative scalar structure, considering that ascalar edge is essentially not a discontinuity. In fact,most shocks cannot produce scalar structures that dissi-pate scalar fluctuations at the same level as shocks forenergy dissipation. To produce such a structure, a shockneeds to reduce the length scale of a scalar structure closeto the diffusion scale. Whether a shock can give rise tothis scale reduction depends on the shock strength, whichdetermines the squeezing factor across the shock. We willshow that strong shocks capable of such a scale reduc- t c , v , d i ss / t dyn M t c t v t diss t c / t v M Fig. 2.—
The scalar dissipation timescale, τ c , the viscous dis-sipation timescale, τ diss , and the timescale for the overall kineticenergy loss, τ v , as functions of the Mach number. The timescalesare normalized to the dynamical time τ dyn , such that the behaviorof the timescales with M depends mainly on the effect of compress-ible modes relative to the solenoidal modes. The inset shows theratio of τ c to τ v , and the filled circle is from the simulation resultby Donzis et al. (2005) for incompressible turbulence. tion are very rare in §
5, and this point is responsiblefor the difference in the roles of shocks in kinetic energydissipation and in scalar dissipation.
Timescales
Kinetic Energy Dissipation
In Fig. 2, we plot the timescales defined in § § τ dyn . The behavior of the normalizedtimescales as a function of M is expected to essentiallyreflect the contribution from compressible modes relativeto solenoidal modes.The dotted line in Fig. 2 shows the normalizedtimescale for the kinetic energy loss by viscous dissipa-tion, τ diss /τ dyn . This timescale steadily decreases as M increases from 0.9 to 3, and shocks become stronger andmore frequent. Shocks are strong dissipative structuresfor kinetic energy, and thus compressible modes providean extra fast channel for the viscous dissipation of kineticenergy, which does not exist in incompressible flows. Thesteady decrease in the τ diss /τ dyn curve with M for M ∼ < M ∼ >
3, corresponding to the slower decreaseof τ diss /τ dyn in that range of M . This suggests that thefaster viscous dissipation at larger M may be primar-ily due to a larger contribution from the compressiblemodes. Overall the density-weighted divergence to vor-ticity variance ratio increases by about 50% as the Machnumber increases from M = 0 . τ v /τ dyn , the normalizedtimescale for the overall kinetic energy loss by both pdVwork and viscous dissipation. This timescale is signifi-0 Mixing in Supersonic Turbulencecantly smaller than τ diss because the pdV work gives aconsiderable contribution to the kinetic energy loss in thechosen range of M .The fraction of kinetic energy that is converted to heatby pdV work is given in the 6th column of Table 1, whichshows that the fraction peaks at the intermediate Machnumber, M ≃
2. At low Mach numbers, increasing M allows more compressible modes to appear and larger |∇ · v | to be available, leading to a rapid increase in therate of pdV work. On the other hand, the importance ofpdV work in converting kinetic energy into heat decreases with M for M ∼ >
2. At these large Mach numbers, theamplitude of the divergence scales linearly with M , cor-responding to the saturation of the fraction of energy incompressible modes (see Table 1) and the increase in theamplitude of pressure fluctuations is slower than ∝ M .Therefore the conversion rate by pdV work does not in-crease faster than ∝ M . Since the viscous dissipationrate increases as ∝ v , this means that the fraction ofenergy loss by pdV work decreases at large M .This trend is also expected from another perspective.The viscous dissipation rate by compressible modes is asecond-order function of the velocity divergence and in-creases with M faster than the pdV work, which scaleslinearly with the divergence. Thus, as M increases above2, the fraction of kinetic energy loss by pdV work startsto decrease as the viscous dissipation rate by compress-ible modes becomes a significant fraction of the total vis-cous dissipation rate.The increase in the fraction of kinetic energy loss bypdV work for M ∼ < τ v /τ dyn curve a lit-tle steeper than the curve for, τ diss /τ dyn , the normalizedtimescale for viscous dissipation only. On the other hand,the τ v /τ dyn curve is flatter than the τ diss /τ dyn curve for M ∼ > M >
3) is faster than in thetransonic case with M = 0 . ≈ . τ v from our simulations is in good agreementwith the results given in their Table 1 (the column for t diss /t f ). Scalar Dissipation
The solid curve in Fig. 2 shows the timescale for thescalar dissipation, t c , focusing on the case in which thesource pattern is identical to the flow driving. The be-havior of the scalar dissipation timescale is completelydifferent from that of the kinetic energy, although all thescalar fields are driven in exactly the same manner asthe velocity fields. Unlike the kinetic energy dissipationtimescale, the normalized timescale for scalar dissipation, τ c /τ dyn , increases by about 20% as M goes from 0 . M (see Table 1).This indicates that the existence of more compress-ible modes reduces the mixing efficiency. There are twoeffects that may contribute to the increase of the mix-ing timescale. First, as shown in more detail below( § M (see § § M , which tends to make the structure functionsteeper ( § M = 3, τ c /τ dyn is almost constant. Again,this is consistent with the fraction of kinetic energy insolenoidal modes, which is constant at 66% in the in-ertial range due to energy equipartition. We note thatthere a slight decrease in the mixing timescale curve inFig. 2 as M goes from 3 to 6.1. This is related to the factthat, although the energy fraction in compressible modesin the inertial range is constant in this range of Machnumbers, the probability of strong compression eventsincreases with M , as can be seen from the increase inthe amplitude of density fluctuations. As discussed in §
5, the compressible modes do contribute to the ampli-fication of the scalar gradient, although the contributionis tiny in comparison to the solenoidal modes. At thecurrent resolution (512 ), the contribution to the gradi-ent amplification by compressible modes is a few percentat M = 6 . § M = 3. This is responsible for the slight decrease in themixing timescale in this range of Mach numbers .Like the timescale for kinetic energy dissipation, themixing timescales in our 256 simulations with 1/3 forc-ing energy in the compressible modes are the same asthose given in Fig. 2 from pure solenoidal driving. Al-though one may suspect the larger overall compressiblefraction from the compressible driving scheme could af-fect the mixing efficiency, that is not true because thescalar cascade is controlled by the velocity fluctuations inthe inertial range. Because the compressible energy frac-tion in the inertial range is independent of the compress-ible fraction in the driving force, the mixing timescaleremains the same.The inset in Fig. 2 shows the ratio of the scalar and In § M goes from 3 to 6 .
1, which may tend to givea slight increase in the mixing timescale. The actual decrease inthe normalized mixing timescales in this range of M means thatthis effect is weaker than the opposite effect due to the increas-ing contribution from compressible modes to gradient amplificationdiscussed here. an & Scannapieco 11 TABLE 2Energy dissipation and mixing timescales
M τ diss /τ dyn τ v /τ dyn τ c /τ dyn ≤ k/ π ≤ . ≤ k/ π ≤ . . ≤ k/ π ≤ . . ≤ k/ π ≤ . kinetic dissipation times, τ c /τ v . Because of the contrarybehavior of these two timescales, their ratio increases byabout a factor of 2.6 from incompressible turbulence to M = 3, and saturates at about 1.1 for large Mach num-bers. The filled circle is taken from the simulation resultsby Donzis et al. (2005) for incompressible turbulence andis consistent with our results.In summary, we find that the role of compressiblemodes is completely different for mixing than for kineticenergy dissipation. By providing an additional channelfor energy dissipation, compressible modes lead to fasterenergy loss at larger Mach numbers. On the other hand,mixing tends to be relatively slower if more kinetic en-ergy is contained in compressible modes, suggesting thatthey are not efficient in producing small scalar struc-tures. However, since the majority of kinetic energy iscontained in solenoidal modes even in flows at very largeMach numbers, the mixing rate changes only weakly atincreasing Mach numbers. Thus, the mixing timescaleonly increases by 20% from M = 0 . M = 3 andalways remains comparable to the dynamical timescale.In fact, although their detailed behaviors as a functionof M are different, the overall timescales for the scalardissipation and for the kinetic energy loss are similar atall Mach numbers. This suggests that the production ofscalar structures at small scales occurs through a cascadeprocess similar to that of kinetic energy. In this case, be-cause the scalar and velocities cascades start at the samephysical scale, their decay timescales are both expectedto be close to the turnover time of the largest turbulenteddies, which is essentially the dynamical timescale. Thisis exactly what we find in our simulations. Scale Dependence of the Mixing Timescale
In Fig. 3, we show the dependence of τ c on the charac-teristic length scale of scalar sources. The three curvesare for three scalars forced in different wave numberranges: 3 . ≤ k/ π ≤ . , . ≤ k/ π ≤ . , and16 . ≤ k/ π ≤ .
5. The three mixing timescales arenormalized to that of the scalar forced in the same wavenumber range as the flow driving. The correspondingmixing timescales normalized to the dynamical time aregiven in Table 2.The results in Fig. 3 suggest that the mixing timescaleis proportional to the turnover time of eddies at thelength scale of the scalar sources. At all Mach numbers,the mixing timescales decrease with decreasing sourcelength scale. This is because the turnover timescalein smaller eddies is shorter, and thus, starting from asmaller source scale, the cascade to the diffusion scaleis faster. Fig. 3 also shows that, with increasing Machnumbers, the scale dependence becomes weaker, whichcorresponds to the steepening of the velocity structure t c / t c ( £ k / p £ ) M £ k /2 p £ £ k /2 p £ £ k /2 p £ Fig. 3.—
The dependence of the mixing timescale, τ c , on thescalar source spectrum. The timescales are normalized to τ c in thecase in which the scalar and the flow are driven at the same scales.The three lines are results from different source spectra, whosecharacteristic length scales, L s , are, respectively, 0.25 (solid), 0.124(dashed), and 0.061 (dotted) times the computation box size. functions (see § L s = 0 .
124 and 0.061, which lie in the inter-tial range. Roughly fitting the mixing timescales forthese two scalars with a power law, τ c ∝ L α s , gives α = 0 . , . , . , . , . , .
52 for Mach numbersfrom 0.9 to 6.1. As we will see below, these numbersare quite close to the scaling exponents of the turnovertime, l/δv ( l ), as a function of the eddy size l , using theinertial-range scaling of δv ( l ) . We thus conclude that themixing timescale is essentially given by the turnover timeof eddies at the scalar injection scale.
Structure Functions and Power Spectra
In Fig. 4 we show the second order structure functionsfor the velocity and scalar fluctuations. Note that theseresults are from volume-weighted averages, although us-ing density weighting would be more consistent with thedissipation timescales for the density-weighted velocityand scalar variances studied above. However, unlike 1-point statistics, the choice of an appropriate density-2 Mixing in Supersonic Turbulence S c ( l ) ll l Mach S ll ( l ) ll l Mach S c ( l ) l -2/3 S ll ( l ) l -2/3 Fig. 4.—
The longitudinal velocity structure function (left panel) and second order scalar structure function (right panel) at differentMach numbers. The structure functions at M = 0 . l = 256, and, for clarity, the curves are shifted upwardfor larger Mach numbers. The insets show the structure functions compensated by the l / scaling for incompressible turbulence, i.e., theObukohov-Corrsin scaling for the scalar structure function and the Kolmogorov scaling for the velocity structure function. -8 -6 -4 -2
1 10 100 E c ( k ) k /2 p k -1.70 k -1.72 M = 6.1 M = 0.910 -8 -6 -4 -2
1 10 100 E ( k ) k /2 p k -2 k -1.72 M = 6.1 M = 0.9 E c ( k ) k E ( k ) k Fig. 5.—
The scalar spectra (left panel) and the velocity spectra (right panel) at different Mach numbers. The spectra at M = 0 . k − / scaling for incompressible turbulence. weighting scheme is not trivial for 2-point statistical mea-sures, and will be pursued in future studies. Neverthe-less, the volume-weighed structure functions and spectrashown here uncover interesting differences between thevelocity and scalar structures.The curves in the left panel of Fig. 4 are the longitu-dinal velocity structure functions, S ll ( l ) ≡ h [( v ( x + l ) − v ( x )) · l l ] i , at different Mach numbers. The separation l is in units of the resolution scale. At M = 0 .
9, anextended inertial range (about a decade) exists alreadyat the resolution of 512 . In this range, the structurefunction is well fit by a power-law with a slope of 0.67,consistent with the Kolmogorov scaling for incompress-ible turbulence.At larger Mach numbers, the structure function is sig-nificantly steeper, as found in previous studies (e.g., Krit-suk et al. 2007). We find that the scale range with goodpower-law fits decreases with M , such that at M = 6 . l from 20-70, farfrom both the dissipation range and the driving scale. Inthis range, the exponents, ζ v , are 0.67, 0.73, 0.83, 0.89,0.92 and 0.95 for Mach numbers from 0.9 to 0.61.The structure functions steepen with M for two rea-sons. First, the number of shocks and the average shockintensity both increase with increasing Mach number,and shocks tend to give rise to a structure function thatscales linearly with separation. This is why at largeMach numbers the scaling exponent approaches 1, whichis sometimes referred to as the Burgers scaling (e.g., Krit-suk et al. 2007). Second, as discussed earlier, the pressureterm preferentially converts kinetic energy into thermalenergy. Thus along the energy cascade some kinetic en-ergy may be lost by the pdV work, resulting in smallervelocities toward small scales, and thus a steeper struc-ture function. This effect may peak at M ≈ S c ( l ) ≡ h ( C ( x + l ) − C ( x )) i , is shown in the rightpanel of Fig. 4. Its behavior as a function of the Machnumber is very different from the velocity structure func-tion. The scalar structure function first flattens as theMach number increases from 0.9 to 2. In this range ofMach numbers, we find an extended inertial-convectiverange for the scalar structure function, and the scalingexponents can be accurately measured. The scaling ex-ponent ζ c is 0.67 at M = 0 .
9, in agreement with theObukohov-Corrsin scaling, and it decreases to 0.62 and0.59 at M = 1 . .
1, respectively.These values are in good agreement with the predic-tion from the cascade picture, represented by eq. (4).As found earlier, the scaling exponents for the velocitystructure function are ζ v = 0 .
73 and 0 .
83 at M = 1 . M = 2 .
1. Therefore, eq. (4) would predict ζ c = 0 .
64 and0.59 at these two Mach numbers, which are very closeto the values measured from the simulation data. Thisconfirms the validity of the cascade picture for scalars inflows with M ∼ <
2. Physically, the flattening of the scalarstructure function with M occurs because, as the scalingof δv ( l ) steepens at larger Mach number, there is rela-tively less velocity power at small scales, and thus thecascade of scalar fluctuations becomes slower, leading toa flatter scalar structure function and scalar spectrum.On the other hand, eq. (4) is no longer valid at Machnumbers larger than 2 where the scalar structure functionstarts to become steeper. The scaling exponents are ζ c =0.60, 0.67, and 0.72 for M = 3, 4 .
6, and 6.1, respectively.As in the velocity case, the inertial-convective range be-comes smaller at larger Mach numbers, and these expo-nents are measured in the range of l ∈ [20-70]. Moreaccurate measurement for these large Mach numbers re-quires higher resolution.The reason for the steepening of the scalar structurefunctions is probably different from that for the velocitystructure functions. First, for scalars there is no coun-terpart for the effect of pdV work in the velocity case,which, as discussed above, can cause kinetic energy lossalong the cascade and make the velocity scaling steeper.Instead, the scalar energy is exactly conserved along thescalar cascade. Second, although at large Mach num-bers scalar “edges” can contribute to the steepening ofthe scalar structure function, it is not clear how muchcontribution these edges can make. Unlike the velocityshocks, they are not discontinuities (see § M ≃ M ∼ > . The topic of density-weighting forstructure functions is out of the scope of the current pa-per and will be explored separately in a future study.To summarize these results, we find that for M ∼ < M > M ∼ >
3, whichis different from simulations with less diffusive numeri-cal codes (e.g., Kritsuk et al. 2007; Lemaster and Stone2009). As mentioned above, strong numerical diffusionat high Mach numbers may be responsible for the reduc-tion in the inertial scale range with good power-law fitsfor the structure functions shown in Fig. 4. However, asdiscussed below, we find evidence that strong numericaldiffusion in the FLASH code at large M does not affectour conclusions on the behavior of the power spectra andstructure functions in the inertial range. For the moment, our speculation that density weighting couldremove the steepening at large Mach number is supported by pre-liminary calculations in which the density-weighting factor is takento be the average of the densities at the two points. With thisweighting, we found that the 2nd order scalar structure functionsdo not significantly steepen at M ∼ >
3. However, applying thesame weighting to the velocity structure functions seems to givestrange results. To find an appropriate density weighting factorfor both the velocity field and the scalar field, further investiga-tion is needed, which should perhaps aim at determining whethera scheme exists in which eq. (4) is valid for all Mach numbers. M , while the scalar spectraflattens at first and starts to be steeper at M ∼ >
3. Alsosimilar to the structure functions, the scalar spectrumis significantly flatter than the velocity spectrum at allMach numbers except the M = 0 . M = 6 to those from simulations byKristuk et al. (2007) using less diffusive codes, who findprominent bottlenecks in their velocity spectra.Kritsuk et al. (2007) found a slope of 0.95 for thevolume-weighted longitudinal structure function in theinertial range of a M = 6 flow. This slope is exactlythe same as that measured in the scale range of [20-70]resolution scales in our simulations at the same M . Thissuggests that the scale range we chose is essentially inde-pendent of bottlenecks. As seen in Fig. 7 of Kritsuk et al.(2007), bottlenecks may contaminate the velocity spectraonly below ∼
30 resolution scales. At those small scales,i.e., below 20-30 times the resolution scale, the velocitystructure function in our M = 6 flow is steeper thanthat in Kristuk et al. (2007), consistent with strongernumerical diffusion in the FLASH code at large M . Inconclusion, above ∼
20 resolution scales, the slopes ofthe structure functions are not affected by bottlenecks,and the steepening trend at large M found in our sim-ulations is realistic. This conclusion also applies to thescalar structure functions since the numerical diffusionfor scalars in the FLASH code is essentially the same asthe numerical viscosity.Furthermore, a comparison of the velocity spectrum inour M = 6 flow to that in Kritusk et al. (2007) showsthat the stronger numerical viscosity in our simulationsmakes the spectrum slightly steeper (by a few percent).However, this uncertainty of a few percent is smaller thatthe slope change from M = 3 to M = 6 . M is also realistic. Scalar PDF
In this subsection, we study the probability distribu-tion function (PDF) of forced scalars in compressible tur- bulence. Here we first calculated the 1-point density-weighted PDF, P ( C, t ), in each snapshot by P ( C, t ) = V R V δ ( C − C ( x , t ))˜ ρ ( x , t ) d x where V the volume of thesimulation box. Then for each Mach number we ob-tained the average PDF, P ( C ), over 10 snapshots cov-ering about 10 dynamical times. The results for Machnumbers 0.9 and 6.1 are shown Fig. 6, where all the PDFsare normalized to have unity variance. Note that thePDFs shown here are not the PDFs of the scalar differ-ences between two points at different separations, whichwill not be considered in the present paper.The four lines in each panel of this figure correspond todifferent forcing schemes for the scalars. We find that,although the scalar source term is set to be Gaussianin all our simulations, the scalar PDF is Gaussian onlyif the scalar is forced at the flow driving scale. Thisagrees with the Gaussian scalar PDF found in Watanabe& Gotoh (2004) for incompressible turbulence where thescalar source was also injected at the flow driving scale.For scalars forced at smaller scales, the PDFs show broadnon-Gaussian tails, which become broader with decreas-ing source length scale at all Mach numbers.This behavior is likely to be related to the PDF of thevelocity difference as a function of separation. The veloc-ity PDF at the driving scale is nearly Gaussian, and thePDF for the scalar forced at this scale is nearly Gaus-sian as well. On the other hand, the PDF of the ve-locity difference across a small separation is known tohave broad tails, which become broader as the separationdecreases, a phenomenon known as intermittency (e.g.,Frisch 1995). Since the scalar cascade is controlled bythe velocity field, such small-scale non-Gaussian veloc-ity structures are likely to leave a signature on the PDFtails of scalars forced at small scales. Furthermore thesmall-scale velocity structures are known to become moreintermittent at larger Mach numbers (see, e.g., Padoanet al. 2004). Again, this may be directly responsible forthe fact that the PDF tails for scalars forced at smallscales broaden with increasing Mach number, as can beseen by comparing the two panels in Fig. 6.We point out that, although the non-Gaussian veloc-ity structures contribute to the broad tails in the 1-pointPDFs of scalars forced at small scales, it is not clearwhether the non-Gaussianity in the velocity structures isthe only cause for non-Gaussian tails in the scalar PDFs.In other words, it is possible that the scalar PDFs wouldexhibit non-Gaussian tails even if the flow velocity wereGaussian at all scales, and that non-Gaussian velocitystructures only enhance an already existing effect. Infact, previous studies of mixing in incompressible tur-bulence show that at small scales the scalar differencePDF is non-Gaussian even if the advecting velocity fieldis exactly Gaussian (see, e.g., Shraiman & Siggia 2000),meaning that the non-Gaussianity in the scalar differencestatistics is an intrinsic feature of mixing itself. Whetherthis is also true for the 1-point PDFs of scalars forced atsmall scales is a question for future studies. Decaying Scalars
Finally, we consider the evolution of decaying scalarfields with no continuous sources (i.e., pollutants addedto the flow only once). As mentioned earlier, we includedfour decaying scalars in our simulations, each of whichhas a different initial spectrum (or length scale). Fig.an & Scannapieco 15 -6 -5 -4 -3 -2 -1 -10 -5 0 5 10 P ( C ) C M = 0.9 £ k /2 p £
3 3.5 £ k /2 p £ £ k /2 p £ £ k /2 p £ -6 -5 -4 -3 -2 -1 -10 -5 0 5 10 P ( C ) C M = 6.1 £ k /2 p £
3 3.5 £ k /2 p £ £ k /2 p £ £ k /2 p £ Fig. 6.—
The concentration PDFs for scalars forced at different length scales. The left panel and the right panel are for M = 0 . M = 6 .
1, respectively. The PDFs are normalized to have unit variance. Different curves correspond to scalars forced at different ranges ofwave numbers. -5 -4 -3 -2 -1
0 1 2 3 Ær ~ C æ t / t dyn M = 0.91 £ k /2 p £
3 3.5 £ k /2 p £ £ k /2 p £ £ k /2 p £ -5 -4 -3 -2 -1
0 1 2 3 4 Ær ~ C æ t / t dyn M = 6.11 £ k /2 p £
3 3.5 £ k /2 p £ £ k /2 p £ £ k /2 p £ Fig. 7.—
The concentration variance as a function of time for decaying scalars in turbulent flows at M = 0 . M = 6 . M = 0 . M = 6 .
1. The results for other Machnumbers are qualitatively similar. The variance decayas a function of time and its dependence on the initialscalar scale is similar to that found by Eswaran and Pope(1988) for the incompressible case (see their Fig. 10).Each curve in Fig. 7 has a short initial period whereit is flat and smooth. During this period, the cascadedevelops scalar structures toward smaller scales, and truemixing by molecular diffusion (numerical diffusion in oursimulations) simply waits for the structures to reach thediffusion scale. Once structures of the size of the diffusionscale appear, the variance starts to decay steadily.We find that the concentration variance can be approx-imately fit by a single exponential function for decayingscalars with initial length scale equal to the flow driv-ing scale (see 1 ≤ k/ π ≤ . − . τ dyn for M ≈ M is similar tothat for the forced cases, i.e., it increases with M at small M and becomes constant at M ∼ >
3. We note that thesetimescales are smaller than those for the correspondingforced cases (see Table 2). The reason for this is that,without continuous sources at large scales, the spectrumof a decaying scalar is flatter than that of a forced scalar,and, with relatively more fluctuations at smaller scales,the mixing process is faster.For scalars with initial scales smaller than the flowdriving scale, the variance evolution cannot be fit by asingle exponential function (with the 3 . ≤ k/ π ≤ . M = 6 being the only exception). The decaygenerally consists of two phases, an early fast phase anda slower phase that starts at a few dynamical timescales.In the fast decay phase, the timescale for the variancedecay decreases with decreasing initial scalar scale, con-sistent with the scale dependence of the mixing timescalein the forced case. For example, roughly fitting exponen-tials to the fast phases of the bottom three curves in the M = 6 . . τ dyn , 0 . τ dyn and 0 . τ dyn . Again, these timescales for the variance de-cay are smaller than the corresponding timescales for theforced cases. In the fast phase, the scalar fluctuations are6 Mixing in Supersonic Turbulencesignificantly homogenized, with the variance reduced bya factor of 100 -1000 in a few dynamical timescales.After that, a slower phase starts, which correspondsto the evolution of the scalar spectrum. We find thatif the initial scalar spectrum is within the inertial rangeof the flow, the spectrum spreads out in both directions,i.e., toward smaller scales, due to the cascade, and alsotoward larger scales. The spread to larger scales is dueto the fact that interactions of scalar fluctuations at agiven scale with velocity fluctuations at larger scales canproduce scalar fluctuations at larger scales. We refer tothis effect as the backward transfer. In a few dynami-cal timescales, the scalar power at small scales is signifi-cantly erased by mixing, and the fluctuation power startsto be dominated by structures at larger scales causedby the backward transfer. This results in an increasein the characteristic length scale for the scalar fluctua-tions, leading to slower mixing rate at later times. Thebackward transfer is expected to stop when the typicalscalar length scale becomes close to the driving scale ofthe flow, L f , beyond which there are no considerable ve-locity power. Therefore, in the later phase, the scalarlength scale is always close to the flow driving scale, andthe decay time is approximately given by the turnovertime of the eddies at the flow driving scale. This is whyat late time the decay rates for all scalars tend to becomesimilar.This picture also explains why the variance of a decay-ing scalar can be approximately fit by a single exponen-tial function if its initial length scale is close to the flowforcing scale. In that case, the characteristic length scaleof the scalar is fixed around that flow driving scale anddoes not change with time.In summary, the timescale for the variance decay issmaller by 30-40 % than that obtained in the correspond-ing forced case. It is about 0.5-0.6 dynamical timescaleif the initial scalar length scale is close to the flow driv-ing scale. The dependence of the decay timescale on theinitial length scale and on the Mach number for decay-ing scalars is generally consistent with that in the forcedcase.Fig. 8 shows the evolution of the density weightedscalar probability distribution, P ( C, t ), for the scalarwhose initial spectrum is in the range 1 ≤ k/ π ≤ DISCUSSION: THE EFFECT OF COMPRESSIBLE MODES
Our calculations for the energy dissipation timescaleand the mixing timescale in § §
4, every shock is a strongly dissipativestructure for kinetic energy, and thus the existence ofshocks significantly contributes to energy dissipation. Onthe other hand, the effect of shocks on mixing is muchweaker. Shocks do not give rise to scalar discontinu-ities, instead they amplify scalar gradients by a finitefactor which is equal to the density jump. For a shock toproduce a scalar structure that can dissipate the scalarfluctuations at the same level as it dissipates energy, ithas to be strong enough to reduce the scalar length scalean & Scannapieco 17 -3 -2 -1 -4 -2 0 2 4 P ( C ) / P m a x CM = 0.91 £ k /2 p £ t = 0 0.8 t dyn t dyn t dyn -3 -2 -1 -4 -2 0 2 4 P ( C ) / P m a x CM = 6.11 £ k /2 p £ t = 0 1.0 t dyn t dyn t dyn Fig. 8.—
The density-weighted concentration PDF as a function of time for a decaying scalar in turbulent flows at M = 0 . M = 6 .
1. The scalar field has an initial spectrum in the 1 ≤ k/ π ≤ to around the diffusion scale, so that molecular diffusioncan efficiently mix right after the compression. Since thediffusion scale is typically much smaller than the scalarinjection scale, only the strongest shocks or successivecompressions by strong shocks from different directionscould possibly produce such structures.We find that such strong shocks are very rare. Con-sider our simulations for example. For the scalar forcedat the flow driving scale, a compression by a factor of ∼
100 is needed to bring a scalar structure at the sourcescale down to the diffusion scale, since the diffusion scalein our simulations is around the resolution scale, and theresolution is 512 . The frequency of such strong shocksmay be characterized by the fraction of regions with den-sity 100 times larger than the average because the scalarscale reduction by shocks follows the density jump. Us-ing the density PDF in our simulated flow at Mach 6, themass fraction of these regions is ∼ − . This is muchsmaller than the mass fraction of regions compressed byintermediate or weak shocks. For example, regions com-pressed by shocks with a jump factor larger than 10 makeup a mass fraction of about 0.1, 3 orders of magnitudelarger than that of regions affected by shocks with a jumpfactor of 100. Therefore, strong shocks that can directlyproduce structures at the diffusion scale are rare. Theprobability is even smaller at lower Mach numbers. Inrealistic environments with much larger Reynolds andPeclet numbers, the diffusion scale relative to the sourcescale is much smaller, while the density jump is proba-bly independent of the viscosity or the Reynolds number.Therefore at realistic Reynolds and Peclet numbers, therewould be essentially no shocks that are strong enoughto reduce the scalar structure directly to the diffusionscale. We thus conclude that shocks generally do notproduce scalar structures that dissipate scalar fluctua-tions as strongly as they dissipate kinetic energy. This isresponsible for the different roles of compressible modesin energy dissipation and scalar dissipation.We next show that the compressible modes are less ef-ficient at enhancing mixing than the solenoidal modes.As mentioned earlier, compressible modes include bothcompressions and expansions, and compressions amplifyscalar gradients, while expansions reduce the scalar gra- dients. Because the impact of compressions and expan-sions on the scalar gradient is proportional to the densitychange, the overall effect of compressible modes on thedensity-weighted gradient variance, h ˜ ρ∂ i C∂ i C i , may beestimated as h ρ i / ¯ ρ (where one factor of ρ/ ¯ ρ is fromdensity-weighting). We find that this factor is about 90in our flow at M = 6 .
1. This suggests that, along withthe development of the density fluctuations, which takesabout a dynamical time, the variance of the scalar gradi-ents may be amplified by a factor of 90 by compressiblemodes in about a dynamical timescale.Although a factor of 90 seems significant, it is smallcompared to the enhancement of the scalar gradientsachieved by solenoidal modes. Based on the mixingtimescale in incompressible flows, incompressible modescan reduce the scalar structures to the diffusive scale inabout a dynamical timescale. Consider our simulationsagain as an example. The diffusion scale in our simu-lations is roughly 1/100 of the source length scale, andthus the incompressible turbulent flow can increase thescale gradient variance by a factor of 10 in about a dy-namical timescale. This is 100 times more efficient thancompressible modes in the M = 6 . ∼
1% inour simulations for M = 6 .
1. It is even smaller at smaller M where the amplitude of density fluctuations is smaller.As pointed out earlier, this may be responsible for theslight decrease of the normalized mixing timescale in therange 3 ≤ M ≤ (e.g.,Lemaster & Stone 2008). The more efficient generationof small structures by solenoidal modes is probably be-cause stretching by shears and vortices can continuously8 Mixing in Supersonic Turbulencereduce the scalar structures to progressively small scales.On the other hand, the degree of compression is limitedby gas pressure, which prevents the same fluid elementfrom being compressed continuously, and the existenceof expansions also tends to counteract the effect of com-pressions. This suggests that the solenoidal modes aremore efficient at amplifying the scalar gradients than thecompressible modes, which is responsible for the decreaseof the mixing efficiency when more kinetic energy is con-tained in compressible modes. CONCLUSIONS
Mixing in compressible turbulence plays a key role indetermining, for example, the metallicity dispersion ofstar forming regions, the abundance scatter in coeval fieldstars, and the transition from primordial to PopII starformation. Yet, there have been surprisingly few theo-retical or numerical studies of this process.We have conducted the first systematic numericalstudy of the physics of mixing in supersonic turbulence.Here we give a summary of main results of our study.1. If the typical length scale of the scalar source isclose to the flow driving scale, the mixing timescaleis similar to the timescale for the kinetic energy dis-sipation at all Mach numbers. Furthermore, bothof these timescales are on the order of the dynami-cal timescale, which suggests that the generation ofsmall-scale scalar fluctuations is through a cascadesimilar to that of the kinetic energy.2. The fraction of kinetic energy contained in com-pressible modes increases with Mach number for M ∼ <
3, and becomes constant at larger M . Forkinetic energy in the inertial range, the fractionis 1/3 for M ∼ >
3, indicating an equipartition be-tween compressible modes and solenoidal modes.Equipartition is established at each wave numberin the inertial range of flows with M ∼ > M ∼ <
3, and is essentially constant at largerMach numbers where the fraction of kinetic energycontained in compressible modes is constant. Com-pressible modes are less efficient at enhancing mix-ing, because they have a much lower efficiency atamplifying scalar gradients and producing small-scale scalar structures than the solenoidal modes,which are the primary “mixer” at all Mach num-bers. Solenoidal modes contain most of kinetic en-ergy even at very high Mach numbers, and the dif-ference in the normalized mixing timescale at allMach numbers is smaller than ≈ . ∼ < M ∼ < .
1. Thefraction of energy loss by pdV work peaks at M ≈ / M ≃ M are consistent with the prediction from the cas-cade picture with the measured exponents for thevelocity structure functions. As the Mach numberincreases further, the scalar structure function be-comes steeper, and the relation between the scalarand the velocity scalings from the cascade pictureis not obeyed. However, this probably does notmean the scalar cascade picture is not valid; insteadit motivates a future study to find an appropriatedensity-weighting scheme to compensate the effectof compressions and expansions on the scalar andvelocity structure functions.9. The PDFs of the forced scalars are generally non-Gaussian even if the scalar sources are set to beGaussian in our simulations. The scalar PDF isfound to be Gaussian only for the scalar forced atthe same length scale as the flow driving scale. ThePDFs show broad tails for scalars forced at smallerscales, and the tail broadens with decreasing sourcelength scale. The tails are also broader at largerMach numbers.10. For decaying scalars, the timescale for the scalarvariance decay is slightly faster than the dissipationtimescale in the corresponding forced case. Thedecay timescale for a scalar with an initial lengthscale close to the flow driving scale is in the range0 . − . τ dyn . The dependence of the variance decaytimescale on the initial scalar length scale and onthe Mach number is similar to that in the forcecases. The PDF of a decaying scalar narrows withtime and displays non-Gaussian tails.While these conclusions shed insight on mixing in as-trophysical environments, they also bring up new ques-tions and point the way for further investigations. In thisstudy, we have analyzed only low-order statistics for thevelocity and the scalar structures in simulations with aan & Scannapieco 19moderate resolution of 512 . Examination of higher or-der statistics for these structures require higher resolu-tion simulations, which will help probe deeper into theimportant physics of mixing in supersonic turbulence.We are grateful to Robert Fisher for his helpful discus-sions. We thank the referee, Paolo Padoan, for commentsand suggestions to improve the paper. We acknowledge the support from NASA theory grant NNX09AD106. Allsimulations were conducted on the “Saguaro” cluster op-erated by the Fulton School of Engineering at ArizonaState University. The results presented here were pro-duced using the FLASH code, a product of the DOEASC/Alliances-funded Center for Astrophysical Ther-monuclear Flashes at the University of Chicago.. Examination of higher or-der statistics for these structures require higher resolu-tion simulations, which will help probe deeper into theimportant physics of mixing in supersonic turbulence.We are grateful to Robert Fisher for his helpful discus-sions. We thank the referee, Paolo Padoan, for commentsand suggestions to improve the paper. We acknowledge the support from NASA theory grant NNX09AD106. Allsimulations were conducted on the “Saguaro” cluster op-erated by the Fulton School of Engineering at ArizonaState University. The results presented here were pro-duced using the FLASH code, a product of the DOEASC/Alliances-funded Center for Astrophysical Ther-monuclear Flashes at the University of Chicago.