© The Authors 2024

1 Introduction

The Milky Way’s Central Molecular Zone (CMZ) is a ring-like ~3 × 10⁷ M_⊙ accumulation of molecular gas within Galactocentric radius R ≃ 200 pc (Morris & Serabyn 1996; Henshaw et al. 2023). It is generated by the Galactic bar, which efficiently transports gas from Galactocentric radius R ≈ 3 kpc down to R ≈ 200 pc (Sormani & Barnes 2019; Hatchfield et al. 2021). The CMZ is the most extreme star-forming environment in the entire Milky Way and has emerged in the last decade as an important astrophysical laboratory to study the physics of the interstellar medium (ISM) and star formation (Henshaw et al. 2023).

The CMZ is permeated by a strong magnetic field that is likely to play an important role in the dynamics of the ISM and in the process of star formation (Ferrière 2009; Morris 2015; Butterfield et al. 2024). However, many aspects of the magnetic field configuration are poorly understood. For example, it is unclear whether the CMZ is immersed in a pervasive |B| ≳ 1 mG field (Morris & Yusef-Zadeh 1989; Morris 2006), or whether the magnetic field drops to |B| ~ 100 µG or less in the more diffuse inter-cloud medium (Tsuboi et al. 1985; Yusef-Zadeh & Morris 1987; Lang et al. 1999a,b; LaRosa et al. 2005; Ferrière 2009; Yusef-Zadeh et al. 2022) while reaching mG strengths only inside very dense gas (Schwarz & Lasenby 1990; Killeen et al. 1992; Plante et al. 1995; Uchida & Guesten 1995; Marshall et al. 1995; Yusef-Zadeh et al. 1999; Pillai et al. 2015).

The large-scale geometry of the magnetic field is also unclear. Observations show that the magnetic field in the dense and cold ISM is oriented predominantly parallel to the Galactic plane (in projection on the plane of the sky), while the magnetic field in the more diffuse ISM, including a population of prominent filamentary structures known as non-thermal filaments (NTFs, Yusef-Zadeh et al. 1984; Tsuboi et al. 1986; Heywood et al. 2022), is predominantly perpendicular to the Galactic plane (Novak et al. 2003; Chuss et al. 2003; Nishiyama et al. 2010; Mangilli et al. 2019; Guan et al. 2021; Hu et al. 2022b,c; Butterfield et al. 2024; Paré et al. 2024). This can be appreciated for example in Fig. 6 of Guan et al. (2021), which shows a change in the observed magnetic field geometry as the fractional contribution of synchrotron radiation (tracing ionised gas) and thermal dust emission (tracing cold neutral gas) varies with frequency, or in Fig. 1 of Nishiyama et al. (2010), which shows that B is prevalently parallel to the Galactic plane at latitudes |b| < 0.4°, where the dense gas dominates their measurements, and becomes perpendicular to the Galactic plane at |b| > 0.4°, where the diffuse ionised gas dominates the measurements. This might point to the presence of two separate magnetic field systems, one predominantly perpendicular and one predominantly parallel to the plane (Morris 2015). It is presently unclear how the two systems relate to each other and what is their three-dimensional geometry.

Since the CMZ is a star-forming nuclear ring similar to those that are commonly found at the centres of barred galaxies (Mazzuca et al. 2008; Comerón et al. 2010; Ma et al. 2018), it is reasonable to expect that its magnetic field system shares many similarities with those of other nuclear rings. Measurements have been performed for a handful of galaxies (for a review, we refer to e.g. Beck 2015), the best studied example probably being NGC 1097 (Beck et al. 2005; Tabatabaei et al. 2018; Lopez-Rodriguez et al. 2021; Hu et al. 2022a). These measurements lead to the following general conclusions: (i) The magnetic field lines inferred from radio polarisation maps of synchrotron-emitting gas in the bar region surrounding the nuclear ring are approximately aligned with the gas streamlines. In particular, the magnetic field in the bar lanes that transport the gas towards the nuclear ring is approximately parallel to the lanes in the frame co-rotating with the bar (e.g. Fig. 2 of Beck et al. 2005). (ii) The magnetic field in the ring spirals towards the centre with a relatively large pitch-angle both in radio and far-infrared polarisation maps (e.g. Fig. 1 of Lopez-Rodriguez et al. 2021). The pitch-angle and general geometry of the magnetic field inside the ring are tracer dependent (Lopez-Rodriguez et al. 2021; Hu et al. 2022a), which is reminiscent of the Milky Way, where as mentioned above the projected orientation of the field near the mid-plane depends on the tracer. Measurements of the magnetic field in dense molecular gas are currently not available for external galaxies (Lopez-Rodriguez et al. 2021). (iii) The equipartition magnetic field in synchrotron-emitting gas in the nuclear ring of NGC 1097 is ~60 µG, which is of the same order of magnitude as similar estimates for the diffuse gas in the CMZ (LaRosa et al. 2005; Morris 2006; Yusef-Zadeh et al. 2022).

The effects of the magnetic fields on the dynamics of the ISM in the CMZ are also poorly understood (e.g. Morris 2006, 2015). Magneto-hydrodynamic (MHD) instabilities are believed to be one of the primary mechanisms for mass and angular momentum transport in astrophysical accretion discs (Balbus & Hawley 1998). A back-of-the-envelope calculation suggests that the magneto-rotational instability (MRI; Balbus & Hawley 1991) should transport gas in the CMZ at a rate given by: $$\dot{M}=3\alpha {\text{M}}_{\odot }{\text{yr}}^{-1}\left(\frac{M}{5\times {10}^7{\text{M}}_{\odot }}\right)\left(\frac{\sigma }{15\text{\hspace{0.17em}km\hspace{0.17em}}{\text{s}}^{-1}}\right)\left(\frac{h}{40\text{pc}}\right){\left(\frac{R}{100\text{pc}}\right)}^{-2},$$(1)

where M is the total gas mass of the CMZ, σ is the gas velocity dispersion, h is the gas scale-height, R is the Galactocentric radius, α is the Shakura & Sunyaev (1973) coefficient which we have assumed to be determined by the MRI, and we have inserted typical CMZ values in the denominators. Assuming α ≃ 0.01, the predicted inflow rate is Ṁ = 0.03 M_⊙ yг⁻¹. This value would be significant, because at this rate the entire circum-nuclear disc (CND; Genzel et al. 1985; Mills et al. 2017; Hsieh et al. 2021), which is the closest large gas reservoir to SgrA* with a mass of M_CND ~ 5 × 10⁴ M_⊙ (Etxaluze et al. 2011; Requena-Torres et al. 2012), would build up on a rather short timescale of M_CND/Ṁ ~ 1.7 Myг. However, simple order-of-magnitude estimates such as these are inherently limited. The MRI-driven transport is traditionally studied in the context of Keplerian, weakly magnetised accretion discs (e.g Balbus 2003), and it is much less understood in the context of non-Keplerian, non-axisymmetric potential of the Galactic centre (Kim & Stone 2012), and in the cold ISM regime with small plasma β that is relevant there (Kim & Ostriker 2000; Kim et al. 2003; Piontek & Ostriker 2007; Jacquemin-Ide et al. 2021). Therefore, it is currently unclear how effective MRI-driven transport is in the Galactic centre, and how important it is compared to other effects such as redistribution of angular momentum due to stellar feedback (Tress et al. 2020).

There have been a number of theoretical studies of magnetised gas flow in barred galaxies, which can broadly be divided in two groups. The first group uses dynamo theory to follow the evolution of magnetic fields using an approximate set of equations in a prescribed velocity field (e.g. Otmianowska-Mazur et al. 2002; Moss et al. 2001, 2007). The velocity field can be taken for example from purely hydrodynamical simulations of gas flow in barred galaxies (e.g. Athanassoula 1992). This approach is computationally fast but requires assumptions on how velocity fluctuations act on the magnetic fields on unresolved scales and ignores the back reaction of the magnetic field on the gas. The second approach uses the full set of MHD equations without approximations. For example, Kim & Stone (2012) performed MHD simulations of barred galaxies, finding that magnetic fields can enhance inflows and that the bar potential plays a role in dynamo action. However, their simulations are two-dimensional and therefore unable to capture the potential effects of poloidal fields and other dynamical processes that may be important in three dimensions. Suzuki et al. (2015) and Kakiuchi et al. (2024) performed three-dimensional MHD simulations of gas flow in the innermost few kiloparsec of the Milky Way, but they assumed an axisymmetric potential and neglected the presence of the Galactic bar. Moon et al. (2023) run semi-global MHD simulations of nuclear rings in an axisymmetric potential subject to a prescribed mass inflow rate, and found that magnetic fields can drive radial flows from the ring inwards and that they can suppress star formation in the ring. However, to the best of our knowledge, a global sub-pc 3D MHD simulation of gas flow in a barred potential is still missing.

In this paper we use global 3D MHD simulations of gas flow in a Milky Way barred potential to address the following open questions:

• How do magnetic fields affect the gas morphology in the CMZ? (Sect. 3).

• What is the geometry of the magnetic field in the CMZ and in the surrounding bar region? (Sect. 4).

• Do magnetic fields drive turbulence? (Sect. 5).

• How are magnetic fields amplified and maintained? (Sect. 6).

• Do magnetic fields enhance the bar-driven inflow from the Galactic disc to the CMZ? Do magnetic fields drive a nuclear inflow from the CMZ towards the central few parsecs? (Sect. 7).

The aim is to study magnetic fields and their effect on the gas dynamics. We therefore deliberately choose not to include the gas self-gravity, nor any type of star formation and stellar feedback. Such additional processes would make it very difficult to isolate the contribution of magnetic fields for example in driving turbulence or changing the probability density function of the gas. Studying these processes and their (possibly non-linear) interaction with the magnetic fields is an important avenue for future work.

This paper is structured as follows. In Sect. 2 we describe our numerical methods. Sections 3-7 are dedicated to addressing the open questions listed above. We sum up our results in Sect. 8.

2 Numerical methods

We run three-dimensional MHD simulations of gas flow in the inner regions of the Milky Way using the moving-mesh code AREPO (Springel 2010; Pakmor et al. 2016; Weinberger et al. 2020). The simulated gas disc covers the entire region within Galactocentric radius R = 5 kpc. We assume ideal MHD which is generally an excellent approximation for the ISM (e.g. Marinacci et al. 2018a) and use the standard MHD implementation contained in AREPO, which has been previously employed for galaxy simulations and tested against a number of standard test problems including the development of the magneto-rotational instability (e.g. Pakmor & Springel 2013; Marinacci et al. 2018b). The equations we solve are: $$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho v)=0,$$(2) $$\frac{\partial (\rho v)}{\partial t}+\nabla \cdot (\rho vv+P\mathbb{I}+\mathbb{T})=-\rho \nabla \mathrm{\Phi },$$(3) $$\frac{\partial (\rho e)}{\partial t}+\nabla [(\rho e)v+(P\mathbb{I}+\mathbb{T})\cdot v]=\rho \frac{\partial \mathrm{\Phi }}{\partial t}-ℒ,$$(4) $$\frac{\partial B}{\partial t}=\nabla \times (v\times B),$$(5)

where ρ is the gas density, v is the velocity, P is the thermal pressure, 𝕀 is the identity matrix, B is the magnetic field, 𝕋 = B²/(8π)𝕀 − BB/(4π) is the Maxwell stress tensor under the approximation of non-relativistic ideal MHD, Φ is the external gravitational potential, ρe = ρe_th + ρv²/2 + B²/(8π) + ρΦ is the total energy per unit volume, which is the sum of the thermal (ρe_th), kinetic (ρv^2;/2), magnetic (B²/(8π)), and gravitational (ρΦ) contributions, and ℒ is the net cooling (or heating) rate per unit volume. Magnetic field divergence errors can arise as a result of the discretization of the MHD equations. In Appendix F we checked that these are always under control and do not dominate the dynamics of the ISM.

The gas is assumed to flow in an externally imposed Milky Way barred potential. The potential is identical to that used in Ridley et al. (2017) and is described in detail in Section 3.2 of that paper (their Fig. 4 shows the rotation curve). The bar rotates rigidly with a pattern speed Ω_p = 40 km s⁻¹ kpc, consistent with recent determinations for the Milky Way (e.g. Sormani et al. 2015b; Portail et al. 2017; Sanders et al. 2019; Li et al. 2022; Clarke & Gerhard 2022). This places the (single) inner Lindblad resonance (ILR) calculated in the epicyclic approximation at R_ILR = 1.1 kpc and the corotation resonance at R_CR = 5.9 kpc. This potential was chosen to allow direct comparison with previous simulations in Ridley et al. (2017) and Sormani et al. (2018) that used the same potential. We did not include gas self-gravity nor the consequent star formation in this paper as we aim to isolate the effects of the magnetic fields from other competing effects on the dynamics of the ISM.

We run simulations using two different thermodynamic setups. The ‘isothermal’ simulations used an isothermal equation of state $$P=c_{\text{s}}^2\rho ,$$(6)

where c_s is a constant that we set to either c_s = 1 km s⁻¹ or c_s = 10 km s⁻¹ (representative of ‘cold’ and ‘warm’ ISM). In the isothermal simulations, ℒ = 0 in Eq. (4). We note that the isothermal approximation does not correspond to a physical situation where there is no heating and cooling. Instead, it means that cooling and heating always exactly balance in such a way that the gas temperature is kept constant (e.g. Klessen & Glover 2016). For example, the energy released in shocks is instantaneously radiated away in the isothermal approximation.

The ‘chemistry’ simulations used the adiabatic equation of state $$P=(\gamma -1)\rho e_{\text{th}},$$(7)

where γ = 5/3 is the adiabatic index. These simulations account for the chemical evolution of the gas using an updated version of the NL97 chemical network from Glover & Clark (2012), which is based on the work of Glover & Mac Low (2007a,b) and Nelson & Langer (1997). This network solves for the non-equilibrium abundances of H, H₂, H⁺, C⁺, O, CO, and free electrons. The heating and cooling contained in the term ℒ in Eq. (4) are calculated on-the-fly by the network based on the instantaneous chemical composition of the gas and taking into account a number of processes, including radiative cooling, heat released by the formation of H₂ on dust grains, and an averaged interstellar radiation field (ISRF) and cosmic ray ionization rate (CRIR). The ISRF is set to the standard value G₀ measured in the solar neighbourhood (Draine 1978) diminished by a local attenuation factor which depends on the amount of gas present within 30 pc of each computational cell. This attenuation factor is introduced to account for the effects of dust extinction and H₂ self-shielding and is calculated using the TREECOL algorithm described in Clark et al. (2012). The value was chosen originally as it was similar to the typical observed separation of OB stars in the Solar neighbourhood. Although in the dense CMZ environment the separation might be smaller, we choose to keep the same value here for consistency with previous simulations (Tress et al. 2020). The cosmic ray ionisation rate (CRIR) is fixed to ζ_H = 3 × 10⁻¹⁷ s⁻¹ (Goldsmith & Langer 1978). Although these values are typical for the Solar neighbourhood and likely too small for the CMZ (Clark et al. 2013; Oka et al. 2019), we expect this to have little effects on the dynamics of the gas discussed in this paper. Indeed, Sormani et al. (2018) (we refer also to the discussion in Section 2.3 of Tress et al. 2020) have shown that the strength of the ISRF and CRIR do not affect the large-scale dynamics of the gas in the Galactic Centre region significantly. The main effect is to modify the amount of gas in different ISM phases since cosmic rays and UV photons dissociate and ionise molecular gas. Increasing the ISRF and CRIR has an effect on the dynamics (and on the MRI and inflow rates) that is similar to increasing the effective sound speed of the gas. To check this, we have run some exploratory simulations with higher ISRF and CRIR, and we found that the dynamical behaviour of these simulations is similar to the high-c_s simulations discussed below in Sect. 3. A more complete study of the effects of varying the ISRF and CRIR is outside the scope of this work.

Finally, we imposed a numerical temperature floor T_floor = 20 K on the simulated ISM. Without this floor, the code would occasionally produce anomalously low temperatures in cells close to the resolution limit undergoing strong adiabatic cooling, causing it to crash. The chemical network is the same as used in Sormani et al. (2018) and Tress et al. (2020), and more details can be found in Section 3.4 of the former or Section 2.3 of the latter. Figure 1 shows the typical density-temperature phase diagram in our chemistry simulations.

The number density is defined as $$n=\frac{\rho }{\mu m_{\text{p}}},$$(8)

where µ is the mean molecular weight and m_p is the proton mass. As a reference, at the assumed solar metallicity the mean molecular weight is µ = 0.67, 1.27, 2.23 for fully ionised, neutral, and fully molecular gas respectively. The temperature of the gas in the chemistry simulations is defined as T = P/(nk_b), where k_b is the Boltzmann constant.

Fig. 1 Gas temperature as a function of total gas density in our fiducial CHEM_MHD model at t = 196 Myr. Blue and red indicate the regions at R > 0.5 kpc and R < 0.5 kpc, respectively. Contours contain [99, 95. medium with a cold (T = 102 K) and warm (T = 104 K) phase. The hot phase (T = 106 K) is absent since our simulations do not include stellar feedback or other processes that can create it.

2.1 Initial conditions

We initialised the density according to the following axisymmetric density distribution: $$\rho (R,ɀ)=\frac{{\Sigma }_0}{4{ɀ}_{\text{d}}}\exp \left(-\frac{R_{\text{m}}}{R}-\frac{R}{R_{\text{d}}}\right){\mathrm{sech}}^2\left(\frac{ɀ}{2{ɀ}_{\text{d}}}\right),$$(9)

where (R, ϕ, ɀ) denote standard cylindrical coordinates, ɀ_d = 85 pc, R_d = 7 kpc, R_m = 1.5 kpc, ∑₀ = 50 M_⊙pc⁻², and we cut the disc so that ρ = 10⁻²⁸ g/cm³ for R ≥ 5 kpc. This profile roughly matches the observed radial distribution of gas in the Galaxy (Kalberla & Dedes 2008; Heyer & Dame 2015) and is identical to the one used in Tress et al. (2020). The total initial gas mass in the simulation is ≃1.5 × 10⁹ M_⊙. The computational box has a total size of 24 × 24 × 24 kpc with periodic boundary conditions. The box is sufficiently large that the outer boundary has a negligible effect on the evolution of the simulated galaxy.

In order to avoid transients, we introduced the bar gradually (e.g. Athanassoula 1992). We started with gas in equilibrium on circular orbits in an axisymmetrised potential and then we turned on the non-axisymmetric part of the potential linearly during the first 146 Myr (approximately one bar rotation) while keeping constant the total mass which generates the underlying external potential (not to be confused with the mass of the gas in the simulation). Therefore, only the simulation at t ≥ 146 Myr, when the bar is fully on, will be considered for the analysis in this paper. The simulations were run until t = 300 Myr.

The initial temperature for the chemistry simulations is T₀ = 10⁴ K everywhere. The precise value does not affect the outcome of the simulation since a new equilibrium is reached within a few megayears (and well before the bar is fully turned on) through the balance of heating and cooling processes.

Unless otherwise specified, we started with a purely poloidal uniform ‘seed’ magnetic field of $B_0=0.02\mu \text{G}\widehat{z}$. We have also experimented with different initial magnetic field strengths and with initial toroidal (rather than poloidal) geometry; the results of these experiments are briefly discussed in Appendix C.

2.2 Summary of simulation runs

Table 1 shows a summary of the main simulations presented in this paper. In addition to these simulations, we have run various tests in which we varied parameters such as the resolution, the initial magnetic field, or where we cut out the CMZ to isolate it from interaction with the large-scale environment. These additional simulations are introduced and discussed when appropriate throughout the paper.

The resolution of our simulations is specified by the target mass of each computational cell. The target mass for all the simulations listed in Table 1 is 100 M_⊙ for R > 500 pc and 10 M_⊙ for R < 500 pc. The resolution is therefore higher in the CMZ region than in the Galactic disc. The system of mass refinement in AREPO splits cells whose mass becomes greater than twice the target mass and merges cells whose mass becomes lower than half the target mass. Because this keeps the mass of the cells approximately constant, our spatial resolution varies as a function of the local gas density. We also implemented a minimum cell volume to prevent excessive refinement and computational slowdown in areas of high density: cells with an effective cell radius less than r_cell = [3V_cell/(4π)]^1/3 = 0.1 pc, where V_cell is the cell volume, were not permitted to divide into smaller cells. The typical number of cells in our simulations is around 25 million, of which approximately 10 million are located in the higher-resolution region at R < 500 pc. Figure 2 shows the resolution as a function of density for our fiducial model CHEM_MHD.

3 Gas morphology

Figures 3 and 4 show the face-on gas surface density and magnetic fields for the six simulations listed in Table 1 as well as the magnetic field for the three magnetised simulations, while Figs. 5 and 6 show the time evolution of our fiducial CHEM_MHD simulation. The gas morphology and general flow pattern in all simulations have the typical characteristics of gas flow in barred potentials, such as the presence of large-scale shocks on the leading side of the bar that transport the gas towards the centre (also known as ‘bar dust lanes’, e.g. Athanassoula 1992), and a central ring-like accumulation of gas, that in the Milky Way corresponds to the CMZ. The general characteristics of this flow have been extensively discussed in numerous papers to which we refer for more in-depth discussions (e.g. Athanassoula 1992; Sellwood & Wilkinson 1993; Fux 1999; Kim et al. 2012; Sormani et al. 2015a, 2018; Tress et al. 2020). Here, we focus only on the differences that appear when a magnetic field is introduced.

The first difference is that the magnetic fields tend to decrease the radius of the CMZ ring-like structure. This is noticeable if we compare the ISO_10_HD to the ISO_10_MHD simulations in the right column of Fig. 3. The CMZ in the magnetised simulation is slightly smaller than in the non-magnetised one. The explanation is likely the following. It is well-known that the radius of the nuclear ring in simulations is strongly dependent on the sound speed (e.g. Englmaier & Gerhard 1997; Patsis & Athanassoula 2000; Li et al. 2015; Sormani et al. 2015a). Sormani et al. (2024) argued that this dependence can be explained in terms of density waves excited by the bar potential. These density waves remove angular momentum, clear out a region around the inner Lindblad resonance, and transport the gas inwards where it accumulates into a ring. When the sound speed is larger, density waves are stronger and can be excited over a more extended region, and transport the gas into a ring of smaller radius. Magnetic fields increase the effective sound speed of the gas by exerting magnetic pressure, and therefore produce smaller rings. The amount by which the effective sound speed is increased by magnetic fields can be roughly estimated by adding in quadrature the Alfvén velocity defined as $$v_A=\frac{B}{{(4\pi \rho )}^{1/2}}.$$(10)

In our simulations, the Alfvén velocity in the dense gas in the CMZ ring is typically of the order of |v_A| ≃ 5 km s⁻¹ (Sect. 4.3), and indeed the effect seen in Fig. 3 is comparable to the effect seen in isothermal unmagnetised simulations when the sound speed is increased by roughly this amount in quadrature (Sormani et al. 2024).

A second difference is the probability density function (PDF). It is well known that magnetic fields can affect the density PDF (e.g. Federrath & Klessen 2013). Figure 7 shows that in the unmagnetised simulations most of the gas mass in the CMZ (thick blue line) lies at the highest densities (n ≳ 10⁶ cm⁻³), because all the gas in the dense ring tends to occupy the same orbit and there is only the thermal pressure preventing further compression. In the magnetised simulation instead, the mass PDF has a peak at a density of n ~ 10³ cm⁻³. This is because the magnetic fields provide pressure support and also drive turbulence, which increases the random motions of the gas and prevents it from accumulating at too high densities (see also Molina et al. 2012). Indeed, Fig. 4 shows that the CMZ ring in the CHEM_HD simulation is very thin and dense¹, while in the CHEM_MHD simulation it is puffed up by turbulence. Turbulence also puffs up the disc in the vertical direction and increases the disc vertical scale-height, which is known to be directly related to the amount of turbulence in galactic discs (e.g. Ostriker & Kim 2022). We will discuss turbulence and the mechanism driving it more in detail in Sect. 5.

A third difference occurs in the region inside the dense ring. Figure 6 shows that the region inside the ring in the CHEM_MHD simulation is devoid of gas at t = 147 Myr, and then gradually gets filled with gas. By contrast, in the CHEM_HD simulation the region inside the ring remains devoid of gas at all times. Figure 4 illustrates this difference in the HD and MHD simulations by comparing snapshots at the same time t = 196 Myr (compare the CHEM_HD panel with the CHEM_MHD panel). The filling up of the ring interior in the magnetised simulation occurs because magnetic fields drive inward accretion from the CMZ towards the central few parsecs. It is interesting to note that supernova feedback can also produce a similar effect of filling up the ring (Fig. 9 in Tress et al. 2020). Thus, it will be important in the future to understand which effect is stronger, and what is the non-linear interaction between the two. We discuss further the inflows driven by the magnetic field in Sect. 7.

Fig. 2 Resolution in our fiducial CHEM_MHD simulation (Table 1) at t = 196 Myr. Top: spatial resolution as a function of total gas density. rcell is the radius of a sphere with the same volume as the cell. The horizontal black dashed line indicates the volume limit of the cells (Sect. 2.2). Middle: mass of cells as a function of density. In both panels, blue contours are for the region at R > 0.5 kpc with a target mass of 100 M⊙, while red contours are for the R < 0.5 kpc with higher resolution at a target mass of 10 M⊙ (Sect. 2.2). Bottom: λMRI/rcell as a function of density for cells in the CMZ (R < 0.5 kpc), where λMRI = 2πυA/Ω. is the characteristic wavelength of the MRI, υA is the Alfven speed, and Ω = υϕ/R is the angular velocity of each gas cell. We see that the MRI is well resolved in the CMZ. The contours contain [100, 99, 90, 75, 50, 25] % of the total number of cells.

Fig. 3 State of the system at t = 196 Myr for our set of simulations. Top and middle: face-on gas surface density for all the simulations listed in Table 1. Bottom: projected mass-weighted magnetic field $\left|{\langle B\rangle }_z\right|={\left({\langle B_x\rangle }_z^2+{\langle B_y\rangle }_z^2+{\langle B_z\rangle }_z^2\right)}^{1/2}$ for the three magnetised simulations in Table 1, where 〈X〉ɀ = (∫ ρX dɀ)/(∫ ρ dɀ) denotes the vertically integrated mass-weighted density of a quantity X. Contours denote the streamlines of the projected field ${\langle B_x\rangle }_z{\widehat{e}}_x+{\langle B_y\rangle }_z{\widehat{e}}_y$. The dashed circles indicate the region with higher numerical resolution at R < 500 pc (Sect. 2.2). Rotation is clockwise. Movies are available online.

Fig. 4 Same as Fig. 3, but zooming-in onto the Central Molecular Zone. Movies are available online.

Fig. 5 Time evolution of the surface gas density in the CHEM_MHD simulation. Rotation is clockwise.

Fig. 6 Time evolution of the gas density in the central regions of the CHEM_MHD simulations. This is the same as Fig. 5, but zooming-in in the central region. The interior of the CMZ gas ring is empty at t = 147 Myr and then is gradually filled with gas due to MRI-driven accretion (Sect. 7).

4 Magnetic fields in the bar and CMZ regions

A general impression of the magnetic field geometry and strength in our simulations can be obtained from the bottom rows in Figs. 3 and 4. An alternative visualisation of the magnetic field in the CMZ for our fiducial CHEM_MHD simulation is shown in Fig. 8. When we consider the magnetic field properties in more detail, we find the following general characteristics, which are explored in the dedicated subsections below:

1. The magnetic field can be understood as the sum of a ‘regular’ time-averaged component and a fluctuating ‘turbulent’ component (Sect. 4.1).

2. The magnetic field is generally aligned with the gas velocity vectors, and the magnetic field geometry changes from toroidal near the ɀ = 0 plane to poloidal at |z| > 0 (Sect. 4.2).

3. The magnetic field strength scales as a function of gas density roughly as B ∝ n^0.33(Sect 4.3).

Fig. 7 Mass-weighted probability density function (PDF) of the CHEM_MHD and CHEM_HD simulations (Table 1) at t = 196 Myr. The thick lines are for the CMZ (R < 500 pc), while the thin lines are for the entire simulation. Each curve is normalised so that ∫(dM/d log n) d log n = 1. In the CHEM_HD simulation, most of the gas in the ring accumulates at the highest densities purely due to orbital convergence. This is not due to the gas self-gravity, since it is switched off in our simulations. The presence of magnetic fields shifts the peak to lower densities by providing pressure support and driving random motions. The total distribution is bimodal, as expected from a two-phase medium (e.g. Seta & Federrath 2022).

4.1 Decomposition into regular and turbulent components

It is instructive to decompose the magnetic field as $$B(R,\varphi ,ɀ,t)={\langle B\rangle }_t(R,\varphi ,ɀ)+\mathrm{\Delta }B(R,\varphi ,ɀ,t)$$(11)

where 〈B〉_t is a time-averaged ‘regular’ component, ΔB is an irregular instantaneous ‘turbulent’ fluctuation, and the time average of a quantity X over an interval ∆t is defined as: $${\langle X\rangle }_t=\frac{{\displaystyle {\int }_t^{t+\mathrm{\Delta }t}X}(R,\varphi ,ɀ,t)\text{d}t}{\mathrm{\Delta }t}.$$(12)

We note that 〈∆B_i〉_t = 0 and ${\langle B_i^2\rangle }_t={\langle B_i\rangle }_t^2+{\langle \mathrm{\Delta }{B}_i^2\rangle }_t$ by definition. Moreover, although the flow reaches an approximate steady-state at t > 150 Myr, some slow changes in the global gas and magnetic field configuration continue to happen as a consequence of the continuous gas inflow towards the centre. However, these changes are slow enough (typical timescale ~100 Myr) that the decomposition into time-averaged and instantaneous components proves to be useful over timescales of tens of megayears (corresponding to a few rotations in the gas ring, where the orbital period is ~10 Myr). Here we considered time averages over t = 200–250 Myr. The conclusions of this subsection are not significantly affected by this choice.

Figure 9 compares the instantaneous field B with the time-averaged 〈B〉_t for our fiducial CHEM_MHD simulation. The time-averaged magnetic field has a very regular structure which resembles the velocity streamlines of gas flowing in a bar potential. Albeit regular, the time-averaged field is far from simple, and exhibits a complex ‘butterfly’ morphology in the xz and yz planes, which we discuss in more detail in Sect. 4.2. The difference ΔB = B – 〈B〉_t is larger where the gas is more turbulent, for example in the nuclear ring, indicating that the turbulence causes magnetic field fluctuations.

We quantified the strength of the regular and irregular magnetic fields using the following mass-weighted and volume-weighted averages: $$B_{\text{tot},\rho }={\displaystyle ⟨}\frac{{\displaystyle {\int }_V|}B|\rho \text{\hspace{0.17em}d}V}{{\displaystyle {\int }_V\rho }\text{\hspace{0.17em}d}V}{\displaystyle ⟩}$$(13) $$B_{\text{reg\hspace{0.17em}},\rho }=\frac{{\displaystyle {\int }_V\left|{\langle B\rangle }_t\right|}{\langle \rho \rangle }_t\text{\hspace{0.17em}d}V}{{\displaystyle {\int }_V{\langle \rho \rangle }_t}\text{\hspace{0.17em}d}V}$$(14) $$B_{\text{trb},\rho }=⟨\frac{{\displaystyle {\int }_V}\left|\mathrm{\Delta }B\right|\rho \text{\hspace{0.17em}d}V}{{\displaystyle {\int }_V\rho }\text{\hspace{0.17em}d}V}⟩$$(15)

and $$B_{\text{tot},V}=⟨\frac{{\displaystyle {\int }_V}\left|B\right|\text{d}V}{{\displaystyle {\int }_V\text{\hspace{0.17em}d}}V}⟩$$(16) $$B_{\text{reg\hspace{0.17em}},V}=\frac{{\displaystyle {\int }_V\left|{\langle B\rangle }_t\right|}\text{d}V}{{\displaystyle {\int }_V\text{\hspace{0.17em}d}}V}$$(17) $$B_{\text{trb},V}=⟨\frac{{\displaystyle {\int }_V}\left|\mathrm{\Delta }B\right|\text{d}V}{{\displaystyle {\int }_V\text{\hspace{0.17em}d}}V}⟩$$(18)

where the integrals are carried out over a volume V. Taking V as the region where R < 500 pc and the time-averaged density is 〈n〉_t > 10² cm⁻³, we find B_tot,ρ ≃ 305 µG, B_reg,ρ ≃ 34 µG, B_trb,ρ ≃ 298 µG, and B_tot,V ≃ 33 µG, B_reg,V ≃ 16 µG, B_trb,V, ≃ 29 µG. The mass-weighted turbulent component is larger than the volume-weighted turbulent component because denser gas (where most of the mass is) is more turbulent than the diffuse gas on average.

Figures 10–12 offer a 3D visualisation of the instantaneous and time-averaged magnetic field, confirming the complex nature of the instantaneous field and the regular structure of the time-average component discussed above.

4.2 Magnetic field geometry and relation with the velocity field

Figures 8 and 9 suggest that the magnetic field in the ɀ = 0 plane is mostly parallel to the plane, and tends to be oriented parallel to the velocity vectors of the gas flowing in the barred potential.

We explore further the alignment between velocity and magnetic fields using their relative angle, defined as $$\cos \alpha =\frac{v\cdot B}{\left|v\right|\left|B\right|}.$$(19)

Figure 13 shows the distribution of |cos α| as a function of density, where each slice at a given density was normalised separately for clarity. When B and v are perfectly aligned, |cos α| = 1. When the orientation between B and v is completely random (i.e. B is uniformly distributed on the surface of a sphere around v), the distribution is uniform in |cos α| in the interval [0,1], with an average value of |cos α| = 0.5. The figure therefore shows that the orientation becomes progressively more random as the gas gets denser. The distribution is shown for a single snapshot, but is qualitatively similar throughout the simulation.

Figure 14 plots an xy map of |cos α| for a slice in the plane ɀ = 0. This shows that in regions of comparable densities, the alignment becomes more random where there is more turbulence. For example, the bottom panel shows significant dis-alignment in the intra-lane region outside the CMZ ring where gas is more turbulent, because no closed orbits exist for the gas to flow on, than in regions where the gas can smoothly flow on x₁/x₂ orbits (Fig. 5 in Sormani et al. 2015a). This suggests that it is the turbulence that tends to disalign the fields.

This point is supported by Fig. 15, where we plotted the distribution of |cos α| as a function of the ratio of the kinetic energy in turbulent motions to the magnetic energy E_K/E_B, which is a measure of the dynamical importance of the turbulence, in 20 pc bins (Sect. 5). In the CMZ, the magnetic field is aligned with the velocity in regions where E_B > E_K, but the distribution becomes progressively more random as the turbulence becomes more dynamically important. This is similar to the finding of Iffrig & Hennebelle (2017), with the difference that in their case the turbulence was driven by supernova feedback, while in our case it is driven by the magnetorotational instability (Sect. 5). The dust lanes behave differently than the CMZ in this respect and were excluded form the distribution shown in Fig. 15. In the dust lanes, the magnetic field stays aligned even in regions with extremely high E_K/E_B ratios. This is likely because these regions are large-scale galactic shocks, and thus have an exceptionally high shear and sharp velocity discontinuities.

Since it is the turbulence that ‘disaligns’ the velocity and magnetic fields, we would expect the regular time-averaged 〈B〉_t, in which turbulent fluctuations are averaged out, to follow even more closely the velocity streamlines. Indeed, Figs. 9 and 11 confirm this. In particular, the regular magnetic field in the CMZ ring is nearly parallel to the ring itself (Figs. 9 and 11), and the magnetic field in the ‘bar lanes’ is roughly parallel to the lanes themselves (Figs. 3, 9, and 14), and therefore to the velocity, since the latter is approximately parallel to the lanes in the frame co-rotating with the bar.

The xz and zy slices of the time-averaged regular field in Fig. 9 illustrate an interesting characteristic of the magnetic field geometry: the field transitions from mostly toroidal (i.e. along ê_ϕ) near the ɀ = 0 plane to mostly poloidal (i.e. along ê_R and ê_ɀ) at |ɀ| > 0. The transition happens through a complex ‘butterfly’ pattern, in which the field wraps around the dense gas in the ring (xz projection of the time-averaged field in Fig. 9). The transition can also be appreciated in the 3D visualisation of Fig. 12, which follows field lines as they change from toroidal just above the midplane to nearly vertical away from the plane. Figure 16 separates the magnetic field geometry into the cold (T < l0³ K) and warm (T > 10⁴ K) phases, illustrating that the vertical field above and below the plane mostly belongs to the warm diffuse phase.

Fig. 8 Visualisation of the magnetic field and column density in the CMZ in our fiducial CHEM_MHD simulation at t = 196 Myr. The со lour-scale represents the total gas column density. The line pattern indicates the orientation of the mass-weighted integrated magnetic field (i.e. 〈Bx〉ɀêx + 〈By〉ɀêy where (X)ɀ = (∫ρXdz)/(∫ ρ dz) for the xy panel, and analogous definitions for the xz and yz panels) obtained with the line integral convolution method. Movies are available online.

Fig. 9 Instantaneous (left) and time-averaged (right, Eq. (11)) density and magnetic fields. Plotted are slices in the ɀ = 0, y = 0, and x = 0 planes (not integrated quantities). Lines show the magnetic vector field. The time-average is taken between 200–250 Myr.

4.2.1 Implications for the Milky Way

The transition from horizontal (parallel to the plane) to vertical field as we move away from the mid-plane that we see in our simulations is reminiscent of the similar transition observed in the Milky Way as we move from diffuse to denser gas mentioned in Sect. 1. However, we must be careful in drawing a comparison. The geometry of the magnetic field in the CMZ is probably affected by the presence of a Galactic outflow (e.g. Ponti et al. 2021; Heywood et al. 2022), which is absent in our simulations due to the lack of stellar feedback and cosmic ray physics (e.g. Girichidis et al. 2024). Nevertheless, it is interesting to note that a transition to a perpendicular field as we leave the plane also happens independent of a Galactic outflow.

Based on our finding that the magnetic field vectors tend to be aligned with gas velocity vectors, especially in the diffuse phase, we might speculate that the vertical magnetic field lines observed in the Milky Way diffuse gas at latitudes |b| > 0.4° are tracing vertical streaming of the gas associated with the multiphase Galactic outflow (Ponti et al. 2021). We might expect that potential de-alignment due to turbulence is not happening in the diffuse gas above and below the plane, as it is sufficiently far from the midplane for turbulence driven by magnetic fields (Sect. 5) and/or stellar feedback (which predominantly occurs in the dense gas) to be ineffective.

4.2.2 Implications for external barred galaxies

Our finding that the magnetic field on large (kiloparsec) scales in the bar region is approximately aligned with the gas velocity streamlines is consistent with observations of polarised radio continuum emission of nearby barred galaxies such as NGC 1097 and NGC 1365 (Moss et al. 2001; Beck et al. 2005). The comparison of the orientation of the magnetic field in the nuclear ring is more tricky. The observed pitch angle of the magnetic field inferred from synchrotron-emitting gas in NGC 1097 is rather large, θ ~ 40°. The pitch angle of the regular 〈B〉_t component in our simulations is much smaller (Fig. 9). However, (i) the pitch angle of the instantaneous magnetic field B often appears much larger due to the presence of fluctuations perpendicular to the ring (Figs. 4 and 8); (ii) it is not clear to what extent our figures, which display the magnetic field in all gas components, are representative of synchrotron-emitting gas. A proper comparison will require synthetic observations of the synchrotron-emitting gas and a more careful comparison with observations, which is outside the scope of this paper.

The fact that the magnetic field is parallel to the bar lanes emerges spontaneously from the global flow in our simulations, and justifies the assumption of Moon et al. (2023), who injected the magnetic field parallel to the velocity vectors into the computational box in their semi-global simulations.

Fig. 10 3D visualisation of the magnetic field lines in our fiducial CHEM_MHD simulation at t = 176 Myr. The field lines are constructed by starting at a given set of points, and following the field lines until they close on themselves or leave the domain to go to infinity (field lines cannot ‘end’ within the domain because ∇ ⋅ B = 0, i.e. they behave like velocity streamlines in an incompressible flow). In the top panel we use a set of starting points distributed on a hexagonal prism centred on the Galactic centre and whose faces are located inside the gas ring, roughly midway between the centre and the CMZ gas ring. In the bottom panel we use a set of starting points located near the end of the ‘bar lanes’, just outside the CMZ ring. The blue solid surface represents an isodensity surface at n = 100 cm−3.

4.3 Magnetic field strength as a function of density

The top panel in Fig. 17 shows the strength of the B field as a function of total gas density in our fiducial CHEM_MHD simulation. We find that the magnetic field scales approximately as $$B=102\mu \text{G}{\left(\frac{n}{{10}^3{\text{cm}}^{-3}}\right)}^k,$$(20)

where k = 0.33. This scaling can be approximately understood as follows. Consider expansion or contraction of gas under the assumption of flux freezing. We can distinguish the following three limiting cases (described for example in Section 3.3.1 in the book of Kulsrud 2005): $$B\propto {\rho }^1\text{\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}the\hspace{0.17em}expansion\hspace{0.17em}or\hspace{0.17em}contraction\hspace{0.17em}is\hspace{0.17em}perpendicular\hspace{0.17em}to\hspace{0.17em}}B\text{;\hspace{0.17em}}$$(21) $$B\propto {\rho }^0\text{\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}the\hspace{0.17em}expansion/contraction\hspace{0.17em}is\hspace{0.17em}parallel\hspace{0.17em}to\hspace{0.17em}}B\text{;\hspace{0.17em}}$$(22) $$B\propto {\rho }^{\frac{2}{3}}\text{\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}the\hspace{0.17em}expansion\hspace{0.17em}or\hspace{0.17em}contraction\hspace{0.17em}is\hspace{0.17em}isotropic}.$$(23)

When В is dynamically dominant (compared to the turbulent motions that cause expansion or contraction), we expect the gas to flow preferentially parallel to B, and therefore we expect k < 2/3. In our simulations, as we will see in Sect. 5, the magnetic energy density is typically 20–40% of the turbulent kinetic energy density, and the Alfvén speed is comparable to the turbulent velocity dispersion. Thus, we are in a trans-Alfvénic non-self gravitating regime, in which magnetic fields play a non-negligible dynamical role. We therefore expect the gas to flow more readily in the direction parallel to В than perpendicular to it, and therefore we expect k < 2/3.

The exponent exhibits a small secular evolution in the simulation, changing from k ≃ 0.4 at t = 150 Myr to k ≃ 0.28 at the end of the simulation, the mean being k ≃ 0.33. This variation is somewhat smaller than what is seen in galaxy-scale simulations including supernova feedback (Konstant’mou et al. 2024), which are likely a source of higher time variation of the exponent.

The magnetic field strength as a function of density in our simulations is consistent with the rather sparse and uncertain measurements of the magnetic field strengths in the literature reported in the top panel of Fig. 17. It is also interesting to note that the power-law index of 0.33 in Eq. (20) is similar to the value of 0.4 reported by Liu et al. (2022), which was obtained by compiling polarised dust emission observations of star forming regions from the literature and computing the magnetic field strength using the Davis-Chandrasekhar-Fermi (DCF) method, although one should note that the estimated power law index has large variations depending on how the magnetic field strength was estimated and the same authors also report a larger value of 0.57 when they estimate the field differently (their Section. 3.2.1 and their Fig. 3). Finally, it is worth noting that we expect the introduction of the gas self-gravity, star formation and of the associated stellar feedback, that are switched off in our simulations, to likely affect the scaling of the magnetic field strength with density (Girichidis et al. 2018).

Fig. 11 Same as Fig. 10, but for the time-averaged magnetic and density fields shown in Fig. 9. This figure clearly illustrates that the streamlines are parallel to the bar dust lanes, and the regularity of the time-averaged field.

Fig. 12 Similar to Fig. 10, but using the time-averaged magnetic field as in Fig. 9 and a set of starting points that are located approximately 30 pc above the plane. The streamlines in this spiral up vertically, illustrating the transition from toroidal to poloidal magnetic field.

Fig. 13 Relative orientation of the gas velocity and magnetic fields as a function of total gas density for our fiducial CHEM_MHD simulation at t = 186 Myr. Colours indicate the angle |cos α| defined by Eq. (19). Each slice at fixed n is normalised separately for improved clarity. The black solid line shows the average value of cos α at the given density while the dotted one shows the median instead. If the orientation of v with respect to В were completely random (in the solid angle), the distribution would be uniform and the black lines would be horizontal at a value |cos α| = 0.5. The plot shows that the velocity and magnetic fields become progressively more random as we move to higher densities.

Fig. 14 Relative orientation of the gas velocity and magnetic fields in a slice at ɀ = 0 in our fiducial CHEM_MHD simulation at t = 186 Myr. A snapshot at early time is chosen here to show the orientation in the dust lanes as well, which are depleted at later times. The colour-scale indicates cos α defined in Eq. (19), where cos α = 1 and cos α = 0 correspond to В and v being parallel and perpendicular respectively. Streamlines with red arrows trace gas velocity in the frame co-rotating with the bar. The plot is separated into two panels for densities n > 20 cm−3 (top) and n < 20 cm−3 (bottom) for clarity. Velocity and magnetic fields are well aligned in the less turbulent regions (cos a close to 1), while their orientation becomes progressively more random in regions where the turbulence is dynamically more important.

5 Turbulence

We have already noted in Sect. 3 that the introduction of magnetic fields causes turbulence that changes the density PDF and ‘puffs up’ the gas in the ring. To quantify the magnetic-driven turbulence in more detail we covered the R < 500 pc region with non-overlapping cubical bins 20 pc on-a-side and calculated the vertical velocity dispersion in each bin, defined as $${\sigma }_{ɀ}^2=\frac{1}{N}{\displaystyle \sum _i\mathrm{\Delta }{v}_{ɀ,i}^2,}$$(24)

where N is the total number of cells in the 20 pc bin, $$\mathrm{\Delta }{v}_i=v_i-v_{\text{CM},i}$$(25)

is the velocity of the i-th cell relative to the centre of mass of the 20 pc bin (we subtracted the centre of mass velocity, since bulk motions do not contribute to turbulent kinetic energy, e.g. Stewart & Federrath 2022), and the sum extends over all cells in the bin. We used the dispersion in the ɀ direction to quantify the turbulence as it is less affected by streaming and rotational motions than the dispersion in other directions. We have also checked that once the streaming motions are taken into account the velocity dispersion is roughly isotropic in our simulations, which we find to be true within a factor of ~1.5².

Figure 18 compares σ_ɀ in the magnetised CHEM_MHD simulation and in its unmagnetised version CHEM_HD (Table 1) at t = 186 Myr. Similarly, in Fig. 19 we show a map of σ_ɀ in the CMZ for the two simulations. It is clear that the introduction of magnetic fields causes a significant increase of the velocity dispersion and of the turbulent kinetic energy. The increment is more significant in the high-density gas. For example, at bin-averaged densities of 〈n〉_bin = 10² cm⁻³ the velocity dispersion increases from σ_z ≃ 1 km s⁻¹ to σ_ɀ ≃ 5 km s⁻¹. Comparing these numbers to the sound and Alfvén speeds (dashed and dotted lines in Fig. 18) shows that the turbulence is supersonic and trans-Alfvénic.

The question arises as to what physical mechanism drives the turbulence in these simulations. We included neither the gas self-gravity nor star formation, so self-gravity and stellar feedback are ruled out as possible sources of the turbulence. The bar-driven inflow onto the CMZ can drive turbulence by converting bulk kinetic energy into turbulent motions (Kruijssen et al. 2014; Sormani & Barnes 2019; Henshaw et al. 2023). However, the unmagnetised CHEM_HD simulation, which also includes the bar-driven inflow, displays a much lower level of turbulence than the magnetised CHEM_MHD simulation. To quantify the relative importance of the bar-driven inflow on the turbulence in the CHEM_MHD simulation, we calculated the velocity dispersion in this simulation at much later times (t > 250 Myr), after the bar-driven inflow effectively shuts down because it runs out of gas (Sect. 7). We find that after the bar inflow shuts down the turbulence decreases until it settles to an intermediate value of σ_ɀ ≃ 3 km s⁻¹ at bin-averaged densities of 〈n〉_bin = 10²cm⁻³. This value is then maintained for very long times, well beyond the turbulence decaying times (or vertical crossing time). This suggests that during the time in which the bar inflow is active (as at t = 186 Myr in Fig. ) the turbulence is driven by a combination of bar inflow and magnetic processes, while the turbulence at later times (after the bar-driven inflow shuts off) is maintained by purely magnetic processes. To further investigate this we have done the following experiment (Appendix D): we have stopped the simulation at t = 166 Myr, removed all the gas at R > 500 pc so that we are left only with the CMZ gas ring, and then restarted the simulation. In this way, we remove the large-scale bar inflow and continue the simulation with only the CMZ ring evolving ‘in isolation’. We find that the turbulence settles to the same intermediate value that we find in the CHEM_MHD simulations at late times, which results in σ_ɀ 3 km s⁻¹ at bin-averaged densities of 〈n〉_bin = 10² cm⁻³. We repeated this test using an axisymmetrised potential after restarting the simulation, to exclude any possible influence of the non-axisymmetric gravitational potential, and find again that turbulence is maintained at the same level. These tests confirm that the turbulence at t = 186 Myr in Fig. 18 is driven by a combination of bar inflow and magnetic processes, while the turbulence at later times (after the bar inflow shuts down) is purely due to magnetic processes.

A well-known and effective mechanism to generate turbulence in astrophysical accretion discs is the magneto-rotational instability (MRI; Velikhov 1959; Balbus & Hawley 1991; Hawley & Balbus 1991; Balbus & Hawley 1998). This instability occurs in every (even weakly) magnetised disc in differential rotation in which the angular speed Ω.(R) decreases as a function of radius, and has been shown to work in the β < 1 limit that is relevant here (for example Kim & Ostriker 2000; Piontek & Ostriker 2007, and Appendix C of Jacquemin-Ide et al. 2021). The MRI generates turbulence by extracting the energy stored in differential rotation and converting it into turbulent fluid motions. Thus, in MRI-driven turbulence the magnetic stresses act as a mediator, allowing the turbulence to tap into the differential rotation that would otherwise not be converted into turbulent motions.

Our simulated CMZ satisfies the conditions for the onset of the MRI, and Fig. 2 shows that the MRI is well-resolved in our simulations. The MRI is therefore the most natural candidate to drive turbulence. The MHD code AREPO that we are using has been already tested to correctly reproduce the linear phase of the MRI (Pakmor & Springel 2013). We therefore conclude that the MRI (in its saturated state) is driving the turbulence in our magnetised simulations at late times (after the bar inflows shuts off).

To compare turbulent and other types of energy, we computed the turbulent kinetic, magnetic, and thermal energies in each 20 pc bin as follows: $$E_{\text{k}}={\displaystyle \sum _i\frac{1}{2}}{m}_i{\left|\mathrm{\Delta }{v}_i\right|}^2,$$(26) $$E_B={\displaystyle \underset{i}{{\sum }^{\text{}}}}{V}_i\frac{|B{|}_i^2}{8\pi },$$(27) $$E_{\text{th}}={\displaystyle \underset{i}{{\sum }^{\text{}}}}{m}_i{e}_{\text{th},i},$$(28)

where m_i is the mass of the i-th cell, B_i is its magnetic field, e_th,i its thermal energy per unit mass, V_i its volume, and the sum extends over all cells in the bin.

Figure 20 plots the energy ratios at t = 186 Myr, when turbulence is driven by a combination of bar-driven inflow and MRI. We find that (i) in high-density gas (bin-averaged density (〈n〉_bin > 1 cm⁻³), the magnetic energy is 20–40% of the turbulent kinetic energy; (ii) in low-density gas (〈n〉_bin < 0.1cm⁻³), the magnetic energy is approximately in equipartition with the thermal energy, and the turbulent energy is small. These ratios are similar to those found in studies of feedback-driven and gravity-driven turbulence, which suggest that in general E_k/E_B ≳ 2 (e.g. Federrath et al. 2011; Rieder & Teyssier 2017; Gent et al. 2021; Higashi et al. 2024). In contrast, at later times, after the bar inflow shuts off, the ratio between the turbulent kinetic energy and the magnetic energy decreases and reaches mass-weighted average values E_k/E_B < 1, which is typical of purely MRI-driven turbulence (e.g. Balbus & Hawley 1998; Kim et al. 2003; Sano et al. 2004; Minoshima et al. 2015). The decrease in ratio is mainly driven by a decrease in E_k, while E_B remains approximately constant (Sect. 6). These findings corroborate the idea that while the bar inflow is active the turbulence is driven by a combination of the inflow and the MRI, while when the bar inflow shuts off the turbulence is purely MRI-driven.

In conclusion, we have found that the combination of bar-driven inflow and MRI turbulence sustains vertical velocity dispersions that on scales of 20 pc are of the order of σ_ɀ ~ 5 km s⁻¹, while the MRI alone sustains σ_ɀ ~ 3 km s⁻¹. Both of these numbers are smaller than the σ_los ~ 10km s⁻¹ line-of-sight dispersion observed in the CMZ on the same scales (e.g. Shetty et al. 2012; Henshaw et al. 2016). This suggests that a further ingredient, likely stellar feedback such as supernovae and/or stellar winds, is necessary to explain the observed levels of turbulence in the CMZ (for example Tassis & Pavlidou (2022) and Section 4.3.4 in the review of Henshaw et al. 2023). However, we note that both the bar-driven turbulence and MRI-driven turbulence are expected to be primarily solenoidal (Gong et al. 2020), so they might be at the origin of the solenoidal driving of turbulence observed in the ‘Brick’ cloud (Federrath et al. 2016).

Fig. 15 Relative orientation of the gas velocity and magnetic fields as a function of Ek/EB averaged over 20 pc cubical bins. The green distribution in the background includes only gas in the CMZ and excludes the gas in the bar lanes. The black lines show the average (solid) and median (dotted) of the distribution. Once the turbulence becomes dynamically important, the relative orientation becomes more random. An exception is given by gas in the dust lanes, for which the average (solid grey line) and median (dotted grey line) are shown.

Fig. 16 Magnetic field in the different ISM phases. Colours represent magnetic field intensity and lines are magnetic field lines obtained with the integral convolution method. Top: all the gas. Middle: warm phase (T > 103 K). Bottom: cold phase (T < 103 K). The magnetic field in the cold gas is mostly parallel to the plane, while the vertical field above and below the plane belongs to the warm diffuse phase.

Fig. 17 Magnetic field properties as a function of density at R < 500 pc in our fiducial CHEM_MHD simulation at t = 186 Myr. Top: magnetic field strength. The В ∝ п2/3 indicates the scaling expected from simple isotropic collapse with flux-freezing (Mestel 1965). The dotted line В ∝ n0.33 is the fit to the average instead. Symbols denote various estimates of the magnetic field in the CMZ from the literature. References are abbreviated as follows: LaRosa et al. (2005, LRO5), Yusef-Zadeh et al. (2013, YZ13), Yusef-Zadeh et al. (2022, YZ22), Schwarz & Lasenby (1990, S&L90), Uchida & Guesten (1995, U&G95), Marshall et al. (1995, M95), Crutcher et al. (1996, C96), Yusef-Zadeh et al. (1996, YZ96), Yusef-Zadeh et al. (1999, YZ99), Mangilli et al. (2019, M19), Inoue et al. (1984,184), Tsuboi et al. (1986, T86), Yusef-Zadeh & Morris (1987, YZ&M87), Gray et al. (1995, G95), Pillai et al. (2015, P15). Middle: Alſvén velocity (Eq. (10)) and thermal sound speed cs = γP/ρ where P is the thermal pressure (Eq. (7)). Bottom: plasma β = P/Pmag, where Pmag = B2/(8π) is the magnetic pressure. In all panels, each density slice is normalised separately for clarity of visualisation, and the solid lines show the average value at the given density.

Fig. 18 Velocity dispersion σɀ in the ɀ direction in 20 pc on-a-side cubical bins for the magnetised CHEM_ MHD (blue) and unmagnetised CHEM_HD (red) simulations at R < 500 pc and t = 186 Myr. The simulations are identical except that the one has magnetic fields and the other does not. On the x axis is the average density within the bin. Red and blue solid lines indicate the average σɀ at the given density. Dashed lines indicate the mass-weighted averaged Alſvén speed, and dotted lines the thermal sound speed. The magnetised simulation is significantly more turbulent. The turbulence is driven by the magneto-rotational instability (Sect. 5).

6 Growth of the magnetic field

Figure 21 plots the time evolution of the volume- and mass-weighted magnetic fields in the region R < 500 pc in our fiducial simulation as a function of time. We started with an initially uniform seed field of $B_0=0.02\mu \text{G}\widehat{z}$. The field grows until it saturates at typical values B_tot ~ 200 µG. We split the plot into different density ranges, the cutoff density chosen at typical values for which the ISM becomes molecular (>10² cm⁻³). We find that saturation is reached more quickly in the dense gas, where turbulence is dynamically more important, and more slowly in the diffuse gas (top panel in Fig. 21). In Appendix В we show that the saturation field strength does not depend on the numerical resolution, while in Appendix С we show that it does not depend on the strength and orientation of the magnetic fields in the initial conditions.

Chandran et al. (2000) proposed that the magnetic field in the CMZ grows by accumulation of magnetic flux that is frozen into the bar-driven inflow and is advected into the CMZ. However, a key assumption of their model is that the magnetic field in the inflow is vertical (i.e. in the ê_ɀ direction), so that magnetic field lines get squeezed together as they move inwards, leading to the В field amplification. This assumption is invalid in our simulations because, as discussed in Sect. 4.2, the magnetic field in the bar lanes that transport the inflow is parallel to the velocity vectors, which mostly lie in the plane ɀ = 0. This implies that the magnetic field does not grow by magnetic flux accumulation via the mechanism envisioned by Chandran et al. (2000) in our simulations. To confirm this, we have performed an experiment similar to the one described in Sect. 5: we stopped the simulation at t = 166 Myr, removed all the gas at R > 500 pc so that we are left only with the CMZ ring, reset the magnetic field to the initial value $B_0=0.02\mu \text{G}\widehat{z}$ everywhere, and then restarted the simulation. In this way, we remove the bar-driven inflow and any related magnetic flux accumulation. We find that the magnetic field still grows and reaches the same saturation level as in the ‘full’ simulation (Appendix D). Repeating the test using an axisymmetrised potential after restarting the simulation leads to the same result. We conclude that the magnetic field in the CMZ does not grow by magnetic flux accumulation advected with the bar-driven inflow.

It is likely that magnetic fields in our simulations grow by dynamo action. Dynamo action can be defined as the process by which motions in the fluid amplify the magnetic fields over time. Differential rotation can amplify a toroidal magnetic field by shearing and stretching the radial field (the so-called Ω effect, for example Parker 1955; Moffatt 1978; Mestel 2012). Turbulent motions can lift the gas upwards in the plane and create and/or amplify a poloidal component from the toroidal component by inducing stretch-twist-fold motions (an example of such motions is the so called α effect, e.g. Parker 1955, 1971, 1992; Childress & Gilbert 1995). The two effects together can produce a cycle that leads to a net increase of the magnetic field intensity over time.

Differential rotation is naturally present in our simulations. Turbulence in our simulations is mostly driven by the MRI as we discussed in Sect. 5. It is therefore likely that the magnetic field in our simulations grows by a combination of Ω-dynamo induced by the differential rotation and an MRI-driven dynamo. Indeed, it is well-known that the MRI can drive dynamo action (e.g. Brandenburg et al. 1995; Stone et al. 1996; Ziegler & Rüdiger 2000; Vishniac 2009; Guan & Gammie 2011; Hawley et al. 2013; Bodo et al. 2014; Dhang et al. 2020). In an MRI-driven dynamo, magnetic fields are not only amplified by the turbulent velocity fluctuations, but they also produce the turbulent velocity field itself via the MRI (Balbus & Hawley 1998). In this aspect, MRI-driven dynamo action is different from dynamos driven by supernova feedback for example (or mean-field dynamos), where the magnetic field responds to velocity perturbations induced by something external (in this case, the supernovae). This is reflected by the fact that ratios between turbulent kinetic energy and magnetic energy are typically lower in MRI-driven dynamos than in stellar feedback-driven dynamos (discussion in Sect. 5).

In summary, dynamo action via differential rotation and MRI-driven turbulent motions is likely responsible for the growth of magnetic fields in our simulation. A more complete investigation of the dynamo action in these simulations is out of the scope of this paper but is a worthwhile direction for future studies. In particular, one could analyse the turbulence in the contexts of mean-field and MRI-driven dynamos to see which framework better describes the simulations. Whether quantities such as kinetic helicity, which is a prerequisite for the mean-field α-Ω mechanism, are the main driver for magnetic field growth here (e.g. Ntormousi et al. 2020) or if small-scale dynamo powered by MRI-driven turbulence are more important. By following high-turbulence regions as they evolve in their orbit around the Galactic centre we could understand how the energy is transferred between turbulent motion and В and vice-versa.

Fig. 19 Vertical velocity dispersion σɀ in the CMZ ring in the CHEM_MHD (top) and CHEM_HD (bottom) simulations at t = l86Myr. The velocity dispersion in the dense gas ring is much larger in the magnetised simulation than in the unmagnetised one due to the MRI-driven turbulence. Only regions with surface density Σ > 100 M⊙pc−2 are shown. The velocity dispersion is calculated using cubical bins 20 pc on-a-side as in Fig. 18.

Fig. 20 Energy ratios in 20 pc on-a-side cubical bins in our CHEM_HD (red) and CHEM_MHD (blue) simulations at R < 500 pc and t = 186 Myr. Blue is the ratio Ek/EB between turbulent kinetic energy (Eq. (26)) and magnetic energy (Eq. (27)). Red is the ratio Eth/EB between thermal energy Eth (Eq. (28)) and magnetic energy. On the x axis is the average density within the respective bin. The solid lines represent the average values at the given density. The average of the turbulent kinetic energy and thermal energy is also shown as the grey dashed line. At high densities (〈n〉bin > 1 cm−3), the kinetic energy is a few times the magnetic energy (Ek/EB ~ 2–5), while the thermal energy is negligible compared to the magnetic energy. At low densities ((〈n〉bin < 1 cm−3) the situation is inverted, with the thermal energy surpassing the magnetic energy and the turbulent kinetic energy becoming negligible at very low densities.

Fig. 21 Time evolution of the magnetic field in our fiducial CHEM_MHD simulation at R < 500 pc and |z| < 100 pc. Top: the full lines are the volume-weighted magnetic fields (Eq. (13)) averaged over the indicated density range. The black dashed line is the mass-weighted magnetic field averaged over all densities (Eq. (13)). The figure shows that the field grows in time starting from the initial seed as a result of dynamo action (Sect. 6). Bottom: integrated magnetic energy in the same density regimes. We note that while the denser part of the ISM has higher magnetic field values, most of the magnetic energy is actually in the low density regime.

Fig. 22 Mass accretion from the disc to the CMZ. Top: mass contained in the cylindrical volume R < 500 pc as a function of time for the CHEM_MHD and CHEM_HD simulations. Bottom: inflow rate into the same volume, which represents the inflow onto the CMZ ring. The black lines represent the running average (over 7 Myr) of the inflow rate, to smooth out the high time variability of Ṁ. The instantaneous (non-averaged) inflow rate is shown as the light grey lines. The inflow rate in the MHD and HD simulation has a similar order of magnitude, but it lasts longer in the MHD simulation (roughly until t ≃ 200 Myr in the HD simulation versus t 240 Myr in the MHD simulation). This is because magnetic fields transport gas radially inwards within the disc at R > 3 kpc and replenish the gas reservoir that supplies the bar-driven inflow at the outskirts of the bar, which instead runs out of gas in the purely HD simulations. Inflow at these scales is mainly driven by the gravitational torques of the Galactic bar (Sect. 7).

7 Inflow

The inflow of gas towards the centre in our simulations can be schematically divided into two regimes operating in different radial ranges which correspond to different physical driving mechanisms. These two regimes are:

1. The bar-driven inflow: from the outskirts of the bar (R ≃ 3 kpc) down to the CMZ gas ring (R ≃ 200–300 pc).

2. The nuclear inflow: from the CMZ gas ring to the central few parsecs.

Figures 22 and 23 quantify the mass inflow rates at two radii that correspond to these two regimes. Figure 22 shows that the bar-driven inflow rate has a similar magnitude (Ṁ ≃ 1 М_⊙ yr⁻¹) in both the MHD and HD simulations during the first ≃30 Myr after the bar is fully on (t = 146 Myr, Sect. 2.1), and then steadily declines. The reason for the decline is primarily that the gas reservoir located at the outskirts of the bar, which supplies the bar-driven inflow, runs out of gas. In other words, once the bar has cleared out all the inner disc region within its reach at R ≲ 3 kpc, the inflow stops. The inflow lasts a bit longer in the MHD simulation than in the purely HD simulation because the MRI-driven turbulence transports some additional gas from the outer disc at R > 3 kpc down to R ≃ 3 kpc where it can be ‘captured’ by the bar. Eventually, the bar-driven inflow runs out of gas in the MHD simulation too. This is expected since our simulations do not include the most efficient mechanisms that are believed to replenish the gas supply at the outskirts of the bar, such as raining of gas with low angular momentum from the circumgalactic medium or interactions between bar and spiral arms (e.g. Lacey & Fall 1985; Bilitewski & Schönrich 2012). Magnetic stress in the bar lanes can also enhance the bar inflow by removing angular momentum (Kim & Stone 2012), but the torques analysis below shows that this is a secondary effect and that gravitational torques dominate over Maxwell torques in this regime. Thus, the bar-driven inflow is only marginally affected by the presence of the magnetic fields.

Figure 23 shows that the nuclear inflow is practically zero in the purely HD simulation. In this simulation, all the mass is accumulating in the ring and no gas is moving further in. In contrast, the MHD simulation has a significant inflow rate of order Ṁ ≃ 0.01–0.1 M⊙ yr⁻¹ that is time-varying with a general trend upwards with time. Thus, in contrast to the bar-driven inflow that is almost unaffected by the presence of magnetic fields, the nuclear inflow changes dramatically when magnetic fields are introduced. As we shall see below, the mass transport in the nuclear inflow is driven by the MRI. The nuclear inflow is one to two orders of magnitude smaller than the bar-driven inflow, so there is a net mass accumulation in the CMZ ring. However, the numerical values depend on the numerical resolution and do not appear to be converged at the maximum resolution we can afford. In Appendix В we perform a resolution study and we show that the nuclear inflow tends to decrease as the resolution is increased. Indeed, it is well known that convergence is hard to achieve in global simulations of MRI-driven accretion discs (Hawley et al. 2013). Therefore, the nuclear inflow rates derived here should be considered upper limits.

We now analyse in more detail the physical mechanisms that drive the inflows. We start by looking at the transport of angular momentum in our simulations. Consider the cylindrical region within R = R₀, with volume V and surface S. Combining Eqs. (2) and (3) the rate of change of the total angular momentum contained within this cylindrical volume can be expressed as (Appendix A): $$\frac{\partial L_{ɀ}}{\partial t}=F_{\text{Rey\hspace{0.17em}}}+F_{\text{Max\hspace{0.17em}}}+F_{\text{Grav\hspace{0.17em}}},$$(29)

where $$L_{ɀ}={\displaystyle {\int }_V\rho }R{v}_{\varphi }\text{\hspace{0.17em}d}V,$$(30)

is the total ɀ angular momentum contained inside the cylinder, and $$F_{\text{Rey\hspace{0.17em}}}=-{\displaystyle {\int }_{S}R}\rho v_{\varphi }{v}_R\text{\hspace{0.17em}d}S,$$(31) $$F_{\text{Max}}=-{\displaystyle {\int }_{S}R}{T}_{R\varphi }\text{\hspace{0.17em}d}S,$$(32) $$F_{\text{Grav\hspace{0.17em}}}=-{\displaystyle {\int }_V\rho }\frac{\partial \mathrm{\Phi }}{\partial \varphi }\text{\hspace{0.17em}d}V,$$(33)

are the fluxes of angular momentum in and out of the cylinder, T_Rϕ = −B_ϕB_R/(4π) is the component of the Maxwell stress tensor defined in Sect. 2, dV denotes the volume element, dS is the surface area element, and (R, ϕ, ɀ) denote Galactocentric cylindrical coordinates

Equation (29) states that the change in the total angular momentum of the gas contained within the cylinder is the sum of three contributions: (i) the Reynolds flux F_Rey due to bulk motions of the gas entering or leaving the cylinder; (ii) the Maxwell flux F_Max due to magnetic forces; and (iii) the gravitational term F_Grav due to gravitational torques from the external bar potential. In Fig. 24 we performed a sanity check by calculating separately the left-hand-side (LHS) and right-hand-side (RHS) of Eq. (30). The two agree well as a function of R, which gives us confidence that our code is working correctly.

To explore the mass transport, first we decomposed the velocity field as $$v=v_0{\widehat{e}}_{\varphi }+u,$$(34)

where $$v_0(R,t)=\frac{{\langle \rho v_{\varphi }\rangle }_{\varphi z}}{\langle \rho \rangle },$$(35)

is the average mass-weighted rotation velocity, u represents the deviations and $${\langle X\rangle }_{\varphi z}=\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }\text{\hspace{0.17em}d}}\varphi {\displaystyle \int X}\text{\hspace{0.17em}d}z,$$(36)

denotes the vertical and azimuthal average of a physical quantity X. We note that u_R = υ_R by definition. We also decomposed the Reynolds flux (Eq. (31)) into an average and turbulent component: $$F_{\text{Rey\hspace{0.17em}}}=F_{\text{Rey\hspace{0.17em}}}^{\text{ave}}+F_{\text{Rey}}^{\text{trb}}$$(37)

where $$F_{\text{Rey\hspace{0.17em}}}^{\text{ave}}=-{\displaystyle {\int }_{S}R}\rho v_0{v}_R\text{\hspace{0.17em}d}S,$$(38) $$F_{\text{Rey\hspace{0.17em}}}^{\text{trb}}=-{\displaystyle {\int }_{S}R}\rho u_{\varphi }{v}_R\text{\hspace{0.17em}d}S.$$(39)

In Appendix A we show that mass accretion rates of a quasi-steady state can be understood as the sum of three contributions: turbulent Reynolds stresses, Maxwell stresses, and gravitational torques from the bar³. By determining which contribution dominates at each radius in our simulation, we can identify the physical mechanism driving the inflow.

Figure 25 shows the turbulent Reynolds, Maxwell, and gravitational stresses as a function of time for two selected radii. The top panel is for R = 3.1 kpc. On these scales, the gravitational torques dominate at all times, demonstrating that they are the ones driving the bar-driven inflow. The gravitational torques display regular oscillations on a timescale of t ≃ 7 Myr, which arise for the following reason. Let us consider a particle on a closed periodic elongated orbit in a non-axisymmetric bar potential, for example an x₁ orbit (Contopoulos & Grosbol 1989). The angular momentum is not conserved in a bar potential, so the angular momentum of the particle will oscillate with a period equal to the orbital period. During half of the orbit, the particle loses angular momentum to the bar, for half of the orbit it gains angular momentum from the bar, while there is no net gain in the long term since the orbit is periodic. The oscillations in Fig. 25 are simply a torque-weighted version of this type of orbital oscillations, averaged over all clouds at R < 3.1 kpc. These oscillations are visible because the gas is not symmetrically distributed with respect to the Galactic centre as shown in Fig. 27 (otherwise contributions on opposite sides would cancel out). These considerations explain the oscillations, but they do not explain why the average value of F_Grav is less than zero nor the average upward trend of the F_Grav. These two are explained by the fact that the gas orbits are not exactly periodic, and there is a net inflow. The upward trend in the F_Grav is because the bar-driven inflow decreases over time (Fig. 22). Eventually the gravitational torques even become slightly positive at t > 220 Myr, so the gas is gaining angular momentum from the bar. This happens because, after the bar-driven inflow shuts down, there are residual large-scale oscillations of the overall gas distribution at the outskirts of the bar from the initial gradual turn-on of the bar, which albeit slow, is not completely ‘adiabatic’. This can be seen in the bottom panel of Fig. 27, which shows that the x₁ ring is tilted so that the tips at x ≃ ±2.5 kpc are in the positive-torques quadrants.

The bottom panel in Fig. 25 is for R = 149 pc, which corresponds to regions inside the CMZ gas ring. Here, the Maxwell stresses dominate at t > 200 Myr, when the nuclear inflow is significant (bottom panel of Fig. 23). Gravitational stresses are negligible at all times at this radius. This shows that the magnetic stresses are the mechanism driving the nuclear inflow in our simulations.

Figure 26 corroborates the same conclusions by showing the stresses as a function of radius. At radii larger than that of the CMZ gas ring, the gravitational torques dominate (when the bar-driven inflow is active). Inside the CMZ ring, the magnetic stresses dominate (when the nuclear inflow is active). The turbulent Reynolds stresses are always negligible.

A natural question to ask is whether MRI-driven accretion would be the dominant inflow mechanism when further processes that are not included in the present simulations are taken into account. Tress et al. (2020) quantified the contribution of supernova feedback in driving a nuclear inflow from the CMZ inwards in simulations that included star formation and self-gravity but not magnetic fields. They found that supernova-driven turbulence can drive a nuclear inflow of approximately Ṁ ~ 0.03 M_⊙ yr⁻¹, which is highly variable in time. This is of the same order of magnitude of the MRI-driven inflow that we found here (Ṁ = 0.01–0.1 M_⊙ yr⁻¹). Thus, it is not obvious which mechanism dominates, or even if there is a single dominant mechanism. Furthermore, supernova-driven feedback and magnetic fields could interact in a non-linear way when they are both present. Moon et al. (2023) run semi-global simulations that included both supernova feedback and magnetic fields, and found that the latter can significantly enhance nuclear accretion flows compared to the supernova-only case. Understanding the dominant mechanism for the nuclear inflow will require a careful comparison that explores all relevant physical processes under a computational setup that allows a clear comparison (same gravitational potential, resolution, code).

In summary, we can clearly distinguish between two regimes in our simulations. The transport of gas from the Galactic disc to the CMZ gas ring (R ≃ 3 kpc → 300 pc) is dominated by the gravitational torques of the Galactic bar, and it is only mildly affected by the presence of magnetic fields. This is what we refer to as the ‘bar-driven inflow’. The transport of gas from the CMZ towards the innermost few parsec is entirely due to MRI-driven turbulence and is mediated by magnetic stresses. We refer to this as the ‘nuclear inflow’.

Fig. 23 Mass accretion from the CMZ inwards. Top: mass contained in the cylindrical volume R < 150 pc as a function of time for the CHEM_MHD and CHEM_HD simulations. Bottom: inflow rate into the same volume, which represents the inflow from the CMZ ring inwards. Similarly to Fig. 22, we show the time averaged inflow rate (black lines) as well as the non-averaged one (grey lines). The inflow is zero in the HD simulation, while it is significant in the MHD simulation. The spikes in the top panel for the HD simulation are caused by orbiting clouds that temporarily enter the R = 150 pc region, but that eventually leave (M at the end of the HD simulation is zero as can be seen in the top panel). Inflow at these scales is mainly driven by the magnetorota-tional instability (Sect. 7).

Fig. 24 Rate of change of the total angular momentum contained within cylindrical radius R (Eq. (29)) for the snapshot t = 196 Myr of the CHEM_MHD simulation. ALɀ/∆t is calculated by summing the angular momentum of all gas cells (Eq. (30)) in two different snapshots separated by ∆t ≃ 1 Myr, taking the difference and dividing by At. The other contributions are calculated using Eqs. (31)–(33).

Fig. 25 Turbulent Reynolds, Maxwell, and gravitational stresses at a given radius as a function of time in the CHEM_MHD simulation. Top: at R = 3.1 kpc, which tracks the entire bar region. The oscillations on timescales t ≃ 7 Myr represent a torque-weighted orbital period of clouds orbiting in the bar region (Fig. 27 and recalling that angular momentum is not conserved in a non-axisymmetric bar potential). The overall upward trend is because angular momentum is removed by the bar at a decreasing rate as the gas reservoir that supplies the bar-driven inflow is depleted and the bar-driven inflow stops (Fig. 22 and Sect. 7). Bottom: R = l49 pc, which tracks what happens inside the CMZ gas ring. Maxwell stresses become dominant at t > 200 Myr, and show that magnetic stresses are responsible for the nuclear inflow from the CMZ ring inwards.

Fig. 26 Turbulent Reynolds, Maxwell, and gravitational stresses at a given time as a function of radius for the CHEM_MHD simulation. Top: for R < 5 kpc. Bottom: for R < 200 pc. At 200 < R < 3 kpc the gravitational stresses dominate, and therefore the bar-driven inflow is driven by the gravitational torques of the Galactic bar (as the name would imply). At 0 < R < 200 pc the Maxwell stress dominate, and therefore the nuclear inflow is driven by magnetic fields. We note that the top and bottom panels are for different snapshots, which are selected to better capture the moments the bar-driven and nuclear inflows take place.

8 Summary and conclusions

We investigated the impact and properties of the magnetic fields in the central regions of the Milky Way using 3D magnetohy-drodynamical simulations of non-self gravitating gas flowing in an externally imposed barred potential. We found the following results:

8.1 Gas morphology

• Magnetic pressure tends to increase the effective sound speed of the gas, and to decrease the radius of the nuclear ring (Sect. 3).

• Magnetically-driven turbulence puffs the gas up and increases the scale-height compared to non-magnetised simulations (Sect. 3).

• Magnetically-driven accretion fills the region inside the CMZ gas ring with gas, which would remain devoid of gas in the absence of magnetic fields or other physical processes not included in the present simulations such as supernova feedback (Sect. 3).

Fig. 27 Mass-weighted gravitational torques map for two selected snapshots of our CHEM_MHD simulation. The bar removes (adds) angular momentum in the two red (green) quadrants.

8.2 Magnetic field properties

• The magnetic field can be conveniently decomposed into a regular time-averaged component and a random turbulent component (Sect. 4.1).

• The regular component is generally well aligned with the gas velocity vectors of gas flowing in a bar potential. In particular, the magnetic field in the bar lanes that transport the gas from the Galactic disc to the CMZ is parallel to the lanes. Turbulence tends to ‘disalign’ the magnetic and velocity fields (Sect. 4.2).

• The field geometry transitions from toroidal near the ɀ = 0 plane to poloidal at |z| > 0. The transition happens through a complex ‘butterfly’ pattern (Sect. 4.2).

• The magnetic field scales as a function of density as В ≃ 102(n/10³cm⁻³)^{0 33} µG. This can be explained by the magnetic field playing a non-negligible dynamical role (Sect. 4.3).

8.3 Turbulence

• The combination of bar inflow and magneto-rotational instability (MRI) drives turbulence in the CMZ and can maintain a velocity dispersion of the order σ ~ 5 km s⁻¹ on a scale of 20 pc. The MRI alone in the absence of bar inflow maintains σ_ɀ ~ 3 km s⁻¹. Both these values are lower than the velocity dispersion observed in the CMZ on the same scale, indicating that magnetic fields and bar driven inflow alone cannot drive the full amount of turbulence observed in the CMZ. Stellar feedback is likely the missing ingredient that is necessary to fully explain the observed turbulence (Sect. 5).

• When turbulence in the CMZ is driven by both the bar inflow and the MRI, the ratio between the turbulent kinetic and magnetic energy is E_k/E_B ≳ 2, similar to the value found in studies of stellar feedback-driven and gravity-driven turbulence. When the turbulence is driven by the MRI alone, this value decreases to E_k/E_B < 1, similar to studies of purely MRI-driven turbulence (Sect. 5).

8.4 Growth of the magnetic field

• Magnetic fields grow in our simulations because of dynamo action driven by a combination of MRI and differential rotation, until they saturate at a mass-weighted value in the CMZ ring of approximately B_tot ≃ 200 µG.

• The saturation value is not sensitive to the initial strength and orientation of the magnetic field, and does not change significantly if we increase the resolution of the simulations (Sect. 6).

8.5 Inflows

• We can clearly distinguish two inflow regimes acting in different radial ranges in our simulation: (1) The bar-driven inflow that transports the gas from the Galactic disc (R ≃ 3 kpc) down to the CMZ gas ring (R = 200–300 pc); (2) The nuclear inflow that transports the gas from the CMZ inwards towards the central few parsec (Sect. 7).

• The bar-driven inflow is driven by the gravitational torques of the Galactic bar and is only marginally influenced by the presence of magnetic fields. The inflow rate is of the order Ṁ ≃ l M_⊙ yr⁻¹ (Sect. 7).

• The nuclear inflow is driven by magnetic stresses in our simulations. The inflow rate is of the order of M = 0.01−0.1 M_⊙ yr⁻¹. This suggests that MRI-driven transport is a viable mechanism to transport gas to the Nuclear Star Cluster (NSC) that will contribute to its in-situ star formation. A resolution study shows that the nuclear inflow rate decreases with increasing numerical resolution, and our simulations do not appear to be converged at the maximum resolution we can afford. The above values should therefore be considered as upper limits (Sect. 7).

Data availability

Movies associated with Figs. 3, 4, 8 are available at https://www.aanda.org

Acknowledgements

The authors acknowledge support by the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant INST 35/1597-1 FUGG and DFG grant INST 35/1134-1 FUGG. M.C.S. acknowledges financial support from the European Research Council under the ERC Starting Grant “GalFlow” (grant 101116226) and from the Royal Society (URF\R1\221118). J.D.H. gratefully acknowledges financial support from the Royal Society (University Research Fellowship; URF/R1/221620)A.S.M. acknowledges support from the RyC2021-032892-I grant funded by MCIN/AEI/10.13039/501100011033 and by the European Union ‘Next GenerationEU’/PRTR, as well as the program Unidad de Excelencia María de Maeztu CEX2020-001058-M. P.G. was supported by the Chinese Academy of Sciences (CAS), through a grant to the CAS South America Center for Astronomy (CASSACA) in Santiago, Chile, and by the “Comisión Nacional de Ciencia y Tecnología (CONICYT)” now ANID via Project FONDECYT de Iniciación 11170551. V.M.R. also acknowledges support from the grant number RYC2020-029387-I funded by MICIU/AEI/10.13039/501100011033 and by “ESF, Investing in your future”, and from the Consejo Superior de Investi-gaciones Científicas (CSIC) and the Centro de Astrobiología (CAB) through the project 20225AT015 (Proyectos intramurales especiales del CSIC). L.C. and V.M.R. acknowledge funding from grants No. PID2019-105552RB-C41 and PID2022-136814NB-I00 by the Spanish Ministry of Science, Innovation and Universities/State Agency of Research MICIU/AEI/10.13039/501100011033 and by ERDF, UE. R.S.K., S.C.O.G. and P.G. gratefully acknowledge support from the ERC via the ERC Synergy Grant “ECOGAL” (grant 855130), from the German Excellence Strategy via the Heidelberg Cluster of Excellence (EXC 2181 – 390900948) “STRUCTURES”, and from the German Ministry for Economic Affairs and Climate Action in project “MAINN” (funding ID 50OO2206). E.S. acknowledges support from the Marie Skłodowska-Curie Grant 101061217. J.M.D.K. gratefully acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme via the ERC Starting Grant MUSTANG (grant agreement number 714907). COOL Research DAO is a Decentralised Autonomous Organisation supporting research in astrophysics aimed at uncovering our cosmic origins. J.G. is a member of the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD) and acknowledges financial support by the European Research Council via the ERC Synergy Grant “ECOGAL” (project ID 855130). A.G. acknowledges support from the NSF under grants AAG 2206511 and CAREER 2142300. C.B. gratefully acknowledges funding from National Science Foundation under Award Nos. 1816715, 2108938, 2206510, and CAREER 2145689, as well as from the National Aeronautics and Space Administration through the Astrophysics Data Analysis Program under Award No. 21-ADAP21-0179 and through the SOFIA archival research program under Award No. 09_0540.

Appendix A Conservation of angular momentum and mass transport

In this appendix we derive equations for the transport of angular momentum and mass. Our treatment mostly follows Balbus & Papaloizou (1999) and Moon et al. (2023).

An equation for the angular momentum transport can be obtained from Eq. (3). Multiplying the azimuthal component of Eq. (3) by R using standard cylindrical coordinates (R, ϕ, z) and rearranging gives: $$\frac{\partial \left(l_z\right)}{\partial t}+\nabla \cdot F_J=-\rho \frac{\partial \mathrm{\Phi }}{\partial \varphi },$$(A.1)

where $$l_z=\rho R{v}_{\varphi },$$(A.2) $$F_J=R\left(\rho v_{\varphi }v+P{\widehat{e}}_{\varphi }-\frac{B_{\varphi }}{4\pi }B+\frac{B^2}{8\pi }{\widehat{e}}_{\varphi }\right).$$(A.3)

The quantity l_ɀ is the angular momentum per unit volume, while F_J is the flux of angular momentum, which is the sum of contributions due to bulk motions of the gas, pressure forces, and magnetic forces. The term ρ∂Φ/∂ϕ is a source term representing the changes in angular momentum due to torques from the external potential. When ρ∂Φ/∂ϕ = 0, Eq. (A.1) implies that the total angular momentum is conserved. Indeed, the only agent that can change the total angular momentum in our simulations is the external bar potential. We note that the angular momentum per unit volume l_ɀ does not contain any contribution due to magnetic fields. This is because, although electromagnetic fields can in general carry angular momentum, this contribution is neglected in ideal MHD since it is of order (υ/c)² with respect to the angular momentum contained in the gas.

The angular momentum flux in the radial direction is: $$F_{JR}=R\left(\rho v_{\varphi }{v}_R+T_{R\varphi }\right),$$(A.4)

where $$T_{R\varphi }=-\frac{B_{\varphi }{B}_R}{4\pi },$$(A.5)

is the component of the Maxwell stress tensor defined in Sect. 2. When (A.4) is integrated over the surface of a cylinder of radius R, it gives the flux of angular momentum through the surface of the cylinder. The term ρυ_ϕυ_R quantifies the angular momentum flux carried by gas that is physically crossing the surface R, while T_Rϕ is the contribution due to magnetic torques. Equation (29) in the main text is obtained by integrating Eq. (A.1) over the volume of a cylinder of radius R₀ and using the divergence theorem.

We now relate these quantities to the mass accretion rate. Expanding Eq. (A.1) in cylindrical coordinates, integrating in both the vertical (ɀ) and azimuthal (ϕ) direction, and assuming that the boundary terms in the vertical direction vanish, we obtain: $$\frac{\partial \langle l_z\rangle }{\partial t}+\frac{1}{R}\frac{\partial \langle R{F}_{JR}\rangle }{\partial R}=-\langle \rho \frac{\partial \mathrm{\Phi }}{\partial \varphi }\rangle ,$$(A.6)

where $$\langle X\rangle =\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }\text{\hspace{0.17em}d}}\varphi {\displaystyle \int X}\text{\hspace{0.17em}d}z,$$(A.7)

denotes the vertical and azimuthal average of a physical quantity X. Integrating and averaging the continuity equation (2) in the same way gives: $$\frac{\partial \langle \rho \rangle }{\partial t}+\frac{1}{R}\frac{\partial }{\partial R}\left[\langle R\rho v_R\rangle \right]=0.$$(A.8)

Next, we decompose the velocity as $$v=v_0{\widehat{e}}_{\varphi }+u,$$(A.9)

where $$v_0(R,t)=\frac{\langle \rho v_{\varphi }\rangle }{\langle \rho \rangle },$$(A.10)

is the average rotation velocity and u represents the fluctuations. Substituting (34) into (A.6), using (A.4), (A.1), and the averaged continuity equation (A.8), and using the fact that 〈ρu_ϕ〉 = 0 by our definition of υ₀, we obtain $$R\langle \rho \rangle \frac{\partial v_0}{\partial t}+\langle \rho u_R\rangle \frac{\partial \left(R{v}_0\right)}{\partial R}+\frac{1}{R}\frac{\partial \langle R^2\rho u_R{u}_{\varphi }\rangle }{\partial R}+\frac{1}{R}\frac{\partial \langle R^2{T}_{R\varphi }\rangle }{\partial R}=-\langle \rho \frac{\partial \mathrm{\Phi }}{\partial \varphi }\rangle ,$$(A.11)

Equation (A.11) can be rewritten in a more illuminating form as: $$\dot{M}={\dot{M}}_{\text{M}}+{\dot{M}}_{\text{R}}+{\dot{M}}_{\text{G}}+{\dot{M}}_t,$$(A.12)

where $$\dot{M}=-2\pi R\langle \rho u_R\rangle ,$$(A.13)

is the total mass accretion rate at radius R, and $${\dot{M}}_{\text{R}}=2\pi {\left[\frac{\partial \left(R{v}_0\right)}{\partial R}\right]}^{-1}\frac{\partial \langle R^2\rho u_R{u}_{\varphi }\rangle }{\partial R},$$(A.14) $${\dot{M}}_{\text{M}}=2\pi {\left[\frac{\partial \left(R{v}_0\right)}{\partial R}\right]}^{-1}\frac{\partial \langle R^2{T}_{R\varphi }\rangle }{\partial R},$$(A.15) $${\dot{M}}_{\text{G}}=2\pi R{\left[\frac{\partial \left(R{v}_0\right)}{\partial R}\right]}^{-1}\langle \rho \frac{\partial \mathrm{\Phi }}{\partial \varphi }\rangle ,$$(A.16) $${\dot{M}}_t=2\pi R^2{\left[\frac{\partial \left(R{v}_0\right)}{\partial R}\right]}^{-1}\langle \rho \rangle \frac{\partial v_0}{\partial t}.$$(A.17)

The first three terms are the contributions due to Reynolds, Maxwell, and gravitational stresses respectively, while the last term is related to changes in the average rotation velocity with time, which are typically small and can be neglected in a quasi-steady state.

Appendix B Resolution study

We perform a resolution study to see how the magnetic field strength at saturation and the inflow rate depend on the numerical resolution. Figure B.1 shows that the saturation magnetic field strengths are relatively insensitive to numerical resolution. In contrast, the MRI-driven inflow rate decreases as the resolution is increased, and does not appear to be converged even at the highest resolution we can afford. We can therefore only put an upper limit to the MRI-driven inflow rate.

Fig. B.1 Impact of resolution on magnetic fields and inflow. Top: volume-weighted average magnetic field intensity in the CMZ (R < 500 pc, |z| < 100 pc) in three different density regimes at different resolutions as a function of time. Bottom: mass accretion as a function of time from the CMZ inwards (R < 100 pc) for the three simulations at different resolutions. The saturation levels of the magnetic field are relatively insensitive to the resolution, while the inflow rate is highly sensitive to it.

Appendix C Impact of initial conditions

We tested the impact of the initial conditions on the magnetic field evolution. Figure C.1 compares the magnetic field strengths and inflow rates in the CMZ as a function of time for simulations that start with a uniform Poloidal seed magnetic field B₀ = 0.02 µGê_ɀ (as in the main text) versus a purely Toroidal magnetic field B₀ = 0.02 µGê_ϕ. We find that the magnetic field configuration and inflow rates are relatively insensitive to the initial orientation of the magnetic field. This is likely because the CMZ builds its gas and magnetic field reservoirs by bar-driven accretion during the course of the simulation, and memory of the initial conditions is erased in the process.

Fig. C.1 Impact of the initial conditions on magnetic fields and inflow. Top: volume-weighted average magnetic field intensity in the CMZ (R < 500 pc, |ɀ| < 100 pc) in three different density regimes for poloidal initial magnetic field and toroidal initial magnetic field. Both tests are run at a resolution of 100 M⊙ per cell. Bottom: mass inflow rate as a function of time from the CMZ inwards (R < 100 pc) in the two simulations. The figure shows that the magnetic field intensity and inflow rates are relatively insensitive to the initial orientation of the magnetic field.

Appendix D Cutout tests

To assess the importance of the bar-driven inflow on the turbulence and magnetic field properties in the CMZ we performed the following numerical experiments. In the first experiment we stop the simulation at t = 166 Myr, remove all the gas at R > 500 pc so that we are left only with the CMZ gas ring, and then restart the simulation. In this way, we remove the large-scale bar inflow and continue the simulation with only the CMZ ring evolving ‘in isolation’. In the second experiment, we do the same and in addition we also axisymmetrised the gravitational potential and reset the magnetic field to the initial seed value B₀ = 0.02 µGê_z. Figure D.1 plots the evolution of the magnetic field and the inflow rates in these experiments. We find that the magnetic field intensity and inflow rates are not affected in the first experiment, and in the second both the magnetic field intensity and inflow rate eventually ramp up to the same values they had in the other simulations. We therefore conclude that the bar-driven inflow has little or no influence on the magnetic field intensity levels at saturation in the CMZ, nor on the inflow rate. These tests are consistent with our suggestion that the MRI is primarily responsible for driving the nuclear inflow within the CMZ ring in our simulations (Sect. 7).

Fig. D.1 Effect of the cutout test on magnetic fields and inflow. Top: volume-weighted average magnetic field intensity in the CMZ in three different density regimes for three simulaitions: (i) the fiducial CHEM_MHD simulation; (ii) the test where we cut out the CMZ; (iii) and the test where we also reset the B field (Appendix D for description). The latter was run until 400 Myr to give time for dynamo processes and MRI to develop again. Bottom: mass accretion as a function of time from the CMZ inwards (R < 100 pc) for the three simulations. The vertical dashed line indicates the time at which we performed the cut-out.

Appendix E Stokes parameters

The interstellar magnetic field is mostly measured indirectly. Dust polarisation observations can constrain the average magnetic field orientation along the line-of-sight, but does not tell us much about the magnetic field strength. Moreover the direction of the polarised emission does not follow simple rules of vector addition along the line-of-sight as the magnetic field does, but instead it is necessary to compute and integrate the Stokes parameters.

To get a projected map of the magnetic field orientation that can be more easily related to observed dust polarisation maps, we computed mass-weighted normalised stokes parameters (for instance eqs. 15, 16, and 17 of Soler et al. 2013): $$Q_{xy}=\frac{B_x^2-B_y^2}{B_x^2+B_y^2},$$(E.1) $$U_{xy}=\frac{2B_x{B}_y}{B_x^2+B_y^2}.$$(E.2)

The mass weighted average magnetic field strength and its orientation are then $${\langle \left|B^{\prime }\right|\rangle }_z=\sqrt{{\langle Q_{xy}\rangle }_z^2+{\langle U_{xy}\rangle }_z^2}$$(E.3)

and $${\mathrm{\Phi }}_{{\langle \left|B^{\prime }\right|\rangle }_z}=\frac{1}{2}\arctan \left(\frac{{\langle U_{xy}\rangle }_z}{{\langle Q_{xy}\rangle }_z}\right),$$(E.4)

where 〈X〉_ɀ = (∫ ρX dɀ) / (∫ ρ dɀ). A Stokes parameter map of our simulations is shown in Fig. E.1. We find only relatively minor differences with the direct B field maps in Fig. 8.

Fig. E.1 Same as Fig. 8 but here we integrate the Stokes parameters along the line-of-sight instead of directly integrating the magnetic field.

Appendix F Divergence of the magnetic field

Discretization of the MHD equations introduces divergence errors of the magnetic field which are obviously spurious. These can quickly grow and dominate the dynamics of the gas. Various methods have been developed to keep them at bay. We use the Powell divergence cleaning technique (Powell et al. 1999), implemented in AREPO by Pakmor & Springel (2013). Here, source terms proportional to ∇ ⋅ B are introduced in the equations such that the generated divergence is advected away, preventing it from growing further. This scheme is robust, local, and does not require additional restrictions on the time step. Moreover, it can be easily implemented on the unstructured moving mesh of the AREPO code. However, since the scheme is not divergence-free by construction, it is recommended to check the impact of these errors in the simulations. The implementation did not produce any obvious artefacts in test problems as well as in more complex and generic astrophysical simulations (Pakmor & Springel 2013; Pakmor et al. 2020) and the relative divergence error was contained provided high enough resolution.

We define the relative error of the divergence as $$\frac{(\nabla \cdot B)r}{\left|B\right|},$$(F.1)

where r is the radius of the cell for which the relative error was computed.

In Fig. F.1 we show the relative divergence error map for a slice through the midplane at t = 196 Myr. In the CMZ the divergence error is always small (of the order of a few percent). Cells with higher divergence errors are present mostly in the low density region surrounding the CMZ or at the interfaces between low and high density regions. Here the magnetic field gradients are highest and greater divergence errors are to be expected.

In the top panel, the divergence error is shown including its sign. There, we see that no large patches of only positive (or negative) divergence errors are present. This means that even in regions with higher divergence errors, the error mostly averages out on larger scales, such that the dynamics of scales including more than just a few cells are properly followed.

In the bottom panel, the absolute value in log scales is shown instead. We see that even for regions of relatively high divergence error, the error never grows excessively.

In Fig. F.2, we show the volume averaged divergence error in the central 500 pc as a function of time. When averaged considering its sign (top panel), the error is always much less than 1%, and is not diverging in time, but oscillates around the zero value.

Even if the average is performed by taking the absolute value of the relative divergence error (bottom panel), it is contained to a few percent.

In Fig. F.2 we also included the values computed for simulations at lower resolutions. As expected, the divergence error reduces with resolution, and even at the lowest resolution the divergence error is never significant.

Fig. F.1 Map of the relative divergence error of the magnetic field for a slice through the ɀ = 0 plane at t = 196 Myr. Top: the divergence error map is shown considering its sign as well. Bottom: map of the absolute value of the divergence error in log scale. The contour corresponds to an isodensity contour of n = 20 cm−3.

Fig. F.2 Relative divergence error of the magnetic field, volume-averaged over the central 500 pc, as a function of time for the different simulations at different resolution. Top: the divergence error is shown considering its sign as well. Bottom: volume-average of the absolute value of the divergence error in log scale.

References

All Tables

All Figures

	Fig. 1 Gas temperature as a function of total gas density in our fiducial CHEM_MHD model at t = 196 Myr. Blue and red indicate the regions at R > 0.5 kpc and R < 0.5 kpc, respectively. Contours contain [99, 95. medium with a cold (T = 102 K) and warm (T = 104 K) phase. The hot phase (T = 106 K) is absent since our simulations do not include stellar feedback or other processes that can create it.
In the text

Fig. 2 Resolution in our fiducial CHEM_MHD simulation (Table 1) at t = 196 Myr. Top: spatial resolution as a function of total gas density. rcell is the radius of a sphere with the same volume as the cell. The horizontal black dashed line indicates the volume limit of the cells (Sect. 2.2). Middle: mass of cells as a function of density. In both panels, blue contours are for the region at R > 0.5 kpc with a target mass of 100 M⊙, while red contours are for the R < 0.5 kpc with higher resolution at a target mass of 10 M⊙ (Sect. 2.2). Bottom: λMRI/rcell as a function of density for cells in the CMZ (R < 0.5 kpc), where λMRI = 2πυA/Ω. is the characteristic wavelength of the MRI, υA is the Alfven speed, and Ω = υϕ/R is the angular velocity of each gas cell. We see that the MRI is well resolved in the CMZ. The contours contain [100, 99, 90, 75, 50, 25] % of the total number of cells.

In the text

Fig. 3 State of the system at t = 196 Myr for our set of simulations. Top and middle: face-on gas surface density for all the simulations listed in Table 1. Bottom: projected mass-weighted magnetic field $\left|{\langle B\rangle }_z\right|={\left({\langle B_x\rangle }_z^2+{\langle B_y\rangle }_z^2+{\langle B_z\rangle }_z^2\right)}^{1/2}$ for the three magnetised simulations in Table 1, where 〈X〉ɀ = (∫ ρX dɀ)/(∫ ρ dɀ) denotes the vertically integrated mass-weighted density of a quantity X. Contours denote the streamlines of the projected field ${\langle B_x\rangle }_z{\widehat{e}}_x+{\langle B_y\rangle }_z{\widehat{e}}_y$. The dashed circles indicate the region with higher numerical resolution at R < 500 pc (Sect. 2.2). Rotation is clockwise. Movies are available online.

In the text

	Fig. 4 Same as Fig. 3, but zooming-in onto the Central Molecular Zone. Movies are available online.
In the text

	Fig. 5 Time evolution of the surface gas density in the CHEM_MHD simulation. Rotation is clockwise.
In the text

	Fig. 6 Time evolution of the gas density in the central regions of the CHEM_MHD simulations. This is the same as Fig. 5, but zooming-in in the central region. The interior of the CMZ gas ring is empty at t = 147 Myr and then is gradually filled with gas due to MRI-driven accretion (Sect. 7).
In the text

Fig. 7 Mass-weighted probability density function (PDF) of the CHEM_MHD and CHEM_HD simulations (Table 1) at t = 196 Myr. The thick lines are for the CMZ (R < 500 pc), while the thin lines are for the entire simulation. Each curve is normalised so that ∫(dM/d log n) d log n = 1. In the CHEM_HD simulation, most of the gas in the ring accumulates at the highest densities purely due to orbital convergence. This is not due to the gas self-gravity, since it is switched off in our simulations. The presence of magnetic fields shifts the peak to lower densities by providing pressure support and driving random motions. The total distribution is bimodal, as expected from a two-phase medium (e.g. Seta & Federrath 2022).

In the text

Fig. 8 Visualisation of the magnetic field and column density in the CMZ in our fiducial CHEM_MHD simulation at t = 196 Myr. The со lour-scale represents the total gas column density. The line pattern indicates the orientation of the mass-weighted integrated magnetic field (i.e. 〈Bx〉ɀêx + 〈By〉ɀêy where (X)ɀ = (∫ρXdz)/(∫ ρ dz) for the xy panel, and analogous definitions for the xz and yz panels) obtained with the line integral convolution method. Movies are available online.

In the text

	Fig. 9 Instantaneous (left) and time-averaged (right, Eq. (11)) density and magnetic fields. Plotted are slices in the ɀ = 0, y = 0, and x = 0 planes (not integrated quantities). Lines show the magnetic vector field. The time-average is taken between 200–250 Myr.
In the text

Fig. 10 3D visualisation of the magnetic field lines in our fiducial CHEM_MHD simulation at t = 176 Myr. The field lines are constructed by starting at a given set of points, and following the field lines until they close on themselves or leave the domain to go to infinity (field lines cannot ‘end’ within the domain because ∇ ⋅ B = 0, i.e. they behave like velocity streamlines in an incompressible flow). In the top panel we use a set of starting points distributed on a hexagonal prism centred on the Galactic centre and whose faces are located inside the gas ring, roughly midway between the centre and the CMZ gas ring. In the bottom panel we use a set of starting points located near the end of the ‘bar lanes’, just outside the CMZ ring. The blue solid surface represents an isodensity surface at n = 100 cm−3.

In the text

	Fig. 11 Same as Fig. 10, but for the time-averaged magnetic and density fields shown in Fig. 9. This figure clearly illustrates that the streamlines are parallel to the bar dust lanes, and the regularity of the time-averaged field.
In the text

	Fig. 12 Similar to Fig. 10, but using the time-averaged magnetic field as in Fig. 9 and a set of starting points that are located approximately 30 pc above the plane. The streamlines in this spiral up vertically, illustrating the transition from toroidal to poloidal magnetic field.
In the text

Fig. 13 Relative orientation of the gas velocity and magnetic fields as a function of total gas density for our fiducial CHEM_MHD simulation at t = 186 Myr. Colours indicate the angle |cos α| defined by Eq. (19). Each slice at fixed n is normalised separately for improved clarity. The black solid line shows the average value of cos α at the given density while the dotted one shows the median instead. If the orientation of v with respect to В were completely random (in the solid angle), the distribution would be uniform and the black lines would be horizontal at a value |cos α| = 0.5. The plot shows that the velocity and magnetic fields become progressively more random as we move to higher densities.

In the text

Fig. 14 Relative orientation of the gas velocity and magnetic fields in a slice at ɀ = 0 in our fiducial CHEM_MHD simulation at t = 186 Myr. A snapshot at early time is chosen here to show the orientation in the dust lanes as well, which are depleted at later times. The colour-scale indicates cos α defined in Eq. (19), where cos α = 1 and cos α = 0 correspond to В and v being parallel and perpendicular respectively. Streamlines with red arrows trace gas velocity in the frame co-rotating with the bar. The plot is separated into two panels for densities n > 20 cm−3 (top) and n < 20 cm−3 (bottom) for clarity. Velocity and magnetic fields are well aligned in the less turbulent regions (cos a close to 1), while their orientation becomes progressively more random in regions where the turbulence is dynamically more important.

In the text

Fig. 15 Relative orientation of the gas velocity and magnetic fields as a function of Ek/EB averaged over 20 pc cubical bins. The green distribution in the background includes only gas in the CMZ and excludes the gas in the bar lanes. The black lines show the average (solid) and median (dotted) of the distribution. Once the turbulence becomes dynamically important, the relative orientation becomes more random. An exception is given by gas in the dust lanes, for which the average (solid grey line) and median (dotted grey line) are shown.

In the text

Fig. 16 Magnetic field in the different ISM phases. Colours represent magnetic field intensity and lines are magnetic field lines obtained with the integral convolution method. Top: all the gas. Middle: warm phase (T > 103 K). Bottom: cold phase (T < 103 K). The magnetic field in the cold gas is mostly parallel to the plane, while the vertical field above and below the plane belongs to the warm diffuse phase.

In the text

Fig. 17 Magnetic field properties as a function of density at R < 500 pc in our fiducial CHEM_MHD simulation at t = 186 Myr. Top: magnetic field strength. The В ∝ п2/3 indicates the scaling expected from simple isotropic collapse with flux-freezing (Mestel 1965). The dotted line В ∝ n0.33 is the fit to the average instead. Symbols denote various estimates of the magnetic field in the CMZ from the literature. References are abbreviated as follows: LaRosa et al. (2005, LRO5), Yusef-Zadeh et al. (2013, YZ13), Yusef-Zadeh et al. (2022, YZ22), Schwarz & Lasenby (1990, S&L90), Uchida & Guesten (1995, U&G95), Marshall et al. (1995, M95), Crutcher et al. (1996, C96), Yusef-Zadeh et al. (1996, YZ96), Yusef-Zadeh et al. (1999, YZ99), Mangilli et al. (2019, M19), Inoue et al. (1984,184), Tsuboi et al. (1986, T86), Yusef-Zadeh & Morris (1987, YZ&M87), Gray et al. (1995, G95), Pillai et al. (2015, P15). Middle: Alſvén velocity (Eq. (10)) and thermal sound speed cs = γP/ρ where P is the thermal pressure (Eq. (7)). Bottom: plasma β = P/Pmag, where Pmag = B2/(8π) is the magnetic pressure. In all panels, each density slice is normalised separately for clarity of visualisation, and the solid lines show the average value at the given density.

In the text

Fig. 18 Velocity dispersion σɀ in the ɀ direction in 20 pc on-a-side cubical bins for the magnetised CHEM_ MHD (blue) and unmagnetised CHEM_HD (red) simulations at R < 500 pc and t = 186 Myr. The simulations are identical except that the one has magnetic fields and the other does not. On the x axis is the average density within the bin. Red and blue solid lines indicate the average σɀ at the given density. Dashed lines indicate the mass-weighted averaged Alſvén speed, and dotted lines the thermal sound speed. The magnetised simulation is significantly more turbulent. The turbulence is driven by the magneto-rotational instability (Sect. 5).

In the text

Fig. 19 Vertical velocity dispersion σɀ in the CMZ ring in the CHEM_MHD (top) and CHEM_HD (bottom) simulations at t = l86Myr. The velocity dispersion in the dense gas ring is much larger in the magnetised simulation than in the unmagnetised one due to the MRI-driven turbulence. Only regions with surface density Σ > 100 M⊙pc−2 are shown. The velocity dispersion is calculated using cubical bins 20 pc on-a-side as in Fig. 18.

In the text

Fig. 20 Energy ratios in 20 pc on-a-side cubical bins in our CHEM_HD (red) and CHEM_MHD (blue) simulations at R < 500 pc and t = 186 Myr. Blue is the ratio Ek/EB between turbulent kinetic energy (Eq. (26)) and magnetic energy (Eq. (27)). Red is the ratio Eth/EB between thermal energy Eth (Eq. (28)) and magnetic energy. On the x axis is the average density within the respective bin. The solid lines represent the average values at the given density. The average of the turbulent kinetic energy and thermal energy is also shown as the grey dashed line. At high densities (〈n〉bin > 1 cm−3), the kinetic energy is a few times the magnetic energy (Ek/EB ~ 2–5), while the thermal energy is negligible compared to the magnetic energy. At low densities ((〈n〉bin < 1 cm−3) the situation is inverted, with the thermal energy surpassing the magnetic energy and the turbulent kinetic energy becoming negligible at very low densities.

In the text

Fig. 21 Time evolution of the magnetic field in our fiducial CHEM_MHD simulation at R < 500 pc and |z| < 100 pc. Top: the full lines are the volume-weighted magnetic fields (Eq. (13)) averaged over the indicated density range. The black dashed line is the mass-weighted magnetic field averaged over all densities (Eq. (13)). The figure shows that the field grows in time starting from the initial seed as a result of dynamo action (Sect. 6). Bottom: integrated magnetic energy in the same density regimes. We note that while the denser part of the ISM has higher magnetic field values, most of the magnetic energy is actually in the low density regime.

In the text

Fig. 22 Mass accretion from the disc to the CMZ. Top: mass contained in the cylindrical volume R < 500 pc as a function of time for the CHEM_MHD and CHEM_HD simulations. Bottom: inflow rate into the same volume, which represents the inflow onto the CMZ ring. The black lines represent the running average (over 7 Myr) of the inflow rate, to smooth out the high time variability of Ṁ. The instantaneous (non-averaged) inflow rate is shown as the light grey lines. The inflow rate in the MHD and HD simulation has a similar order of magnitude, but it lasts longer in the MHD simulation (roughly until t ≃ 200 Myr in the HD simulation versus t 240 Myr in the MHD simulation). This is because magnetic fields transport gas radially inwards within the disc at R > 3 kpc and replenish the gas reservoir that supplies the bar-driven inflow at the outskirts of the bar, which instead runs out of gas in the purely HD simulations. Inflow at these scales is mainly driven by the gravitational torques of the Galactic bar (Sect. 7).

In the text

Fig. 23 Mass accretion from the CMZ inwards. Top: mass contained in the cylindrical volume R < 150 pc as a function of time for the CHEM_MHD and CHEM_HD simulations. Bottom: inflow rate into the same volume, which represents the inflow from the CMZ ring inwards. Similarly to Fig. 22, we show the time averaged inflow rate (black lines) as well as the non-averaged one (grey lines). The inflow is zero in the HD simulation, while it is significant in the MHD simulation. The spikes in the top panel for the HD simulation are caused by orbiting clouds that temporarily enter the R = 150 pc region, but that eventually leave (M at the end of the HD simulation is zero as can be seen in the top panel). Inflow at these scales is mainly driven by the magnetorota-tional instability (Sect. 7).

In the text

	Fig. 24 Rate of change of the total angular momentum contained within cylindrical radius R (Eq. (29)) for the snapshot t = 196 Myr of the CHEM_MHD simulation. ALɀ/∆t is calculated by summing the angular momentum of all gas cells (Eq. (30)) in two different snapshots separated by ∆t ≃ 1 Myr, taking the difference and dividing by At. The other contributions are calculated using Eqs. (31)–(33).
In the text

Fig. 25 Turbulent Reynolds, Maxwell, and gravitational stresses at a given radius as a function of time in the CHEM_MHD simulation. Top: at R = 3.1 kpc, which tracks the entire bar region. The oscillations on timescales t ≃ 7 Myr represent a torque-weighted orbital period of clouds orbiting in the bar region (Fig. 27 and recalling that angular momentum is not conserved in a non-axisymmetric bar potential). The overall upward trend is because angular momentum is removed by the bar at a decreasing rate as the gas reservoir that supplies the bar-driven inflow is depleted and the bar-driven inflow stops (Fig. 22 and Sect. 7). Bottom: R = l49 pc, which tracks what happens inside the CMZ gas ring. Maxwell stresses become dominant at t > 200 Myr, and show that magnetic stresses are responsible for the nuclear inflow from the CMZ ring inwards.

In the text

Fig. 26 Turbulent Reynolds, Maxwell, and gravitational stresses at a given time as a function of radius for the CHEM_MHD simulation. Top: for R < 5 kpc. Bottom: for R < 200 pc. At 200 < R < 3 kpc the gravitational stresses dominate, and therefore the bar-driven inflow is driven by the gravitational torques of the Galactic bar (as the name would imply). At 0 < R < 200 pc the Maxwell stress dominate, and therefore the nuclear inflow is driven by magnetic fields. We note that the top and bottom panels are for different snapshots, which are selected to better capture the moments the bar-driven and nuclear inflows take place.

In the text

	Fig. 27 Mass-weighted gravitational torques map for two selected snapshots of our CHEM_MHD simulation. The bar removes (adds) angular momentum in the two red (green) quadrants.
In the text

Fig. B.1 Impact of resolution on magnetic fields and inflow. Top: volume-weighted average magnetic field intensity in the CMZ (R < 500 pc, |z| < 100 pc) in three different density regimes at different resolutions as a function of time. Bottom: mass accretion as a function of time from the CMZ inwards (R < 100 pc) for the three simulations at different resolutions. The saturation levels of the magnetic field are relatively insensitive to the resolution, while the inflow rate is highly sensitive to it.

In the text

Fig. C.1 Impact of the initial conditions on magnetic fields and inflow. Top: volume-weighted average magnetic field intensity in the CMZ (R < 500 pc, |ɀ| < 100 pc) in three different density regimes for poloidal initial magnetic field and toroidal initial magnetic field. Both tests are run at a resolution of 100 M⊙ per cell. Bottom: mass inflow rate as a function of time from the CMZ inwards (R < 100 pc) in the two simulations. The figure shows that the magnetic field intensity and inflow rates are relatively insensitive to the initial orientation of the magnetic field.

In the text

Fig. D.1 Effect of the cutout test on magnetic fields and inflow. Top: volume-weighted average magnetic field intensity in the CMZ in three different density regimes for three simulaitions: (i) the fiducial CHEM_MHD simulation; (ii) the test where we cut out the CMZ; (iii) and the test where we also reset the B field (Appendix D for description). The latter was run until 400 Myr to give time for dynamo processes and MRI to develop again. Bottom: mass accretion as a function of time from the CMZ inwards (R < 100 pc) for the three simulations. The vertical dashed line indicates the time at which we performed the cut-out.

In the text

	Fig. E.1 Same as Fig. 8 but here we integrate the Stokes parameters along the line-of-sight instead of directly integrating the magnetic field.
In the text

	Fig. F.1 Map of the relative divergence error of the magnetic field for a slice through the ɀ = 0 plane at t = 196 Myr. Top: the divergence error map is shown considering its sign as well. Bottom: map of the absolute value of the divergence error in log scale. The contour corresponds to an isodensity contour of n = 20 cm−3.
In the text

	Fig. F.2 Relative divergence error of the magnetic field, volume-averaged over the central 500 pc, as a function of time for the different simulations at different resolution. Top: the divergence error is shown considering its sign as well. Bottom: volume-average of the absolute value of the divergence error in log scale.
In the text