Modeling Turbulence in Galactic Centers
Jesus M. Salas Mark R. Morris Smadar Naoz





Jesus M. Salas
jesusms@astro.ucla.edu



Dept. of Physics & Astronomy, University of California, Los Angeles, CA, 90095, USA


Dept. of Physics & Astronomy, University of California, Los Angeles, CA, 90095, USA


Dept. of Physics & Astronomy, University of California, Los Angeles, CA, 90095, USA
Mani L. Bhaumik Institute for Theoretical Physics, Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA


Turbulence is a prevalent phenomenon in the interstellar medium, and in particular, the environment at the centers of galaxies. For example, detailed observations of the Milky Way's Central Molecular Zone (CMZ) revealed that it has a complex and turbulent structure. Turbulence on galactic scales is often modeled using star formation and feedback. However, these effects do not appear to be sufficient for explaining the high-velocity dispersion observed in the CMZ, indicating that additional gas-stirring processes are likely to be operating. Here we introduce a proof-of-concept method to drive turbulence in gas that orbits under the influence of a galactic potential. Instead of relying on a particular physical mechanism, we have adopted a Fourier forcing module and have applied it using a Smoothed Particle Hydrodynamics code. To test our method, we performed simulations of a simplistic model of the CMZ. 
Our turbulence injection method is capable of balancing the self-gravity of the gas, which allows us to run the simulations for long timescales and thereby follow the evolution of the CMZ. Our results show that turbulence induces a flocculent spiral pattern in our model, analogous to that found in galactic-scale simulations. Furthermore, we find that our turbulence injection method induces inward migration of gas, a result consistent with previous numerical simulations. We submit that this injection method is a promising new tool to simulate turbulence in galactic centers. 






§ INTRODUCTION


Turbulence is one of the major processes that governs
the structure and evolution of the interstellar medium <cit.>. For example, supersonic  turbulence  is  known  to be  a dominant  process in regulating (both inhibiting  and  fostering) star formation in molecular clouds <cit.>. However, galactic
centers show a number of interesting deviations from the star formation behavior seen at larger galactic radii. Molecular clouds in our own Galactic Center (GC) show systematic differences in their properties from
disk clouds: they appear to have thermal, turbulent and magnetic pressures much higher than those present in the large-scale Galactic disk <cit.>. 

The main gaseous feature of the GC, the Central Molecular Zone (CMZ), has a rich and complex structure that extends over a galactocentric radius of ∼300 pc and contains a mass of M∼3-77  <cit.>. 
It is largely composed of relatively dense (n∼10^3-10^5 ) and warm gas (∼70-100 K on average, e.g., ), mostly condensed into Giant Molecular Clouds (GMCs) or dense tidal streams of molecular gas. These relatively warm gas temperatures are one of the key properties of CMZ clouds, and there is evidence showing that the gas is kept warm by the dissipation of turbulence <cit.>. It has also been suggested that turbulence plays a role in the suppression of star formation in this region (e.g., ). However, the driving mechanism for the turbulence in CMZ clouds has not been conclusively identified (see  for a discussion of possible sources of turbulence). Furthermore, the large turbulent velocity dispersion within the CMZ must be responsible for supporting the gas against gravitational collapse, since the thermal pressure of the gas would be insufficient. 

Interstellar turbulence decays quite rapidly, on timescales of the order of the free-fall time of the system <cit.>. Therefore, energy must be injected into the system in order to maintain the turbulence. Simulations of turbulence-driven gas are often 
employed in studies of the interstellar medium and star formation <cit.>. Typically, this is achieved by a Fourier forcing module, which can be modeled with a spatially static pattern in which the amplitude is adjusted in time <cit.>. Other studies employ a forcing module that can vary both in time and space <cit.>.

In the case of galactic-scale simulations, driven turbulence is mimicked by injecting energy due to supernovae (SN). For example, simulations by <cit.>, and <cit.> have modeled turbulence by using star formation and SN feedback models. In general, these models depend on underlying assumptions regarding star formation rates, SN energies and injection rates. Furthermore, recent studies <cit.> have demonstrated that the different choices of SN feedback model (including the underlying physical processes driving the feedback) produce significant differences in morphology, density, etc, of the simulated galaxies. 

Here we introduce a proof-of-concept method to drive turbulence in gas that orbits under the influence of a galactic potential. Instead of relying on a particular physical mechanism, we adopt a Fourier forcing module, which has the advantage of being independent of the source of turbulence. Our turbulence treatment is based on the method described by <cit.>, in which a turbulent velocity field is drawn from a spatially static pattern having a power spectrum P(k)∝ k^-n, where k is the wavenumber. We apply our method to a smoothed particle hydrodynamics (SPH) code, and we test its effectiveness using a simplistic model of the CMZ. Our simulations consider self-gravity (i.e., the mutual gravitational interactions between the SPH particles) and the effects of pressure from a surrounding medium.

This paper is organized as follows: Section <ref> summarizes the numerical methods, with further details on our turbulence method in Appendix <ref>. Section <ref> and Appendix <ref> describe the tests performed to verify the effectiveness of our turbulence method. We present our main results in Section <ref>, and conclude in Section <ref>.




§ NUMERICAL METHODS


We used the N-body/SPH code  <cit.>, which is based on the tree-Particle Mesh method for computing gravitational forces and on the SPH method for solving the Euler equations of hydrodynamics. The smoothing length of each particle in the gas is fully adaptive down to a set minimum of 0.001 pc.  employs an entropy formulation of SPH, as outlined in <cit.>, with the smoothing lengths defined to ensure a fixed mass (i.e., fixed number of particles) within the smoothing kernel volume (set for N_neigh = 64). The code adopts the Monaghan-Balsara form of artificial viscosity <cit.>, which is regulated by the parameter α_MB, set to 0.75.

We modified the standard version of  to include turbulence driving, the gravitational potential of a Milky Way-like galaxy, and the effects of pressure by a surrounding medium. We describe these modifications below. 




 §.§ External Pressure


The interstellar medium of the GC is modeled via an external pressure term to approximate a constant pressure boundary.  
Following <cit.>, we modify 's momentum equation <cit.>:

    d v_i/dt = - ∑_j m_j [  f_i P_i/ρ^2_i∇_i W_ij(h_i) + f_j P_j/ρ^2_j∇_i W_ij(h_j)  ]
 
where v_i is the velocity of particle i, m_j is the mass of particle j, P_i is the pressure, ρ_i is the density, W_ij(h_i) is the kernel function which depends on the smoothing length h_i, and f is a unitless coefficient that depends on ρ_i and h_i. We replace P_i and P_j with P_i - P_ext and P_j - P_ext, respectively, where P_ext is the external pressure. The pair-wise nature of the force summation over the SPH neighbors ensures that P_ext cancels for particles that are surrounded by other particles. At the boundary, where the P_ext term does not disappear, it mimics the pressure contribution from a surrounding medium <cit.>. We set P_ext equal to 10^-10 ergs , an approximate value for the GC <cit.>. 



 §.§ The galactic potential
 

The gravitational potential we use is adopted from <cit.>, which is a modified version of the prolate bar potential introduced by <cit.>. This potential has the form: 

    Φ(r,θ,ϕ) = 4π G ρ_0 r_0^2 (r/r_0)^α P(θ,ϕ)  ,

where (r,θ,ϕ) are spherical coordinates fixed on the rotating bar (r = √(x^2+y^2+z^2), where x,y,z are the standard Cartesian coordinates. The supermassive black hole would be at r = 0, the bar's major axis is aligned with the x-axis, and the z axis represents the vertical direction, with the galactic plane at z=0), and P(θ,ϕ) is the associated Legendre function, which can be written as:

    P(θ,ϕ) = 1/α (1+α) - Y(θ,ϕ)/(2-α)(3+α) .

Y is a linear combination of spherical harmonic functions of the
l= 2, m= 0,2 modes:

    Y(θ,ϕ) = -b_20 P_20 (cosθ) + b_22P_22(cosθ)cos2ϕ .

The parameter b_20 determines the degree of oblateness/prolateness
while b_22 determines the degree of non-axisymmetry. Motivated by the previous work of <cit.>, and more recently of <cit.>[We note that there is a negative sign misprint in <cit.> (their Equation 2) and in <cit.> (their associated Legendre function).], we use the parameters: α = 0.25, b_20=0.3, b_22 =0.1, ρ_0 =40  pc^-3 and r_0 = 100 pc. Given these parameters, a bar with axis ratios of [1: 0.74: 0.65] is obtained for the isodensity surface that intersects points [x = 0, y = ±200 pc, z = 0].
Enclosed masses inside 200 pc and 1000 pc are 10^9  and 7×10^9 , respectively.


In addition to the gravitational force due to the potential
above, we introduced the rotation of the bar by adding centrifugal and Coriolis forces. We used the most recent estimate for the Galactic bar's pattern speed (Ω_bar = 40 km s^-1 kpc^-1; e.g., ).




 §.§ Turbulence Driving


As noted above, driven turbulence is often modeled using one of two types of Fourier forcing modules. Both methods require Fourier transforms on a cubic grid (or lattice) with N^3 points (or a square lattice with N^2 on 2D simulations, where typical values for N = 128, 256, 1024, etc). By imposing this cubic lattice onto a simulation box with a physical size of L per side, we can use the separation between adjacent lattice points (L/N) as a proxy for the resolution of the turbulence. Hence, for a large-scale simulation such as the CMZ environment, e.g., L = 500 pc, and simultaneously resolving turbulence on small scales, e.g., 0.01 pc, the turbulence cubic lattice would have to contain N^3 = 5000^3 points, which would require massive computational resources. To circumvent this limitation, we instead use many smaller turbulence grids to fill the volume of our larger simulation box, as follows:

First, we create a library of 10 files, which our modified version of the  code reads in at the start of the simulation. Each file contains a unique realization of a turbulent velocity field (in the form of a 3D matrix) with power spectrum P(k) ∝ k^-4 (suitable for compressible gas; e.g., ). Each of these 3D matrices of turbulence is generated using the methods described in <cit.> and <cit.>: via fast Fourier transforms inside a 128^3 box. 



These 3D turbulence matrices can be visualized as lattice cubes (or grids) with equally spaced lattice points, containing 128×128×128 points. The parameter that sets the physical size of these cubes is named L_cube. In the test model used in this paper, we make the simplifying assumption that turbulence is injected on scales similar to the average size of CMZ clouds, i.e., L_cube = 8 pc per side. Thus, the separation between two adjacent lattice points along one axis is 8pc/128 = 0.0625 pc. 

In some studies, the driving module only contains power on the larger scales (e.g., ). This type of driving models the kinetic energy input from large-scale turbulent fluctuations, which then break up into smaller structures as the kinetic energy cascades down to scales smaller than the turbulence injection scale. However, in SPH, the artificial viscosity may damp this energy cascade and prevent it from reaching the smaller scales. Consequently, to create the different realizations of turbulent velocity fields, we use a discrete range of k values from k_min = 2 to k_max = 128, thus effectively injecting energy on scales between L_cube/2 = 4 pc (for k=2) and L_cube/128 = 0.06 pc (for k=128). To create the initial (turbulent) velocity field of an individual cloud (see Section <ref>), we use trilinear interpolation to calculate the velocity components for each SPH particle, based on that particle's position on a turbulence cube. This interpolation method results in a turbulent velocity function I⃗(x,y,z).  

Next, we use 64×64×64 cubes of turbulence to fill up the volume of our large simulation box. This gives us a simulation box of size L_global = 8pc×64 = 512 pc per side. Thus, the spatial resolution of the turbulence in our large simulation box is the same as the resolution of an individual turbulence lattice cube. Each of the turbulence lattice cubes that fill up the large simulation box is randomly chosen from our library of 10 files, thereby avoiding a velocity field that is coherent over scales greater than L_cube. A 2D graphic representation of the method we describe here is shown in Figure <ref>. This method raises a concern regarding the interface between the turbulence cubes, namely, that there will be discontinuities in the turbulent velocity kicks throughout the gas at the cubes' interfaces. However, the turbulent velocity added to the gas is quite small in comparison with the orbital velocity, and indeed, we find in practice that there are no obvious shocks or discontinuities at the interfaces that exceed those induced by the injected velocity increments (see Section <ref>). 

To drive the turbulence, we follow a method similar to that described by <cit.>: every N_t timesteps (we fixed the timestep in all our simulations to Δ t_s = 1000 yrs) we add a velocity increment to every SPH particle, given by:

    Δv⃗(x,y,z) = A √(G ρ)I⃗(x,y,z) ,

where I⃗(x,y,z) is the turbulent velocity interpolated from the turbulence field of the cube that contains the particle in question. The amplitude A is chosen to maintain a constant kinetic energy input rate Ė_̇i̇ṅ = Δ E_in/(N_t Δ t_s), and the term √(G ρ) is added to counteract gravitational collapse by ensuring that higher-density regions receive more kinetic energy[This density factor stems from the assumption that molecular clouds are supported against gravitational collapse by turbulence, and therefore the velocity kicks should be inversely proportional to the free-fall time,  t_ff∝ (Gρ)^-1/2.]. Any particle outside the (512 pc)^3 simulation box does not receive any turbulent energy. 

For compressible gas with a time-dependent density distribution, maintaining a constant kinetic energy input rate requires solving a quadratic equation in the amplitude A every time the driving is performed <cit.>. For N_sph particles of mass m_p, each with density ρ_i, A is derived from (see Appendix <ref>):

    Δ E_in = 1/2 m_p ∑_i=1^N_sph[      A^2  G ρ_i I⃗_⃗i⃗·I⃗_⃗i⃗  + 2 A √(G ρ_i)I⃗_⃗i⃗·v⃗_⃗1⃗,⃗i⃗] .

Following <cit.>, we take the larger root as the value for A. 

Therefore, our methodology introduces two adjustable parameters: , the energy input per injection, and , the number of timesteps between injections.

Finally, to mimic the random nature of turbulence, we change the turbulent velocity field of each of the 64^3 cubes every time the driving is performed. This is done by replacing  each of the cubes with a different one, chosen randomly from the 10 files in our library. See Appendix <ref> for further details.




 §.§ Initial conditions


To initiate the computations, we use a fiducial model of the central region of the Milky Way, consisting of a collection of 100 initially isolated, spherical GMCs. Each cloud contains N_p=10^4 particles, with mass m_p = 30  per particle. The clouds are distributed randomly in an annular disk of inner radius 30 pc, outer radius 200 pc, and a Gaussian scale height of 30 pc. Each cloud has a radius of 4 pc, and an initial turbulent velocity field such that |E_turb/E_pot| = 0.5 (i.e., they are initially in virial equilibrium). The clouds' initial center of mass velocities were set so that they move on circular orbits and are parallel to the galactic plane, with their magnitudes (v_ϕ) calculated using the potential described in Section <ref>. In order to also give the system a broader vertical structure, the initial v_z components of each cloud were set such that v_z = 0.5v_ϕ, with the v_z vector always pointing towards the galactic plane. Figure <ref> shows the initial snapshot of our CMZ model.

All simulations were run using an isothermal equation of state with T = 100 K. This assumption of isothermal gas is somewhat crude, but may still provide an adequate physical approximation to the real thermodynamics in dense molecular gas <cit.>.  

During our testing phase, we also considered an initially uniform disk as our initial condition. Except for the time to reach a steady-state (see Section <ref>), the end result was qualitatively identical to the model described in this paper. We therefore conclude that the specific features of the initial conditions are unimportant for our purposes as long as the particles are initially distributed over the same domain. We opted for a collection of clouds as our initial conditions due to the flexibility in setting the position, size, mass, initial turbulent velocity field, etc, of each individual cloud.





§ CONVERGENCE AND CONSISTENCY TESTS
   
When the self-gravity of the gas is included, combined with a relatively low temperature, our CMZ simulations without injected turbulence experience runaway gravitational collapse, causing the simulations to fail within a dynamical timescale (∼0.1-0.3 Myrs).  This effect can be alleviated by increasing the softening length (which in our simulations is set equal to the smoothing length, i.e., 0.001 pc). In such a case, the simulation can run for long timescales, however the gas concentrates to unphysically large densities, causing the simulation to slow down to an impractical pace, and is thus an expensive use of computing resources. This result thereby emphasizes the importance of the turbulence injection method introduced in this paper. However, turbulence injected too infrequently leads to the same result as if there was no turbulence injection: gas collapses locally to unphysically large densities (see Appendix <ref>). Therefore, it is important to inject turbulence relatively often. 

In practice, the turbulence driving method described in Section <ref>, adds an additional velocity “kick” to each particle every  timesteps. The energy  injected per  timestep is kept constant. With these two free parameters we conducted several tests using our CMZ model to find the optimal range of values that give rise to reasonable densities over long timescales. We varied the energy input  from 10^46 to 10^50 ergs, and N_t from 1 to 5. We describe in detail all performed tests in Appendix <ref>. Our tests led us to the choice of parameters  = 2,  = 10^47 ergs (e.g., see Figure <ref>).  




§ RESULTS
 






We ran the system for 250 Myrs, using the initial conditions described in Section <ref>, the turbulence parameters described in Section <ref>, and including self-gravity. The system reaches a semi-steady state after 50-100 Myrs regardless of the turbulence parameters (as described in Appendix <ref>),  and thus the choice of 250 Myrs allows us to capture the relevant dynamics (see below). We also ran a simulation without turbulence (and self-gravity turned off) for comparison. We show the face-on and edge-on views of the resulting column density evolution of our simulations in Figures <ref> and <ref>, respectively. The results shown here indicate that our turbulence injection module is capable of balancing the self-gravity of gas concentrations, which allows us to run the simulation for long timescales. Furthermore, discontinuities due to the grids' interfaces are unnoticeable, as anticipated. 

In both runs, the clouds are tidally stretched relatively quickly, and the gas settles into a disk after ∼50 Myrs. In the simulation without turbulence, the clouds are stretched into gas streams which comprise the disk, and reaches steady state by ∼50 Myrs. The inner inner cavity of radius ∼30 pc, which is a product of the initial conditions, remains unchanged for the entire simulation. 

However, there are two major differences between the two simulations. First, the resulting streams in the simulation with turbulence coalesce into a disk with a flocculent spiral pattern. Second, the inner cavity fills in slowly over time. We address each of these differences below.



 §.§ Spiral structure

After running for ∼50 Myrs, the gas in the run with turbulence settles into a quasi-steady state with a flocculent spiral pattern (because of the constantly injected turbulence, a perfect steady state cannot occur). These spiral segments are attributable to the dynamical response of our self-gravitating, shearing disk to local density perturbations (e.g., ), which in this case are caused by the forced turbulence. It is a common result that turbulence promotes the development of high-density regions due to convergent flows (e.g., ). 

A similar result was found on galactic scales by <cit.>. The spiral patterns in their simulations of self-gravitating disks of stars are not global as predicted by classical static density wave theory, but locally they appear to fluctuate in amplitude with time. Their spirals are actually segments produced by sheared local under-dense and over-dense regions. These under-dense and over-dense regions act as gravitational perturbers, maintaining the local spiral morphology. In our case, it is the injected turbulence that acts as the local perturber in the gas.

Furthermore, observations of galactic centers have revealed intricate dust structures that are often organized in a clear spiral pattern. For example, the survey studies by <cit.> and <cit.> indicate that ∼50-80% of galaxies in their samples possess nuclear spirals, regardless of their nuclear activity.  The  pattern  of  some of the observed nuclear spirals is highly organized, similar to grand-design spirals in main galactic discs. Others display a more chaotic, or flocculent, spiral pattern. A direct comparison between our simulations and observations is difficult due to the limited spatial resolution of the observations, however, our results suggest that turbulence could be (at least partially) responsible for the spiral pattern observed in the gas layer in the centers of galaxies.  



  §.§.§ Effects of external pressure

As discussed in Section <ref>, we modified Gadget2 to model an external pressure boundary, as opposed to the vacuum or periodic boundary conditions that are the only choices in the standard version of the code.  We used an observationally motivated value of P_ext = 10^-10 erg cm^-3 for the Galactic Center. However, it is worth investigating how varying this pressure term affects the gas morphology. We ran two additional simulations with P_ext = 0 (standard vacuum boundary conditions) and P_ext = 10^-9 erg cm^-3, which is an order of magnitude greater than the value from our fiducial simulations. We show the column density maps (both face-on and edge-on) at t = 100 Myr in Figure <ref>. As depicted, without the external pressure term (see left panels in Figure <ref>), the gas disk has a smooth boundary at the edges, as well as material above and below the plane. When the external pressure term is introduced with our fiducial value, it pushes the low density gas at the edge towards the disk, and it increases the density contrast of the flocculent spiral pattern. Similarly, when the pressure term is large (P_ext = 10^-9 erg cm^-3, right panels in Figure <ref>), the density contrast of the spirals is increased further. Furthermore, the material above and below the plane is pushed towards the disk. Interestingly, this behaviour suggests an observational test for the value of external pressure, which may be used for systems in which the pressure is undetermined. 





 §.§ Inward migration




The inner cavity closes slowly over time, as gas migrates inward. This might imply a loss of angular momentum in the system due to the turbulence. In Figure <ref> we show the total (specific) angular momentum over time for both simulations. Despite our strategy for avoiding adding a net angular momentum to the simulations (see Appendix <ref>), we see that, in the run with turbulence, angular momentum slightly decreases, from  ∼1.45×10^10 km s^-1 pc at 50 Myrs to ∼1.4×10^10 km s^-1 pc at 250 Myrs, a difference of ∼3%. Clearly, the sign switching strategy we describe in Appendix <ref> does not completely ensure angular momentum conservation. Furthermore, it has been shown that turbulence can impose a nonzero angular momentum <cit.>. However, this change is quite small over long timescales, and thus we do not consider this to be the reason for the inward migration[ We verified that the inward migration is indeed a feature of the injected turbulence by considering the possibility that the interaction between turbulence and the rotating, non-axisymmetric potential could have led to this angular momentum decrease. We tested this idea by running a simulation with a rotating, axisymmetric potential (by setting b_22 = 0 and Ω_bar = 40 km s^-1 kpc^-1, see Section <ref>), and another with a non-rotating, axisymmetric potential (by setting b_22 = 0 and Ω_bar = 0). However, in both cases, we recover the same negative slope as in our fiducial model (red line in Figure <ref>). Thus, we conclude that this small angular momentum reduction is a feature of the injected turbulence. ].

Conversely, it has been shown that supersonic turbulence inside accretion disks (e.g., ) can promote accretion onto SMBHs by enhancing angular momentum transfer (e.g., ). In particular, <cit.> demonstrated using numerical simulations of supermassive black hole accretion that turbulence can broaden the angular momentum distribution, setting some gas on low angular momentum orbits. We find a similar result in Figure <ref>, where we show the mass fraction as a function of specific angular momentum at different times for both simulations. We find that in the simulation without turbulence (left panel of Figure <ref>) the gas settles into a triple-peaked distribution by 100 Myrs (red curve), and the location of these peaks, as well as the overall distribution, does not change over time. In contrast, the simulation with turbulence exhibits an angular momentum distribution with a single peak (right panel of Figure <ref>). This peak then slowly moves towards lower angular momentum values over time. 

This transport of angular momentum can be explained in terms of the elementary theory of accretion disks <cit.>: viscous torques between adjacent annuli of gas in an accretion disk provokes mass to flow inwards. While Gadget2 does contain an artificial viscosity, this treatment is dedicated to capturing shocks in SPH, and has no effect on the transfer of angular momentum, as seen in our turbulence-free simulation (Figure <ref>). The transport of angular momentum is due to “turbulent” viscosity induced by our driving method. This explains the slow inward migration of gas in our simulation. 

To estimate the viscosity, ν_ turb, induced by our turbulence method, we use the accretion rate due to α-viscosity <cit.>:

    Ṁ≈ 3πν_ turbΣ    ,

where Ṁ is the mass inflow rate and Σ is the surface density. 





To estimate Ṁ, we calculate the mass inside a radius of 30 pc over time (see Figure <ref>). We approximate the rate from 150 Myrs to 250 Myrs as a straight line, and find a slope of Ṁ = 0.004 /yr. Similarly, we adopt the value of Σ at r = 30 pc to be 300 /pc^2, which is the average value between 150 and 250 Myrs (see Figure <ref>). Thus, we evaluate the viscosity to be ν_ turb=4.2×10^23 cm^2/s. 

Alternatively, we can use the definition of α-viscosity ν_ turb = α c_s H, where α≤ 1 is a parameter that adjusts the strength of the viscosity, c_s is the sound speed, H = c_s/Ω, and Ω is the angular velocity due to the gravitational potential. Assuming the gas in our simulations is mainly composed of molecular hydrogen (H_2), c_s = 0.64 km/s. The value of H at r = 30 pc is ∼0.2 pc (see Figure <ref>). This results in a value for ν=3.7×10^22 cm^2/s (assuming α=1). However, the actual thickness of the disk at r = 30 pc is much larger, ∼2 pc (see Figure <ref>). Using this value, we obtain ν_turb=3.7×10^23 cm^2/s, which is comparable to the value we calculated above using the mass inflow rate.

To understand the implications of this calculated turbulent viscosity, we make a comparison similar to that described in <cit.>, as follows: to significantly affect the dynamics of a gaseous disk, the viscous timescale, t_ν∼ R^2/ν (where R is the radius of the disk), must be shorter than the Hubble time (t_H = 14 Gyrs). This condition gives a lower limit to the value of the viscosity, ν_min≈ R^2/t_H. Using R=200 pc (the radius of our simulated CMZ disk) gives a minimum viscosity of ν_min = 8×10^23 cm^2/s. This value  of ν_min is roughly a factor of 2 higher than the value ν_turb≈ 4×10^23 cm^2/s calculated from our simulations. This result is consistent with the findings by <cit.>, where they found that the viscosity used in their simulations of galactic nuclear rings, ν=3×10^24 cm^2/s (which was sufficient to significantly affect the morphology of their simulated rings) was a factor of ∼10 smaller than the corresponding minimum viscosity for nuclear rings, ν_min = 8×10^25 cm^2/s (for a R=1 kpc ring). Thus, our results support <cit.> conclusion that viscosity may be more effective in influencing the dynamics of gaseous systems than implied by generic estimates such as the one described above.

We note that our viscosity calculations are only approximations, since the values of Ṁ, Σ, and the thickness of the disk, vary over radius and time. Thus, the viscosity induced by our turbulence method will also vary with radius and time. A more careful analysis of the turbulent viscosity imposed by our driving method is reserved for a future paper.




§ SUMMARY AND DISCUSSION


Supersonic turbulence occurs over a wide range of length scales in the interstellar medium, especially within molecular clouds. The importance of turbulence in modulating star formation in the interstellar medium was further highlighted recently by a combination of numerical and analytical studies  <cit.>. It has also been suggested that turbulence plays a key role in forming the very first star-clusters and perhaps even globular cluster progenitors <cit.>. Furthermore, turbulence in the centers of galaxies, particularly in our own CMZ, seems to greatly influence its thermal structure and star formation rate <cit.>. 








Numerical simulations have shown that turbulence decays quickly, within a few dynamical timescales <cit.>. Thus, in order to properly study the dynamics of the CMZ environment with numerical simulations, it is necessary to include a mechanism for turbulence driving. Typically, simulating turbulence in gas on galactic scales is achieved by modeling star formation and feedback. However, we note that despite the high densities and the large amount of available gas, there is about an order of magnitude less active star  formation in the CMZ than might be expected from the quantity and surface density of molecular gas <cit.>. Therefore, other physical processes may supplement this methodology to explain the high-velocity dispersions observed in the CMZ. 

This prompted us to develop a more general mechanism for driving turbulence. Consequently, we adapted the Fourier forcing module described by <cit.> in order to simulate not just the CMZ, but this method can be applied in general to gas that orbits under the influence of a gravitational potential. We implemented this turbulence method to the SPH code Gadget2.

Rather than depending on a single grid of turbulence, as is typical with this method, we create instead many turbulent grids and use them to fill the volume of a larger simulation box. We then use trilinear interpolation to calculate and add a velocity kick to each SPH particle. The amplitude of these velocity increments is adjusted every injection time to maintain a constant energy input. As shown in Section <ref>, discontinuities in the added velocity field to the interfaces between grids are indistinguishable in the simulations. 

Using a simplistic model of the CMZ consisting of a collection of Giant Molecular Clouds, our turbulence driving method allows us to study the dynamics of the gas and turbulence over long timescales. One of the main results in our simulations is that turbulence induces inward migration of gas, a result that is consistent with previous numerical simulations <cit.>. 

This paper focusses on testing this new method in the context of a single physical scenario. However, our turbulence module is flexible and can be applied to study different physical scales. There are a number of parameters that can be adjusted:

    
  * Power spectrum index: For incompressible turbulence, the Kolmogorov power spectrum in three dimensions is ∝ k^-11/3 <cit.> (for a two-dimensional distribution, P ∝ k^-8/3). For compressible turbulence, the power law index has been shown to be slightly steeper <cit.>. Milky Way observations have found power law slopes of -2.8 to -3.2 in 2D maps <cit.>. In our tests we used a power spectrum index of -4, but the power spectrum index can be modified to match observations in order to create more applicable simulations. 
    
  * The number of distinct realizations of turbulence grids: in order to avoid coherence in scales larger than L_cube, each of the turbulence grids has a different realization of velocity fields, drawn randomly from a library of 10 realizations. This library can be expanded with more realizations, which could potentially improve the “randomness” of the turbulence. However, this library will be limited by the memory constraints of the computing resources.
    
  * The size of the turbulence grids (L_cube): in our code, this parameter also represents the largest scales on which turbulence is injected. Thus, modifying this allows flexibility in studying different turbulence injection scales depending on the physical environments to be simulated. For example, in <cit.>, we use our driving method to study how turbulence affects the formation of galactic nuclear rings. There, we discuss galactic scale (L_global = 4 kpc) simulations using larger grids of turbulence (L_cube = 64 pc).
    
  * The number of turbulence grids per side: modifying this parameter (in tandem with the previous one, L_cube), allows for adjusting the size of the overall simulation domain, L_global.
    
  * The total energy of injection, , and time interval between injections, : these parameters need to be tuned depending on the physical environment to be simulated. In Appendix <ref> we vary these parameters in order to find the optimal values to use. However, different gas configurations will require different parameters than the ones used here. For example, we expect that smaller values of  and larger values of  may be sufficient to balance self-gravity in lower density gas.




J.M.S would like to thank Paul Clark for his help on adapting the external pressure term into . Furthermore, J.M.S would like to thank Blakesley Burkhart and Diederik Kruijssen for very helpful discussions about turbulence, and finally, to Sofia G. Gallego and Sungsoo S. Kim for their help in implementing the gravitational potential of the inner Galaxy. This material is based upon work supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE-1144087. SN acknowledges the partial support of NASA grant No. 80NSSC20K0500 and thanks Howard and Astrid Preston for their generous support. This work used computational and storage services associated with the Hoffman2 Shared Cluster provided by UCLA Institute for Digital Research and Education's Research Technology Group. This work also used the Extreme Science and Engineering Discovery Environment (XSEDE) Comet at the San Diego Supercomputer Center at UC San Diego through allocations AST170039 and AST180051. XSEDE is supported by National Science Foundation grant number ACI-1548562.

Software:
Figure <ref>, <ref> and <ref> were done using the SPH visualization software  <cit.>. We used Gadget2 <cit.> to build our turbulence method. The version of the code that includes our turbulence routine can be found at https://github.com/jesusms007/turbulencehttps://github.com/jesusms007/turbulence.








§ CONSTANT ENERGY INPUT RATE



In this appendix we derive Equation <ref>. 

Our algorithm adds a velocity “kick” to each particle every N_t timesteps while  maintaining a constant energy input rate Ė_̇i̇ṅ = Δ E_in / (N_t Δ t), where Δ t is the simulation timestep (fixed to be equal to 1000 yrs), and:

    Δ E_in = E_2 - E_1 = 1/2 m_p ∑_i=1^N_sphv⃗_⃗2⃗,⃗i⃗·v⃗_⃗2⃗,⃗i⃗ - 1/2 m_p ∑_i=1^N_sphv⃗_⃗1⃗,⃗i⃗·v⃗_⃗1⃗,⃗i⃗ ,

where v⃗_⃗1⃗,⃗i⃗ is the velocity vector of a particle at time t_1 (before the kick), and

    v⃗_⃗2⃗,⃗i⃗ = v⃗_⃗1⃗,⃗i⃗ + AΔv⃗_⃗i⃗ = v⃗_⃗1⃗,⃗i⃗ + A √(G ρ_i)I⃗_⃗i⃗(x,y,z)

is the velocity of the particle at time t_2 (after the kick). Therefore, t_2 - t_1 = N_tΔ t_s. The function I⃗(x,y,z) is the interpolation function that represents the turbulent velocity increment based on the particle's position (hereafter called I⃗), and A is the target variable.

Equation <ref> then becomes:

    Δ E_in  = 1/2 m_p [ ∑_i=1^N_sph(v⃗_⃗1⃗,⃗i⃗ + A √(G ρ_i)I⃗_⃗i⃗) ·(v⃗_⃗1⃗,⃗i⃗ + A √(G ρ_i)I⃗_⃗i⃗) -  ∑_i=1^Nv⃗_⃗1⃗,⃗i⃗·v⃗_⃗1⃗,⃗i⃗]


By simplifying Equation <ref>, the result is:

    Δ E_in = 1/2 m_p ∑_i=1^N_sph[      A^2  G ρ_i I⃗_⃗i⃗·I⃗_⃗i⃗  + 2 A √(G ρ_i)I⃗_⃗i⃗·v⃗_⃗1⃗,⃗i⃗]


This equation always gives both positive and negative values for A and our code always chooses the positive value. However, to ensure that the injection of turbulence does not violate the conservation of total angular momentum of the gas, we simply multiply A by a factor of -1 or +1, alternating between these two factors every time the driving is performed. This ensures that, over time, the net angular momentum added to the gas particles is approximately zero while keeping the same increase in energy Δ E_in. 

Furthermore, we change each of the 64^3 cubes of turbulence every time the driving is performed. This is to mimic the random nature of turbulence. Regardless of the source of turbulence, we expect that a parcel of gas will experience a coherent turbulent driving force during a given time interval. We can justify changing each grid every N_t timesteps if we consider that the crossing time of a parcel of gas traveling at an orbital speed of v≈150 km s^-1 through the average scale of a turbulence cube, i.e., L_ave≈1 pc, is t_cross≈6000 yrs. Hence, we conclude that there is little need for the velocity field to be coherent on timescales longer than 6 timesteps, thus justifying our replacement of each turbulence cube as long as N_t is more than a few.



§ DEPENDENCE ON TURBULENCE PARAMETERS


Here we describe the tests we carried out to study the performance of our turbulence algorithm in order to choose  and  values for our CMZ model. While these parameters are better represented in form of a rate, the results of this section indicate that the time interval at which the turbulence is injected affects the subsequent evolution of the gas. 

We tested driving the turbulence every 1 to 5 timesteps for a total of 25 tests: for each , we used  = 10^46 to 10^50 ergs, in a factor of 10 increments. We ran each simulation for 100 Myrs to allow the system to reach a quasi-steady state. 

As proxy for the evolution of the systems, we plot the average number density (n_ave) as a function of time, as shown in Figures <ref> to <ref>. For the simulations with = 2 and 3, the time evolution of n_ave is very similar (except for the case with = 10^50 ergs, whose n_ave evolution diverges from all other runs). The average density oscillates due to the interplay between turbulence and self-gravity, potentially reaching a steady-state by ∼50 Myrs.

We also plot in Figures <ref> to <ref> the RMS number density (n_RMS) from t = 50-100 Myrs, in order to better discern differences between runs. The runs with parameter = 4 and 5 show higher average densities, as expected, because in these cases, turbulence is injected less frequently, allowing gravity more time to compress the gas to higher densities. However, given their high density RMS peaks, these runs are less consistent and more chaotic than those with lower values of . Also, the runs with  = 5 were only run for 50 Myrs due to the fact that the tests with  = 10^46 to 10^48 ergs exhibited high density clumps which slowed down the computation time. These clumps originate because in those cases, the turbulence is not injected often enough to support the gas against local gravitational collapse. Particles pile on top of each other, and due to the nature of the kernel used by , once the distance between particles approaches the smoothing length, the pressure gradient is no longer correct and the particles stick together, creating very high density clumps.



Additionally, in Figures <ref> to <ref>, we plot the distribution of mass fraction as a function of density (density PDF) at t = 50 and 100 Myrs for each run. Inside a molecular cloud, the density PDF is shaped by the complex interaction between turbulence, self-gravity, magnetic fields and stellar feedback. As a result, it is an effective tool to determine the dynamical state of the gas (e.g. ). In most of our test runs, the calculated density PDF resembles a lognormal distribution, a result that several groups have found for isothermal gas (e.g.,  and others). This gives us confidence that our new turbulence forcing module is consistent with previous studies of turbulence. Thus, we can fit the density PDFs to Gaussian functions of x = ln(n/n_0) with mean μ and dispersion σ:

    f_M = C  exp[ (x - |μ|)^2/σ^2] .

The mass-weighted median number density (half of the mass is at densities above and below this value) is proportional to e^μ. We indicate this parameter with vertical bars in both figures. For each value of  = 1, 2 and 3, the value mass-weighted median number density is consistent for all energies, except  = 10^50.

Finally, it is a standard result in the literature that high-density regions are created by turbulence due to supersonic turbulent convergent flows (e.g., ). The increasing values of n_ave with  seen in Figures <ref> to  <ref> agrees with this result.

The tests with parameter  = 2 and 3 show better consistency in the range of densities in the simulations (n∼ 10^3-4 ) than the other test runs. Based on their density evolution, we notice that the exact choice of of energies in this range of injection intervals has little effect on the final results. Thus we opt for the parameters  = 2 and  = 10^47 ergs as the standard choice for our CMZ model simulated in Section <ref>. 

As a final note, we remind the reader that we are only testing this parameter-space of injection energies and intervals in a specific physical scenario, namely, that of a CMZ-like environment. We expect that similar choice of parameters will be useful in different environments, but further testing is needed. However, our goal here is to show that our proof-of-concept method gives consistent results and agrees with previous results in the literature. A more comprehensive analysis of the turbulence method that we introduce in this paper is left for future work. 











aasjournal


