Simulations of the Milky Way's central molecular zone -- II. Star formation
Mattia C. Sormani, Robin G. Tress, Simon C.O. Glover, Ralf S. Klessen, Cara D. Battersby, Paul C. Clark, H Perry Hatchfield, Rowan J. Smith
MMNRAS , 000–000 (0000) Preprint 3 July 2020 Compiled using MNRAS L A TEX style file v3.0
Simulations of the Milky Way’s central molecular zone - II. Starformation
Mattia C. Sormani , Robin G. Tress , Simon C.O. Glover , Ralf S. Klessen , ,Cara D. Battersby , Paul C. Clark , H Perry Hatchfield and Rowan J. Smith Universität Heidelberg, Zentrum für Astronomie, Institut für theoretische Astrophysik, Albert-Ueberle-Str. 2, 69120 Heidelberg, Germany Universität Heidelberg, Interdiszipliäres Zentrum für Wissenschaftliches Rechnen, Im Neuenheimer Feld 205, 69120 Heidelberg, Germany University of Connecticut, Department of Physics, 196 Auditorium Road, Unit 3046, Storrs, CT 06269, USA School of Physics and Astronomy, Queen’s Buildings, The Parade, Cardiff University, Cardiff, CF24 3AA, UK Jodrell Bank Centre for Astrophysics, School of Physics and Astronomy, University of Manchester, Oxford Road, Manchester M13 9PL, UK
ABSTRACT
The Milky Way’s central molecular zone (CMZ) has emerged in recent years as a uniquelaboratory for the study of star formation. Here we use the simulations presented in Tress et al.2020 to investigate star formation in the CMZ. These simulations resolve the structure of theinterstellar medium at sub-parsec resolution while also including the large-scale flow in whichthe CMZ is embedded. Our main findings are as follows. (1) While most of the star formationhappens in the CMZ ring at R (cid:38)
100 pc, a significant amount also occurs closer to SgrA* at R (cid:46)
10 pc. (2) Most of the star formation in the CMZ happens downstream of the apocentres,consistent with the “pearls-on-a-string” scenario, and in contrast to the notion that an absoluteevolutionary timeline of star formation is triggered by pericentre passage. (3) Within thetimescale of our simulations ( ∼
100 Myr), the depletion time of the CMZ is constant within afactor of ∼
2. This suggests that variations in the star formation rate are primarily driven byvariations in the mass of the CMZ, caused for example by AGN feedback or externally-inducedchanges in the bar-driven inflow rate, and not by variations in the depletion time. (4) We studythe trajectories of newly born stars in our simulations. We find several examples that have ageand 3D velocity compatible with those of the Arches and Quintuplet clusters. Our simulationssuggest that these prominent clusters originated near the collision sites where the bar-driveninflow accretes onto the CMZ, at symmetrical locations with respect to the Galactic centre,and that they have already decoupled from the gas in which they were born.
Key words:
Galaxy: centre - Galaxy: kinematics and dynamics - ISM: kinematics anddynamics - ISM: clouds - ISM: evolution - stars: formation
The central molecular zone (CMZ, R (cid:46)
200 pc) is the Milky Way’scounterpart of the star-forming nuclear rings that are commonlyfound in the central regions of external barred galaxies such asNGC 1300 (see for example the atlas of nuclear rings of Comerónet al. 2010). Being a hundred times closer than the nucleus of thenext comparable galaxy, Andromeda, it offers us the possibility tostudy a nuclear ring in unique detail.The CMZ has emerged in the last decade as a unique laboratoryfor the study of star formation (e.g. Molinari et al. 2011; Kruijssenet al. 2014; Armillotta et al. 2019). The main reason is that theenvironmental conditions in which stars are born are more extremethan anywhere else in the Galaxy. Indeed, the physical propertiesof the interstellar medium (ISM) in the CMZ are substantially dif-ferent from those in the Galactic disc: average gas volume densities(Guesten & Henkel 1983; Walmsley et al. 1986; Longmore et al. 2017; Mills et al. 2018), temperatures (Immer et al. 2016; Ginsburget al. 2016; Krieger et al. 2017; Oka et al. 2019), velocity disper-sions (Shetty et al. 2012; Federrath et al. 2016), and magnetic fieldstrengths (Morris 2015; Mangilli et al. 2019) are all much higherthan in the disc. The interstellar radiation field and higher cosmicray ionisation rate (Clark et al. 2013; Ginsburg et al. 2016; Oka et al.2019) are also much stronger. In addition, the CMZ region is char-acterised by the presence of Galactic outflows (Ponti et al. 2019), bythe widespread presence of radio-emitting magnetised non-thermalfilaments (Heywood et al. 2019), and by a strong hydrodynamicalinteraction with the larger-scale gas inflow driven by the Galacticbar (Sormani et al. 2018a). The star formation process, which isdetermined by the complex interplay of all these physical agents, istherefore expected to proceed differently in the CMZ. Observationsconfirm this, by showing that the CMZ does not obey some star for-mation relations that are valid in the disc (Longmore et al. 2013a; © 0000 The Authors a r X i v : . [ a s t r o - ph . GA ] J u l Sormani et al.
Kauffmann et al. 2017a). Hence, understanding star formation inthe CMZ is important for understanding the star formation processin extreme environments, as well as in general by probing a peculiarcorner of parameter space.In a companion paper (Tress et al. 2020, hereafter Paper I) wehave presented sub-parsec resolution hydrodynamical simulationsand have used them to study the gas dynamics in the CMZ. In thispaper, we use the same simulations to investigate star formation.Open questions that we address in the current work include:(i) What is the temporal distribution of star formation in theCMZ? (Section 3.1)(ii) What is spatial distribution of star formation in the CMZ?(Section 3.2)(iii) What is the impact of the orbital dynamics on star formation?Can we identify an absolute evolutionary timeline of star formationas suggested by Longmore et al. (2013b) and Kruijssen et al. (2015)?(Section 4.2)(iv) What drives the time variability of star formation in theCMZ? (Section 4.1)(v) Are the Arches and Quintuplet cluster on a common orbitwith gas in the CMZ ring (Kruijssen et al. 2015) or are they on othertypes of orbits (Stolte et al. 2008)? (Section 4.4)The paper is structured as follows. In Section 2 we give a briefsummary of our numerical simulations. In Section 3 we study thetemporal and spatial distribution of star formation and the trajecto-ries of newly born stars. In Section 4 we discuss the implications ofour results for some of the open questions raised above. We sum upin Section 5.
Our simulations have been presented in detail in Paper I. Hence wegive here only a very brief overview, and refer to that paper for moredetails.
The simulations are similar to those we previously discussed inSormani et al. (2018a, 2019), with the following differences: (i)inclusion of gas self-gravity; and (ii) inclusion of a sub-grid pre-scription for star formation and stellar feedback. In particular, weemploy exactly the same externally-imposed rotating barred poten-tial, the same chemical/thermal treatment of the gas, and the sameinitial conditions as in Sormani et al. (2019).We use the moving-mesh code arepo (Springel 2010; Wein-berger et al. 2019). The simulations are three-dimensional and un-magnetised, and include a live chemical network that keeps trackof hydrogen and carbon chemistry. The simulations comprise inter-stellar gas in the whole inner disc ( R ≤ Φ ext ( x , t ) which is constructed to fit the propertiesof the Milky Way. The bar component rotates with a pattern speed Ω p =
40 km s − kpc − , consistent with the most recent determina-tions (e.g. Sormani et al. 2015; Portail et al. 2017; Sanders et al.2019), which places the (only) inner Lindblad resonance (ILR) cal-culated in the epicyclic approximation at R ILR = . R CR = . ρ c = − g cm − is removed from the simulation and replacedwith a non-gaseous sink particle, provided that it is unambiguouslygravitationally bound and not within the accretion radius of an exist-ing sink particle. The sink particle does not represent an individualstar, but rather a small cluster which contains both gas and stars.(ii) Once a sink is created, a stellar population is assigned to itby drawing from an initial mass function (IMF) according to thePoisson stochastic method described in Sormani et al. (2017).(iii) Sink particles are allowed to accrete mass at later times,provided that the gas is within the sink accretion radius r acc = M ≥ (cid:12) ) assigned to the sink, weproduce a supernovae (SNe) event with a time delay which dependson the stellar mass. Each SNe event injects energy and/or momentuminto the ISM and gives back to the environment part of the gas“locked-up” in the sink. Energy is injected only if the local resolutionof the Voronoi mesh is high enough to resolve the supernova remnantat the end of its Sedov-Taylor phase; otherwise, an appropriateamount of momentum is injected instead. SNe feedback is the onlytype of feedback included in the simulation.(v) When all the SNe associated with a sink have exploded and allof its gas content has been given back to the environment, the sink isconverted into a collisionless N-body star particle with a mass equalto the stellar mass of the sink. This N-body particle continues to existindefinitely in the simulation and affects it through its gravitationalpotential, but, unlike a sink, it can no longer accrete new gas or formnew stars.When making projections onto the plane of the Sky, we assumean angle between the Sun-Galactic centre line and the bar major axisof φ = ◦ , a Sun-Galactic centre distance of 8 . v (cid:12) =
235 km s − (Schönrich et al. 2010; Reid et al.2019), as in Paper I. As in Paper I, we subdivide our simulation into three spatial regionsin order to facilitate the subsequent analysis (see Figure 1): • The
CMZ is defined as the region within cylindrical radius R ≤
250 pc. • The dust lane region (DLR) is the elongated transition regionbetween the CMZ and the Galactic disc, where highly non-circulargas motions caused by the bar are present. • The disc is defined as everything outside the DLR.
Figure 2 shows the star formation rate (SFR) as a function of timein our simulation, calculated as a running average over the last
MNRAS , 000–000 (0000) tar formation in the CMZ − − − x [ kpc ] − − y [ kp c ] t = . CMZDLR Disc
Figure 1.
Definition of the three regions (CMZ, DLR, disc) into which wesubdivide our simulated Galaxy for subsequent analysis. See Section 2.2 formore details. . The thinblue line shows the total SFR in the entire simulation box (CMZ +DLR + disc). This is roughly constant at a value of approximately ∼ (cid:12) yr − , consistent with typically reported values of the MWtotal SFR derived from observations ( ∼ (cid:12) yr − , e.g. Kennicutt& Evans 2012) when we take into account that our simulated disconly extends to R (cid:39) (cid:39) . × M (cid:12) ) is only ∼ / R ≤
250 pc, see Figure 1). The insert pan-els correlate the SFR with the CMZ gas morphology at differenttimes. At t =
146 Myr (when the bar potential is fully turned on,see Section 2.7 in Paper I), the SFR in the CMZ has a value of ∼ . (cid:12) yr − , consistent with observational estimates (see Yusef-Zadeh et al. 2009; Immer et al. 2012; Longmore et al. 2013a; Barneset al. 2017 and Section 4.1), and the total gas mass of the CMZ is ∼ × M (cid:12) , which also agrees well with observational values( ∼ × M (cid:12) , Dahmen et al. 1998; Longmore et al. 2013a).At later times ( t >
146 Myr) the SFR of the CMZ slowlybut steadily increases with time, with small fluctuations on shorttimescales ( ∼ ∼ τ depl = M / SFR), is shownby the blue dashed line in Figure 3. It is approximately constant intime. Therefore, the SFR in the CMZ in our simulation is roughlyproportional to its total mass, and variations in the value of the SFRare determined by variations in the total mass.Figure 3 also shows that the depletion time in the disc (yellowdashed line) is a factor of ∼ Observationally determined rates are more often averaged over longertimescales ( ∼
10 Myr). We will briefly discuss the distribution of older starsin Section 3.2, while we defer a more observationally oriented approach andsynthetic observations to future work. when considering different portions of the Galaxy. The variationsin the depletion times can be explained by the different stellar grav-itational potential, whose vertical gradient is stronger in the CMZthan in the disc. This can be seen as follows. For a medium in whichthe turbulence is driven by supernovae feedback and assuming thatthe vertical force of the gravitational potential is balanced by theturbulent pressure (both conditions that are approximately verifiedin our simulations) the analytical model of Ostriker & Shetty (2011)predicts that (see their Equation 13): Σ SFR ∝ ( + χ ) Σ , (1)where Σ SFR is the SFR surface density, Σ is the total gas surfacedensity, χ = C /( + √ + C ) , C = ζ d ρ b σ z /( π G Σ ) , σ z isthe vertical velocity dispersion, ρ b is the stellar midplane density(which is proportional to the strength of the gravitational potential), G is the gravitational constant and ζ d (cid:39) / C (cid:29) χ (cid:39) ( C ) / and therefore Σ SFR ∝ ρ / b Σ . The depletion time is then τ depl = Σ / Σ SFR ∝ ρ − / b . For the potential employed in our simulations, wefind [ ρ b ( R =
150 pc )/ ρ b ( R = )] − / (cid:39)
6, in good agreementwith the results in Figure 3 considering the uncertainties presentboth in the simulations and in the simplifying assumptions on whichthe theory of Ostriker & Shetty (2011) is based. The agreementbetween our simulation and the theory of Ostriker & Shetty (2011)is consistent with the fact that the integrated properties of the CMZfollow well star formation relations based on the total or moleculargas surface density, such as the Schmidt-Kennicutt or the Bigielet al. (2008) relation, and only become peculiar when consideringthe very dense gas (see Section 4.3).Another factor that is likely to contribute to lowering the deple-tion time in the CMZ, and which is not accounted for in the verticalequilibrium theories of Ostriker et al. (2010) and Ostriker & Shetty(2011), is the increased number of shocks due to the large-scale barflow, which cause local compressions and therefore enhanced starformation (Mac Low & Klessen 2004; Klessen & Glover 2016).How would the CMZ mass/SFR evolve if we continue oursimulation beyond the maximum time shown in Figure 2? Assumingthat the depletion time remains constant at the value τ depl , CMZ (cid:39) × yr inferred from Figure 3, we might extrapolate that themass of the CMZ would keep increasing until the SFR matchesthe bar-driven inflow rate. For an inflow rate of (cid:219) M (cid:39) (cid:12) yr − (see Paper I), the equilibrium CMZ mass would be (cid:219) M τ depl , CMZ (cid:39) × M (cid:12) . However, there are several factors that might invalidatethis extrapolation: (i) at a mass (cid:39) × M (cid:12) , the gravitationalpotential of the gas would become comparable to that of the stars,which would affect the depletion time (see discussion immediatelyafter Equation 1 above); (ii) at a SFR of (cid:39) (cid:12) yr − , the increasedSN feedback rate might also change the depletion time; (iii) the bar-driven inflow rate will decrease once the reservoir at R (cid:38) MNRAS000
6, in good agreementwith the results in Figure 3 considering the uncertainties presentboth in the simulations and in the simplifying assumptions on whichthe theory of Ostriker & Shetty (2011) is based. The agreementbetween our simulation and the theory of Ostriker & Shetty (2011)is consistent with the fact that the integrated properties of the CMZfollow well star formation relations based on the total or moleculargas surface density, such as the Schmidt-Kennicutt or the Bigielet al. (2008) relation, and only become peculiar when consideringthe very dense gas (see Section 4.3).Another factor that is likely to contribute to lowering the deple-tion time in the CMZ, and which is not accounted for in the verticalequilibrium theories of Ostriker et al. (2010) and Ostriker & Shetty(2011), is the increased number of shocks due to the large-scale barflow, which cause local compressions and therefore enhanced starformation (Mac Low & Klessen 2004; Klessen & Glover 2016).How would the CMZ mass/SFR evolve if we continue oursimulation beyond the maximum time shown in Figure 2? Assumingthat the depletion time remains constant at the value τ depl , CMZ (cid:39) × yr inferred from Figure 3, we might extrapolate that themass of the CMZ would keep increasing until the SFR matchesthe bar-driven inflow rate. For an inflow rate of (cid:219) M (cid:39) (cid:12) yr − (see Paper I), the equilibrium CMZ mass would be (cid:219) M τ depl , CMZ (cid:39) × M (cid:12) . However, there are several factors that might invalidatethis extrapolation: (i) at a mass (cid:39) × M (cid:12) , the gravitationalpotential of the gas would become comparable to that of the stars,which would affect the depletion time (see discussion immediatelyafter Equation 1 above); (ii) at a SFR of (cid:39) (cid:12) yr − , the increasedSN feedback rate might also change the depletion time; (iii) the bar-driven inflow rate will decrease once the reservoir at R (cid:38) MNRAS000 , 000–000 (0000)
Sormani et al.
Figures 4 and 5 show the spatial distribution of the SFR density in atypical simulation snapshot. As before, the SFR is calculated as therunning average over the last 0 . material crashes into the dust lanes,which in Sormani et al. (2019) we have interpreted as producing theobserved extended velocity features (EVF), are sites of enhancedstar formation. However, by the time this star formation is visible,these regions will have moved at high speed ( ∼
200 km s − ) inwardstowards the CMZ. The time delay between sink formation in ourmodel and the star formation actually becoming visible will dependon our choice of star formation rate tracer, but we would expect itto be at least ∼ . α that are highlysensitive to dust obscuration. This is consistent with observationsof Bania Clump 2 (one of the most prominent EVF), which despitecontaining dozen of 1.1 mm clumps, has been found to be deficientin near- and mid-infrared emission in the Spitzer images and hasbeen suggested to be in a pre-stellar stage of cloud evolution byBally et al. (2010). Our simulations therefore support the idea thatBania Clump 2 will shortly begin to form massive stars.A noteworthy feature of the averaged as well as of the instan-taneous SFR density distributions (bottom panel of Figure 6 andFigure 5) is that there is a site of star formation inside the CMZring radius, after a radial gap. Indeed, we noted in Section 3.4 ofPaper I that gas can be found inside the CMZ radius in these simu-lations (in contrast to our previous non self-gravitating simulationsin Sormani et al. 2019, in which there was no gas inside the CMZring). This star formation might be associated with star formationoccurring near SgrA* ( R ≤
10 pc). This would be consistent withclaims of observational evidence for ongoing star formation in thisregion (Yusef-Zadeh et al. 2008, 2015), although we note that theseclaims are controversial at the moment (Mills et al. 2017). Suchstar formation might also be related to the formation of the nuclearstellar cluster (NSC, see for example Genzel et al. 2010; Schödelet al. 2014; Gallego-Cano et al. 2020) by providing in-situ newlyborn stars and, since such stars are rotating, it might contribute toits observed rotation (Feldmeier et al. 2014; Feldmeier-Krause et al.2015; Chatzopoulos et al. 2015; Tsatsi et al. 2017; Neumayer et al.2020). We use the term “overshooting” to denote material that, after plungingtowards the CMZ along one of the dust lanes, passes close to the CMZ butdoes not stop and continues towards the dust lane on the opposite side. Seefor example Figure 4 in Sormani et al. (2019).
Figure 7 analyses the radial distribution of Σ gas , Σ SFR and ofthe depletion time. The lines show the time-averaged values, whilethe shaded regions show the scatter. This figure indicates that both Σ gas and Σ SFR increase considerably in the centre, while the ratiobetween the two, the depletion time, decreases by a factor of ∼ R (cid:39)
500 pc,in the terminal part of the dust lanes. This is where gas reaches thehighest bulk speeds (and observed line of sight velocities) over theentire MW disc, and may indicate that star formation is suppressedat these sites due to the very high shear, in line with the argumentspresented in Renaud et al. (2015) and Emsellem et al. (2015). Inorder to check this, we plot in Figure 9 the quantity τ = (cid:34)(cid:18) ∂ V x ∂ y + ∂ V y ∂ x (cid:19) + (cid:18) ∂ V x ∂ x − ∂ V y ∂ y (cid:19) (cid:35) / , (2)where V i = ∫ ∞∞ ρ v i d z / ∫ ∞∞ ρ d z is the density-weighted projectedvelocity. The quantity τ is a good indication of shear for a 2D flow,and has the desirable property of being invariant under rotationsof the coordinates since it is the magnitude of the eigenvalues ofthe traceless shear tensor D ij = (cid:2) ∂ j V i + ∂ i V j − δ ij (∇ · V ) (cid:3) / ∂ i V j using finitedifferences with a resolution ∆ x = ∼ l (cid:39) . ◦ and l = − ◦ (lower panel),roughly consistent with observations which have peaks at the posi-tion of SgrB2 and SgrC (see for example Figure A1 of Barnes et al.2017). Again, fluctuations of the instantaneous distribution aroundthe averaged distribution can be quite large, and the peaks can be MNRAS , 000–000 (0000) tar formation in the CMZ more or less evident in the instantaneous distributions dependingon the particular snapshot chosen. Once a sink particle is formed, it typically follows a different tra-jectory than the gas. As already noted in Section 3.2, this can beseen for example from Figure 8, which show how gas and stars inthe CMZ quickly decouple and have achieved significantly differentdistribution within ∼ − ), and after eachpassage reach several kpc out from the centre. The second panelshows stars that formed downstream along the dust lanes, where thegas is accreting onto the CMZ. These stars will overshoot a little bitand typically have elongated orbits which are a factor of 2-3 largerthan the CMZ ring. Typical orbital speeds of these stars are larger( ∼
150 km s − ) than gas in the CMZ ring ( ∼ − ). Thethird panel shows stars formed within the CMZ ring. These starswill stay within the ring and have typical orbital velocities compa-rable to the gas in the ring ( ∼ − ), but after a few Myrthey will decouple from the gas. The accumulation of stars similarto those shown in the second and third panel is what forms the nu-clear stellar disc over time (NSD, see for example Launhardt et al.2002; Nishiyama et al. 2013; Schönrich et al. 2015; Baba & Kawata2020). Finally, the last panel shows stars that have formed from gasinside the CMZ ring. These typically follow roughly circular orbitswith moderate speeds ( ∼
80 km s − ), so they will remain inside theCMZ ring. As noted in Section 3.2, such star formation might alsobe related to the formation of the nuclear stellar cluster (NSC, seefor example Genzel et al. 2010; Schödel et al. 2014; Gallego-Canoet al. 2020).The trajectories of the sink particles in our simulation canbe compared with the kinematics of star clusters and HII regions.In Section 4.4 we compare them with the Arches and Quintupletclusters. In a an upcoming paper (Anderson et al., in preparation)we will compare them with HII regions in the SgrE complex.Finally, it is worth mentioning a limitation of our simulations.In the code, the gravitational force is calculated using a softeninglength, which for the gas is adaptive and depends on the cell sizewith a lower limit set at 0 . The current global SFR in the CMZ (intended here as the regionwithin R (cid:46)
200 pc, or | l | (cid:46) . ◦ assuming a distance to the Galacticcentre of 8 . (cid:39) . (cid:12) yr − (e.g. Yusef-Zadeh et al. 2009; Immer et al. 2012; Longmore et al. 2013a; Barnes et al. 2017). Thisnumber is obtained by combining different independent methods,including direct counting of young stellar objects and integratedlight measurements. All these methods agree with each other withina factor of two (see Table 3 in Barnes et al. 2017), and also agreewith the number obtained from counts of supernovae remnants(see Section 8.9 in Ponti et al. 2015). Since these various methodstrace SF over different timescales in the range 0 . ∼ . . (cid:12) yr − , i.e. a factor of a few higher thanthe rate averaged over the last 5 Myr. They also found that the SFRhas been variable during the past Gyr, with periods of more intenseactivity ( ∼ . (cid:12) yr − ). This suggests that the SFR in the CMZ isnot constant, but varies in time. Evidence for time variability in thestar formation activity has also been found by Sarzi et al. (2007) forexternal galactic nuclei by analysing the star formation history of asample of nuclear rings. It is therefore natural to ask: what drivesthe time variability in the SFR of the CMZ?A possible explanation is that the CMZ goes through episodicstarbursts driven by feedback instabilities (Krumholz & Kruijssen2015; Krumholz et al. 2017; Torrey et al. 2017; Armillotta et al.2019). In this scenario, the CMZ has a roughly constant gas mass butorder-of-magnitude level variations in the SFR. The depletion timeis not constant, but has large variations over time. The large scatter( ∼ . . (cid:12) yr − ∼
30 Myr ago (Nogueras-Lara et al. 2020) to the value ∼ . (cid:12) yr − inferred for the last 5 Myr (Barnes et al. 2017). This would becompatible with the currently estimated ages of the Fermi Bubbles(see for example Mou et al. 2018 and references therein). The massof the CMZ could also change due to variations in the accretionrate, induced for example by an external perturbation such as amerger. We note that at the current estimated mass inflow rateof ∼ (cid:12) yr − (Sormani & Barnes 2019), the entire current gasmass of the CMZ ( (cid:39) × M (cid:12) ) can be accumulated in just50 Myr, so a change in this rate could potentially induce mass andSFR variability within the timescales required by observations. We MNRAS000
30 Myr ago (Nogueras-Lara et al. 2020) to the value ∼ . (cid:12) yr − inferred for the last 5 Myr (Barnes et al. 2017). This would becompatible with the currently estimated ages of the Fermi Bubbles(see for example Mou et al. 2018 and references therein). The massof the CMZ could also change due to variations in the accretionrate, induced for example by an external perturbation such as amerger. We note that at the current estimated mass inflow rateof ∼ (cid:12) yr − (Sormani & Barnes 2019), the entire current gasmass of the CMZ ( (cid:39) × M (cid:12) ) can be accumulated in just50 Myr, so a change in this rate could potentially induce mass andSFR variability within the timescales required by observations. We MNRAS000 , 000–000 (0000)
Sormani et al.
80 100 120 140 160 180 200 220 240 t [ Myr ] − − SF R [ M (cid:12) y r − ] SFR(CMZ)SFR(DLR)SFR(Disc)SFR totalObserved CMZ t = . t = . t = . t = . t = . t = . N [ cm − ] N [ cm − ] Figure 2.
Star formation rate as a function of time in our simulation. The thick blue, thin pink and thin yellow lines are the SFR in the three different spatialregions (CMZ, DLR, disc) in which we have subdivided our simulation (see Figure 1). The thin black line is the total SFR (CMZ+DLR+disc). The insert panelsshow total gas surface density maps that allow us to correlate the SFR with the instantaneous CMZ morphology. The blue shaded horizontal region indicatesthe observed current SFR of the CMZ, taken to be in the conservative range 0 . . (cid:12) yr − (see references in Section 4.1) The grey shaded area indicatesthe times when the bar potential is still gradually turning on (see Section 2.7 in Paper I), which are excluded from the analysis. also note that much higher accretion rates seem to be possible inbarred galaxies: for example Elmegreen et al. (2009) reports a bar-driven inflow rate of 40 M (cid:12) yr − in NGC 1365. Our scenario is alsosupported by the work of Seo et al. (2019), who run hydrodynamicalsimulations of gas flowing in a live N-body barred potential and findthat the SFR correlates well with the bar-inflow rate. In our scenario,the large scatter in the depletion times observed in the centre ofexternal barred galaxies (Leroy et al. 2013; Utomo et al. 2017)would be explained as due to different environmental conditionsrather than to high time variability. For example, different strengthsof the stellar gravitational potential might contribute to the scatter in the depletion times (see Equation 1 and related discussion). Notealso that some of the scatter in these values may be driven bydifferences in the size of the CMZ-like region in different galaxies,since this region is typically not resolved in the kpc-scale moleculargas maps considered in Leroy et al. (2013) and Utomo et al. (2017).What causes the differences between the results presented hereand those in Armillotta et al. (2019)? There are several factors thatcould contribute to this and it is difficult to point to which oneis most important. First, the two papers use significantly differenttreatments of ISM cooling. Armillotta et al. (2019) treat gas coolingusing equilibrium cooling curves provided by the grackle astro- MNRAS , 000–000 (0000) tar formation in the CMZ TOTAL
CMZDLRDiscObserved CMZwithout gas in sinkswith gas in sinks τ d e p [ y r ] H
100 125 150 175 200 225 250 t [ Myr ] HI Figure 3.
Depletion time ( τ dep ) as a function of simulation time ( t ) for thevarious regions defined in Figure 1. The blue shaded region indicates theobserved depletion time of the CMZ, 0 . − × M (cid:12) , see refer-ences in Section 3.1) by the observed SFR of the CMZ shown in Figure 2(0 . . (cid:12) yr − ). The grey shaded area indicates where the bar potentialis gradually turning on, which are excluded from the analysis. chemistry and cooling package (Smith et al. 2017), which potentiallyyield differences in behaviour compared to the fully non-equilibriumtreatment we use here. In addition, they treat photoelectric heatingas a uniform heating process and do not account for variations in theheating rate due to changes in the fractional ionisation of the gas orits degree of dust shielding. Although the two treatments result inISM phase diagrams that are qualitatively similar in many aspects(compare their Figure 5 with Figure 11 in Paper I), there are clearquantitative differences that may have some impact on the predictedstar formation rates.Second, the star formation prescription used in Armillotta et al. (2019) is also quite different from that used in our code. In theirapproach star particles are stochastically formed in gas denser than10 cm − , provided that it is gravitationally bound, cold and self-shielded. Compared to our scheme, the main differences are theirchoice of density threshold and the fact that in their scheme, sig-nificant quantities of dense gas can accumulate above the densitythreshold, something which is impossible by design in our scheme.Third, Armillotta et al. (2019) include the effects of photoionisationfeedback as well as supernova feedback, while we concentrate heresolely on the latter.Finally, there is a substantial difference in the mass resolutionachieved in dense gas in the two simulations. In our simulation,gas at densities around 10 cm − is typically resolved with Voronoicells with a mass of around 2 M (cid:12) (see Figure 3 in Paper I). Incontrast, the default particle mass in Armillotta et al. (2019) is2000 M (cid:12) , a factor of 1000 worse than we achieve here. Armillottaet al. (2019, 2020) also present results from a “high resolution”run with a particle mass of 200 M (cid:12) , which they carried out for amuch shorter period than their main run, but even this has a muchworse resolution than our simulation. An important consequence ofthis difference in resolution is that in the Armillotta et al. (2019)simulation, the Sedov-Taylor phase following a supernova explosionis resolved only for supernovae exploding in low density gas with n < − , whereas in our simulation it remains well-resolvedeven for supernovae exploding in gas with a density close to oursink creation threshold. Therefore, Armillotta et al. (2019) primar-ily inject momentum with their supernovae, since the associatedthermal energy is rapidly radiated away, whereas we are able tofollow the injection of both thermal energy and momentum in amore self-consistent fashion. This results in a clear difference inthe morphology of the supernova-affected gas: in our simulation,supernova explosions produce large holes in the gas distribution,while corresponding features are rarely seen in the Armillotta et al.(2019) simulation.In view of these significant differences in numerical approach,together with the fact that the results Armillotta et al. (2019) obtainfor the star formation rate of the CMZ are clearly not numericallyconverged (see their Figure A2), and that the simulations span aquite different period in the life of the CMZ ( ∼
100 Myr after barformation in our run vs. 500 Myr in their simulation), it is difficultto assess the reasons for the difference in results regarding the timevariability of the SFR in the CMZ. This is an issue that we hope toaddress further in future work.
Longmore et al. (2013b) and Kruijssen et al. (2015, 2019) sug-gested that star formation follows an evolutionary timeline as thegas clouds orbit the CMZ ring. In this scenario, star formation istriggered when the clouds are compressed during pericentre pas-sage, i.e. when the clouds pass closest to the Galactic centre. Thisscenario is at variance with the two scenarios for star formationin nuclear rings that are more commonly discussed in the extra-galactic context, namely the “popcorn” and the “pearls on a string”scenarios, which are schematically depicted in Figure 7 of of Bökeret al. (2008). In the “pearls on a string” scenario, star formation oc-curs prevalently at the contact point between the dust lanes and thegas ring, which typically coincides with the ring apocentre ratherthan with the pericentre. In the “popcorn” scenario, star formationoccurs uniformly along the ring. The observational evidence for aclear evolutionary sequence as implied by the pericentre passagescenario is mixed (Kauffmann et al. 2017a; Krieger et al. 2017),
MNRAS000
MNRAS000 , 000–000 (0000)
Sormani et al. − − − − − − − − − − − . . . . . . x [ kpc ] . . . . . . y [ kp c ] − − Σ SFR [ M (cid:12) yr − kpc − ] Figure 4.
SFR density for various snapshots in our simulation. Shown is the very recent (0 . surfacedensity. Compare with the time-averaged SFR density shown in the bottom panel of Figure 6. while the peals on a string scenario has obtained some mild supportfrom observations of nearby galaxies (see for example Section 4.1in Böker et al. 2008, see also Mazzuca et al. 2008).These three scenarios make different predictions that can betested with our simulations. The pericentre scenario predicts thatstar formation occurs predominantly after the pericentre passage.According to this scenario, very young stars should be found shortlyafter the passage, while stars of increasing age should be found fur-ther downstream of the pericentre. The pearls on a string scenariopredicts that star formation happens predominantly downstream ofthe contact point between the dust lanes and the CMZ ring, i.e.downstream of the apocentre. The popcorn scenario predicts thatstar formation is distributed uniformly along the ring, without pre-ferred locations.In order to test these predictions, we look at the time-averageddistribution of very young stars (age t ≤ .
25 Myr), which is shownin the top row of Figure 8. These stars trace where the star formation is being triggered. The right panel in the top row shows the azimuthaldistribution of stars in the CMZ. The apocentres of the CMZ ring areat θ = θ = ◦ , and coincide with the contact points betweendust lanes, while the pericentres are at θ = ◦ and θ = ◦ . Thispanel shows that the distribution of very young stars has a bi-periodicstructure with two strong peaks at the apocentres, consistent withthe prediction of the pearls on a string scenario.The above analysis considers all the star formation within R ≤
250 pc, including some that strictly speaking is outside the “ring”structure. In order to investigate this aspect in more detail, we focusspecifically on the ring in Figure 12. The right panel shows thetime-averaged surface density of very young stars ( Σ YS ) within By plotting the surface density rather than a histogram of the mass distri-bution as a function of azimuth, we avoid any potential bias due to geometriceffects caused by the area within the ellipse not being constant in each angu-lar range. For example, if the surface density were constant along the ring,MNRAS , 000–000 (0000) tar formation in the CMZ − . − . . . . − . − . . . . − . − . . . . − . . . − . − . . . . − . . . − . . . . . . . . . x [ kpc ] . . . . . . y [ kp c ] − Σ SFR [ M (cid:12) yr − kpc − ] Figure 5.
Same as Figure 4, but zooming onto the CMZ. Compare with the time-averaged SFR density shown in the bottom panel of Figure 6. the elliptical ring shown in the middle panel. It can be seen thatmost of the star formation occurs downstream of the apocentres, butbefore the pericentre passage. This is consistent with the predictionof the pearls on a string scenario, but not with the pericentre passagescenario. Note however that the maxima are quite broad, and starformation away from these maxima is not zero. Thus, while themaxima constitute a region of more intense star formation, they arenot the only regions where stellar birth takes place.Let us now consider the distribution of older stars, which canbe seen from Figure 8. For ages 1 < t < the azimuthal distribution of mass would not be constant, although the 2Dface-on maps would look perfectly uniform. ( t >
10 Myr), the bipolar structure precesses as a consequence ofthe precession of the apocentres of the stellar orbit (which at theirformation prevalently coincide with the contact point between ringand the dust lanes, but change at later times). The bipolar structurealso becomes less pronounced, and the distribution more uniform,as stars mix in phase spaceWe remark that all our conclusions above come from analysingthe time-averaged distributions. As discussed in Section 3.2, theinstantaneous star formation distribution fluctuates strongly aroundthe average (compare the left panels in Figure 8 and the various pan-els in Figure 5 with the middle-right panels of Figure 8). Because ofthese fluctuations, it is much harder to tell whether our simulationsare consistent with the pearls on a string scenario by looking justat a single snapshot. Moreover, while the time-averaged distribu-tions favour the pearls on a string scenario, there is also significantstar formation throughout the ring and away from the apocentres.These complications should be taken into account when analysingobservations, which only constitute individual snapshots.
MNRAS000
MNRAS000 , 000–000 (0000) Sormani et al. − − t = . t = . − − y [ kp c ] − − − x [ kpc ] − − − − Σ HI [ M (cid:12) pc − ] − − Σ H [ M (cid:12) pc − ] − − Σ SFR [ M (cid:12) yr − kpc − ] Figure 6.
Time-averaged plot of top: H surface density; middle : Hi surfacedensity; bottom : star formation rate (SFR) density. The average is calculatedover the range t = . . x orbits. From a physical point of view, there are two reasons whyenhanced star formation should be expected at the apocentres: (i)they are collisions sites where the gas from the dust lanes crashesinto the ring (see e.g. the left panel in Figure 12); (ii) gas slowsdown at the apocentre of an orbit, causing it to pile up and becomemore dense. Our simulations suggest that these effects are dominantover the tidal compression at the pericentre proposed by Kruijssenet al. (2015). Even neglecting these two dominant effects, thereis evidence that the pericentre passage only has a minor role intriggering star formation events. We note that in the simulations ofDale et al. (2019) the pericentre has a rather weak effect in enhancingthe SFR (compare the circular and non-circular orbits in Figures 3and 9 of Dale et al. 2019). Jeffreson et al. (2018) also estimates thatonly a small fraction ( ∼ Σ G a s [ M (cid:12) / p c ] Total gasH HIGas in sinks − − − − − Σ SF R [ M (cid:12) / y r / p c ] SFR . . . . . . . . . R [ kpc ] τ d e p [ y r ] Total gasH2HIH2 + gas in sinks . . . . . . . . . . . . − − − − − . . . . . . Figure 7.
Top : time-averaged radial distribution of gas surface density.
Mid-dle : SFR density.
Bottom : depletion times. Plots are averaged over time in therange t = ( . , . ) Myr. Shaded areas show the 1-sigma scatter. Thezoom-in inlays show the time-averaged quantities in the innermost 0 . dust lanes and gas slowing down at the apocentre, so it is likely thatthe actual number is significantly lower than this. Finally, Kruijssenet al. (2019) also acknowledge that star formation might be triggeredby accretion, similarly to the pearls-on-a-string scenario. However,in their discussion the accumulation of gas in the CMZ takes placewithin the context of the Krumholz & Kruijssen (2015) model ratherthan from direct accretion from the dust lanes. As we have arguedin Section 6.2 of Paper I, the theoretical framework of Krumholz& Kruijssen (2015) and Krumholz et al. (2017) does not capturewell the physics of the CMZ since it predicts the existence of a MNRAS , 000–000 (0000) tar formation in the CMZ Figure 8.
Left column : instantaneous spatial distribution of stars with age in the given range, for a typical snapshot of our simulation. This is a scatter plot.
Middle-left column : same as the left panel, but binned in a 2D histogram weighted by mass.
Middle-right column : time-averaged spatial distribution of stars byage, i.e. obtained by time-averaging the middle-left column. The time averaged is taken over t = ( , ) . Right column : azimuthal distribution of stars withage in the given range in the CMZ ( R ≤
250 pc), obtained by looking at the azimuthal distribution of the histograms in the middle-right column. This showsthat the stars are not distributed uniformly through azimuth, but have distinct peaks whose azimuthal position depends on the age range of the stars considered. quasi-axisymmetric outer CMZ extending out to R (cid:39)
450 pc, whichis not supported by either observations or simulations. Moreover,Sormani & Li (2020) have shown that the acoustic instability onwhich these models are based is a spurious result which cannotdrive turbulence and mass transport in the interstellar medium.We conclude that our simulation supports a scenario which is amixture of the pearls on a string and of the popcorn scenarios. Mostof the star formation happens downstream of the apocentres, buta significant amount of star formation also takes place distributedalong the ring. Our results do not support the pericentre passagescenario.
Star formation relations are empirical correlations between the SFRand properties of the interstellar medium (ISM) from which starsare born. It has been extensively discussed in the literature that theCMZ follows some star formation relations but not all of them (e.g.Yusef-Zadeh et al. 2009; Longmore et al. 2013a; Kruijssen et al.2014; Kauffmann et al. 2017b,a). In particular, it has been shownthat the global
SFR of the CMZ is consistent with the Schmidt-Kennicutt density relation (Schmidt 1959; Kennicutt 1998), withthe Bigiel et al. (2008) molecular gas relation, and with the Bac-
MNRAS000
MNRAS000 , 000–000 (0000) Sormani et al. − . − . . . . x [ kpc ] − . . . y [ kp c ] ρ High τ τ [ km / s / kpc ] Figure 9.
Shear map for the snapshot at t = . τ defined in Equation (2), which is a good indication of shearfor a 2D flow. The white contour indicates where the total surface density ofthe gas is N = × cm − . The red arrow indicates the terminal part ofthe dust lane, a region of high density and high shear, where the maximumof the depletion time as a function of Galactocentric radius is reached (seebottom panel in Figure 7 and discussion in Section 3.2). The map also showsthat high shear occurs predominantly in the dust lanes and in expanding SNshells. − . − . − . − . . . . . . . . . . . . . SF R p e r un it l [ M (cid:12) y r − d e g − ] Averaged over ∆ t = . t = . − − l [ deg ] . . . . . . . . SF R p e r un it l [ M (cid:12) y r − d e g − ] Figure 10.
Star formation rate as a function of Galactic longitude in oursimulation. The magenta line shows the time-averaged distribution, whilethe grey line shows the instantaneous distribution at t = . t = ( . , . ) Myr. − . − . − . . . . . − . − . − . . . . . y [ kp c ] − . − . − . − . . . . . . x [ kpc ] − . − . − . . . . . − . − . − . − . . . . . . − . − . . . . y [ kp c ] Figure 11.
Typical trajectories of newly born stars in our simulations.
Toppanel : for stars formed upstream along the dust lane.
Middle-top panel : forstars formed downstream along the dust lanes.
Middle-bottom panel : forstars formed from gas orbiting in the CMZ ring.
Bottom panel : for starsformed at radii smaller than the CMZ ring. Star markers indicate the currentposition of the stars. Round markers indicate the location where they formed.Full lines indicate the past trajectories, while dashed lines indicate the futuretrajectories. Cross markers log the position of the stars at equal time intervalsof ∆ t = surface density at current time.MNRAS , 000–000 (0000) tar formation in the CMZ − . . . x [ kpc ] − . . . y [ kp c ] Inflow from DLRCloudCollisions − . . . x [ kpc ] ◦ ◦ ◦ ◦ High SF Low SF θ [ deg ] . . . . . . . Σ Y S [ M (cid:12) / p c ] Pericentre PericentreApocentreApocentre − − N H2 [ cm − ] Σ YS [ M (cid:12) / pc ] Figure 12.
Left : time-averaged H surface density. Middle : time-averaged surface density of the very young stars (age t ≤ .
25 Myr) formed in our simulation.
Right : time-averaged surface density of very young stars (age t ≤ .
25 Myr) as a function of azimuth along the elliptical ring shown in the middle panel. Thetime averages are calculated over the time range t = ( . , . ) Myr. chini et al. (2019a,b) volumetric star formation relation (Bacchini,private communication). However, the global SFR of the CMZ isnot consistent with the SFR-dense gas relation observed by e.g. Gao& Solomon (2004), Wu et al. (2005) and Lada et al. (2010, 2012).This is a linear relation between the quantity of dense gas (as tracedby HCN emission or high dust extinction) and the SFR. It has beenshown to work well both for the total (integrated) properties of ex-ternal galaxies, and for local molecular clouds in the MW, whichmade it apparently valid over an impressive 9 orders of magnitude(although with a gap in the middle, see Figure 2 in Lada et al.2012). This generated the expectation that the same law should bevalid for the CMZ, but the data shows that it is not (see Longmoreet al. 2013a; Kruijssen et al. 2014 and Figure 1 in Kauffmann et al.2017a). This expectation, and the universality of the SFR-dense gasrelation, is also challenged by observations that suggest that thecentres of nearby galaxies lie on average below the Lada et al. 2012relation (see Gallagher et al. 2018; Jiménez-Donaire et al. 2019 andin particular Figure 13 in the latter).Figure 13 shows that the CMZ in our simulation follows theSchmidt-Kennicutt relation (Schmidt 1959; Kennicutt 1998) and theBigiel et al. (2008) molecular gas relation, consistent with observa-tional findings. This reassures us that our numerical star formationsubgrid model is working correctly. Unfortunately, our simulationsdo not have the resolution to probe the dense-gas star formationrelations, which the CMZ has been shown to be not consistent with(e.g. Longmore et al. 2013a; Kruijssen et al. 2014; Kauffmann et al.2017b,a). To do that, we would need to increase the sink formationdensity threshold ρ c (see Appendix A) to densities of n (cid:39) cm − ,which is the dense gas formation threshold in the CMZ estimatedby Kruijssen et al. (2014) and Kauffmann et al. (2017a). This isimpractical with our current simulations owing to the very highcomputational expense, but is a worthwhile direction for future in-vestigations. The Arches and Quintuplet clusters are two young massive ( M (cid:38) M (cid:12) ) clusters found close to the Galactic centre ( (cid:39)
30 pc inprojected distance). They have estimated ages of 3 . ± . Figure 13.
Schmidt-Kennicutt plot for our simulation. We bin face-on H and SFR surface densities with a grid size of 100 pc. Each point in this graphrepresents one such bin. The points are coloured based on the position of thecentre of the bin. The underlying distribution is obtained by Gaussian kerneldensity estimation of the points associated to the disc. To increase statis-tics especially for the CMZ we include surface densities of eight differentconsecutive snapshots i.e. over a time of approximately 2 Myr. The CMZapproximately follows the SK as well as the Bigiel et al. (2008) relation, asfound in observations (e.g. Figure 2 in Kruijssen et al. 2014). and 4 . ± . v los = ± − (Figeret al. 2002) and a proper motion velocity of v pm = ±
15 km s − (Clarkson et al. 2012), which yields a 3D orbital velocity of v = ±
17 km s − in the direction of increasing longitude (Clarksonet al. 2012). The Quintuplet cluster has a line-of-sight velocityof v los = ± − and a proper motion velocity of v pm = ±
15 km s − , which yields a 3D orbital velocity of v = ±
15 km s − , also in the direction of increasing longitude (Stolte et al.2014).The observed motions of the Arches and Quintuplet clusters MNRAS000
15 km s − , also in the direction of increasing longitude (Stolte et al.2014).The observed motions of the Arches and Quintuplet clusters MNRAS000 , 000–000 (0000) Sormani et al. l = − ◦ l = ◦ l = ◦ Arches cluster t = .
100 pc l = − ◦ l = ◦ l = ◦ t = . l = − ◦ l = ◦ l = ◦ t = . l = − ◦ l = ◦ l = ◦ Quintuplet cluster t = .
100 pc l = − ◦ l = ◦ l = ◦ t = . l = − ◦ l = ◦ l = ◦ t = . Figure 14.
Sink particles in our simulations with properties (age, line-of-sight velocity, proper motion velocity) within the observational constraints ofthe Arches (left panels) and the Quintuplet (right panels) clusters. Red/violettriangles denote the present day position, while solid and dotted lines showthe past (from the birth site to the current position) and the future trajectories(for the next 5 Myr) respectively. The crosses log the position of the clusterat equal time intervals of 1 Myr. The background shows the gas total densitydistribution at the time when the clusters are at their present day position. can be compared with the trajectories of our sink particles dis-cussed in Section 3.3. We have searched in our simulations for sinkparticles that are within 30 pc of the Galactic centre (in projecteddistance) on the positive longitude side and that have age, line-of-sight velocity and proper motion velocities compatible with thoseof the observed clusters within the observational uncertainties givenabove. Figure 14 shows trajectories for a sample of sinks which arefound according to this procedure. This figure suggests that (i) theArches cluster (left panels) formed from gas that is colliding intothe far side dust lane at negative longitudes while orbiting in theCMZ. All the clusters in our simulation compatible with the aboveobservational constraints are consistent with this picture. (ii) TheQuintuplet cluster formed either in a similar scenario as the Archescluster, but from gas colliding onto the near side dust lane (top-rightpanel), or more probably by gas in the terminal part of the dust laneswhich is just entering the CMZ (middle- and lower-right panels).Occurrences of the second type are more frequent (roughly by afactor of ∼ x and x orbits, since this transition happens atthe contact point between the dust lanes (compare Figure 14 withFigure 12 in Stolte et al. 2014). Kruijssen et al. (2015) have pro-posed that the clusters originated on the same orbit that they use tofit dense gas data. However, we find that sink particles with prop-erties compatible with the observed kinematics of the clusters havetypically decoupled from the gas in which they are born by the timethe clusters have reached their present age. Moreover, the gas orbit- ing in the CMZ ring has typically lower absolute 3D velocities thanthose of the clusters. Therefore, the scenario proposed by Kruijssenet al. (2015) seems to be inconsistent with the result of the presentsimulation. We have used the high-resolution hydrodynamical simulations pre-sented in Paper I to study star formation in the CMZ. These includea realistic Milky Way external barred potential, a time-dependentchemical network that keeps track of hydrogen and carbon chem-istry, a physically motivated model for the formation of new stars us-ing sink particles, and supernovae feedback. The simulations reachsub-parsec resolution in the dense regions and allow us to resolve in-dividual molecular clouds which are formed self-consistently fromthe large-scale flow.Our main conclusions are as follows: • We have studied the temporal distribution of star formation. Wefind that the depletion time in the CMZ is approximately constantin time. This implies that variations in the SFR of the CMZ areprimarily driven by variations in its mass, caused for example bychanges in the bar-driven inflow rate, AGN events or other externalfactors, while the observed scatter in the depletion time of externalgalactic centres is interpreted as variations in the environmentalfactors (e.g. the stellar surface density, Jiménez-Donaire et al. 2019).Contrary to the findings of Armillotta et al. (2019), we do not findthat the depletion time in the CMZ goes through strong oscillatorycycles, at least within the timescale of our simulation ( ∼
100 Myr,see Sections 3.1 and 4.1). • We have studied the spatial distribution of star formation. Mostof the star formation happens in the CMZ ring at R (cid:38)
100 pc, buta significant amount of star formation also occurs closer to SgrA*( R ≤
10 pc, see Section 3.2 and Figure 5). While the time-averagedspatial distribution of the SFR is typically smooth, the instantaneousdistribution can have complex and transient fluctuations which devi-ate significantly from the average morphology (compare the bottompanel in Figure 6 with Figures 4 and 5). Molecular clouds formedself-consistently from the large-scale flow, and their embedded starformation, exhibit complicated filamentary morphologies and donot resemble the idealised “spherical clouds” that are often used asa model to understand star formation. We have also investigated howthe spatial distribution changes when we consider stars in differentage ranges, and found that a bi-polar structure persists even for starswith age 10-20 Myr (see Section 3.2 and Figure 8). • We tested the predictions of the three main scenarios that havebeen put forward to explain the spatial and temporal distribution ofstar formation in the centre of barred galaxies, namely the “pearlson a string", the “popcorn" and the “pericentre passage" scenarios.We found that our simulations are consistent with a mixture of thepearls on a string and popcorn scenarios, while they are inconsistentwith the pericentre passage scenario (see Section 4.2). • We have studied the trajectories of newly born stars (see Figure11). We find that gas and stars typically decouple within at most2-3 Myr (see Sections 3.2 and 3.3). • We have used the trajectories of newly born stars to provide adetailed analysis of the origin of the Arches and Quintuplet clusters.Our simulation favour a scenario in which the Arches cluster isformed from gas that crashed into the far side dust lane at negativelongitudes while orbiting in the CMZ, while the Quintuplet clusteris either formed in a similar event but with the roles of the near/farsides the Galaxy reversed, or more likely by gas in the terminal part
MNRAS , 000–000 (0000) tar formation in the CMZ of the near side dust lane which was just entering the CMZ (seeFigure 14 and Section 4.4). ACKNOWLEDGEMENTS
We thank the anonymous referee for constructive comments thatimproved the paper. MCS thanks Ashley Barnes, Adam Gins-burg, Jonathan Henshaw, Zhi Li, Diederik Kruijssen, AlessandraMastrobuono-Battisti, Elisabeth Mills, Nadine Neumayer, FranciscoNogueras-Lara and Ralph Schoenrich for useful comments and dis-cussions. MCS and RGT both acknowledge contributing equallyto this paper. MCS, RGT, SCOG, and RSK acknowledge financialsupport from the German Research Foundation (DFG) via the col-laborative research center (SFB 881, Project-ID 138713538) “TheMilky Way System” (subprojects A1, B1, B2, and B8). CDB andHPH gratefully acknowledges support for this work from the Na-tional Science Foundation under Grant No. (1816715). PCC ac-knowledges support from the Science and Technology FacilitiesCouncil (under grant ST/K00926/1) and StarFormMapper, a projectthat has received funding from the European Union’s Horizon 2020Research and Innovation Programme, under grant agreement no.687528. HPH thanks the LSSTC Data Science Fellowship Pro-gram, which is funded by LSSTC, NSF Cybertraining Grant No.1829740, the Brinson Foundation, and the Moore Foundation; hisparticipation in the program has benefited this work. RJS grate-fully acknowledges an STFC Ernest Rutherford fellowship (grantST/N00485X/1) and HPC from the Durham DiRAC supercomput-ing facility (grants ST/P002293/1, ST/R002371/1, ST/S002502/1,and ST/R000832/1). The authors acknowledge support by the stateof Baden-Württemberg through bwHPC and the German ResearchFoundation (DFG) through grant INST 35/1134-1 FUGG. The au-thors gratefully acknowledge the data storage service SDS@hd sup-ported by the Ministry of Science, Research and the Arts Baden-Württemberg (MWK) and the German Research Foundation (DFG)through grant INST 35/1314-1 FUGG. The data underlying thisarticle will be shared on reasonable request to the correspondingauthor.
REFERENCES
Armillotta L., Krumholz M. R., Di Teodoro E. M., McClure-Griffiths N. M.,2019, MNRAS, 490, 4401Armillotta L., Krumholz M. R., Di Teodoro E. M., 2020, MNRAS, 493,5273Baba J., Kawata D., 2020, MNRAS, 492, 4500Bacchini C., Fraternali F., Iorio G., Pezzulli G., 2019a, A&A, 622, A64Bacchini C., Fraternali F., Pezzulli G., Marasco A., Iorio G., Nipoti C.,2019b, A&A, 632, A127Bally J., et al., 2010, ApJ, 721, 137Barnes A. T., Longmore S. N., Battersby C., Bally J., Kruijssen J. M. D.,Henshaw J. D., Walker D. L., 2017, MNRAS, 469, 2263Bigiel F., Leroy A., Walter F., Brinks E., de Blok W. J. G., Madore B.,Thornley M. D., 2008, AJ, 136, 2846Böker T., Falcón-Barroso J., Schinnerer E., Knapen J. H., Ryder S., 2008,AJ, 135, 479Chatzopoulos S., Fritz T. K., Gerhard O., Gillessen S., Wegg C., Genzel R.,Pfuhl O., 2015, MNRAS, 447, 948Clark P. C., Glover S. C. O., Ragan S. E., Shetty R., Klessen R. S., 2013,ApJ, 768, L34Clarkson W. I., Ghez A. M., Morris M. R., Lu J. R., Stolte A., McCrady N.,Do T., Yelda S., 2012, ApJ, 751, 132 Comerón S., Knapen J. H., Beckman J. E., Laurikainen E., Salo H., Martínez-Valpuesta I., Buta R. J., 2010, MNRAS, 402, 2462Dahmen G., Huttemeister S., Wilson T. L., Mauersberger R., 1998, A&A,331, 959Dale J. E., Kruijssen J. M. D., Longmore S. N., 2019, MNRAS, 486, 3307Elmegreen B. G., Galliano E., Alloin D., 2009, ApJ, 703, 1297Emsellem E., Renaud F., Bournaud F., Elmegreen B., Combes F., GaborJ. M., 2015, MNRAS, 446, 2468Federrath C., et al., 2016, ApJ, 832, 143Feldmeier-Krause A., et al., 2015, A&A, 584, A2Feldmeier A., et al., 2014, A&A, 570, A2Figer D. F., et al., 2002, ApJ, 581, 258Gallagher M. J., et al., 2018, ApJ, 858, 90Gallego-Cano E., Schödel R., Nogueras-Lara F., Dong H., ShahzamanianB., Fritz T. K., Gallego-Calvente A. T., Neumayer N., 2020, A&A, 634,A71Gao Y., Solomon P. M., 2004, ApJ, 606, 271Genzel R., Eisenhauer F., Gillessen S., 2010, Reviews of Modern Physics,82, 3121Ginsburg A., et al., 2016, A&A, 586, A50Gravity Collaboration et al., 2019, A&A, 625, L10Guesten R., Henkel C., 1983, A&A, 125, 136Heywood I., et al., 2019, Nature, 573, 235Immer K., Schuller F., Omont A., Menten K. M., 2012, A&A, 537, A121Immer K., Kauffmann J., Pillai T., Ginsburg A., Menten K. M., 2016, A&A,595, A94Jeffreson S. M. R., Kruijssen J. M. D., Krumholz M. R., Longmore S. N.,2018, MNRAS, 478, 3380Jiménez-Donaire M. J., et al., 2019, ApJ, 880, 127Kauffmann J., Pillai T., Zhang Q., Menten K. M., Goldsmith P. F., Lu X.,Guzmán A. E., Schmiedeke A., 2017a, A&A, 603, A90Kauffmann J., Goldsmith P. F., Melnick G., Tolls V., Guzman A., MentenK. M., 2017b, A&A, 605, L5Kennicutt Robert C. J., 1998, ApJ, 498, 541Kennicutt R. C., Evans N. J., 2012, ARA&A, 50, 531Klessen R. S., Glover S. C. O., 2016, in: Star Formation in Galaxy Evolution:Connecting Numerical Models to Reality, Saas-Fee Advanced Course,43, 85Krieger N., et al., 2017, ApJ, 850, 77Kruijssen J. M. D., Longmore S. N., Elmegreen B. G., Murray N., Bally J.,Testi L., Kennicutt R. C., 2014, MNRAS, 440, 3370Kruijssen J. M. D., Dale J. E., Longmore S. N., 2015, MNRAS, 447, 1059Kruijssen J. M. D., et al., 2019, MNRAS, 484, 5734Krumholz M. R., Kruijssen J. M. D., 2015, MNRAS, 453, 739Krumholz M. R., Kruijssen J. M. D., Crocker R. M., 2017, MNRAS, 466,1213Lada C. J., Lombardi M., Alves J. F., 2010, ApJ, 724, 687Lada C. J., Forbrich J., Lombardi M., Alves J. F., 2012, ApJ, 745, 190Launhardt R., Zylka R., Mezger P. G., 2002, A&A, 384, 112Leroy A. K., et al., 2013, AJ, 146, 19Li Z., Shen J., Kim W.-T., 2015, ApJ, 806, 150Longmore S. N., et al., 2013a, MNRAS, 429, 987Longmore S. N., et al., 2013b, MNRAS, 433, L15Longmore S. N., et al., 2017, MNRAS, 470, 1462Mac Low M.-M., Klessen R. S., 2004, Reviews of Modern Physics, 76, 125Maciejewski W., 2008, in Bureau M., Athanassoula E., Barbuy B., eds, IAUSymposium Vol. 245, Formation and Evolution of Galaxy Bulges. pp161–164 ( arXiv:0709.3959 ), doi:10.1017/S1743921308017547Mangilli A., et al., 2019, A&A, 630, A74Mazzuca L. M., Knapen J. H., Veilleux S., Regan M. W., 2008, ApJS, 174,337Mills E. A. C., Togi A., Kaufman M., 2017, ApJ, 850, 192Mills E. A. C., Ginsburg A., Immer K., Barnes J. M., Wiesenfeld L., FaureA., Morris M. R., Requena-Torres M. A., 2018, ApJ, 868, 7Molinari S., et al., 2011, ApJ, 735, L33Morris M. R., 2015, Manifestations of the Galactic Center MagneticField. Springer International Publishing, p. 391, doi:10.1007/978-3-319-10614-4_32MNRAS000
Armillotta L., Krumholz M. R., Di Teodoro E. M., McClure-Griffiths N. M.,2019, MNRAS, 490, 4401Armillotta L., Krumholz M. R., Di Teodoro E. M., 2020, MNRAS, 493,5273Baba J., Kawata D., 2020, MNRAS, 492, 4500Bacchini C., Fraternali F., Iorio G., Pezzulli G., 2019a, A&A, 622, A64Bacchini C., Fraternali F., Pezzulli G., Marasco A., Iorio G., Nipoti C.,2019b, A&A, 632, A127Bally J., et al., 2010, ApJ, 721, 137Barnes A. T., Longmore S. N., Battersby C., Bally J., Kruijssen J. M. D.,Henshaw J. D., Walker D. L., 2017, MNRAS, 469, 2263Bigiel F., Leroy A., Walter F., Brinks E., de Blok W. J. G., Madore B.,Thornley M. D., 2008, AJ, 136, 2846Böker T., Falcón-Barroso J., Schinnerer E., Knapen J. H., Ryder S., 2008,AJ, 135, 479Chatzopoulos S., Fritz T. K., Gerhard O., Gillessen S., Wegg C., Genzel R.,Pfuhl O., 2015, MNRAS, 447, 948Clark P. C., Glover S. C. O., Ragan S. E., Shetty R., Klessen R. S., 2013,ApJ, 768, L34Clarkson W. I., Ghez A. M., Morris M. R., Lu J. R., Stolte A., McCrady N.,Do T., Yelda S., 2012, ApJ, 751, 132 Comerón S., Knapen J. H., Beckman J. E., Laurikainen E., Salo H., Martínez-Valpuesta I., Buta R. J., 2010, MNRAS, 402, 2462Dahmen G., Huttemeister S., Wilson T. L., Mauersberger R., 1998, A&A,331, 959Dale J. E., Kruijssen J. M. D., Longmore S. N., 2019, MNRAS, 486, 3307Elmegreen B. G., Galliano E., Alloin D., 2009, ApJ, 703, 1297Emsellem E., Renaud F., Bournaud F., Elmegreen B., Combes F., GaborJ. M., 2015, MNRAS, 446, 2468Federrath C., et al., 2016, ApJ, 832, 143Feldmeier-Krause A., et al., 2015, A&A, 584, A2Feldmeier A., et al., 2014, A&A, 570, A2Figer D. F., et al., 2002, ApJ, 581, 258Gallagher M. J., et al., 2018, ApJ, 858, 90Gallego-Cano E., Schödel R., Nogueras-Lara F., Dong H., ShahzamanianB., Fritz T. K., Gallego-Calvente A. T., Neumayer N., 2020, A&A, 634,A71Gao Y., Solomon P. M., 2004, ApJ, 606, 271Genzel R., Eisenhauer F., Gillessen S., 2010, Reviews of Modern Physics,82, 3121Ginsburg A., et al., 2016, A&A, 586, A50Gravity Collaboration et al., 2019, A&A, 625, L10Guesten R., Henkel C., 1983, A&A, 125, 136Heywood I., et al., 2019, Nature, 573, 235Immer K., Schuller F., Omont A., Menten K. M., 2012, A&A, 537, A121Immer K., Kauffmann J., Pillai T., Ginsburg A., Menten K. M., 2016, A&A,595, A94Jeffreson S. M. R., Kruijssen J. M. D., Krumholz M. R., Longmore S. N.,2018, MNRAS, 478, 3380Jiménez-Donaire M. J., et al., 2019, ApJ, 880, 127Kauffmann J., Pillai T., Zhang Q., Menten K. M., Goldsmith P. F., Lu X.,Guzmán A. E., Schmiedeke A., 2017a, A&A, 603, A90Kauffmann J., Goldsmith P. F., Melnick G., Tolls V., Guzman A., MentenK. M., 2017b, A&A, 605, L5Kennicutt Robert C. J., 1998, ApJ, 498, 541Kennicutt R. C., Evans N. J., 2012, ARA&A, 50, 531Klessen R. S., Glover S. C. O., 2016, in: Star Formation in Galaxy Evolution:Connecting Numerical Models to Reality, Saas-Fee Advanced Course,43, 85Krieger N., et al., 2017, ApJ, 850, 77Kruijssen J. M. D., Longmore S. N., Elmegreen B. G., Murray N., Bally J.,Testi L., Kennicutt R. C., 2014, MNRAS, 440, 3370Kruijssen J. M. D., Dale J. E., Longmore S. N., 2015, MNRAS, 447, 1059Kruijssen J. M. D., et al., 2019, MNRAS, 484, 5734Krumholz M. R., Kruijssen J. M. D., 2015, MNRAS, 453, 739Krumholz M. R., Kruijssen J. M. D., Crocker R. M., 2017, MNRAS, 466,1213Lada C. J., Lombardi M., Alves J. F., 2010, ApJ, 724, 687Lada C. J., Forbrich J., Lombardi M., Alves J. F., 2012, ApJ, 745, 190Launhardt R., Zylka R., Mezger P. G., 2002, A&A, 384, 112Leroy A. K., et al., 2013, AJ, 146, 19Li Z., Shen J., Kim W.-T., 2015, ApJ, 806, 150Longmore S. N., et al., 2013a, MNRAS, 429, 987Longmore S. N., et al., 2013b, MNRAS, 433, L15Longmore S. N., et al., 2017, MNRAS, 470, 1462Mac Low M.-M., Klessen R. S., 2004, Reviews of Modern Physics, 76, 125Maciejewski W., 2008, in Bureau M., Athanassoula E., Barbuy B., eds, IAUSymposium Vol. 245, Formation and Evolution of Galaxy Bulges. pp161–164 ( arXiv:0709.3959 ), doi:10.1017/S1743921308017547Mangilli A., et al., 2019, A&A, 630, A74Mazzuca L. M., Knapen J. H., Veilleux S., Regan M. W., 2008, ApJS, 174,337Mills E. A. C., Togi A., Kaufman M., 2017, ApJ, 850, 192Mills E. A. C., Ginsburg A., Immer K., Barnes J. M., Wiesenfeld L., FaureA., Morris M. R., Requena-Torres M. A., 2018, ApJ, 868, 7Molinari S., et al., 2011, ApJ, 735, L33Morris M. R., 2015, Manifestations of the Galactic Center MagneticField. Springer International Publishing, p. 391, doi:10.1007/978-3-319-10614-4_32MNRAS000 , 000–000 (0000) Sormani et al.
Mou G., Sun D., Xie F., 2018, ApJ, 869, L20Neumayer N., Seth A., Boeker T., 2020, arXiv e-prints, p. arXiv:2001.03626Nishiyama S., et al., 2013, ApJ, 769, L28Nogueras-Lara F., et al., 2020, Nature Astronomy, 4, 377Oka T., Geballe T. R., Goto M., Usuda T., Benjamin McCall J., Indriolo N.,2019, ApJ, 883, 54Ostriker E. C., Shetty R., 2011, ApJ, 731, 41Ostriker E. C., McKee C. F., Leroy A. K., 2010, ApJ, 721, 975Ponti G., et al., 2015, MNRAS, 453, 172Ponti G., et al., 2019, Nature, 567, 347Portail M., Gerhard O., Wegg C., Ness M., 2017, MNRAS, 465, 1621Reid M. J., et al., 2019, ApJ, 885, 131Renaud F., et al., 2013, MNRAS, 436, 1836Renaud F., et al., 2015, MNRAS, 454, 3299Sanders J. L., Smith L., Evans N. W., 2019, arXiv e-prints 1903.02009,Sarzi M., Allard E. L., Knapen J. H., Mazzuca L. M., 2007, MNRAS, 380,949Schmidt M., 1959, ApJ, 129, 243Schneider F. R. N., et al., 2014, ApJ, 780, 117Schödel R., Feldmeier A., Kunneriath D., Stolovy S., Neumayer N., Amaro-Seoane P., Nishiyama S., 2014, A&A, 566, A47Schönrich R., Binney J., Dehnen W., 2010, MNRAS, 403, 1829Schönrich R., Aumer M., Sale S. E., 2015, ApJ, 812, L21Seo W.-Y., Kim W.-T., Kwak S., Hsieh P.-Y., Han C., Hopkins P. F., 2019,ApJ, 872, 5Shetty R., Beaumont C. N., Burton M. G., Kelly B. C., Klessen R. S., 2012,MNRAS, 425, 720Smith B. D., et al., 2017, MNRAS, 466, 2217Sormani M. C., Barnes A. T., 2019, MNRAS, 484, 1213Sormani M. C., Li Z., 2020, arXiv e-prints 2002.10559,Sormani M. C., Binney J., Magorrian J., 2015, MNRAS, 454, 1818Sormani M. C., Treß R. G., Klessen R. S., Glover S. C. O., 2017, MNRAS,466, 407Sormani M. C., Treß R. G., Ridley M., Glover S. C. O., Klessen R. S.,Binney J., Magorrian J., Smith R., 2018a, MNRAS, 475, 2383Sormani M. C., Sobacchi E., Fragkoudi F., Ridley M., Treß R. G., GloverS. C. O., Klessen R. S., 2018b, MNRAS, 481, 2Sormani M. C., et al., 2019, MNRAS, 488, 4663Springel V., 2010, MNRAS, 401, 791Stolte A., Ghez A. M., Morris M., Lu J. R., Brand ner W., Matthews K.,2008, ApJ, 675, 1278Stolte A., et al., 2014, ApJ, 789, 115Su M., Slatyer T. R., Finkbeiner D. P., 2010, ApJ, 724, 1044Torrey P., Hopkins P. F., Faucher-Giguère C.-A., Vogelsberger M., QuataertE., Kereš D., Murray N., 2017, MNRAS, 467, 2301Tress R. G., Sormani M. C., Glover S. C. O., Klessen R. S., Battersby C. D.,Clark P. C., Hatchfield H. P., Smith R. J., 2020, Paper I (submitted),Tsatsi A., Mastrobuono-Battisti A., van de Ven G., Perets H. B., BianchiniP., Neumayer N., 2017, MNRAS, 464, 3720Utomo D., et al., 2017, ApJ, 849, 26Walmsley C. M., Guesten R., Angerhofer P., Churchwell E., Mundy L.,1986, A&A, 155, 129Weinberger R., Springel V., Pakmor R., 2019, arXiv e-prints;arXiv:1909.04667,Wu J., Evans Neal J. I., Gao Y., Solomon P. M., Shirley Y. L., Vanden BoutP. A., 2005, ApJ, 635, L173Yusef-Zadeh F., Braatz J., Wardle M., Roberts D., 2008, ApJ, 683, L147Yusef-Zadeh F., et al., 2009, ApJ, 702, 178Yusef-Zadeh F., et al., 2015, ApJ, 808, 97
APPENDIX A: RESOLUTION STUDY
In this appendix we show the results of a resolution study that wehave conducted in order to assess the impact of varying the reso-lution and the sink particle creation threshold ρ c (see Section 2.4 name M base [ M (cid:12) ] ρ c [ g cm − ] m1000densc1e2 1000 10 − m300densc1e3 300 10 − m100densc1e3 100 10 − m100densc1e4 (fiducial) 100 10 − Table A1.
Summary of the simulations considered in the resolution study. M base is the base target cell mass. No cells in the simulations are allowedto fall below this resolution (i.e., no cells can have mass higher than M base ). ρ c is the sink particle formation threshold (see Section 2.4 of Paper I).m100densc1e4 is the fiducial simulation considered for analysis in the maintext of this paper and of Paper I. − − − ρ [ g / cm ] − M ce ll [ M (cid:12) ] m100densc1e4m100densc1e3m300densc1e3m1000densc1e2 − − − ρ [ g / cm ] Figure A1.
Mass resolution of the four simulation considered in the reso-lution study. The histogram shows the distribution of number of cells in the ( M cell , ρ ) plane. of Paper I). We consider 4 simulations, whose properties are sum-marised in Table A1.The simulations differ for two parameters: the base target cellmass M , and the sink particle formation density threshold ρ c . Thefiducial simulation (m100densc1e4) has both the smallest M (high-est resolution) and the highest ρ c . The simulation m100densc1e3has the same M but lower ρ c . This allows us to assess the impact ofhaving lower resolution in the high density regions where the grav-itational collapse is happening, which is important in the context ofstar formation. Then we consider a simulation with a higher M (i.e.,lower resolution), m300densc1e3, in order to assess how a differentbase mass resolution affects the various phases of the ISM. Finally,we consider a very low resolution simulation, m1000densc1e2, asa general benchmark. Figure A1 shows the mass resolution of thefour simulations as a function of density.Figure A2 shows the behaviour of various quantities as a func-tion of time for the four simulations considered in the resolutionstudy. From this figure we see that the largest difference betweenthe different simulations is seen in the chemical mass fraction: athigher resolution there is roughly a factor of 2 more gas in molecularform (H ) than in lower resolution simulations (see second panelfrom top to bottom). This induces a similar difference in the H depletion times. The SFR and the total gas depletion times does notappear to change substantially between the different simulations.While this is encouraging and gives us confidence in the results ofour main simulation, we caution against drawing too many conclu-sions about convergence from this. We cannot rule out that a furtherincrease in resolution may show major differences, since the starformation process is not resolved in our simulations. MNRAS , 000–000 (0000) tar formation in the CMZ − − S FR [ M (cid:12) y r − ] SFR
CMZGlobalm1000densc1e2m300densc1e3m100densc1e3m100densc1e4 . . . . . . M i / M t o t Chem mass fractions H HIH + τ d e p [ y r ] Total gas depletion time without gas in sinkswith gas in sinks τ d e p [ y r ] H depletion time
80 100 120 140 160 180 200 220 t [ Myr ] τ d e p [ y r ] HI depletion time
Figure A2.
Various quantities as a function of time for the four simulationsconsidered in our resolution study. Different colours indicate different sim-ulation (see Table A1). Thick lines indicate the CMZ, defined as the regionwithin R ≤
250 pc, while thin lines indicate all the simulated box.MNRAS000