© The Authors 2026

1 Introduction

Galactic bars significantly affect the gas dynamics in the centers of barred galaxies and fundamentally change the picture of how the interstellar medium (ISM) and star formation (SF) evolve in the galactic environment. Bars efficiently drive large amounts of gas into inflows from the disk to the inner galaxy, where the gas generally settles into a central ring or disk-like structure (Mazzuca et al. 2008; Comerón et al. 2010). In the Milky Way, this structure is known as the Central Molecular Zone (CMZ) (Morris & Serabyn 1996; Henshaw et al. 2023) and has a diameter of approximately 600 pc.

Recent multifrequency observations by the PHANGS (Physics at High Angular resolution in Nearby GalaxieS) collaboration (Leroy et al. 2021 ; Emsellem et al. 2022; Lee et al. 2023) have started to uncover the physical conditions in the CMZs of many nearby spiral galaxies (e.g. Schinnerer et al. 2023; Sun et al. 2024), emphasizing the complexity and variety in terms of sizes, star formation, and morphology of these regions (for a review, see Schinnerer & Leroy 2024). Our own Milky Way Galaxy is barred as well, and exhibits a similar configuration, although its 3D morphology and structure are much more challenging to discern (Armillotta et al. 2019; Henshaw et al. 2023; Battersby et al. 2025a,b; Walker et al. 2025; Lipman et al. 2025; Sofue et al. 2025).

The surface density in the CMZ is higher by up to 10³ times than typical values in the disk. This is close to conditions experienced in the early Universe. More subtly, the bar induces a variety of other extreme and unique conditions, so that the CMZ is much more than a high surface density version of the solar neighborhood. For instance, high shear, short orbital and dynamical times, and frequent cloud collisions due to bar-induced large-scale gas inflows all occur on timescales relevant to the evolution of giant molecular clouds (GMCs) and SF (Regan et al. 1997; Laine et al. 1999; Elmegreen et al. 2009; Shimizu et al. 2019; Seo et al. 2019; Hatchfield et al. 2021; Sormani et al. 2023). It is a highly turbulent environment (Rathborne et al. 2014, 2015; Federrath et al. 2016), exhibiting an unusually steep size-linewidth relation (Kauffmann et al. 2017; Tassis & Pavlidou 2022). Finally, observations of the magnetic field in the CMZ suggest that the region is highly magnetized, with a complex magnetic field structure (see, for example, Crutcher et al. 1996; Lu et al. 2024; Butterfield et al. 2024a,b). The abundance of prominent nonthermal filaments observed in the Galactic center are likely associated with a strong pervasive magnetic field (Heywood et al. 2022; Yusef-Zadeh et al. 2022). In the dense molecular clouds, intensities of up to several milligauss have been measured, suggesting that magnetic field contributions to the dynamics of the ISM in the region are important (Chuss et al. 2003; Pillai et al. 2015). Strong magnetic fields can suppress star formation, and magnetic instabilities can induce turbulence and drive inflows (see, for example, Hennebelle & Inutsuka 2019; Moon et al. 2023; Tress et al. 2024).

It is still unclear whether these effects drive a significant deviation from the classic SF scaling relations (eg. Sun et al. 2023) and regulating mechanisms (e.g. Ostriker & Kim 2022). Because of these extreme conditions, the Galactic center is indeed a unique laboratory for testing the universality of the star formation process. The Milky Way center seems to follow the scaling law between the gas column density and the surface density of the star formation rate given by the Kennicutt-Schmidt relation (Kennicutt & Evans 2012). However, it underproduces stars by about one order of magnitude relative to what is expected by its dense (>10⁴ cm⁻³) gas content (Longmore et al. 2013), and some molecular clouds seem to be essentially devoid of massive star formation despite their high masses and densities (Rodríguez & Zapata 2013; Kauffmann et al. 2017; Walker et al. 2021). Even so, on a statistical level, the departure from these scaling relations mostly affects the more extended dense gas, while on the scale of star-forming molecular clouds, the agreement is better instead (Battersby et al. 2025b). Moreover, molecular gas properties in barred galaxy centers seem distinctly offset compared to other systems with higher molecular gas surface densities and velocity dispersions (Sun et al. 2020).

These special environmental conditions make it clear that isolated cloud simulations, while useful, are insufficient for capturing the complex feedback loops that are present in real CMZs, even if they are exposed to CMZ-like radiation fields and cosmic-ray fluxes (Clark et al. 2013; Bertram et al. 2015; Cusack et al. 2025). While achieving excellent resolution close to the pro-toplanetary disk scales, they fail to include the larger context, and they generally fail to reproduce observed properties such as star formation efficiencies. An improvement can be achieved by evolving similar isolated clouds in the galactic potential of the Galactic center, through which effects such as galactic shear and tidal compressions can be captured (Dale et al. 2019; Kruijssen et al. 2019; Petkova et al. 2023), but it still lacks all the interactions and complexity arising from a rich barred ecosystem, such as gas inflow and cloud collisions.

Galactic-scale simulations, on the other hand, often lack the necessary resolution to accurately capture star formation and GMC structure, and/or they miss important physics that are needed to account for a realistic matter cycle (Seo et al. 2019; Tress et al. 2020; Hatchfield et al. 2021; Tress et al. 2024; Verwilghen et al. 2024, 2025). Idealized setups have been developed to mimic Galactic center conditions and inflow by the bar (Moon et al. 2021, 2022, 2023). These models still overlook the self-consistent interplay between all the important components acting nonlinearly in unison to shape the gas flows in the barred region. For instance, instabilities arise in the bar lanes that cause the intermittent inflow toward the CMZ (Wada & Koda 2004; Sormani et al. 2017, 2018), torques between the gas and the bar regulate inflow rates, and SF in the bar lanes can modulate the accretion and the young stars can provide peculiar high-angular momentum feedback sources in the CMZ. Moreover, a coherent and systematic comparison of ISM properties and SF in the CMZ with respect to the disk using the same numerical framework is still lacking. However, it would be extremely useful for highlighting differences and understanding the regulating mechanisms in the Galactic Center.

The aim of this paper is to use a powerful new zoom simulation that combines the Galactic-scale flows with detailed high-resolution physics to fill this gap. This allows us to study how the complex environment of the CMZ affects the evolution of a typical young star-forming region. In particular, we looked at how the effect of the dense environment, high turbulence, frequent cloud collisions, strong galactic shear, and short orbital times affect how stellar feedback couples to the surrounding ISM compared to the solar neighborhood. The paper is constructed as follows: we first present our numerical setup in Section 2 and then compare the properties of the ISM in the solar neighborhood and the CMZ in Section 3. We then focus on the surrounding environment of young massive stars and the effect of their feedback on the ISM (Section 4). The results of that analysis are discussed in Section 5, where we introduce a new framework of how feedback couples to the ISM in the CMZ. We discus some caveats in Section 6 and present our conclusions in Section 7.

2 Methods

We designed radiation magnetohydrodynamic (MHD) simulations of the ISM of a Milky Way-like barred disk galaxy. We solve the ideal MHD equations using the code AREPO (Springel 2010, Sect. 2.1). The gas evolves in a background barred potential fine-tuned to the Milky Way (Sect. 2.2) and high-resolution regions were installed in the two environments we want to compare (Galactic Center vs solar neighborhood) (Sect. 2.3). The code is coupled to the relevant chemistry (Sect. 2.4), star formation (Sect. 2.5) and stellar feedback modules such as radiative feedback and supernovae (Sect. 2.6), necessary to capture GMC and SF physics. After a build-up phase (Sect. 2.7), we evolve the system for ∼110 Myr. This model builds upon the Rhea simulations (Göller et al. 2025) and can therefore be referred to as Rhea-RT.

2.1 Solving the MHD equations

The magnetized ISM evolves by obeying the ideal MHD equations, which we solved numerically using the moving-mesh code AREPO (Springel 2010; Pakmor et al. 2016; Weinberger et al. 2020). For the divergence control of the magnetic field, we used the Powell method, which has proven to be a robust and efficient implementation, despite not enforcing the divergence-free constraint (∇ · B = 0) to machine precision. Saturation levels of the magnetic field have been found to be relatively robust with this method, although timescales over which the magnetic field grows through dynamo action are probably not (Pakmor et al. 2016).

Conserved quantities are discretized on an unstructured grid, and fluxes between cells are computed by solving the Riemann problem in the rest frame of the faces of the cells using the HLLD Riemann solver (Miyoshi & Kusano 2005). The grid is constructed on each timestep by building the Voronoi mesh of a set of mesh-generating points, which move with the direction of the flow, minimizing advection errors, improving Galilean invariance, and achieving automatic adaptive resolution with density. Cells can, however, still be de/refined by customized refinement schemes based on the problem (see Sect. 2.3). To close the set of fluid equations, we used an adiabatic equation of state for the gas with adiabatic index γ = 5/3.

2.2 Galactic potential

We modeled the gravity of the different mass components of a Milky Way type galaxy by an externally imposed Galactic potential, which defines the gravitational forces felt by the gas. We used the potential introduced by Hunter et al. (2024), which is fine-tuned to reproduce the Milky Way case and is coupled to AREPO through the galaxy modeling software AGAMA (Vasiliev 2019). The external gravitational potential is defined by the following mass components: $$\begin{array}{rl}{\mathrm{\Phi }}_{\mathrm{t}\mathrm{o}\mathrm{t}}=& {\mathrm{\Phi }}_{\mathrm{N}\mathrm{S}\mathrm{C}}+{\mathrm{\Phi }}_{\mathrm{N}\mathrm{S}\mathrm{D}}+{\mathrm{\Phi }}_{\mathrm{B}\mathrm{a}\mathrm{r}}\\ & +{\mathrm{\Phi }}_{\mathrm{D}\mathrm{i}\mathrm{s}\mathrm{k}}+\mathrm{\Delta }{\mathrm{\Phi }}_{\mathrm{S}\mathrm{p}\mathrm{i}\mathrm{r}\mathrm{a}\mathrm{l}}+{\mathrm{\Phi }}_{\mathrm{D}\mathrm{M}},\end{array}$$(1)

where Φ_NSC is the nuclear star cluster; Φ_NSD the nuclear stellar disk; Φ_Bar the Galactic bar which is divided into the X-shaped boxy-peanut component and the long bar; Φ_Disk the Galactic stellar disk, which in turn is divided into the thin and thick disk; ∆Φ_Spiral a rotating spiral arm perturbation; and Φ_DM an Einasto (Einasto 1965; Merritt et al. 2006) dark matter halo. We refer to Hunter et al. (2024) for the details of these components and their functional dependency and parameters.

Improvements compared to our previously employed potential (Tress et al. 2020, 2024) are mainly in the mass distribution of the inner Galaxy, which now is in much better agreement with observations. In particular, the inner ∼300 pc are now modeled in greater detail. We include now the contribution of the nuclear star cluster, which contains ∼2.5 × 10⁷ M_⊙ (e.g. Launhardt et al. 2002; Schödel et al. 2014a; Feldmeier et al. 2014; Chatzopoulos et al. 2015; Fritz et al. 2016; Feldmeier-Krause et al. 2017) and has an effective radius of about 4-5 pc (e.g. Schödel et al. 2014b; Gallego-Cano et al. 2020), dominating the gravity within the central ∼20 pc. We also include the nuclear stellar disk, which dominates the dynamics up to 300 pc instead. This is a more extended and flattened structure with a stellar mass near 10⁹ M_⊙ (e.g. Launhardt et al. 2002; Nogueras-Lara et al. 2020; Sormani et al. 2020, 2022b), a scale length of ∼100 pc, and a scale height of ∼40 pc (e.g. Launhardt et al. 2002; Gallego-Cano et al. 2020; Sormani et al. 2022b. Also the Galactic Bar model has improved substantially, which now follows the made-to-measure model from Portail et al. (2017) and adapted by Sormani et al. (2022a).

The gravity and dynamics of the gaseous component and the newly formed stars are followed self-consistently instead by solving the Poisson equation using the N-body gravity solver of the AREPO code. The gravitational softening of the gas is adaptive and set to twice the size of a cell. Star particles, on the other hand, have a gravitational softening of 10 pc.

The central supermassive black hole (SMBH) of the Milky Way dominates the potential only within the central parsec and is modeled in our simulations as an accreting sink particle following the approach by Petersson et al. (2025) in the Noctua suite of simulations implemented in the arepoNoctua branch of AREPO. The inclusion of this component is necessary to avoid accumulation of gas in the center which would result in spurious star formation and feedback events. The central sink particle can also be used to study gas accretion onto Sgr A* and its resulting active galactic nucleus (AGN) activity, although these are outside of the scope of our current study. The accretion radius of the central sink is fixed at 3 pc. Its mass is 4 × 10⁶ M_⊙ which is reinitialized at the end of the pre-processing phase (Sect. 2.7). This is slightly below the estimated mass of Sgr A* of 4.3 × 10⁶ M_⊙ (GRAVITY Collaboration 2022). As the black hole mass will grow during the simulation, it will ideally reach the observed mass in the middle of the simulation time. The gravitational softening of the SMBH particle is set to 2 pc, within the accretion radius.

2.3 Resolution

Our aim was to compare ISM conditions, star formation, early cluster evolution, and feedback coupling in the CMZ with a solar neighborhood environment. To achieve this, we set up two highresolution regions in the simulation, one region at 7.5 kpc < R < 8.5 kpc, which throughout the text we refer to as the solar circle throughout, and another region in the central R < 4 kpc where the dynamics are controlled by the Galactic bar, which we refer to as the barred region. The central R < 500 pc we define as the CMZ region, and our analysis focused on the ISM and SF in the CMZ region by comparing it to the solar circle. We still simulated the entire barred region at higher resolution, as the bar dynamically connects the inner part of the disk with the CMZ on timescales relevant for GMC evolution, and it is comparable to orbital timescales of gas in the CMZ. The different regions are highlighted in Fig. 1. We set a mass resolution of M_cell = 20 M_⊙ in these regions, and M_cell = 500 M_⊙ for the rest of the galaxy.

In the region R < 11 kpc and |z| < 1 kpc, we set a maximum volume of V_cell,max < 10⁵ pc³, which corresponds to a spherical radius of r_cell ≃ 30 pc. This is to improve the resolution in the low density region, and avoid un-physically large SN injection regions, when clustered SNe carve out low density cavities. Moreover, we restricted the volume ratio of neighboring cells to a maximum of 8.

When a density regime is reached where with the target mass resolution the code is unable to even resolve a single Jeans mass, it is pointless to follow the collapse further. Once this regime is reached, the gas cells were not allowed to further shrink in size as this would lead to spurious fragmentation of the unresolved collapse. To this end, we set a minimum volume for the resolution elements of 1.7 pc³ which is given by M_cell = M_J ie. $V_{\mathrm{m}\mathrm{i}\mathrm{n}}\propto M_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}^3$ and corresponds to a spherical radius of ≃0.2 pc. Moreover this gas is converted into stars with an enhanced efficiency of e_ff = 1 (see Section 2.5 below).

The mass and spatial resolution as a function of density is summarized in Fig. 2. In the high-resolution region we resolve the Jeans mass with at least 8 resolution elements up to a density of n = 2 × 10² cm⁻³ and with at least one resolution element up to a density of n = 1.3 × 10⁴ cm⁻³. At those densities the cold ISM has average temperatures which are above our imposed temperature floor, and so we expect the dynamics and fragmentation properties of the Jeans-resolved ISM to be properly captured. Gas at even higher densities is still present in our simulations but should be considered as unresolved and be analyzed with this caveat in mind. The dotted lines in Fig. 2 indicate in which density regime HII regions around young O stars are resolved. The strong dependence of the Strömgren radius with density makes it prohibitively expensive to properly resolve the early most embedded phase of every HII region. We discuss this in further details in the caveats section (Sect. 6).

Fig. 1 Face-on view of the simulated Milky Way analog. The grayscale shows the total gas column density. The highlighted areas correspond to the enhanced-resolution regions of the simulation.

2.4 Chemistry and thermal treatment

We explicitly follow hydrogen (Glover & Mac Low 2007a,b) and CO (Nelson & Langer 1997) chemistry by using a nonequilibrium time-dependent chemical network. For details of the network see Glover & Clark (2012, their network NL97). We opted for this rather simplified CO treatment, sacrificing accuracy in CO abundances for better performance. This is a valid compromise, as even in our high-resolution regions, the scales necessary to have converged CO abundances are mostly unresolved (Joshi et al. 2019), rendering the use of a more accurate network pointless. Still, even with this simplified network the thermodynamic behavior of the gas is excellently captured, and only minor differences are seen in the dynamics and properties of the ISM and SF at these scales compared to more expensive and accurate networks.

Radiative heating and cooling of the gas is followed using a detailed atomic and molecular cooling function, based originally on work by Glover et al. (2010) and Glover & Clark (2012) and described most recently in Mackey et al. (2019). At low temperatures, the cooling is dominated by the ground-state fine structure transitions of C⁺ and O and by rotational emission from CO molecules. The contributions of these processes are computed using the nonequilibrium abundances of these species provided by our chemical network. At higher temperatures (T > 10⁴ K), permitted atomic transitions dominated. We include the effects of atomic hydrogen cooling using the nonequilibrium abundance of H⁺ provided by our chemical network, but compute the cooling from He and metals in this temperature regime using the values tabulated by Gnat & Ferland (2012) that assume collisional ionization equilibrium.

We assume the ISM to be at solar metallicity, and no chemical enrichment by feedback is considered throughout the simulation. For numerical stability reasons we set a temperature floor of T_floor = 10 K. The same network has been employed in similar context for previous studies (Sormani et al. 2018; Tress et al. 2020, 2024).

Fig. 2 Resolution as a function of density. Distribution of the spherical radius (top panel) and mass (bottom panel) of the cells as a function of density. The cells in the high- and low-resolution regions are shown in red and blue, respectively. The contours contain 100, 99, 90, 75, 50, and 25% of the total number of cells. The total number of cells in each region is indicated (bottom left corner in the bottom panel). We highlight the target mass in the two regions, as well as the maximum and minimum volume limits we imposed. In the bottom panel, we show with a gray overlay the region for which gas cells are more massive than one and eight times the Jeans mass. This gas is Jeans-unresolved and potentially star forming. The dotted lines indicate the demarcation lines for where the Strömgren radius of a 15 and 120 M⊙ star (which is the mass range of O-type stars) is resolved.

2.5 Star formation

To model star formation, we employ a stochastic star particle approach, widely used in galaxy evolution research (e.g. Kretschmer & Teyssier 2020; Smith et al. 2021). Here, gas particles are stochastically converted into collisionless stellar particles based on the local properties of the gas. Contrary to sink particles, which are more commonly used in problems where local centers of collapse within the dense ISM are well resolved, these star particles interact only through gravity with the gas, and cannot accrete after they have formed (see Kang et al. 2025, for a comparison between sink and star particles). Given our resolution (see Sect. 2.3), the simulations are on the verge of a regime where sink particles are a valid tool. However, we still opted for the star particle approach, which is better suited for the low-resolution regions and provides a good approximation in the higher resolution ones, enabling a seamless and coherent transition between the two.

We used the implementation described in Göller et al. (2025). Star formation is initiated when the local collapse cannot be comfortably resolved anymore based on the resolution of the simulation, i.e. when M_cell > M_J/8 (Truelove et al. 1997, for the equivalent MHD resolution criterion, see e.g. Heitsch et al. 2001 or Federrath et al. 2011). At the current fiducial resolution we are able to resolve centers of collapse within GMCs, therefore we are more stringent in our ISM requirements to differentiate between simply dense gas and actual collapsing regions. We require that ∇ · v < 0 and that the gas cell has to be at a local minimum of the gravitational potential.

Gas cells satisfying these criteria are stochastically converted into star particles with an efficiency per free-fall time of e_ff = 0.1. As stated in Sect. 2.3, at densities where M_cell > M_J we use e_ff = 1 instead.

The star particles formed by this method do not represent individual stars; rather, they are considered to be small stellar associations with masses ranging from 10-1000M_⊙. We therefore use the mass of the star particle to perform an Initial Mass Function (IMF) sampling and populate the particle with stars on a subgrid level. We assume a Kroupa (2001) IMF throughout the simulation domain. In the CMZ there are observational and theoretical hints for a top-heavy IMF (e.g. Hosek et al. 2019; Chabrier & Dumond 2024), so an interesting future research direction would be to investigate whether the resulting difference in feedback would have a noticeable effect on gas dynamics in the Galactic Center. We use the IMF sampling scheme of Sormani et al. (2017), which by design ensures that the IMF is always reproduced, regardless of the mass distribution of the particles to be sampled. This leads to the case where sometimes the drawn stellar mass is larger than the mass of the star particle itself. However, the drawn stars are only used to determine the feedback output from the particle, such that this prescription does not compromise mass conservation. Moreover, this only happens for a small fraction of stellar particles (Fig. 3) (see Hirai et al. 2021, for an alternative treatment to prevent such cases). We only keep information about the mass of the O and B stars drawn from the IMF (i.e. stars with M > 8 M_⊙), since those are the ones that are relevant for the photoionization and supernovae (see Section 2.6 below). These are the most important feedback processes at the scales we study (Gatto et al. 2017; Rathjen et al. 2021).

Throughout the paper, whenever we use the term stellar or star particle, we refer to the AREPO particle, which does not necessarily represent an individual star. The term star is instead used to indicate the individual stars drawn by the IMF sampling and assigned to each stellar particle.

2.6 Stellar feedback

Star particles are sources of feedback based on the number and mass of massive stars assigned to each of them.

2.6.1 Ionizing radiation feedback

Each star particle is a potential source of ionizing radiation feedback. Each OB star drawn from the IMF and assigned to the star particle contributes to the radiation flux of the particle for the duration of its main sequence lifetime. The radiation flux is determined by assuming each star to be a black body of a given radius and surface temperature, by taking the stellar Zero Age Main Sequence (ZAMS) values from Ekström et al. (2012). The radiation for each star is kept constant for the lifetime of the star.

Under these assumptions, in Fig. 4 we show how much stars of a given mass contribute to the total photon budget of ionizing radiation given our assumed IMF. Less than 10% of ionizing photons come from B stars despite their longer lifetimes and the fact that more than 60% of massive stars are B stars. We still consider the contribution of B stars to the ionizing radiation feedback. However, throughout the analysis we only focus on O stars and their surroundings.

The radiation is then transported through the domain by solving the radiative transfer equations using SWEEP (Peter et al. 2023), which is a discrete ordinates method coupled to AREPO, assuming infinite speed of light approximation under steady-state-conditions. Here the equations are discretized in time, position, frequency and angular direction. For each angular direction, the cells are sorted to determine a dependency graph which is then solved using a task-based parallelization approach. This is swiftly achieved by exploiting the interconnectivity of AREPO cells through the Delaunay triangulation, which makes looping over neighboring cells trivial. Apart for a post-processing step necessary to the analysis in Section 3.3, we followed a single frequency bin for photons with energy above 13.6 eV.

Transport sweeps are performed on every AREPO synchronization timestep, when the Delaunay triangulation for the entire domain is reconstructed. This is when photon fluxes for each cell are updated. These are then passed to the chemistry routines which compute temperature and ionization state of each gas particle. In between synchronization timesteps photon fluxes are kept constant. We forced the code to perform one synchronization timestep at least once every 4 kyr such that the radiation is updated on a high enough cadence, and we were so able to capture variation on timescales relevant to the evolution of HII regions (e.g. Maillard et al. 2021).

2.6.2 Diffuse interstellar radiation field

The background interstellar radiation field at lower frequencies (<13.6 eV), which is relevant for instance for molecular hydrogen photo-dissociation and photoelectric heating, was assumed to have a constant background value depending on the galactocentric radius. Locally, for each gas cell, we computed anisotropic attenuations based on the surrounding (<30 pc¹) gas by computing column densities along 48 different directions with the TreeCol algorithm (Clark et al. 2012). To estimate the radial dependence, we assumed that its profile is shaped by the younger stellar populations. Similar to what was done by Smith et al. (2023), we computed the time-averaged star formation rate surface density as a function of galactocentric radius of a low-resolution simulation of the system, and used a rescaled version of this for the Interstellar Radiation Field (ISRF). The profile used is shown in Fig. 5. The same radial dependence is then also used for the cosmic-ray ionization rate.

Fig. 3 Properties of the star particles in the high-resolution regions. We show the distribution of the mass of the star particles (first column), the density of the gas particle that generated the star particle (second column), the ionizing radiation flux assigned to the stellar particle (third column), the total mass of stars drawn from the IMF (fourth column) and the number of stars drawn (last column). Often, no massive star is drawn from the IMF sampling. The mass and parent density distribution of those star particles is shown in the upper right.

Fig. 4 Contribution to the total ionizing photon budget of massive (>8 M⊙) stars given their mass weighted with the IMF. The dotted line shows the cumulative fraction of massive stars, the solid line the cumulative fraction of ionizing photons at birth of the star, while the dashed line shows the cumulative fraction of the ionizing photons generated over the entire life of the star as a function of mass.

2.6.3 Supernovae

We included supernova (SN) feedback from both type II and type Ia. The former are coupled to stellar particles, and inject energy at their current location into the ISM at the end of the main sequence lifetime of their associated massive stars. These naturally generate clustered SN feedback, as stars form in associations. Type Ia supernovae instead are associated with older stellar populations and therefore distributed more randomly. We assigned them the same distribution as the density profiles of the stellar components used to generate the background potential in which the gas evolves. They were injected with a uniform distribution in time, with a rate of 1 /250 yr⁻¹.

The energy is inserted into the nearest neighbors to the cell closest to the location of the SN event, which are accessed by looping through the Delaunay connections of this cell (Smith et al. 2018). If the adiabatic Sedov-Taylor phase is resolved (Blondin et al. 1998), a total of 10⁵¹ erg in the form of thermal energy is injected. The hydro solver is then able to self-consistently evolve this energy into a blast wave, converting it partly to kinetic energy which will be deposited into the ISM, while still contributing to the hot T ≳ 10⁶ K ISM phase. On the other hand, if the Sedov-Taylor phase is not resolved, we estimate the energy of the SN remnant at the end of the momentum conserving snow-plow phase and inject this energy in the form of kinetic energy (Gatto et al. 2017). In this case, the SN event will not contribute to the hot phase, but it will deposit the right amount of kinetic energy into the ISM. This energy would have been un-physically radiated away if it would have been injected in the form of thermal energy (e.g. Kim & Ostriker 2015). Most of the SNe in the high-resolution regions (Section 2.3) are resolved, and ≲10% require momentum injection.

To ensure that the energy is always injected into active cells, we modified the time-step criteria such that cells within a 50 pc radius of the star particle have to be active when a SN event is scheduled. For stability reasons, during injection we limit the maximum temperature to 5 × 10⁷ K. Except during the setup phase (Section 2.7) where we want to avoid gas depletion, each SN event only injects energy, and no mass.

Fig. 5 Background interstellar UV radiation field (photons with energy <13.6 eV) as a function of galactocentric radius. The value is scaled to the solar neighborhood value G0 (Draine 1978). The cosmic-ray ionization rate follows the same profile, scaled to the solar neighborhood value.

2.7 Initial conditions

We used the same initial density profile as in the simulations of Tress et al. (2024). The disk was truncated at 10 kpc as we are not interested in the outer Galaxy in this study. The magnetic field was initialized to a purely vertical component of intensity |B₀| = B_z = 0.02 μG. Turbulence and galactic rotation will bring this initial field closer to an observed configuration with saturation levels around 10 μG at the solar circle. To generate turbulence in the disk and avoid violent and coherent star formation cycles induced by the initial collapse of the smooth disk, we performed a pre-processing step. Here we reduced the lifetime of stars by a factor of five such that the feedback coupling to the ISM happens on shorter timescales and an equilibrium state can be achieved more smoothly. During this phase, the only source of feedback is in the form of SNe, and stars reinject their mass fully into the ISM at the end of their lifetime, so that the gas is not depleted during this stage. This is performed at a low mass resolution of 500 M_⊙ over a period of 100 Myr. During this time, the background potential is kept in its axisymmetric form. The next 150 Myr were devoted to the smooth introduction of the Galactic bar. The non-axisymmetric version of the barred potential is introduced linearly over this duration by keeping the total mass constant, leading to a steady formation of the bar lanes and CMZ (Athanassoula 1992). Stellar lifetimes were slowly returned to their nominal values over the first 100 Myr of this phase. Finally, 20 Myr before the Bar is fully on, we increased the resolution to the reference values described in Sect. 2.3 and introduced the ionization feedback. After the Bar is fully on, we turned off the 100% gas re-injection, and reset the mass of the central SMBH to its value defined in Section 2.2 as it could have experienced excessive accretion events during this unphysical buildup phase. The simulation was then advanced in a production phase for ∼110 Myr. This is enough integration time to have several independent star formation cycles and multiple rotation periods of the gas on x2 orbits in the CMZ, such that our results are statistically significant. We denote the time when the bar is fully on as t = 0 Myr and restrict our analysis to t > 0 Myr. During the production phase, we produced one full snapshot every 0.25 Myr. This high cadence is necessary to follow the fast evolution of the early star formation phase, in particular in the dynamic environment of the CMZ.

3 Comparison of the ISM properties of the CMZ and the solar circle

In this section, we explore in detail the properties of the simulated ISM. We emphasize the differences between the two high-resolution regions to understand and quantify the effect of the special environment of the CMZ on the ISM and SF.

3.1 Morphology

Compared to our previous modeling of the CMZ (Tress et al. 2020, 2024), this simulation is more advanced in terms of ISM and feedback processes included. The general large scale morphology, however, is not radically affected by this and exhibits the same features seen in simulations with simpler prescriptions. The Galactic bar dominates the gas dynamics in the central ∼3 kpc, inducing strong non-axisymmetric motions. As the bar redistributes angular momentum, the gas is driven inward and flows almost radially along so-called bar lanes toward the center. The gas accumulates in a ring-like structure about 300 pc in radius, the CMZ of this simulation. Here, torques from the bar are unable to further remove angular momentum (see Sormani et al. 2024, for an explanation of what sets the size of the CMZ), and other mechanisms are at play to drive gas further inward (Shlosman et al. 1989; Tress et al. 2024). The size of our simulated CMZ appears slightly larger than what is observed (estimated around 100 pc in radius Kruijssen et al. 2015), although there is some ambiguity in where exactly the CMZ ring ends (Section 4.3.2.2 of Henshaw et al. 2023). This is a difficult parameter to fine-tune, as it depends on the included ISM physics and feedback, as well as on resolution, so it can only be measured a posteriori once the high-resolution simulation is completed. Future iterations will account for this fact and adjust for it. The background potential also includes a fourarmed spiral pattern (Section 2.2), and in the disk we indeed see spiral structures develop throughout the simulation, induced by the combined effect of this spiral potential perturbation and by the rigidly rotating bar (see Göller et al. 2025, where the impact of the external potential on the morphology is studied in more details).

In Fig. 6, we show a closeup of the CMZ (at the top) and a region of the same size of 2 kpc across in the disk, at the solar circle (bottom) at t = 15 Myr. Without the loss of generality, we chose this simulation time as our fiducial snapshot. Many of the following figures refer to this snapshot. The ISM has settled in an equilibrium state through several SF cycles. Inspecting the bottom panel, we see dense molecular clouds forming, embedded in cold atomic envelopes. These are the formation sites of new SF events, assembling stellar clusters of different sizes and masses. The ionizing radiation associated with these young stellar clusters generate HII regions of different sizes and degrees of compactness. The combined feedback of these clusters is disrupting the molecular complexes from within, carving out hot expanding bubbles into the disk. The oldest clusters are found at the center of the largest bubbles, powering their expansion through SN feedback. This mechanism drives and maintains turbulence in the ISM, sustains the extended scale-height, and drives atomic and ionized outflows, which will fall back onto the disk in a fountain flow, powering future SF cycles.

The morphology of the CMZ here confirms the general ringlike structure identified in previous simulations. Tress et al. (2024) and Sormani et al. (2022a) found that both magnetic instabilities and supernova feedback are able to transport gas from the x₂ ring toward the center. This model now includes both magnetic fields and stellar feedback, and indeed the same steady inflow is observed throughout the simulation. In this sense, the inner region of the CMZ is relatively devoid of gas in the initial stage of the simulation, and is gradually filled with gas at later times, similar to what is shown in Tress et al. (2024, their figure 6). It is interesting to investigate how the two processes act in unison to transport gas toward Sgr A*, whether the two effects compound on each other or mutually cancel each other out. We defer this analysis to some future study.

The CMZ structure is much more disturbed during the initial phases of the simulation, when the accretion from the bar lanes is still prominent and highly variable (top panel of Fig. 7). This perturbs the CMZ gas on x₂ orbits, and might cause its asymmetry (Sormani et al. 2018). At later times in the simulation, the inflow is much smaller and steadier, resulting in a more symmetric appearance of the CMZ. This explains the more disk-like structure of the CMZ toward the end of the simulation. The real CMZ exhibits a highly asymmetric gas distribution (Bally et al. 1988; Jenkins & Binney 1994; Sormani et al. 2018; Battersby et al. 2025a), and is likely to still experience strong accretion events from the bar lanes (Sormani & Barnes 2019; Su et al. 2025; although lower values are reported by Gramze et al. 2023, their Table 2), such that an early snapshot is probably closer to what is observed in the Milky Way in terms of gas distribution.

By a visual comparison of the panels of the two environments in Fig. 6, we can anticipate that the CMZ has much more extreme conditions. High surface densities in the CMZ ring lead to intense star formation in a much smaller region. In the disk, for the same area, we see a handful of star-forming regions and young stellar clusters. In the CMZ, on the other hand, most of the stars form on these x₂ orbits, generating a ring of young stars and high stellar densities. Most of the clusters are sheared apart in the dynamic environment of the inner Galaxy, leading to streams of coeval stars (Age_StarP panel of Fig. 6). Only few clusters survive and stay compact (Σ_starp panel). The resulting HII regions are generally much smaller in the CMZ, but the high density of young stars also substantially increases the diffuse contribution of ionized gas compared to the disk. At the solar circle, stellar clusters are able to drive extended superbubbles that have sizes comparable to the entire CMZ. Instead, in the Galactic Center we do see feedback-driven bubbles in the ISM, but of considerably smaller size, and which are not able to significantly affect the overall structure of the CMZ.

In the CMZ a lot of complexity is happening in a very small region. It is extremely important for a simulation to be able to capture these complex dynamics. Inability to do so might result in widely over estimation of inflow rates (e.g. Appendix B of Tress et al. 2024), BH accretion, star formation rates and major differences in size and morphology of the CMZ.

3.2 Star formation

In Fig. 7 (middle panel), we show the star formation rate (SFR) in the CMZ region as a function of time. Compared to the simulations of Tress et al. (2020) where only SN feedback was considered, the SFR is more variable in time (Sormani et al. 2020, their Figure 2). Still, the variation is contained within one order of magnitude in a time span of 100 Myr and the period of intense star formation (≳0.3 M_⊙ yr⁻¹) is correlated to the initial phase of intenser inflow rates (top panel of Fig. 7), suggesting that in our simulation, the SFR of the CMZ is regulated by the inflow rate onto the central region (like in the simulations of Moon et al. 2022). In contrast with the simulation of Armillotta et al. (2019), we do not see episodes where the feedback is able to significantly disrupt the gas in the CMZ after a period of intense star formation, completely shutting star formation. Our simulations, therefore, favor the inflow regulated star formation scenario over the episodic feedback regulated one. The second part of our study (Sections 4 and 5) is dedicated to understand how feedback is coupled to the ISM in the CMZ and provides a possible explanation of why star formation is inflow-regulated there. At later times in the simulation the SFR is comparable to observed values. Most of the star formation is contained in the CMZ ring, and only very few stars form in the region between 0.5-2.5 kpc where the gas dynamics are dominated by the bar lanes (Fig. 8).

In the bottom panel of Fig. 7 we compare depletion times (τ = M_gas/SFR) in the two regions of interest. The HI depletion times is higher by an order of magnitude in the solar neighborhood than in the CMZ, while the H₂ depletion times are similar in both regions. This suggests that the widely different conditions in the CMZ makes the transition from atomic to molecular gas much more efficient. At first sight, it seems that the star formation process in the dense gas may proceed in a similar fashion in both region. However, subtle differences exist. This is emphasized in Fig. 9 where we show star formation rate surface densities (Σ_SFR) as a function of molecular gas surface densities (Σ_H2) in 100 pc bins. CMZ bins have generally higher Σ_H2 values, but have lower values of Σ_SFR compared to bins at the solar circle at similar surface densities. This indicates that in the CMZ the dense molecular gas is less efficiently converted into stars as if it would have been located at the solar circle with the same density.

3.3 ISM phases

Before comparing the different ISM phases in the two regions of interest, we first did a post-processing radiative transfer step to improve the ionization state and temperature of HII regions and better identify the photo-ionized gas phase. During this stage we only allow the radiation transfer and the thermochemistry to evolve in time performing 10 radiation sweeps over 10 kyr. Contrary to the production run, in this step we include 3 frequency bins for the ionizing radiation ((13.6-15.4) eV, (15.4-24.6) eV, (24.6-∞) eV)² instead of a single >13.6 eV one.

Phases other than the photoionized gas are largely unaffected by this.

We show the mass weighted density-temperature phase space distribution in Fig. 10 and the mass weighted probability density function (pdf) of the density (Fig. 11) and the temperature (Fig. 12) in the two regions. The density pdf in both cases roughly has a log-normal profile with a high-density power-law tail in molecular gas (e.g. Brunt 2015), and is shifted to higher densities in the case of the CMZ compared to the solar circle. The total density pdf in the CMZ peaks at n ≃ 1.3 × 10² cm⁻³, and at the solar circle at n ≃ 6 × 10⁻¹ cm⁻³. The shift toward higher density is present even for the distribution of ionized, atomic and molecular gas taken individually.

The ISM is then divided into the following phases (Draine 2011):

• Cold Molecular Medium (CMM) when x_H2 > 0.25;

• Cold Neutral Medium (CNM) when x_H > 0.5 and T < 500 K;

• Unstable Neutral Medium (UNM) when x_H > 0.5 and 500K < T < 6000 K;

• Warm Neutral Medium (WNM) for x_H > 0.5 and T > 6000 K;

• Warm Ionized Medium (WIM) when T < 3.5 × 10⁴ K and x_H⁺ > 0.5;

• Hot Ionized Medium (HIM) when T > 3.5 × 10⁴ K.

In Fig. 13, we show the mass fraction in these different phases for the two regions for a snapshot of the simulation at t = 15 Myr. A larger fraction of the ISM is molecular in the CMZ compared to the solar circle, and the majority of the gas is at low temperatures, either in the CNM or CMM. Still, a significant amount of gas is actually in the form of atomic hydrogen. The WIM is gas which is photo-ionized, either in the form of diffuse ionized gas, or in HII regions, while the HIM is mostly generated by SNe. Both the CMZ and the solar circle have comparable fractions of gas in the ionized phases (WIM and HIM) despite the much higher star formation density in the CMZ which accounts for higher feedback.

Fig. 6 Projections of the CMZ (top) and a region of the disk (bottom). The panels from left to right and top to bottom depict the total, atomic, molecular, ionized gas column densities, mass-weighted thermal pressure, temperature and magnetic field intensity averages, star formation rate surface density, stellar surface density and star particles colored by their age. These last three quantities are overlaid on the total gas column density in gray.

Fig. 7 Gas inflow/outflow to the CMZ (top), star formation rate in the CMZ (middle), and depletion times for the CMZ and solar circle (bottom) as a function of simulation time. Inflows and outflows are shown separately and the different gas phases stacked in different color to show the contribution of each phase to the total gas inflow/outflow shown as the solid purple line. Depletion times for the atomic, molecular and total gas are shown separately. The observed SFR value for the CMZ and its uncertainty is shown as the dashed line and shaded region in the middle panel.

Fig. 8 Cumulative star formation as a function of galactocentric radius at t = 15 Myr.

Fig. 9 Star formation rate surface density as a function of molecular gas surface density. The values are computed over a 100 pc sized bins and shown separately in different colors for the CMZ (teal) and the solar circle (orange). The values are taken for several snapshots throughout the simulation. The colored area shows the kernel density estimation of the distributions and the black lines show lines of constant depletion time.

Fig. 10 Density-temperature phase diagram of the gas in the simulation at t = 15 Myr. The colors indicate the gas mass fraction in the given (n – T) bin normalized by the bin size (dMgas(n, T)/Mgas,tot/d(log10 T)/d(log10 n)). The gas in the CMZ is shown in the left panel, the gas at the solar circle in the middle panel, and the difference on the right.

Fig. 11 Number density probability distribution function for the CMZ (top) and the solar circle (bottom) at t = 15 Myr. The contribution of the different ISM phases is shown in different colors. Dotted vertical lines indicate average density values for the different phases.

Fig. 12 Same as Fig. 11, but for the temperature.

3.4 Magnetic field

Due to flux freezing, magnetic fields are generally higher in regions of higher density. As the average density in the Galactic Center is much larger than the solar circle, it is natural to expect higher magnetic field averages as well. We measured magnetic fields in the high-resolution regions averaged over the simulation time when the Bar is fully formed. The volume-weighted average magnetic field in the CMZ is 11.9 ± 1.4 μG compared to 2.7 ± 0.3 μG for the solar circle. If we compute the mass-weighted average instead, we get values of 131.4 ± 33.6 μG for the CMZ, and 10.2 ± 2.2 μG for the solar circle.

But magnetic fields in the CMZ are not only higher due to a larger fraction of gas being at higher densities. Even in the same density regime, magnetic fields seem to be higher in the CMZ environment. This can be seen in Fig. 14, where the average magnetic field intensity as a function of density is shown for the solar circle and the CMZ in the simulation. As the gas in the Galactic Center is more turbulent compared to the disk (see Sect. 3.5), and orbital times are much shorter, we can have a more efficient action of small-scale dynamo and large-scale mean field dynamo, increasing the magnetic field. This can also be visualized in terms of energy equipartition, which implies higher magnetic fields for regions of higher kinetic energy in turbulent motions. Both environments largely follow a relation B ∝ n^0.5, consistent with what is reported by Whitworth et al. (2025) for observations and simulations (see also Crutcher 2012). The difference in the exponent between the solar circle and the CMZ for intermediate density regimes might be explained by different strengths and modality of the turbulent driving mechanism (Sur et al. 2012) in the two regions.

The difference in B as a function of density for the two regions might not seem significant, but folding the information of Fig. 14 with the density distribution in Fig. 11, we notice that most of the gas in the CMZ sits at densities where the difference is largest. Moreover, the value of B in the CMZ consistently lies above the solar circle value throughout the different density regimes, which accounts for a significant total energy difference in magnetic fields for the two regions.

Our results of a stronger field in the CMZ than in the solar neighborhood are in line with observational findings by Akshaya & Hoang (2023, 2024). Using multiwavelength dust polarization studies, they found plane-of-sky strengths of 2-7 mG in the circumnuclear disk and mini-spiral, with pronounced spatial variations. Their finding of particularly weak fields along the Eastern Arm of the mini-spiral, where the Alfven Mach number exceeds unity, resonates with our picture of local variations in feedback coupling and turbulence driving.

Fig. 13 Mass fraction in different ISM phases for the two highresolution regions at t = 15 Myr.

Fig. 14 Average (dashed line) and median (solid line) magnetic field intensity as a function of density for gas in the CMZ (teal) and at the solar circle (orange) for the snapshot at t = 15 Myr. The shaded area around the mean indicates ±σ dispersion.

3.5 Turbulence

To estimate the turbulence in our simulated ISM, we examined the velocity dispersion of the gas. In particular, in the Galactic Center, the velocity dispersion can be quite asymmetric as the in-plane motions are dominated by the strong shear induced by the differential galactic rotation. To avoid accounting for these coherent shearing motions, we only computed the z component of the velocity dispersion, which instead is not directly affected by this.

In Fig. 15, we show a σ_z map of the CMZ and an equally sized region at the solar circle at t = 15 Myr. The velocity dispersion is computed in bins 20 pc in size. The higher degree of turbulence becomes evident. To quantify this further, we provide in Fig. 16 the (spatially-) averaged value of σ_z as a function of density for the two regions of interest. Especially for higher density regimes the gas in the CMZ is significantly more turbulent with velocity dispersions averaging around ∼9 km/s compared to values around 2 km/s for the solar circle. The reason for this is more energetic stellar feedback due to higher star formation rate surface densities in the Galactic center, stronger galactic shear cascading down to turbulent motions, more frequent cloud collisions in the dense and dynamic environment, violent gas accretion from the dust lanes and magnetic instabilities whose nonlinear regime manifests as turbulence. In addition, the feedback of wandering young massive stars could also be an important source of turbulence in the CMZ (Tassis & Pavlidou 2022).

3.6 Summary: Differences between the solar circle and the CMZ

In Table 1, we report properties of the ISM in the CMZ and solar circle regions, time-averaged over 0 < t < 110 Myr.

To retrieve total, atomic, molecular, ionized gas column densities and SFR surface densities, we first computed the profile of these properties as a function of galactocentric radius for each snapshot (using a 30 pc radial bin). The solar circle value for each snapshot is then taken by averaging these profiles over 7.5 kpc < R < 8.5 kpc. In the case of the CMZ, we took the maximum value of the profile within R < 0.5 kpc instead. This is to get a better estimate of the actual column density encountered in the CMZ ring, as its size varies slightly over the course of the simulation, and averaging over the entire R < 0.5 kpc would arbitrarily dilute these values.

Velocity dispersion and magnetic field intensities for each snapshot are obtained by computing the mass-weighted average of these values over each region.

To infer the SN-rate surface density, we logged the position of each SN over 5 Myr bins. We account for the galactic rotation by rotating the SN positions by θ = Ω(R)(t – t_0,bin), where Ω(R) = v_rot/R is the angular rotation velocity of the galaxy. For the SNe in the solar circle we take Ω(R = 8 kpc) and Ω(R = 0.2 kpc) for the CMZ. This ensures that the SN-rate surface density does not artificially decline simply because the parent cluster moves along its Galactic orbit. However, shearing forces, disrupting the cluster, will decrease this values. The reported value in Table 1 is then obtained by time-averaging all the values obtained over these 5 Myr time intervals.

Finally the mean ionization parameter was obtained by computing a photon number weighted average over each region of the ionization parameter U_i = φ_i/(n_H,tot,ic) where φ is the photon flux and n_H,tot the nucleon number density of each cell. Again, the reported value is then the time-average of these values.

Fig. 15 Velocity dispersion (σz) map for the CMZ and a section of the disk at t = 15 Myr. Values are obtained by inferring the mass-weighted velocity dispersion of cells in 20 pc bins. The contours show column densities of [10,50] cm−3. Lower density bins are excluded for clarity.

Fig. 16 Average (dashed line) and median (solid line) z component of the velocity dispersion as a function of density for gas in the CMZ (teal) and at the solar circle (orange) at t = 15 Myr. The shaded area around the mean indicates the standard deviation.

4 Environment around O stars

After highlighting the differences between the two environments, we focus now on the question whether these differences can affect stellar feedback coupling with the surrounding gas. To investigate this, we more closely inspect the environment around young O stars. These are the main source of feedback and are central for the early disruption of the parental clouds. How effective are they in doing so in our simulation and in particular in the two different environments?

Fig. 17 Schematic representation of the assembly and disruption of a GMC. Stars form as the cloud is assembling (panels a-c), until their combined feedback is strong enough to disrupt the cloud from within (panel d).

4.1 Life and death of a GMC at the solar circle

A classical star formation site in a dynamic environment would follow a time evolution as illustrated schematically in Fig. 17 (e.g. Geen et al. 2016; Grudić et al. 2021 ; Kim et al. 2021 ; Brucy et al. 2025). Stars form in the dense parts of a GMC, which itself forms due to several instabilities in the gas-rich disk of a Milky Way-type galaxy. As soon as stars start to form, massive stars will influence the surrounding ISM by strong ionization feedback and stellar winds³. However, a single star is unlikely to be able to disperse the entire dense surrounding molecular gas. Moreover, in a dynamic and interacting environment, more gas might be flowing toward the star formation site and more molecular gas might be condensing from the atomic phase, which leads to more star formation. Stars, therefore, never form alone, and the formation phase in a GMC complex is not instantaneous, but might be happening over timescales of a few Myr. In this initial phase, stars might be partially or totally embedded in the dense cloud and the HII region is ionization bounded (no ionizing radiation is leaving this region). Sooner or later, however, enough stars might have formed such that the combined feedback is able to disperse the surrounding gas, and shut further star formation off in the region. The cluster is now fully exposed, and the HII region becomes density bounded, allowing a large fraction of the ionizing photons to escape.

In Fig. 18 (top), we show the evolution of an example star-forming region at the solar circle as experienced from the position of an O star associated with the region. It essentially follows the story sketched out in Fig. 17, as the GMC structure is slowly being disrupted from within by stellar feedback. More quantitatively, in Fig. 19 (top left panel), we look at the spherically-averaged density profile around this O star as a function of time. The progressive effect of internal feedback disrupting the cloud around the star and expanding the ionized bubble becomes evident. The left panel of the second row of Fig. 19 shows the structure of the HII bubble around the star with time. We see that the new-born star is fully embedded and the HII region is initially small, with size ≲10 pc. It remains ionization bounded until 1-2 Myr when the radiation is able to carve out small windows and a small fraction of photons is able to escape to large distances. Slowly, the HII region tends to widen, as the windows’ size increases and a progressively larger fraction of photons escape, until the stars are fully exposed. Overall, the evolution is very monotonous. At a given radius from the star, e.g. 10 pc, the density is smoothly decreasing with time, and the fraction of ionized gas is increasing (left column of Fig. 19).

When averaging over all O stars at the solar circle (middle panel of the two first rows of Fig. 19), the same progressive evolution emerges. Only less massive O stars (M ≲ 20 M_⊙) survive until 10 Myr, and those have a much lower energy contribution (Fig. 4), which is why the average ionization fraction is less than one at 10 Myr. Those low mass O stars are struggling to fully disrupt their high density surrounding.

4.2 The case of the CMZ

In order to know whether the extreme environment of the CMZ still allows for the classical GMC evolution sketched in Fig. 17, we follow O stars in the CMZ and explore their surrounding.

The CMZ being a much higher surface-density environment, we expect ionizing radiation there to be less efficient at carving out a low density cavity, disrupting the parental cloud from within. In Fig. 18 (bottom) we have selected an O star in the CMZ, and an observer at its position indeed experiences higher surface densities, a generally smaller HII region, and smaller windows where the ionizing radiation can escape.

On the other hand, due to the higher star formation rate surface density, we also expect more massive young clusters, whose combined feedback action could disrupt the high-density environment. In practice, the short dynamical timescale in the CMZ is comparable to the feedback timescale, meaning that fresh dense molecular gas can easily flow into the star formation site, re-embedding the O star. This is the picture that we see by following the early lifetime of an O star in the CMZ, shown in Fig. 18 (bottom). Initially, the strong ionizing radiation from the cluster is able to expose the star, and carve out an HII region. Yet, after some time, owing to the dense and dynamic nature of the surroundings, new molecular material is flowing in, re-embedding the star and closing the windows for radiation to escape. At later stages yet, the star can become re-exposed once more.

If we look at the density profile and structure of the HII region surrounding this randomly selected O star (Fig. 19, left panel of the two bottom rows), we can spot this episodic nature as well. The star is born in a high-density environment, is able to disrupt its surrounding within 10-20 pc over their first ∼4 Myr lifetime, only to be re-embedded fully at around 5-7 Myr. At an age of 10 Myr the stellar particle is exposed again.

Looking at the averaged density and ionization profiles of all stellar particles holding at least one O star in the CMZ (Fig. 19, middle column), we are not able to see this episodic behavior as the time evolution of each star is different and will average out (see the right column of the same figure). However, we learn that on average, the background density at 50 pc of HII regions in the CMZ is around ∼30 cm⁻³, much higher than at the solar circle. Moreover, on average the HII region is able to ionize only a small fraction of their surrounding, ionizing at most ∼30% of a surface with r = 50 pc. This is due to the higher surface density environment, but also because of the large fraction of stars undergoing these episodic high-density embedding events. We see that at least on average the star is still able to clear its high density birth environment over its initial 4-6 Myr, but densities increase again and HII regions become more embedded at later times. The reason can be found in the higher cadence of re-embedding events for less massive O stars at later ages in a dynamic high-surface density environment, as the more massive and energetically radiating O stars died off earlier.

To get a more quantitative comparison, we classify stars based on the evolution of their surrounding environment. We use the ionized surface fraction f_Sion(r) = S_ion(r)/(4πr²) at a distance of r = 10 pc, where S_ion(r) is the surface around the star at a distance r which the star contributes at keeping ionized, i.e. when a ray from the star only passes through ionized gas (x_H⁺ > 0.5). The star is categorized as always embedded if the feedback is unable to expose the star, and f_Sion(10pc) < 0.5 throughout their lifetime up to a maximum of t = 10 Myr. On the other hand, we classify star particles as exposed if at some point fs_ion(10pc) > 0.5 and it remains that way for the rest of their lifetime up to 10 Myr. Finally, we label stellar particles as re-embedded if at some point, after being exposed with f_Sion(10pc) > 0.5, the star enters a high-density environment and is strongly embedded again, with f_Sion(10pc) < 0.2.

In Fig. 20, we show the fraction of star particles in each category with at least one O star, for the CMZ and the solar circle separately. We see that most stars in the solar circle end up being exposed without experiencing any re-embedding event, which is consistent with the scenario sketched of Fig. 17⁴. In contrast, most of the stars in the CMZ experience at least one reembedding event. This implies that the classical GMC evolution scenario does not apply in the CMZ and a new model is needed. Only a small fraction of stars remain fully embedded throughout their lifetime in both environments, which is a confirmation that early ionization feedback is effective in generating low-density HII regions.

Fig. 18 Time evolution of the environment surrounding an O star at the solar circle (top three rows) and in the CMZ (bottom three rows). In each panel, the first row shows the H2 density slice surrounding the selected star, which is centered in the box. Iso-density contours are shown at a density of 100 cm−3. The inlay panels shows the H2 column density of the larger-scale context. Surrounding O stars are shown as white dots in the inlay and turquoise stars in the main panel. We highlight the age of the selected star particle, the total number of OB stars assigned to it, and the mass of the most massive one among those. The magenta area is the region visible by the ionizing radiation of the star and which the star contributes to keeping ionized, i.e. the region of space for which rays connecting to the star only pass through ionized gas (xH+ > 0.5). The second and third rows show the H2 column density and size of the HII region, as seen by an observer at the position of the star. The second row is the Mollweide H2 column density projection. The third row shows the size of the HII region surrounding the star. As defined here, the size is the distance in each direction at which the fractional ionization of the gas first drops below 50% (i.e. xH+ < 0.5). Note that the panels are not equally spaced in time, but instead a few salient snapshots in the evolution of the star are selected.

Fig. 19 Density and ionization profile around O stars. The top two rows are the case for the solar circle, and the bottom two rows for the CMZ. For each row, the first column is an example for one selected star, while the second column is the average over all stellar particles in the region (solar circle or CMZ) that contain at least one O star. For each region, the first row shows the density as a function of distance from the O stars, averaged over spherical shells. The second row depicts the fraction of the surface that the O stars are helping to keep ionized, i.e. the fraction of a surface of a given radius centered on the star for which a ray passes only through ionized gas (xH+ > 0.5). These quantities are shown at different ages of the star particles. The last column shows the time evolution of the environment around a selection of seven O stars at a radius of 10 pc.

Fig. 20 Fraction of stars classified based on their evolution (clarified in the text) for the solar circle and the CMZ. Most of the O stars in the CMZ will undergo at least one event where the star after being exposed, will be re-embedded in a high density environment, in stark contrast to what O stars at the solar circle experience.

5 Toward a new model for feedback coupling in the CMZ

5.1 Dynamical decoupling of stars and gas

The task now is to understand and describe the causes of the differences in behavior of how these stars affect their surroundings in the two different environments. Just referring to the more extreme conditions of the CMZ is unsatisfactory as an explanation. These certainly contribute in making it harder for the ionizing feedback to efficiently disrupt the surrounding molecular gas, but what makes the CMZ so special is that the gas there follows the very specific dynamics shaped by the barred potential. The gas circles around the Galactic Center on x₂ orbits with orbital times of the order of 5-10 Myr, comparable to GMC lifetimes, so this has to be factored in in any description of star formation and cloud evolution in the region.

One major effect is that gas and stars undergo different physical processes, and while gas follows the equation of hydrodynamics and responds to pressure forces such as gas collisions and shocks, stars are collisionless and only feel gravitational forces. This leads to gas and star decoupling on times comparable to the orbital timescale which in the CMZ are short enough to matter for star formation and GMC evolution. In essence, the gas in the CMZ can only follow nearly closed x₂ orbits. Otherwise, it would collide with its own self-intersecting orbit, dissipating energy and falling back on the nearly closed orbit. Stars, on the other hand, will generally follow open Rosetta-type orbits (Nieuwmunster et al. 2024), and therefore decouple from the gas on an orbital timescale (Fig. 21). The high-density and turbulent environment exacerbates this behavior, and any cloud collision and shock will enhance the velocity and spatial difference between star and surrounding material from which the star formed. Tassis & Pavlidou (2022) estimate the stellar velocity dispersion to be large enough such that even young O stars can cross entire GMCs within their lifetime. Moreover, small differences in velocity and/or position, can quickly grow on short timescales in the environment governed by strong potential gradients, which leads to high orbital angular velocities and strong shear. Finally, the CMZ experiences violent and frequent gas accretion events through radial gas-inflow from the bar lanes. These inflow events happen on timescales relevant to star formation and GMC evolution as well, and can play an important role in shaping the dynamics of the CMZ. The new inflowing material can cover up freshly exposed stars and clusters, and colliding with the gas in the CMZ, can enhance gas-star decoupling.

In periods of more violent accretion, the CMZ can be disrupted and become asymmetric. As a consequence, stellar orbits can be more heated in phase-space, leading to more extreme and frequent decoupling between gas and star. This can be seen in the first panel of Fig. 21, where a stellar orbit is shown with respect to the averaged gas distribution at an early time of the simulation, when the accretion from the bar lanes is quite violent. The star exhibits a noticeable open orbit. At later times, when the inflow is lower and steadier, the decoupling is still present, but less extreme (last panel of Fig. 21). For clarity, we only plot one orbit, but most stars follow a similar behavior. The decoupling between gas and stars is taken to the extreme when the star formation is associated with the bar lanes as can be seen in the middle panel of Fig. 21. The gas will collide with CMZ gas, while the star will continue to follow the high-angular momentum x₁ orbits, a family of orbits elongated parallel to the bar major axis and much more elongated than the x₂ orbits. The HII regions of SgrE (Anderson et al. 2020) are an observational example of this behavior.

In essence, in the CMZ, the dynamical decoupling between star and gas happens on a very short timescale, when the star is still in its feedback stage. Because of this decoupling, the behavior of a typical young star-forming region in the CMZ is very different than the classical scenario in the solar circle illustrated in Fig. 17. We have sketched out an adapted scenario of the early evolution of a star-forming region in the CMZ in Fig. 22, accounting for the decoupling between star and gas which, in turn, triggers the re-embedding events described in Section 4.

5.2 Cluster shearing

Another major effect on early evolution of star-forming regions and stars is given by the strong shear in the CMZ region compared to the Disk. This has the effect that small initial velocity and position differences, will quickly increase in amplitude on orbital timescales. Coeval stellar associations which formed from the same cloud complex, will be sheared apart much more efficiently and on smaller timescales compared to the solar circle (Fig. 23). The dissolution of star clusters due to the Galactic tidal field has been estimated to happen on timescales of the order of just a few megayears in the Galactic Center (Portegies Zwart et al. 2002; Kruijssen et al. 2014). Indeed, in the CMZ only two young clusters are known (Arches and Quintuplet) whose masses alone cannot account for all the recent star formation in the region. This is known as the missing clusters problem. Moreover, evidence for the presence of young dispersed stellar associations in SgrB1 (∼10 Myr old, Nogueras-Lara et al. 2022) and Sgr C (∼20 Myr old, Nogueras-Lara 2024) have been reported.

In the simulation, this effect generates stellar streams of coeval stars along the simulated CMZ (Fig. 6, panel of the face-on view of the CMZ where star particles are colored by age). As in the observations, we do see a few surviving clusters of stars in the projected stellar surface density, but the bulk of the young stars are dispersed throughout the CMZ.

Fig. 21 Orbits of a few selected stellar particles containing O stars at birth. The orbit is colored by the age of the stellar particle, and each SN event in its lifetime associated with the particle is marked along its orbit. The birthplace of the star is highlighted with the blue x-marker. In the background in grayscale we show the averaged H2 density over a period of 20 Myr starting from the snapshot where the star was born. The black contour shows an average density of 50 cm−3. In the first panel we show an early snapshot when the inflow from the bar lanes is still quite copious. In the middle panel we show the orbit of a star which formed in a GMC associated with the bar lanes and inherited their velocity. On the right panel a snapshot at a later time was chosen, where the CMZ is more regular, and the inflow smaller and more steady. Stars quickly decouple from the gas on orbital timescales, which are short in the CMZ.

Fig. 22 Schematic representation of the early evolution of a starforming region in the CMZ and the subsequent evolution of the formed stars. Compared to Fig. 17, the situation is more complex in the dynamic environment of the Galactic Center. Stellar-gas decoupling (panel b) due to cloud collisions and short orbital times, leads to re-embedding events in the early evolution of the stars (panel d). This has an effect on how feedback can couple to the dense gas in the CMZ.

5.3 Implications on the efficiency of stellar feedback

Both the dynamical decoupling of the gas and stars (Section 5.1) and the shearing of the stellar clusters (Section 5.2) have a noticeable impact on how the feedback couples back to the surrounding ISM and affect how efficiently clouds are disrupted.

Classically, in the disk, ionizing radiation and SNe will disrupt the GMC from which the star was formed. However, in the CMZ, stellar decoupling and cluster shearing have the consequence that the feedback is injected into gas which is dynamically de-correlated with the star. The dense gas, left without ionizing stars, may continue its collapse and form new stars, such that these processes have an impact on the star formation efficiency of GMCs. By contrast, dense gas in other parts of the CMZ may host wandering O stars and be disrupted before being able to collapse. Tassis & Pavlidou (2022) estimate that such wandering O stars could explain the high level of turbulence as well as account for the steep size-linewidth relation in CMZ clouds.

Additionally, cluster shearing will preclude the compounding effect of clustered feedback, which is an important aspect for cloud disruption and feedback coupling (e.g. Walch et al. 2015; Gentry et al. 2019). Indeed, unless stars are formed in a sufficiently massive and compact bound cluster, more extended stellar associations will not survive long in their x₂ orbit. As a consequence, the feedback from the stars is less concentrated, and instead more spread-out over the entire CMZ and thus less efficient in disrupting the surrounding molecular gas. These stars will act as a more evenly distributed background feedback source, driving turbulence, without being able to create large superbubbles. In isolation, they are then much more easily re-embedded in the turbulent flow of the CMZ.

These special CMZ conditions will have an effect on the environment into which supernovae are injected. In the Disk, they are predominantly injected at the center of superbubbles, in progressively lower density. In the CMZ, instead, due to the previously mentioned effects, they will be more randomly distributed, such that SNe in higher density gas are more frequent (Fig. 24). This might change the way in which stellar feedback impacts the ISM in the CMZ as a whole, for instance in driving turbulence (Iffrig & Hennebelle 2015; Walch & Naab 2015; Ohlin et al. 2019).

This could be another puzzle piece in explaining how SF is regulated in the Galactic Center. Contrasting theories, each supported by simulations, state that either the SF in the CMZ is highly episodic (Krumholz et al. 2017; Armillotta et al. 2019; Nogueras-Lara et al. 2020), with phases of intense SF followed by phases where the feedback disrupts the CMZ entirely and shuts SF off. Or conversely, SF is rather constant in time, and regulated by the inflow. The solution to the conundrum reconciling the two hypotheses could lie in how stars are formed in the CMZ. In one case stars form in more dense and compact clusters, surviving the strong shearing environment and the resulting feedback-driven superbubbles are much more disruptive to the gas in the CMZ. In the other case, less compact clusters form and get easily sheared apart. The stellar feedback then provides a background source of turbulence (Tassis & Pavlidou 2022), but is unable to fully shut SF off, as the high-density gas is inefficiently being disrupted, and instead the amount of dense gas is regulated primarily by the inflow.

Stellar-gas dynamical decoupling, and cluster shearing strongly depend on the properties of the local potential, orbital times and size of the CMZ. As the CMZ is believed to be even smaller in observations compared to our simulations, it is plausible that these effects are even more strongly pronounced in the real CMZ. Observations of external CMZs (e.g. Schinnerer et al. 2023; Sun et al. 2024) could be used to corroborate this hypothesis.

Fig. 23 Young stellar association that is quickly sheared apart in the CMZ environment. We highlight new-born (<3 Myr) O stars at a given time in the simulation (left, yellow stars) and select a compact association of coeval O stars in red. We follow these stars in their evolution and show their position with respect to the background molecular gas (in grayscale) after t = 3 Myr (middle) and t = 6 Myr (right), which correspond to orbital fractions of forbit = 0.3 and forbit = 0.6 respectively.

Fig. 24 Distribution of the average gas density of the injection regions (defined in Sect. 2.6.3) of supernova events. The distribution of SNe (solid lines) and the cumulative distributions (dotted lines) are shown. SNe in the CMZ and the solar circle are depicted in different colors. The dashed vertical lines indicate averages of the distributions. In the bottom left we also note the fraction of resolved SN events, for which we can inject thermal energy as opposed to kinetic energy.

6 Caveats

6.1 Resolution of the HII regions

There are a few limitations of these simulations which could quantitatively affect our results. The main limitation is that we are not able to resolve the most compact and embedded initial phase of the evolution of HII regions around O stars.

Indeed, the dotted lines in the top panel of Fig. 2 mark the size of a Strömgren sphere as a function of density and are demarcation lines for when the radiation feedback is resolved, based on the mass of the star. Note the steepness of these lines, which shows how expensive it would be to require the Strömgren radius to always be resolved. On the other hand, at later stages of the evolution of young massive stars, when the density progressively decreases, the interaction between the stellar radiation and the surrounding ISM is captured with increasing accuracy. The situation is even more hopeless if we want to resolve the radiation coming from B stars as well. Luckily most of the energy in ionizing radiation comes from the more massive O stars (see Fig. 4).

Since star particles are considered a small association of stars, more than one OB star might be contributing to the radiation feedback of a single particle, helping us resolving the surrounding Strömgren sphere. All this considered, we can see from Fig. 25 that while the early stages of the HII region are not resolved throughout the simulation, after just 0.5-1 Myr the stellar particles enter density regimes where, on average, excellent resolutions are achieved. Moreover, at any point in time, only ∼ 10% of O stars in the domain are in a regime where the Strömgren sphere is unresolved. Better capturing this initial evolution could lead to more efficient disruption of the surrounding gas, affecting timescales of the evolution.

Fig. 25 Distribution of the ratio of the ionizing photon flux S* of star particles containing at least one O star over the recombination rate of its closest gas cell Nrec. This is a measure of how well the Strömgren radius is resolved. In the top panel, we show the (cumulative) distribution of all stellar particles containing at least one O star for a given snapshot. For the bottom panel, we follow the evolution of individual star particles and show S*/Nrec averaged over all star particles as a function of their age. The shaded region shows ±σ around that average.

6.2 Softening of the gravity solver

Decoupling times and stellar dynamics would change if an unsoftened gravity solver was used and all the complexity of stellar interactions in associations and clusters would have been considered (such as binaries, close encounters etc...). In particular, stellar clusters are more easily dynamically dissolved in less resolved N-body simulations (Dehnen 2001 ; Heggie & Hut 2003). On the other hand, if the formed cluster is tightly bound, gravitational softening does not allow for strong gravity kicks from close stellar encounters, and the cluster stays compact (e.g. Lahén et al. 2025, their Figure 10). However, given our coherent comparison of the two environments within the same simulation framework, we can robustly highlight qualitative differences arising from the different conditions experienced by the stars and gas.

7 Summary and conclusions

We described a radiation-MHD simulation of a truncated Milky Way galaxy. We focused on and compared two high-resolution regions, a ring at the solar radius and the inner Galaxy, to understand and describe how the special environment of the CMZ affects ISM properties, star formation and early evolution, and feedback coupling.

Our main findings are summarized below:

• The CMZ is a much harsher environment than the disk. This affects the ISM properties and SF. The gas in the CMZ is denser on average, and a larger mass fraction is cold and molecular. The magnetic fields in the CMZ are enhanced in every density regime by up to one order of magnitude. The ISM is more turbulent with σ_z ~ 9 km/s against σ_z ~ 2 km/s for the solar circle. The star formation rate surface densities are higher, but the molecular depletion times are similar. However, at comparable molecular gas surface densities, the CMZ seems to be less efficient in converting gas into stars than in the solar circle;

• The CMZ is not just a high surface density analog of the disk (inflow, short orbital times, etc.). This leads to a difference in the interaction of young massive stars with their environment. In the disk, O stars disrupt clouds from within, following a progressive evolution, from dense and highly embedded to progressively more diffuse and exposed. Instead, O stars in the CMZ mostly follow cycles where the star once being in an exposed state, will be re-embedded at some later point;

• The difference in behavior can be explained by the high surface density environment, by the highly turbulent and dynamic ISM, by frequent cloud collisions between GMCs within the CMZ, and with inflowing gas from the bar lanes, but most noticeably, it can be explained by quick decoupling of stars and gas in the CMZ and by the efficient dynamical disruption of stellar associations through shear;

• These two effects strongly determine the coupling of feedback to the (dense) ISM in the CMZ and the regulation of SF. Dynamical decoupling and cluster shearing cause the feedback to be deposited far away from the GMC structure in which the stars formed. Moreover, the shearing deprives the ISM of the combined effect of clustered feedback and its disrupting power through strong superbubbles. Therefore, feedback is less disruptive for the global CMZ structure, but instead acts as a background source of turbulence.

Data availability

Movies of the simulation can be found at https://doi.org/10.7910/DVN/ULCTSZ, which can provide a better understanding of what is described here.

Acknowledgements

We thank the anonymous referee for constructive feedback on the manuscript. The authors are grateful for computing resources provided by the Ministry of Science, Research and the Arts (MWK) of the State of Baden-Württemberg through bwHPC and the German Science Foundation (DFG) through grants INST 35/1134-1 FUGG and 35/1597-1 FUGG, and also for data storage at SDS@hd funded through grants INST 35/1314-1 FUGG and INST 35/1503-1 FUGG. VMP, NB, RSK, SCOG, SR, JS, JG, and PG gratefully acknowledge financial support from the European Research Council via the ERC Synergy Grant “ECOGAL” (project ID 855130). FNL gratefully acknowledges financial support from grant PID2024-162148NA-I00, funded by MCIN/AEI/10.13039/501100011033 and the European Regional Development Fund (ERDF) “A way of making Europe”, from the Ramón y Cajal programme (RYC2023-044924-I) funded by MCIN/AEI/10.13039/501100011033 and FSE+, and from the Severo Ochoa grant CEX2021-001131-S, funded by MCIN/AEI/10.13039/501100011033. MCS, MD and ZXF acknowledge financial support from the European Research Council under the ERC Starting Grant “GalFlow” (grant 101116226). MCS further acknowledges support from Fondazione Cariplo under the grant ERC attrattività n. 2023-3014. NB acknowledges support from the ANR BRIDGES grant (ANR-23-CE31-0005). RSK and TP also acknowledge support from the German Excellence Strategy via the Heidelberg Cluster “STRUCTURES” (EXC 2181 - 390900948). In addition RSK is grateful for funding from the German Ministry for Economic Affairs and Climate Action in project “MAINN” (funding ID 50OO2206), and from DFG and ANR for project “STARCLUSTERS” (funding ID KL 1358/22-1). MH, JP and AP acknowledge funding from the Swiss National Science Foundation (SNF) via a PRIMA grant PR00P2 193577 ‘From cosmic dawn to high noon: the role of BHs for young galaxies’. CB gratefully acknowledges funding from National Science Foundation under Award Nos. 2108938, 2206510, 2414862, and CAREER 2145689, as well as from the National Aeronautics and Space Administration through the Astrophysics Data Analysis Program under Award “3-D MC: Mapping Circumnuclear Molecular Clouds from X-ray to Radio,” Grant No. 80NSSC22K1125 as well as participation in the PRIMA project under Grant No. 80NSSC25K7944.

References

All Tables

All Figures

	Fig. 1 Face-on view of the simulated Milky Way analog. The grayscale shows the total gas column density. The highlighted areas correspond to the enhanced-resolution regions of the simulation.
In the text

Fig. 2 Resolution as a function of density. Distribution of the spherical radius (top panel) and mass (bottom panel) of the cells as a function of density. The cells in the high- and low-resolution regions are shown in red and blue, respectively. The contours contain 100, 99, 90, 75, 50, and 25% of the total number of cells. The total number of cells in each region is indicated (bottom left corner in the bottom panel). We highlight the target mass in the two regions, as well as the maximum and minimum volume limits we imposed. In the bottom panel, we show with a gray overlay the region for which gas cells are more massive than one and eight times the Jeans mass. This gas is Jeans-unresolved and potentially star forming. The dotted lines indicate the demarcation lines for where the Strömgren radius of a 15 and 120 M⊙ star (which is the mass range of O-type stars) is resolved.

In the text

Fig. 3 Properties of the star particles in the high-resolution regions. We show the distribution of the mass of the star particles (first column), the density of the gas particle that generated the star particle (second column), the ionizing radiation flux assigned to the stellar particle (third column), the total mass of stars drawn from the IMF (fourth column) and the number of stars drawn (last column). Often, no massive star is drawn from the IMF sampling. The mass and parent density distribution of those star particles is shown in the upper right.

In the text

Fig. 4 Contribution to the total ionizing photon budget of massive (>8 M⊙) stars given their mass weighted with the IMF. The dotted line shows the cumulative fraction of massive stars, the solid line the cumulative fraction of ionizing photons at birth of the star, while the dashed line shows the cumulative fraction of the ionizing photons generated over the entire life of the star as a function of mass.

In the text

	Fig. 5 Background interstellar UV radiation field (photons with energy <13.6 eV) as a function of galactocentric radius. The value is scaled to the solar neighborhood value G0 (Draine 1978). The cosmic-ray ionization rate follows the same profile, scaled to the solar neighborhood value.
In the text

Fig. 6 Projections of the CMZ (top) and a region of the disk (bottom). The panels from left to right and top to bottom depict the total, atomic, molecular, ionized gas column densities, mass-weighted thermal pressure, temperature and magnetic field intensity averages, star formation rate surface density, stellar surface density and star particles colored by their age. These last three quantities are overlaid on the total gas column density in gray.

In the text

Fig. 7 Gas inflow/outflow to the CMZ (top), star formation rate in the CMZ (middle), and depletion times for the CMZ and solar circle (bottom) as a function of simulation time. Inflows and outflows are shown separately and the different gas phases stacked in different color to show the contribution of each phase to the total gas inflow/outflow shown as the solid purple line. Depletion times for the atomic, molecular and total gas are shown separately. The observed SFR value for the CMZ and its uncertainty is shown as the dashed line and shaded region in the middle panel.

In the text

	Fig. 8 Cumulative star formation as a function of galactocentric radius at t = 15 Myr.
In the text

Fig. 9 Star formation rate surface density as a function of molecular gas surface density. The values are computed over a 100 pc sized bins and shown separately in different colors for the CMZ (teal) and the solar circle (orange). The values are taken for several snapshots throughout the simulation. The colored area shows the kernel density estimation of the distributions and the black lines show lines of constant depletion time.

In the text

	Fig. 10 Density-temperature phase diagram of the gas in the simulation at t = 15 Myr. The colors indicate the gas mass fraction in the given (n – T) bin normalized by the bin size (dMgas(n, T)/Mgas,tot/d(log10 T)/d(log10 n)). The gas in the CMZ is shown in the left panel, the gas at the solar circle in the middle panel, and the difference on the right.
In the text

	Fig. 11 Number density probability distribution function for the CMZ (top) and the solar circle (bottom) at t = 15 Myr. The contribution of the different ISM phases is shown in different colors. Dotted vertical lines indicate average density values for the different phases.
In the text

	Fig. 12 Same as Fig. 11, but for the temperature.
In the text

	Fig. 13 Mass fraction in different ISM phases for the two highresolution regions at t = 15 Myr.
In the text

	Fig. 14 Average (dashed line) and median (solid line) magnetic field intensity as a function of density for gas in the CMZ (teal) and at the solar circle (orange) for the snapshot at t = 15 Myr. The shaded area around the mean indicates ±σ dispersion.
In the text

	Fig. 15 Velocity dispersion (σz) map for the CMZ and a section of the disk at t = 15 Myr. Values are obtained by inferring the mass-weighted velocity dispersion of cells in 20 pc bins. The contours show column densities of [10,50] cm−3. Lower density bins are excluded for clarity.
In the text

	Fig. 16 Average (dashed line) and median (solid line) z component of the velocity dispersion as a function of density for gas in the CMZ (teal) and at the solar circle (orange) at t = 15 Myr. The shaded area around the mean indicates the standard deviation.
In the text

	Fig. 17 Schematic representation of the assembly and disruption of a GMC. Stars form as the cloud is assembling (panels a-c), until their combined feedback is strong enough to disrupt the cloud from within (panel d).
In the text

Fig. 18 Time evolution of the environment surrounding an O star at the solar circle (top three rows) and in the CMZ (bottom three rows). In each panel, the first row shows the H2 density slice surrounding the selected star, which is centered in the box. Iso-density contours are shown at a density of 100 cm−3. The inlay panels shows the H2 column density of the larger-scale context. Surrounding O stars are shown as white dots in the inlay and turquoise stars in the main panel. We highlight the age of the selected star particle, the total number of OB stars assigned to it, and the mass of the most massive one among those. The magenta area is the region visible by the ionizing radiation of the star and which the star contributes to keeping ionized, i.e. the region of space for which rays connecting to the star only pass through ionized gas (xH+ > 0.5). The second and third rows show the H2 column density and size of the HII region, as seen by an observer at the position of the star. The second row is the Mollweide H2 column density projection. The third row shows the size of the HII region surrounding the star. As defined here, the size is the distance in each direction at which the fractional ionization of the gas first drops below 50% (i.e. xH+ < 0.5). Note that the panels are not equally spaced in time, but instead a few salient snapshots in the evolution of the star are selected.

In the text

Fig. 19 Density and ionization profile around O stars. The top two rows are the case for the solar circle, and the bottom two rows for the CMZ. For each row, the first column is an example for one selected star, while the second column is the average over all stellar particles in the region (solar circle or CMZ) that contain at least one O star. For each region, the first row shows the density as a function of distance from the O stars, averaged over spherical shells. The second row depicts the fraction of the surface that the O stars are helping to keep ionized, i.e. the fraction of a surface of a given radius centered on the star for which a ray passes only through ionized gas (xH+ > 0.5). These quantities are shown at different ages of the star particles. The last column shows the time evolution of the environment around a selection of seven O stars at a radius of 10 pc.

In the text

	Fig. 20 Fraction of stars classified based on their evolution (clarified in the text) for the solar circle and the CMZ. Most of the O stars in the CMZ will undergo at least one event where the star after being exposed, will be re-embedded in a high density environment, in stark contrast to what O stars at the solar circle experience.
In the text

Fig. 21 Orbits of a few selected stellar particles containing O stars at birth. The orbit is colored by the age of the stellar particle, and each SN event in its lifetime associated with the particle is marked along its orbit. The birthplace of the star is highlighted with the blue x-marker. In the background in grayscale we show the averaged H2 density over a period of 20 Myr starting from the snapshot where the star was born. The black contour shows an average density of 50 cm−3. In the first panel we show an early snapshot when the inflow from the bar lanes is still quite copious. In the middle panel we show the orbit of a star which formed in a GMC associated with the bar lanes and inherited their velocity. On the right panel a snapshot at a later time was chosen, where the CMZ is more regular, and the inflow smaller and more steady. Stars quickly decouple from the gas on orbital timescales, which are short in the CMZ.

In the text

Fig. 22 Schematic representation of the early evolution of a starforming region in the CMZ and the subsequent evolution of the formed stars. Compared to Fig. 17, the situation is more complex in the dynamic environment of the Galactic Center. Stellar-gas decoupling (panel b) due to cloud collisions and short orbital times, leads to re-embedding events in the early evolution of the stars (panel d). This has an effect on how feedback can couple to the dense gas in the CMZ.

In the text

Fig. 23 Young stellar association that is quickly sheared apart in the CMZ environment. We highlight new-born (<3 Myr) O stars at a given time in the simulation (left, yellow stars) and select a compact association of coeval O stars in red. We follow these stars in their evolution and show their position with respect to the background molecular gas (in grayscale) after t = 3 Myr (middle) and t = 6 Myr (right), which correspond to orbital fractions of forbit = 0.3 and forbit = 0.6 respectively.

In the text

Fig. 24 Distribution of the average gas density of the injection regions (defined in Sect. 2.6.3) of supernova events. The distribution of SNe (solid lines) and the cumulative distributions (dotted lines) are shown. SNe in the CMZ and the solar circle are depicted in different colors. The dashed vertical lines indicate averages of the distributions. In the bottom left we also note the fraction of resolved SN events, for which we can inject thermal energy as opposed to kinetic energy.

In the text

Fig. 25 Distribution of the ratio of the ionizing photon flux S* of star particles containing at least one O star over the recombination rate of its closest gas cell Nrec. This is a measure of how well the Strömgren radius is resolved. In the top panel, we show the (cumulative) distribution of all stellar particles containing at least one O star for a given snapshot. For the bottom panel, we follow the evolution of individual star particles and show S*/Nrec averaged over all star particles as a function of their age. The shaded region shows ±σ around that average.

In the text