The turbulence driving parameter of molecular clouds in disc galaxies
MMon. Not. R. Astron. Soc. , 000–000 (0000) Printed 7 July 2020 (MN L A TEX style file v2.2)
The turbulence driving parameter of molecular clouds in discgalaxies
Bastian Körtgen
Hamburger Sternwarte, Universität Hamburg, Gojenbergsweg 112, D-21029 Hamburg, Germany
Released 2020
ABSTRACT
Supersonic turbulence plays a pivotal role during the formation of molecular cloudsand stars in galaxies. However, little is known about how the fraction of compressiveand solenoidal modes in the velocity field evolves over time and how it depends onproperties of the molecular cloud or the galactic environment. In this work, we carryout magnetohydrodynamical simulations of disc galaxies and study the time evolu-tion of the turbulence driving parameter for an ensemble of clouds. We find that thetime-averaged turbulence driving parameter is insensitive to the position of the cloudwithin the galaxy. The ensemble-averaged driving parameter is found to be rathercompressive with b ∼ . − . , indicating almost time-independent global star for-mation properties. However, each individual cloud shows a highly fluctuating drivingparameter, which would strongly affect the cloud’s star formation rate. We find thatthe mode of turbulence driving can rapidly change within only a few Myr, both fromsolenoidal to compressive and vice versa. We attribute these changes to cloud collisionsand to tidal interactions with clouds or overdensities in the environment. Last, we findno significant differences in the average driving parameter between hydrodynamic andinitially strongly magnetised galaxies. However, the magnetic field tends to reduce theoverall fluctuation of the driving parameter. The average driving as well as its uncer-tainty are seen to be in agreement with recent constraints on the turbulence drivingmode for solar neighbourhood clouds. Key words: galaxies: ISM; ISM: kinematics and dynamics; ISM: magnetic fields;ISM: clouds; stars: formation
Stars are born within self-gravitating, magnetised and tur-bulent molecular clouds. The complex dynamical interplaybetween gravity, turbulence, magnetic fields and (at moreevolved stages) stellar feedback strongly shapes the den-sity distribution of the parental cloud and thus determinesthe rate and efficiency at which stars form (Mac Low &Klessen 2004; Elmegreen & Scalo 2004; Scalo & Elmegreen2004; Federrath & Klessen 2012). For supersonically tur-bulent, isothermal (and magnetised) gas the density prob-ability distribution function (PDF) resembles a lognormaldistribution, which means that the logarithm of the den-sity (contrast) is Gaussian distributed (Vazquez-Semadeni1994; Passot & Vázquez-Semadeni 1998; Federrath et al.2008). However, the interstellar medium (ISM) of galaxies isa multiphase gas, strongly influenced by gravitational forces.Hence, the density PDF will take several forms that indi-cate, which physical process is dominant or at least relevant(Kainulainen et al. 2009; Ballesteros-Paredes et al. 2011;Schneider et al. 2015). For a thermally bistable gas, the den- sity PDF becomes double-peaked. Each peak corresponds toa single thermal phase of the gas, and each phase can be fit-ted by a lognormal distribution (Seifried et al. 2011). Whenthe gas has become sufficiently dense, gravity starts to be-come dominant. The subsequent flow of gas towards smallerscales and the resultant formation of power-law density pro-files induces the formation of a power-law tail in the den-sity PDF at high densities (Kritsuk et al. 2011). Over time,the transition point from the log-normal to the power-lawpart in the PDF will shift towards lower densities, as gravitystarts to impact the more diffuse gas in the cloud (Burkhartet al. 2017). In addition, the power-law tail flattens withincreasing star formation activity of the cloud (Federrath& Klessen 2013). However, in case the bistable gas is sub-jected to stellar feedback or strong turbulent phase mixing,a power-law tail can also form in a non-gravitating medium,thereby complicating the interpretation of the density PDF(Seifried et al. 2011; Tremblin et al. 2014).An essential ingredient for current theoretical models of starformation, based on turbulent fragmentation of the gas, is c (cid:13) a r X i v : . [ a s t r o - ph . GA ] J u l B. Körtgen the width of the density distribution (Krumholz & McKee2005; Padoan & Nordlund 2011; Hennebelle & Chabrier2011; Federrath & Klessen 2012; Zamora-Avilés et al. 2012;Konstandin et al. 2016; Völschow et al. 2017). Theoreticalstudies have shown that there exists a relation between thevariance of the density PDF and gas-dynamical quantities,such as the turbulent Mach number, M , the ratio of thermalto magnetic pressure, β and the fraction of compressive orsolenoidal modes in the turbulent velocity field, typically de-noted by b . For a fully compressive driving of the turbulentfield b = 1 , whereas purely solenoidal driving implies b = 1 / (see Federrath et al. 2010, for a first detailled investigation).Therefore, the gas dynamics, which affect the density vari-ance, implicitly affect the star formation properties of thecloud. For example, Federrath & Klessen (2012) have shownthat the star formation rate of a turbulent molecular cloudwith equal Mach number is increased by a factor of (cid:38) ,when the turbulent velocity field is changed from solenoidalto compressive.Rather than treating the driving parameter, b , as an inputparameter, several recent studies focused on retrieving itfrom the self-consistent gas dynamics in the analysed envi-ronment (Körtgen et al. 2017; Jin et al. 2017; Menon et al.2020). While Körtgen et al. (2017) and Jin et al. (2017)find b to vary between b ∼ . and b ∼ . in their dif-ferent simulations of molecular cloud formation, the studyby Menon et al. (2020) finds a rather stable b ∼ . ina simulated region, which is affected by radiative feedback.Furthermore, Pan et al. (2016) argue that supernova feed-back induces mostly solenoidal motions due to the onset ofthe baroclinic instability in the post-shock region of the su-pernova remnant. These results all emphasise the impact ofthe various feedback mechanisms and involved physics onthe turbulent velocity field.Deriving the turbulence driving parameter from observa-tions of molecular clouds is even more challenging. This ismostly due to the restriction of the data to at most twodimensions. However, it is nevertheless useful and neces-sary to confront theory with reality. For molecular cloudsin the solar neighbourhood, Ginsburg et al. (2013) deriveda lower limit of b ∼ . . Orkisz et al. (2017) analysed thesolenoidal fraction of the velocity field in the Orion B cloudand showed that the fraction varies both as a function ofposition across the complex and radial distance to a givencenter (e.g. the position of the peak density). In terms of thedriving parameter, this latter study deduced that b variesfrom b > . to b (cid:46) , i.e. the gas motion is neither fullysolenoidal nor entirely compressive. Their most likely valuewas about b ∼ . , which indicates equipartition between thetwo modes (Federrath et al. 2010). Kainulainen & Federrath(2017) studied a set of 13 solar neighbourhood clouds andargued that b ∈ [0 . , . provides a strict range of values. Asthese authors point out, the driving parameter was derivedfrom the density variance - Mach number relation, despitethe fact that their data did not necessarily follow such arelation. In contrast to clouds in the solar neighbourhood,the Galactic center may provide a significantly different en-vironment. As Federrath et al. (2016) show in their study ofthe central molecular zone cloud G0.253+0.016, the drivingof turbulence in this region might be almost fully solenoidalwith b ∼ . ± . due to the increased shear. Hence, un-derstanding the driving of turbulence in different environ- ments of the Galaxy and within molecular clouds is crucialfor the understanding of the local gas dynamics and the sub-sequent star formation process.This work strives to add further insights into the driving ofturbulence by analysing the time evolution of the turbulencedriving parameter in a set of molecular clouds formed in discgalaxies. The study is structured as follows: In Sect. 2 webriefly discuss the used numerical tools and initial condi-tions. Sect. 3 discusses the findings of this study, before it isclosed with a summary in Sect. 4. We use the finite volume code flash (version 4.2.2, Dubeyet al. 2008). The equations of ideal magnetohydrodynam-ics are solved each timestep on the numerical grid witha five-wave Riemann solver (Bouchut et al. 2009; Waaganet al. 2011). To ensure the solenoidal constraint of the mag-netic field, we use a hyperbolic cleaning scheme (Dedneret al. 2002). Poisson’s equation for the self-gravity of thegas is solved with a Barnes-Hut tree solver (see e.g. Lukat& Banerjee 2016, for a GPU version). The thermodynamicevolution of the gas is treated by optically thin heatingand cooling, where we follow the approach by Koyama &Inutsuka (2002, with modifications by Vázquez-Semadeniet al. (2007)), i.e. the heating rate is constant and thetemperature-dependent cooling rate is provided in tabulatedform.The cubic simulation domain has an edge length of L = 40 kpc and the root grid is at a grid resolution of ∆ x root = 625 pc . Over the course of the evolution, the gridis adaptively refined (Berger & Oliger 1984) once the localJeans length is resolved with less than 32 grid cells (as rec-ommended based on the results in Federrath et al. 2011).De-refinement is initiated when the Jeans length instead isresolved by more than 64 grid cells. In total, we allow for amaximum of 11 refinement levels, which gives a peak spatialresolution of ∆ x peak = 19 . . The gas scale height is thusresolved with (cid:38) − grid cells in the major part of thedisc (Körtgen et al. 2019). Any fragmentation at later timesfrom R ∼ − on can thus assumed to be physical (seealso Truelove et al. 1997). To extract gas dynamics withinthe formed clouds solely driven by gravity, neither stellarfeedback nor sink particles are included in the current studyand we introduce an artificial pressure term on the highestlevel of refinement. An in depth description of the initial conditions is given inKörtgen et al. (2019, see also Körtgen et al. (2018)), so weprovide only a brief overview here.We set up a thin disc, following the approach by Tasker& Tan (2009), with a radially increasing gas scale height.The gas density profile is adjusted in such way that theToomre stability parameter is initially constant at Q = 2 in the major part of the disc. In the inner most part it is c (cid:13) , 000–000 urbulence in MCs set to Q = 20 to increase the stability of the less well re-solved region. The final mass of the disc is M disc ≈ M (cid:12) ,similar to the mass of the LMC. We use a fixed external log-arithmic gravitational potential to account for the effect ofdark matter and old stars, which also provides a flat rotationcurve. The majority of the performed simulations starts witha purely toroidal magnetic field, where the field strength iscoupled to the gas density, i.e. B ∝ (cid:37) / , which results in aconstant ratio of thermal to magnetic pressure. We do notimpose any turbulent perturbation in order to simplify theinitial conditions. Hence, any local velocity fluctuations atlater stages are the result of gravitational fragmentation andself-consistently generated dynamics (i.e. cloud-cloud inter-actions). Fig. 1 presents column density maps of the hydrodynamic(top panels) and the magnetised (bottom) discs, separatedin time by ∆ t = 200 Myr , which corresponds to almost anentire disc rotation at R = 8 kpc .With time, the initially smooth density profile is distorted,because the discs fragment into individual overdensities.While the hydrodynamic disc breaks up into several ringsdue to the onset of the Toomre instability, the magnetisedgalaxy builds up filamentary structures, which extend intothe radial and azimuthal directions. The latter fragmenta-tion pattern is a consequence of the Parker instability asthe dominant mode of fragmentation (Mouschovias et al.2009; Körtgen et al. 2018, 2019). After one additional discrotation the shape of the galaxies is markedly different. Thehydrodynamic disc has broken up into many almost spheri-cal objects. The magnetised galaxy shows more filamentarystructures. Here, the magnetic field restricts the motion ofthe diffuse gas. Hence, gas is not arbitrarily absorbed by theformed clouds. We emphasise the presence of spiral featuresaround some clouds in both galaxies due to the cloud’s bulkrotation. At even later times, the global disc morphologylooks very similar. This indicates that, at this stage, themagnetic field does not have a pronounced impact on theglobal dynamics of the disc. Filamentary structures are nowalso seen in the hydrodynamic galaxy as a consequence ofcloud-cloud interactions.In the top row of Fig. 2 we show the divergence of the ve-locity field at a time close to the end of the simulation forboth discs (dark: diverging, bright: converging). While thisoperation on the residual velocity field generates filamen-tary structures, there still appear pronounced differences.For instance, the hydrodynamic galaxy reveals a much moreconverging velocity field throughout the disc. In contrast,the magnetised galaxy reveals large patches of diverging gasin between kpc-scale filaments with a primarily convergingfield. In addition, the hydrodynamic disc allows for the for-mation of small scale bundles in the divergence field (seee.g. the pattern near x = +6 kpc and y = +2 kpc ), whoseformation is suppressed in the magnetised galaxy.In the bottom panel of Fig. 2, we show the vorticity. Here,again, filamentary structures can be identified. In contrastto the divergence field, these are not as thick, indicating that the vorticity throughout the disc seems to show a quite com-mon orientation. In general, regions of pronounced changesin the divergence or the vorticity correlate well with gasstructures emerging from dynamical interactions, such ascloud-cloud collisions. However, both discs do not show anysystematic trends e.g. with respect to the distance from thecenter of the galaxy. In this section we focus on the time evolution of dynami-cal quantities of the found clouds. Before we commence, wedescribe the cloud tracking algorithm.
At a given time, a domain of size L x × L y × L z = 20 × × centered on the galaxyis searched for connected regions. The regions are definedby a minimum density, n min = 100 cm − , which is aboutthe average number density of molecular clouds in theMilky Way (Blitz et al. 2007; Dobbs et al. 2014). A gridcell belongs to a certain structure if its density is similar(provided some uncertainty) to the one of its nearestneighbours and further has a spatial connection to them.The search starts in the cell with the highest density. Once,all objects with n min , obj (cid:62) n min have been found, theircurrent center of mass position and bulk velocity are usedto follow their evolution in time. At the next timestep –we use ∆ t = 2 Myr –, a volume of V = 1 kpc centeredon the predicted center of mass is searched for clouds andthe cloud with its center of mass closest to the predictedone is chosen (see also Tasker & Tan 2009). We note herethat, in the current version and as a difference to Tasker& Tan (2009) and Jin et al. (2017), the cloud trackingalgorithm does not search for newly formed clouds and thustraces only the evolution of the firstly identified ensembleof clouds. If not stated differently, the left panel in the following discus-sion corresponds to the hydrodynamic disc, while the rightpanel indicates results of the magnetised galaxy. In all fig-ures, the grey lines indicate the evolution of each individualcloud, the red thick line is the cloud-mass weighted average,and the blue dashed lines represents the evolution of an ex-ample cloud .In Fig. 3 we show the time evolution of the average num-ber density of the clouds. The large difference in densitybetween the two scenarios is initially due to the higher den-sity of the magnetised galaxy, which is necessary to keepthe Toomre parameter at Q = 2 . Apart from this initialdifference, the long term evolution of the cloud’s averagedensity is very similar. There appear, however, subtle vari-ations, i.e. the magnetised clouds reach higher densities ofup to (cid:104) n (cid:105) ∼ cm − , whereas there is only one cloud ata single time in the hydrodynamic disc, which reaches suchhigh average densities. The average instead is less affected, To be more precise, it is the same cloud in all plots.c (cid:13)000
At a given time, a domain of size L x × L y × L z = 20 × × centered on the galaxyis searched for connected regions. The regions are definedby a minimum density, n min = 100 cm − , which is aboutthe average number density of molecular clouds in theMilky Way (Blitz et al. 2007; Dobbs et al. 2014). A gridcell belongs to a certain structure if its density is similar(provided some uncertainty) to the one of its nearestneighbours and further has a spatial connection to them.The search starts in the cell with the highest density. Once,all objects with n min , obj (cid:62) n min have been found, theircurrent center of mass position and bulk velocity are usedto follow their evolution in time. At the next timestep –we use ∆ t = 2 Myr –, a volume of V = 1 kpc centeredon the predicted center of mass is searched for clouds andthe cloud with its center of mass closest to the predictedone is chosen (see also Tasker & Tan 2009). We note herethat, in the current version and as a difference to Tasker& Tan (2009) and Jin et al. (2017), the cloud trackingalgorithm does not search for newly formed clouds and thustraces only the evolution of the firstly identified ensembleof clouds. If not stated differently, the left panel in the following discus-sion corresponds to the hydrodynamic disc, while the rightpanel indicates results of the magnetised galaxy. In all fig-ures, the grey lines indicate the evolution of each individualcloud, the red thick line is the cloud-mass weighted average,and the blue dashed lines represents the evolution of an ex-ample cloud .In Fig. 3 we show the time evolution of the average num-ber density of the clouds. The large difference in densitybetween the two scenarios is initially due to the higher den-sity of the magnetised galaxy, which is necessary to keepthe Toomre parameter at Q = 2 . Apart from this initialdifference, the long term evolution of the cloud’s averagedensity is very similar. There appear, however, subtle vari-ations, i.e. the magnetised clouds reach higher densities ofup to (cid:104) n (cid:105) ∼ cm − , whereas there is only one cloud ata single time in the hydrodynamic disc, which reaches suchhigh average densities. The average instead is less affected, To be more precise, it is the same cloud in all plots.c (cid:13)000 , 000–000
B. Körtgen
Figure 1.
Surface density maps of the hydrodynamic (top) and magnetised disc (bottom) at three times, starting from the point wherethe first diffuse ( n min (cid:62)
10 cm − ) clouds have been identified. Note the different morphology of the fragmented disc and the increaseddensity in the diffuse gas for the MHD disc at later times. with a difference of about a factor of two. The overall timeevolution of the average is more pronounced in the hydro-dynamic disc and a plateau is reached after about 200 Myrof evolution (i.e. one disc rotation). As can be seen from theevolution of the example clouds, both ensembles are affectedby cloud-cloud interactions with actual collisions resultingin a sudden increase of the average density . Interestingly,this effect is seen for almost every cloud. However, colli-sions/mergers induce only a temporarily limited increase inthe average density and it takes several interactions for thecloud to become denser on average with time.The time evolution of the velocity dispersion is shown inFig. 4. As expected for unimpeded gravitational collapse,the velocity dispersion increases over time. The retrievedvelocity dispersion evolves very similar in both galaxies, butis shifted to larger values in the MHD disc. The scatter,in contrast, is larger in the hydrodynamic case, which rep-resents a larger variety of cloud properties. While it takesabout one disc rotation to reach saturation in the magne-tised disc, the hydrodynamic clouds still show a net increaseover the course of about one and a half disc rotations andhave not reached a saturated stage yet. Given the fact thatthe velocity dispersion is calculated from the size-linewidthrelation σ v ∼ . R p , with p = 0 . and R being the cloud-size (Larson 1981), the saturation indicates a phase in whichthe clouds do not grow anymore. The typical cloud tempera-ture is about T ∼
30 K , which corresponds to a sound speed We point out that the impact of the artificial pressure term,which results in an adiabatic behaviour in the innermost regionof the core, does not affect the resulting (average) cloud radii anddensities Derived from a size-linewidth relation due to the limited nu-merical resolution and as a guideline for the internal dynamics. of about c s ∼ . / s . Therefore, sonic Mach numbers, asderived from Fig. 4, are in the range of M s ∼ − . Fig. 5 emphasises the time evolution of the cloud’s AlfvénMach number and plasma- β . Similar to the sonic Mach num-ber, the average Alfvén Mach number appears to reach satu-ration after one disc rotation. However, as inferred from theexample cloud, individual clouds can show dramatic changesin the dynamic interplay of turbulence and magnetic fields.In general, all clouds are seen to be super-Alfvénic. This isalso a result of matter accretion along the field lines, whichresults in a net decrease of the Alfvén speed when the fieldis not amplified sufficiently. Please note that the identifiedclouds start out trans- to slightly super-Alfvénic and reachvalues of M A (cid:38) after ∼
40 Myr . The ratio of thermal tomagnetic pressure of the clouds varies between β ∼ . − for the majority of the clouds, with a rather stable averagevalue of (cid:104) β (cid:105) ∼ . − . , which is also in agreement with theconstraints discussed by Kainulainen & Federrath (2017) fora set of solar neighbourhood clouds. Furthermore, the aver-age value also fits the global disc average (Körtgen et al.2019). Having the time evolution of the relevant quantities at hand,one can investigate whether there emerges a relation be-tween these and the density-variance. For the logarithmic We expect the initial conditions to be sub- to trans-Alfvénic forincreasing resolution. c (cid:13) , 000–000 urbulence in MCs Figure 2.
Top:
Divergence of the projected residual (rotation curve subtracted) velocity field for the hydrodynamic (left) and magnetised(right) discs shortly before the end of the simulation. Both discs reveal filamentary structures of the divergence field. However, thereappear more converging regions in the hydrodynamic disc. Note the large patches of diverging gas in the MHD disc.
Bottom:
Vorticity.The colour bar is arranged that converging regions are shown in bright colour, while diverging regions are shown as the dark patches. Inthe bottom panel, the same colour coding holds for negative and positive vorticity.
Figure 3.
Average number density of the tracked molecular clouds as a function of time. In grey we show the evolution of each individualcloud, while the thick red line denotes the average over all clouds. The blue dashed lined highlights the evolution of a single cloud. Leftfor the hydrodynamic case. Right for the magnetised case.c (cid:13)000
Average number density of the tracked molecular clouds as a function of time. In grey we show the evolution of each individualcloud, while the thick red line denotes the average over all clouds. The blue dashed lined highlights the evolution of a single cloud. Leftfor the hydrodynamic case. Right for the magnetised case.c (cid:13)000 , 000–000
B. Körtgen
Figure 4.
Velocity dispersion of the tracked molecular clouds as a function of time. In grey we show the evolution of each individualcloud, while the thick red line denotes the average over all clouds. The blue dashed lined highlights the evolution of a single cloud. Leftfor the hydrodynamic case; right for the magnetised case.
Figure 5.
Left:
Alfvén Mach number of the tracked molecular clouds as a function of time.
Right:
Plasma- β as a function of time. Theline colors and styles have the usual meaning. In the right sub-figure, the black line denotes β = 1 and the green dashed lines indicateupper and lower limits for the magnetisation of solar neighbourhood clouds as constrained in Kainulainen & Federrath (2017). density contrast this relation is given by (e.g Federrath &Klessen 2012, and references therein) σ (cid:37)/(cid:37) ) = ln (cid:18) b M ββ + 1 (cid:19) . (1)We stress here that we remove systematic motions in the ve-locity field of the cloud and correct for a gradient in density,which would affect the density variance (Federrath et al.2011). Fig. 6 shows the resulting density variance - Machnumber relation for the hydrodynamic and magnetised discs.We show all ensemble clouds at all times, where the time iscolour coded. As is observed, a relation emerges in the hy-drodynamic case, which is, however, not described by eq. 1.At first, the density variance is too low for the estimatedMach numbers and, secondly, the relation appears rather lin-ear than logarithmic. Interestingly, in the magnetised coun-terpart, a logarithmic relation is slightly recognised, whichindicates that the overall relation is changed quite signifi-cantly by the magnetic field. The spread across the cloudsis larger, but the overall density variances come closer tothe theoretical values. However, for individual snapshots in time, no significant correlation would be observed. This isespecially clear for the MHD disc, where the majority of datapoints at late times (yellow) show a rather flat distribution,while the hydrodynamic clouds set up a point cloud between M ∼ − without any clear trend. In any case, inboth discs the derived driving parameters would be far toolow. This has recently been discussed in Jin et al. (2017).The authors showed that the derived b is only a lower limitand it was seen to increase with increasing numerical reso-lution. We stress here that also no correct relation emerges,when using the velocity dispersion from the size-linewidthrelation. Since Fig. 6 indicates the lack of a relation between the den-sity variance and Mach number as provided by analyticalmodels of turbulent fragmentation (Molina et al. 2012), werefrain from using this relation. We show in the appendix c (cid:13) , 000–000 urbulence in MCs Figure 6.
Density variance – Mach number relation at different times (blue: early, yellow: late). In the right panel, we show themagnetised version, which accounts for the plasma- β . Although the magnetised disc is closer to the theoretical prediction, no relationbetween the density variance and the Mach number term arises. In contrast, a clear relation is seen in the hydrodynamic case. However,the non-converged velocity field and the resulting inaccurate Mach numbers yield a relation that is far off the theoretical ones (purplelines. Dotted: b = 1 , solid: b = 0 . , dashed: b = 1 / ). below that the density PDFs are dominated by a power-lawtail. Hence, the application of the log-normal density vari-ance should be taken with caution. However, using insteadthe density variance for the linear density contrast, (cid:37)/(cid:37) , σ (cid:37)/(cid:37) = b M ββ + 1 (2)gives b -values, which are still too low. We re-iterate thatthe reason for this discrepancy is the lack of numericalresolution in our simulations, as previously found in Jinet al. (2017). We thus concentrate in the following on thedriving parameter as derived from the compressive ratio χ = (cid:10) v (cid:11) / (cid:10) v (cid:11) , where v comp and v sol are the compres-sive and solenoidal components of the velocity field. Thesecan be retrieved via a Helmholtz decomposition in Fourierspace. It can then be shown (Pan et al. 2016) that the com-pressive ratio and the turbulence driving parameter are re-lated via b χ = (cid:114) χ χ . (3)We point out that this relation provides only an approxi-mation to the true value of b , which can be retrieved fromequations 1 and 2. We further note here that we remove anysystematic motion in the velocity field, i.e. bulk motions inspherical shells around the center of mass. Thus, the com-pressive ratio in this study is comparable to the turbulentcompressive ratio given in Pan et al. (2016, see also (Jinet al. 2017; Mandal et al. 2020; Menon et al. 2020)).The resulting time evolution of the cloud-mass weightedaverage driving parameter, as well as of its standard de-viation, is provided in Fig. 7. Surprisingly, there is almostno pronounced difference between the hydrodynamic (left)and magnetised clouds (right). The magnetised driving pa-rameter appears to be slightly higher, but is still around b χ ∼ . . A closer look at the early evolution reveals thatthe magnetised parameter starts fluctuating earlier, after Note that we use b χ and b driv interchangeably in the figures. ∼
100 Myr , while the hydrodynamic driving parameter keepsits initial value for over half a rotation of the galaxy (againestimated at R = 8 kpc ). Over the course of the evolution,the average b χ starts to fluctuate quite strongly, therebytaking values from fully solenoidal to strongly compressive.However, despite the large fluctuations, both average drivingparameters are within the bounds estimated by Kainulainen& Federrath (2017, green dashed lines). Given the fact thatwe only include gravity here, the tendency towards com-pressive driving is not surprising. Please note that, as sum-marised in Federrath et al. (2017), stellar feedback is alsolikely to drive compressive turbulent motions.The thin blue line underlines the time evolution of thedriving parameter for an individual cloud. Here, large fluc-tuations and deviations from the average are evident. Forexample, the hydrodynamic example cloud shows a stagenear t ∼
370 Myr , where its turbulence would be describedas fully compressive. In contrast, not even 30 Myr later,around t ∼
400 Myr , its driving parameter has dropped to b χ ∼ . , which indicates purely solenoidal driving. Similarshort-period variations are observed for the example cloud inthe magnetised disc. Both example clouds furthermore showa long-periodic modulation, with the magnetised modula-tion being less clear. We point out that not all clouds showsuch clear modulation patterns. However, the short- andlong-period features can be interpreted as local and global(i.e. galactic-scale) processes affecting the velocity field inthe clouds.Fig. 8 analyses the behaviour of the turbulence driving modeas a function of distance to the galactic center. The variouslines denote time averages over a period of 50 Myr or overthe full tracking period. Although the spatial distributionof clouds in the hydrodynamic disc is narrower, there is nodifference to the magnetised ensemble of clouds. In addi-tion, we do not find any systematic trends with positionin the galaxy, although the large fluctuations hint towardssmall temporal variations, which could indicate a trend atthe given time. Contrary to the missing radial trend, a fewsmall differences can be recognised. Initially, the clouds in c (cid:13)000
400 Myr , its driving parameter has dropped to b χ ∼ . , which indicates purely solenoidal driving. Similarshort-period variations are observed for the example cloud inthe magnetised disc. Both example clouds furthermore showa long-periodic modulation, with the magnetised modula-tion being less clear. We point out that not all clouds showsuch clear modulation patterns. However, the short- andlong-period features can be interpreted as local and global(i.e. galactic-scale) processes affecting the velocity field inthe clouds.Fig. 8 analyses the behaviour of the turbulence driving modeas a function of distance to the galactic center. The variouslines denote time averages over a period of 50 Myr or overthe full tracking period. Although the spatial distributionof clouds in the hydrodynamic disc is narrower, there is nodifference to the magnetised ensemble of clouds. In addi-tion, we do not find any systematic trends with positionin the galaxy, although the large fluctuations hint towardssmall temporal variations, which could indicate a trend atthe given time. Contrary to the missing radial trend, a fewsmall differences can be recognised. Initially, the clouds in c (cid:13)000 , 000–000 B. Körtgen
Figure 7.
Turbulence driving parameter b as derived via the compression ratio χ as a function of time for the ensemble of clouds identifiedin the hydrodynamic (left) and MHD (right) disc. Although strong variations per cloud from purely solenoidal to largely compressiveare seen, the mass-weighted average driving parameter is about b ∼ . , indicating slightly compressive driving. In blue we show theevolution of b for an exemplary cloud. There are times with little variation and events, where b changes dramatically within only a fewMyr. The stronger fluctuations are due to cloud collisions and tidal interactions. The horizontal black lines denote natural ( b = 0 . ) andpurely solenoidal driving ( b = 1 / ), while the green dashed lines highlight the possible upper/lower limit for solar neighbourhood cloudsdiscussed in Kainulainen & Federrath (2017, i.e. b l = 0 . and b u = 0 . ). the magnetised disc show a slightly lower driving parame-ter . < b (cid:46) . , whereas the hydrodynamic clouds reveal b ∼ . . The fluctuations in b χ are larger in the hydrody-namic ensemble, which was also recognised in Fig. 7. Thisimplies that the magnetic field tends to (slightly) decreasethe overall variation in b χ . Keeping in mind the absence of stellar feedback, the lack of acorrelation between the driving parameter and the positionof the cloud in the galaxy seems to be in contradiction withrecent estimates of the turbulence driving mode in Galacticmolecular clouds, which show mixed to mildly compressiveturbulence in the solar neighbourhood (Ginsburg et al. 2013;Kainulainen & Federrath 2017) and solenoidal forcing nearthe Galactic center (Federrath et al. 2016). Hence, theremust be other mechanisms that induce the derived changesin the forcing over time. In Fig. 9 we show the time evolutionof b χ for the example cloud. We extract a certain intervalof about 30 Myr and present surface density maps centeredon the cloud’s center of mass. The sequence of maps clearlyshows that the cloud undergoes a collision/merger with adifferent, slightly less dense object. The corresponding timesshown in the surface density maps are highlighted as verticallines in the time evolution of b χ . The driving mode is initiallycompressive with b χ ∼ . . Interestingly, the appearance ofthe collider reduces the forcing parameter to b χ ∼ . . Thisimplies that tidal forces induce shear across the cloud, giventhat the net tidal field is disruptive (see e.g. discussion in Jog2013). The later stages, which resemble an in-spiral phaseand the actual collision raise the value to b χ ∼ . . This iswhat should be expected from a collision. After the colli-sion event, however, the driving parameter decreases againto solenoidal values. This is, because the collision was nothead-on and thus increased the shear within the cloud. Wepoint out that this change of b χ from ∼ . to ∼ . and back to solenoidal happens on timescales of less than 10 Myr(in agreement with the findings by Körtgen et al. 2017). Ex-trapolating this to all clouds, we conclude that the natureof the collision event is vital for the short-term temporalvariation of the driving parameter. Finally, we briefly study the influence of the choice of thelower threshold density for the identification of overdense ob-jects. Fig. 10 shows the time evolution of the average drivingparameter of about 16 clouds. The left panel shows the stan-dard threshold density of n = 100 cm − . In the right panelwe have increased the lower threshold by a factor of ten, thuscorresponding to molecular clump or core densities. As ex-pected, the driving mode becomes more compressive for thehigher threshold density since this material is more tightlybound by gravity, although the increase can temporarily berather small (see also Orkisz et al. 2017). The average driv-ing parameter over this period is b cloud ± ∆ b cloud = 0 . ± . for the fiducial threshold and b dense ± ∆ b dense = 0 . ± . in case of the increased lower density bound. Please notethe time evolution of the driving parameter for the exampleclouds. These nicely show that the driving mode can changesignificantly over the course of 25 Myr and also does notrepresent the ensemble well. We present a study on the time evolution of an ensemble ofclouds formed in disc galaxies. We focus on a hydrodynamicdisc and a strongly magnetised disc with an initial plasma- β = 0 . . The clouds are identified shortly after the discshave fragmented on large scales either due to the classicalToomre or the Parker instability. We find that the dynam-ical quantities, relevant for the turbulent properties of the c (cid:13) , 000–000 urbulence in MCs b χ Galactic Distance [kpc] b χ Galactic Distance [kpc]
Figure 8.
Turbulence driving parameter b as derived via the compressive ratio χ as a function of distance from the galactic center atvarious times and averaged over an interval of ∆ t av = 50 Myr . Two facts are observed: 1) The driving parameter fluctuates less in themagnetised case and 2) it is generally slightly smaller than the corresponding hydrodynamic case. There is also only small variation in b over the course of a disc rotation. Note that the magnetised clouds are initially less compressive. clouds, do not vary significantly when including a magneticfield. In a next step, we determined whether there arises a re-lation between the density variance and the turbulent sonicMach number. Both galaxies show some kind of a relation.This, however, does not match the theoretical expectationfor isothermal, turbulent gas, primarily because of the lackof numerical resolution in our simulations. As a consequence,we do not determine the driving parameter from the densityvariance - Mach number relation, but rather via the com-pressive ratio. Our main findings concerning the propertiesand evolution of the driving parameter obtained this wayare: • The derived driving parameter varies between fullysolenoidal ( b χ ∼ / ) and entirely compressive ( b χ ∼ )driving. • The average driving parameter, b χ , does not signifi-cantly vary over the evolution of about ∼ Myr. • In contrast, for individual clouds, we find large fluctu-ations as well as times of little to almost no variation. • The largest fluctuations of b χ can be associated withexternal distortions such as cloud-cloud mergers. • Collisions between clouds induce a rapid change of b χ within only a few Myr. • The tidal forces exerted onto each cloud by its environ-ment can reduce the driving parameter as they increase theshear. • In our current framework, there appear no variationsacross the disc. • Taking the magnetic field into account slightly reducesthe spread/the uncertainty in b χ .We conclude that the merger history and/or the environ-ment play a significant role in shaping the turbulent velocityfield of molecular clouds in galaxies. ACKNOWLEDGEMENT
BK thanks the anonymous referee for her/his comments,which helped to clarify the findings of this study. BK ac-knowledges Paris-Saclay University’s Institut Pascal pro-gram "The Self-Organized Star Formation Process" and theInterstellar Institute for hosting discussions that nourishedthe development of the ideas behind this work. BK furtheracknowledges discussions with Robi Banerjee, Urs Schäfer,Wolfram Schmidt, Simon Selg and Pranjal Trivedi (in alpha-betical order) and funding via the Australia-Germany JointResearch Cooperation Scheme (UA-DAAD) and from DFGgrant BA3706/15-1. The simulations were run on HLRN-IIIunder project grant hhp00043. The flash code was in partdeveloped by the DOE-supported ASC/Alliance Center forAstrophysical Thermonuclear Flashes at the University ofChicago. Figs. 1 and 9 were generated with
GUFY , a graph-ical user interface for the analysis of flash data with yt,developed by F. Balzer (University of Hamburg).
APPENDIX A: DENSITY PROBABILITYDISTRIBUTION FUNCTION
A way to study the dynamics of molecular clouds is byanalysing their (column-) density distribution. The resultsare shown in Fig. A1 for different times. The PDF is almostflat at the earliest time, when the cloud has just formed viacompression of gas flows due to the Parker instability. Overtime, the PDF develops a power-law shape. At some times,a clear log-normal part is observed. In any case, it is clearthat the PDF is more likely to be a power-law than a combi-nation of lognormal and power-law, probably due to the lackof resolution. Hence, deriving the turbulence driving param-eter from the logarithmic density variance - Mach numberrelation should be taken with caution. c (cid:13)000
A way to study the dynamics of molecular clouds is byanalysing their (column-) density distribution. The resultsare shown in Fig. A1 for different times. The PDF is almostflat at the earliest time, when the cloud has just formed viacompression of gas flows due to the Parker instability. Overtime, the PDF develops a power-law shape. At some times,a clear log-normal part is observed. In any case, it is clearthat the PDF is more likely to be a power-law than a combi-nation of lognormal and power-law, probably due to the lackof resolution. Hence, deriving the turbulence driving param-eter from the logarithmic density variance - Mach numberrelation should be taken with caution. c (cid:13)000 , 000–000 B. Körtgen b d r i v Time [Myr] b d r i v Time [Myr]
Figure 9.
Top two panels:
Time evolution of the turbulence driving parameter for the example cloud, where a specific time period ishighlighted.
Bottom panels:
Surface density maps of a . × . area centered on the cloud’s center of mass (indicated by the reddot). Each panel is highlighted by a vertical line in the time evolution plot (sequence: top to bottom and left to right). The sequenceof maps shows a cloud merger and emphasise that this process initiates at first a decrease in the driving parameter due to tidal forcesinducing shear. The first merger process is highly compressive with b ∼ . . Note the inspiral-phase after the closest approach, whichresults in strong shearing within the cloud and subsequently in low values of b . APPENDIX B: RESOLUTION STUDY
In Fig. B1 we study the influence of the numerical resolu-tion on the calculation of the average turbulence drivingparameter, b driv . As investigated in Jin et al. (2017, seealso Körtgen et al. (2017)), full numerical convergence is achieved only with sub-parsec resolution. We show two res-olutions, namely our fiducial resolution with ∆ x = 19 . and a slightly higher one with ∆ x = 9 . . The time axis isgiven in time relative to the time when clouds are identifiedfor the first time, t .The driving parameter in the higher resolution run seems c (cid:13) , 000–000 urbulence in MCs Figure 10.
Dependence of the driving parameter b on changing the density threshold for defining clouds in the simulated galaxies.Evolution of identified clouds at late times in the magnetised disc. The left panel shows clouds with a minimum density of n min =100 cm − , while the right shows clouds with n min = 10 cm − . As expected, the clouds with the higher threshold density show a slightlymore compressive driving parameter, since these are more gravitationally bound/contracting. This is in agreement with Orkisz et al.(2017), who find a smaller solenoidal fraction on smaller scales. However, the time averages differ only little with b cloud ± ∆ b cloud =0 . ± . and b dense ± ∆ b dense = 0 . ± . . N b i n ln( ρ / ρ ) t=260 Myrt=360 Myrt=460 Myrt=560 Myrt=610 Myrt=660 Myr Figure A1.
Probability distribution function (PDF) of the loga-rithmic density contrast for the example cloud at different times.The PDF is initially flat, because the cloud has just condensedout of the diffuse gas and has not developed any substructure, yet.With time, the PDF is transformed into a mixture of a log-normaland a power-law part. However, the power-law part dominates thePDF at most times, which is attributed to the rather coarse nu-merical resolution, which does not allow the cloud to fragmentany further. to fluctuate slightly more, but the very good agreement ofthe averages is convincing. Our explanation for this goodcorrespondence is that we derive the driving parameter viathe compressive ratio χ = (cid:10) v (cid:11) / (cid:10) v (cid:11) , where resolutioneffects nearly cancel out. APPENDIX C: NOTES ON THE FOURIERDECOMPOSITION
The extracted clouds in our simulations resemble regionswith non-periodic boundary conditions. This can lead tosignificant errors in the resulting spectra due to aliasing ef-fects. In such case, one commonly applies a window function,which enforces the field of interest to go smoothly to zero sothat the box can be thought of as being periodic. The resultof this windowing is shown in Fig. C1. The modificationsto the final driving parameter are only minor. This is dueto the fact that it is obtained from the compressive ratio.Thus, possible errors will cancel out as they appear in bothof the decomposed fields.
DATA AVAILABILITY
The data underlying this article will be shared on reasonablerequest to the corresponding author.
REFERENCES
Ballesteros-Paredes J., Hartmann L. W., Vázquez-Semadeni E., Heitsch F., Zamora-Avilés M. A., 2011, MN-RAS, 411, 65Berger M. J., Oliger J., 1984, Journal of ComputationalPhysics, 53, 484Blitz L., Fukui Y., Kawamura A., Leroy A., Mizuno N.,Rosolowsky E., 2007, Protostars and Planets V, pp 81–96Bouchut F., Klingenberg C., Waagan K., 2009, NumerischeMathematikBurkhart B., Stalpes K., Collins D. C., 2017, ApJ, 834, L1Dedner A., Kemm F., Kröner D., Munz C. D., Schnitzer T.,Wesenberg M., 2002, Journal of Computational Physics,175, 645Dobbs C. L., Krumholz M. R., Ballesteros-Paredes J., Bo-latto A. D., Fukui Y., Heyer M., Low M.-M. M., Ostriker c (cid:13)000
Ballesteros-Paredes J., Hartmann L. W., Vázquez-Semadeni E., Heitsch F., Zamora-Avilés M. A., 2011, MN-RAS, 411, 65Berger M. J., Oliger J., 1984, Journal of ComputationalPhysics, 53, 484Blitz L., Fukui Y., Kawamura A., Leroy A., Mizuno N.,Rosolowsky E., 2007, Protostars and Planets V, pp 81–96Bouchut F., Klingenberg C., Waagan K., 2009, NumerischeMathematikBurkhart B., Stalpes K., Collins D. C., 2017, ApJ, 834, L1Dedner A., Kemm F., Kröner D., Munz C. D., Schnitzer T.,Wesenberg M., 2002, Journal of Computational Physics,175, 645Dobbs C. L., Krumholz M. R., Ballesteros-Paredes J., Bo-latto A. D., Fukui Y., Heyer M., Low M.-M. M., Ostriker c (cid:13)000 , 000–000 B. Körtgen t-t [Myr] b d r i v dx=9.70 pcdx=19.5 pc Figure B1.
Resolution study of the evolution of the turbulence driving parameter for a period of 180 Myr starting at the identificationof clouds ( t ). The error bars provide the standard deviation ∆ b driv . For the two resolutions, the driving parameter obtained from theHelmholtz decomposition appears to be largely converged. The grey horizontal lines denote purely solenoidal as well as natural driving. b d r i v Time [Myr]
Hanning WindowFiducial
Figure C1.
Comparison of the resulting driving parameter forthe early phase of evolution of the example cloud. The lines showthe fiducial method (red) and the result after applying a Han-ning window to the turbulent velocity field. The agreement of thetwo methods ensures that no numerical errors emerge from e.g.aliasing effects.
E. C., Vázquez-Semadeni E., 2014, Protostars and PlanetsVI, pp 3–26Dubey A., Fisher R., Graziani C., Jordan IV G. C., LambD. Q., Reid L. B., Rich P., Sheeler D., Townsley D.,Weide K., 2008, in Pogorelov N. V., Audit E., Zank G. P.,eds, Numerical Modeling of Space Plasma Flows Vol. 385of Astronomical Society of the Pacific Conference Series,Challenges of Extreme Computing using the FLASH code. pp 145–+Elmegreen B. G., Scalo J., 2004, ARA&A, 42, 211Federrath C., Klessen R. S., 2012, ApJ, 761, 156Federrath C., Klessen R. S., 2013, ApJ, 763, 51Federrath C., Klessen R. S., Schmidt W., 2008, ApJ, 688,L79Federrath C., Rathborne J. M., Longmore S. N., KruijssenJ. M. D., Bally J., Contreras Y., Crocker R. M., GarayG., Jackson J. M., Testi L., Walsh A. J., 2016, ApJ, 832,143Federrath C., Rathborne J. M., Longmore S. N., Krui-jssen J. M. D., Bally J., Contreras Y., Crocker R. M.,Garay G., Jackson J. M., Testi L., Walsh A. J., 2017,in Crocker R. M., Longmore S. N., Bicknell G. V., eds,IAU Symposium Vol. 322 of IAU Symposium, The linkbetween solenoidal turbulence and slow star formation inG0.253+0.016. pp 123–128Federrath C., Roman-Duval J., Klessen R. S., Schmidt W.,Mac Low M.-M., 2010, A&A, 512, A81Federrath C., Sur S., Schleicher D. R. G., Banerjee R.,Klessen R. S., 2011, ApJ, 731, 62Ginsburg A., Federrath C., Darling J., 2013, ApJ, 779, 50Hennebelle P., Chabrier G., 2011, ApJ, 743, L29Jin K., Salim D. M., Federrath C., Tasker E. J., Habe A.,Kainulainen J. T., 2017, MNRAS, 469, 383Jog C. J., 2013, MNRAS, 434, L56Kainulainen J., Beuther H., Henning T., Plume R., 2009,A&A, 508, L35Kainulainen J., Federrath C., 2017, A&A, 608, L3Konstandin L., Schmidt W., Girichidis P., Peters T., ShettyR., Klessen R. S., 2016, MNRAS, 460, 4483Körtgen B., Banerjee R., Pudritz R. E., Schmidt W., 2018,MNRAS, 479, L40Körtgen B., Banerjee R., Pudritz R. E., Schmidt W., 2019,MNRAS, 489, 5004Körtgen B., Federrath C., Banerjee R., 2017, MNRAS, 472, c (cid:13) , 000–000 urbulence in MCs c (cid:13)000