<figure><img src="image_1.png" /></figure>

# When Gas Dynamics Decouples from Galactic Rotation: Characterizing ISM Circulation in Disk Galaxies

José Utreras<fnref n="1" />, Guillermo A. Blanc<fnref n="1" />,<fnref n="2" />, Andrés Escala<fnref n="1" />, Sharon Meidt<fnref n="3" />, Eric Emsellem<fnref n="4" />,<fnref n="5" />, Frank Bigiel<fnref n="6" />, Simon C. O. Glover<fnref n="7" />, Jonathan Henshaw<fnref n="8" />, Alex Hygate<fnref n="8" />,<fnref n="9" />, J. M. Diederik Kruijssen<fnref n="9" />, Erik Rosolowsky<fnref n="10" />, Eva Schinnerer<fnref n="8" />, and Andreas Schruba<fnref n="11" />

Received 2019 December 9; revised 2020 February 25; accepted 2020 February 25; published 2020 April 1

## Abstract

In galactic disks, galactic rotation sets the bulk motion of gas, and its energy and momentum can be transferred toward small scales. Additionally, in the interstellar medium, random and noncircular motions arise from stellar feedback, cloud–cloud interactions, and instabilities, among other processes. Our aim is to comprehend to what extent small-scale gas dynamics is decoupled from galactic rotation. We study the relative contributions of galactic rotation and local noncircular motions to the circulation of gas, $\Gamma$, a macroscopic measure of local rotation, defined as the line integral of the velocity field around a closed path. We measure the circulation distribution as a function of spatial scale in a set of simulated disk galaxies and model the velocity field as the sum of galactic rotation and a Gaussian random field. The random field is parameterized by a broken power law in Fourier space, with a break at the scale $\lambda_c$. We define the spatial scale $\lambda_{eq}$ at which galactic rotation and noncircular motions contribute equally to $\Gamma$. For our simulated galaxies, the gas dynamics at the scale of molecular clouds is usually dominated by noncircular motions, but in the center of galactic disks galactic rotation is still relevant. Our model shows that the transfer of rotation from large scales breaks at the scale $\lambda_c$, and this transition is necessary to reproduce the circulation distribution. We find that $\lambda_{eq}$, and therefore the structure of the gas velocity field, is set by the local conditions of gravitational stability and stellar feedback.

Unified Astronomy Thesaurus concepts: Hydrodynamical simulations (767); Hydrodynamics (1963); Astronomical simulations (1857); Astrophysical fluid dynamics (101); Galaxy dynamics (591); Interstellar dynamics (839); Galaxy rotation (618); Astrostatistics strategies (1885); Posterior distribution (1926)

## 1. Introduction

The structure of the gas velocity field is crucial for understanding how galaxies and molecular clouds evolve. The dynamical state of gas is one of the key elements in star formation theories, e.g., invoking turbulence at the scale of molecular clouds (Padoan et al. 2012; Semenov et al. 2016) or galactic rotation as a particular parameter controlling star formation at galactic scales (Elmegreen 1997; Silk 1997; Kennicutt 1998; Tan 2000; Krumholz et al. 2012; Utreras et al. 2016; Jeffreson & Kruijssen 2018; Meidt et al. 2018). The common picture of star formation involves self-gravity and sources of energy acting against self-gravity. Galactic rotation is one of those energy sources, acting at the largest spatial scales, where ordered motions make up the bulk of the kinetic energy (Utreras et al. 2016; Colling et al. 2018; Meidt et al. 2020). However, while its importance is evident on large scales, it is not clear down to which spatial scales galactic rotation remains dynamically relevant. At the scales of molecular clouds or stellar cores, gas can be dynamically less coupled with galactic rotation, and local noncircular motions start to dominate.

These nonordered or noncircular motions are originated by gravitational instabilities, hydrodynamical instabilities (Matsumoto & Seki 2010; Renaud et al. 2013; Sormani et al. 2018), cloud–cloud interactions, torques from nonaxisymmetric potentials, gas accretion, and stellar feedback (Goldbaum et al. 2015; Krumholz et al. 2018). These energy sources inject turbulence and induce noncircular motions that cascade toward small and large scales (Kraichnan 1967; Bournaud et al. 2010). The scale at which gas motion goes from being dominated by ordered rotation to these noncircular dynamical regimes depends on the importance of these other processes relative to ordered rotation. While large-scale dynamics are set by the galactic angular velocity $\Omega(R) = V(R)/R$, small-scale noncircular motions are more difficult to model.

At galactic scales many studies have focused on the role of galactic rotation in the stability of gaseous rotating disks described by the Toomre parameter $Q$ (Toomre 1963). If $Q < 1$, the disk is gravitationally unstable to radial perturbations. The classical form of this parameter involves the stability of a razor-thin disk of gas with $Q = \frac{\kappa c_s}{\pi G \Sigma_{\text{gas}}}$, where $c_s$ is the gas sound speed, $\kappa$ is the epicyclic frequency given by $\kappa^2 = 4\Omega^2\left(1 + \frac{1}{2}\frac{\partial \ln \Omega}{\partial \ln R}\right)$, $\Sigma_{\text{gas}}$ is the gas surface density, and $G$

<footnote n="1">¹Departamento de Astronomía, Universidad de Chile, Casilla 36-D, Santiago, Chile; jutreras@ug.uchile.cl</footnote>

<footnote n="2">²Observatories of the Carnegie Institution for Science, 813 Santa Barbara Street, Pasadena, CA 91101, USA</footnote>

<footnote n="3">³Sterrenkundig Observatorium, Universiteit Gent, Krijgslaan 281 S9, B-9000 Gent, Belgium</footnote>

<footnote n="4">⁴European Southern Observatory, Karl-Schwarzschild-Str. 2, D-85748 Garching, Germany</footnote>

<footnote n="5">⁵Université Lyon 1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230 Saint-Genis-Laval, France</footnote>

<footnote n="6">⁶Argelander-Institut für Astronomie, Universität Bonn, Auf dem Hügel 71, D-53121 Bonn, Germany</footnote>

<footnote n="7">⁷Zentrum für Astronomie, Universität Heidelberg, Institut für Theoretische Astrophysik, Albert-Ueberle-Str. 2, D-69120 Heidelberg, Germany</footnote>

<footnote n="8">⁸Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="9">⁹Astronomisches Rechen-Institut, Zentrum für Astronomie der Universität Heidelberg, Mönchhofstraße 12-14, D-69120 Heidelberg, Germany</footnote>

<footnote n="10">¹⁰Department of Physics, 4-181 CCIS, University of Alberta, Edmonton, AB T6G 2E1, Canada</footnote>

<footnote n="11">¹¹Max-Planck Institut für Extraterrestrische Physik, Giessenbachstraße 1, D-85748 Garching, Germany</footnote>
is the gravitational constant. This ideal case illustrates how galactic rotation delivers support against collapse, in particular to perturbations of size $\lambda > \lambda_{\text{rot}} \equiv 4\pi^2 G \Sigma_{\text{gas}} / \kappa^2$, setting a maximum size of collapsing fragments (Escala & Larson 2008). Following this body of work, it is natural to expect that rotation plays a significant role in the dynamics of sufficiently large molecular clouds and ultimately in the process of star formation on galactic scales. In particular, the works of Padoan et al. (2012) and Utreras et al. (2016) show a difference in the efficiency of star formation at these two different scales. Padoan et al. (2012) found that in a turbulent cloud the efficiency is proportional to $\exp(-t_{\text{ff}}/t_{\text{cr}})$, while at galactic scales Utreras et al. (2016) found an efficiency proportional to $\exp(-t_{\text{ff}}/t_{\text{orb}})$, where $t_{\text{ff}}$ is the initial freefall time, $t_{\text{cr}}$ is the cloud crossing time, and $t_{\text{orb}} = 2\pi/\Omega$ is the orbital time. These works invoke different processes as being important to control the star formation process for two different spatial regimes. It is expected that dynamics are linked from large to small scales; however, most observational studies and theories have neglected the multiscale nature of gas dynamics.

As we move our analysis toward the scale of molecular clouds, noncircular motions start to become relevant. A significant body of observational and theoretical research has been devoted to studying the balance between gravitational potential energy $W$ and kinetic energy $K$, commonly described by the virial parameter $\alpha_{\text{vir}} = 2K/W$ (Leroy et al. 2017; Padoan et al. 2017; Sun et al. 2018). CO measurements in the PHANGS-ALMA<fnref n="12" /> survey made by Sun et al. (2018) show that $\alpha_{\text{vir}}$ varies weakly from cloud to cloud, with $\alpha_{\text{vir}} \sim 1.5-3.0$, expected values for marginally bound clouds or free-falling gas. Simulations of turbulent molecular clouds from Padoan et al. (2012) have shown that the efficiency of star formation per freefall time, $\epsilon_{\text{ff}}$, is sensitive to the strength of self-gravity and decreases exponentially with $\alpha_{\text{vir}}$ (see also Federrath & Klessen 2012; Semenov et al. 2016). This anticorrelation between $\alpha_{\text{vir}}$ and $\epsilon_{\text{ff}}$ has been observed in M51 (Leroy et al. 2017) and in low-pressure atomic-dominated regions in nearby galaxies (Schruba et al. 2019).

A common assumption is that at the scale of molecular clouds, most of the kinetic energy $K$ comes from noncircular turbulent motions. This assumption is supported by measurements of velocity gradients of molecular clouds in nearby galaxies. Studying our Galaxy, Koda et al. (2006) estimated that the fractions of clouds with prograde or retrograde rotation with respect to the Galaxy’s spin are similar. Another studied galaxy in this subject is M33: by measuring velocity gradients, Rosolowsky et al. (2003) found that, if clouds do rotate, nearly 40% of molecular clouds are counterrotating with respect to the galaxy. More recently, Braine et al. (2018) found that in M33 molecular clouds do rotate and that their rotation is low, contributing little to the support of the cloud against gravity. These results are expected for clouds dominated by noncircular motions that have randomly aligned spins. In the field of simulations, Tasker & Tan (2009) found similar fractions of prograde and retrograde clouds in a simulated Milky Way–like galaxy, even in the absence of stellar feedback. Tasker & Tan (2009) argued that as time progresses, cloud–cloud interactions inject turbulence at the scale of these interactions.

However, the relevance of galactic rotation versus non-circular motions might depend on the local environment or

position in the galactic disk. To compare galactic rotation and noncircular motions, our interest focuses on in-plane motions. Meidt et al. (2018) argue that radial variations in the galactic potential are able to influence the dynamics of molecular clouds. Particularly, Meidt et al. (2018) argue that the internal dynamics of a cloud depend on the velocity field and the cloud size $R_c$ relative to the epicyclic frequency $\kappa$. For example, if $R_c > \sigma_v/\kappa$, where $\sigma_v$ is the velocity dispersion of gas, the Coriolis force is still relevant in the dynamics of molecular clouds. In other words, for gas structures larger than $\sigma_v/\kappa$ galactic rotation is still relevant. Since $\kappa$ and $\sigma_v$ vary with galactocentric radius, we might expect that the dynamics of molecular clouds change across a galaxy. Simulations of the galactic center (Kruijssen et al. 2019) show that molecular clouds are dominated by strong shear and tidal deformations. Moreover, shear motions from galactic rotation might set the cloud lifetimes in certain conditions (Jeffreson & Kruijssen 2018).

One way to estimate the role of galactic rotation is to measure its impact in the local rotation of gas, i.e., the rotation measured with respect to an inertial reference frame. The local rotation of gas is influenced by the large-scale motion of the galaxy and by noncircular motions acting on multiple scales. Our aim is to create a framework that allows us to obtain the contributions from these two types of motion to the local rotation. Since nonordered motions have multiple sources, we need to adopt a statistical approach. A first-order approximation is to consider noncircular motions as a Gaussian random field (GRF), described by a generating function in Fourier space $V(k)$, where $k$ is the wavenumber. If we know $V(k)$, we also know the magnitude of noncircular motions as a function of spatial scale. Ultimately, $V(k)$ let us know at which scales galactic rotation is still relevant.

We will employ a quantitative measure of the local rotation of gas, the circulation of a fluid $\Gamma$, which is defined as a line integral of the velocity field along a closed path and corresponds to a macroscopic measure of rotation. We define a two-component model for gas motions with a smooth function for large scales and a generating function $V(k)$ to model the noncircular motions. The velocity field arising from $V(k)$ behaves as a GRF. We compare the contributions from each component to the total measured circulation. In this framework, on galactic scales the contribution of noncircular motions to the circulation is negligible compared to the large-scale ordered rotation. On the smallest scales the circulation field is given mostly by $V(k)$. In other words, changes in the behavior of the observed distribution of $\Gamma$ at different scales illustrate how the dynamics transitions from circular to noncircular motions. With this in mind, we can define a spatial scale $\lambda_{\text{eq}}$ at which large-scale rotation and noncircular motions contribute equally to the measured circulation of gas.

To test whether circulation is a useful tool to find the transition scale between galactic rotation and noncircular motions, we use hydrodynamical simulations of galactic disks with different initial conditions. Numerical simulations are an excellent test bed for the study of circulation since they provide the full velocity field and allow us to look for observable signature by changing different physical parameters, such as rotation or self-gravity. In future work we plan to expand these methods to make them applicable to high-resolution observations of gas velocity fields, like those being produced by the PHANGS project.

<footnote n="12">¹² Physics at High Angular-resolution in Nearby GalaxieS with ALMA: www.phangs.org.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Illustration of the local rotation of fluid elements for different azimuthal velocity fields. Blue squares represent the fluid elements in two different positions in their path. The positions on each panel are the same and hence do not correspond to a specific instant in time. Left: solid-body rotation, $v \propto R$. The local rotation of any fluid element (blue squares) is equal to the galactic angular velocity $\Omega$. Middle: flat velocity curve, $v = v_0$. The local rotation of fluid elements $\Omega/2$. Right: irrotational flow, $v \propto 1/R$. For this flow the local rotation of a fluid element is zero.</figcaption></figure>

The paper is organized as follows. In Section 2 we introduce the main quantities analyzed in this work, vorticity and circulation, and we describe our technique to study the circulation in galaxies. In Section 3 we describe the simulations and the application of our technique in those objects. We discuss our results in Section 4 and list possible caveats of our simulations in Section 5.

## 2. Method

### 2.1. Vorticity and Circulation in Gas Dynamics

One of the most useful notions in fluid dynamics is the vorticity vector, $\omega$. In simple words, vorticity is a measure of the local rotation, and its direction is parallel to the spin of a fluid element. An infinitesimal fluid element experiences a rotation of $2\omega$ with respect to a local inertial reference frame. The vorticity is given by the curl of the velocity field

$$ \omega = \nabla \times \mathbf{v}. \quad (1) $$

To make this concept clearer, let us imagine a fluid with a circular velocity field $\mathbf{v} = v(R)\hat{\phi}$, where $\hat{\phi}$ is the azimuthal unit vector in cylindrical coordinates. In this scenario, the vorticity is

$$ \omega = \frac{v(R)}{R} \left( 1 + \frac{\partial \ln v(R)}{\partial \ln R} \right) \hat{z}. \quad (2) $$

We illustrate different velocity fields in Figure 1: solid-body rotation ($v \propto R$), a flat velocity curve ($v = v_0$), and an irrotational fluid ($v \propto 1/R$). We show three fluid elements in two arbitrary positions along their orbits. On each panel we show the fluid elements at the same azimuthal angle. The aim is to compare how much a fluid element has rotated once it passes by the same position on the disk. In the case of solid-body rotation $\omega = 2\Omega\hat{z}$, where $\Omega(R) = v_0/R$, and patches of gas experience a local rotation of $\Omega$ with respect to a local inertial reference frame. For a flat velocity curve $\omega = \Omega\hat{z}$, and the local rotation is half the galactic rotation $\Omega/2$. We can see in the middle panel of Figure 1 that each fluid element has completed half the rotation of the left panel. We notice from Equation (2) that there is a critical case when $v(R) \propto 1/R$. For such a velocity field, $\omega = 0$ and fluid elements experience no local rotation with respect to an inertial reference frame. This kind of fluid is called irrotational, illustrated in the right panel of Figure 1.

There are two useful relations between vorticity and two quantities that are very helpful to have in mind. First, for a fluid with a circular velocity field $\mathbf{v} = v(R)\hat{\phi}$, the vorticity is

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Illustration of the circulation in a closed region for different velocity fields. Left: constant velocity field. Middle: shear across the x-axis, $v_x = 0$ and $v_y \propto x$. Right: solid-body rotation with angular velocity $\Omega_0$. The length of the arrows represents the magnitude of the velocity field.</figcaption></figure>

proportional to the local angular momentum $L \propto \omega$, as demonstrated in Appendix A. Second, for a circular velocity field the vorticity is related to $\Omega$ and $\kappa$ by $\omega = \frac{\kappa^2}{2\Omega}\hat{z}$. This implies that for an irrotational fluid $\kappa = 0$ and the Toomre parameter $Q = 0$. Any perturbation larger than the thermal Jeans scale, $\lambda_J \equiv c_s^2/\pi G\Sigma_{\text{gas}}$, is gravitationally unstable. It is noteworthy that in an irrotational fluid its angular velocity $\Omega \neq 0$ and its shear $\frac{\partial\Omega}{\partial r} \neq 0$. Given the relation between vorticity and local angular momentum, this is not a surprising implication. A parcel with no angular momentum does not have rotational support to halt gravitational collapse.

Unfortunately, vorticity is a local quantity, defined for an infinitesimal fluid element. For finite regions of space it is better to compute the fluid circulation $\Gamma$, which corresponds to a macroscopic measure of rotation (Pedlosky 1992). $\Gamma$ is defined as a line integral of the velocity field along a closed path,

$$ \Gamma = \oint_{\partial S} \mathbf{v} \cdot d\mathbf{l}. \quad (3) $$

For a continuous velocity field we can apply Stokes's theorem, which relates the line integral along a closed path to the surface integral over the area enclosed by it. This allows us to make the connection between $\Gamma$ and $\omega$:

$$ \Gamma = \oint_{\partial S} \mathbf{v} \cdot d\mathbf{l} = \int_S \nabla \times \mathbf{v} \cdot d\mathbf{S} = \int_S \omega \cdot d\mathbf{S}. \quad (4) $$

We can think of circulation as the area-weighted integral of the vorticity field. Figure 2 illustrates the circulation around a circular path for different velocity fields. A fluid with a constant velocity $v_0$ has $\Gamma = 0$, as it is constant in both magnitude and direction, and the line integral cancels out owing to the change of direction of the path with respect to the velocity field. In other words, bulk displacements make no contribution to $\Gamma$. A shear velocity field of the form $\mathbf{v} = (0, \Omega_0 x)$ has $\omega = \Omega_0$ and $\Gamma = \pi r^2 \Omega_0$. The last example shows solid-body rotation, $\mathbf{v} = (-\Omega_0 y, \Omega_0 x)$, with $\omega = 2\Omega_0$ and $\Gamma = 2\pi r^2 \Omega_0$.

However, realistic velocity fields are not completely smooth and behave differently at different scales. For instance, we might find that over a region $S$ the circulation is $\Gamma = 0$. But this does not imply that $\omega = 0$ over the whole region. Since the circulation is defined as a sum, if we divide $S$ into $N$ small subregions $s_i$, and $\Gamma = 0$, then

$$ \Gamma_S = \int_S \omega \cdot d\mathbf{S} = \int_{s_1 + \dots + s_N} \omega \cdot d\mathbf{S} = \sum_{i=1}^{N} \Gamma_{s_i} = 0. \quad (5) $$

The sum of all $\Gamma_i$ gives zero. There are infinite ways to distribute the values of $\Gamma_i$ to get $\Gamma = 0$. The exact distribution of the circulation at this smaller scale will depend on the nature
of the velocity field. This implies that only a multiscale measurement of the circulation can characterize the velocity field.

To have a full picture of the rotation of a fluid, we need to compute the circulation of gas at each point in the fluid on regions of different sizes. In this way we can create distributions of circulation at each spatial scale. To compare $\Gamma$ at different scales, let us define the normalized circulation $\gamma$:

$$ \gamma = \frac{\int_S \omega \cdot dS}{\int_S dS} = \frac{\Gamma}{A}, \qquad (6) $$

where A is the area of S. For solid-body rotation with angular velocity $\Omega_0$, $\gamma = 2\Omega_0$ for any fluid patch, and the distribution of $\gamma$ will be a Dirac delta function $\delta(\omega - 2\Omega_0)$. In the case of a rotating fluid with added random motions the distribution of $\gamma$ will be broader at small scales and will get narrower as we increase the size of the region in question, since we are adding random numbers and then dividing by a larger area. Hereafter and for simplicity we will refer to the normalized circulation simply as the circulation unless explicitly stated.

We can extend these ideas to the case of galactic dynamics. Imagine that the velocity field is composed by an ordered and smooth circular velocity field, $v_{\text{rot}}$, and a noncircular, random field, $v_{\text{nc}}$, i.e., $v = v_{\text{rot}} + v_{\text{nc}}$. This gets translated into two components of the circulation field $\gamma = \gamma_{\text{rot}} + \gamma_{\text{nc}}$. At large scales, $\gamma \approx \gamma_{\text{rot}}$, since the major contribution comes from galactic rotation and random or noncircular components cancel each other. At small scales, the distribution of $\gamma_{\text{nc}}$ gets broader, while the distribution of $\gamma_{\text{rot}}$ converges to $\omega_{\text{rot}}$ given by Equation (2).

In brief, the probability density function (pdf) of $\gamma$ follows pdf($\gamma$) $\simeq$ pdf($\gamma_{\text{rot}}$) at galactic scales and pdf($\gamma$) $\simeq$ pdf($\gamma_{\text{nc}}$) at small, parsec scales. This implies that at a particular scale $\lambda_{\text{eq}}$, $\gamma_{\text{nc}} = \gamma_{\text{rot}}$, i.e., the velocity fields $v_{\text{rot}}$ and $v_{\text{nc}}$ contribute equally to the measured circulation. Since $\lambda_{\text{eq}}$ depends on the local properties of $v_{\text{rot}}$ and $v_{\text{nc}}$, this scale has different values in different regions of a galaxy. For example, near the center of a galactic disk, $\omega_{\text{rot}}$ is higher and the transition to nonordered motions occurs at smaller scales, i.e., smaller $\lambda_{\text{eq}}$. In Section 2.2.5 we define how to compute $\lambda_{\text{eq}}$. We analyze in detail the behavior of $\lambda_{\text{eq}}$ for our set of simulations in Section 4.

## 2.2. The Method: Circulation as a Diagnostic of Gas Dynamics

We model the velocity field as the sum of two fields with different properties: the first being an axisymmetric and smooth velocity field, $v_{\text{rot}} = R\Omega(R)$, which is given by galactic rotation. The second field corresponds to a GRF, $v_{\text{nc}}$. GRFs are fields that follow a Gaussian distribution. In this paper, we will use a continuous GRF that is defined by a generating function in Fourier space that specifies the contribution from each spatial scale to the random velocity field (Lang & Potthoff 2011). These kinds of random fields are widely used in cosmology to model the primordial perturbations of the density field (Pranav et al. 2019). In $v_{\text{nc}}$ we are including any source of noncircular large-scale motions, such as outflows, collapse, turbulence, and other large-scale coherent motions like those induced by spiral arms and bars. While many of these motions are not expected to be random or Gaussian, this is a good first-order approximation for a statistical description. In future work

we can study better models for each component of the velocity field.

The aim is then to obtain the three fields, $\omega$, $\omega_{\text{rot}}$, and $\omega_{\text{nc}}$, from our simulations. We can compute $\omega$ directly from the simulations. To get $\omega_{\text{rot}}$, we need to choose how to model the smooth profile of the velocity field (i.e., the rotation curve). Finally, we will model the random component by means of a function in Fourier space on the spatial coordinates. We have to point out that since we are computing the vorticity field for a discrete grid, this field is also $\gamma$ at the resolution level.

### 2.2.1. Vorticity Field

We calculate $\omega(x, y)$ as follows. First, we compute the two-dimensional velocity field averaging along the z-axis:

$$ V(x, y) = \frac{\int_{-z_0}^{z_0} v(x, y, z) \rho(x, y, z) dz}{\int_{-z_0}^{z_0} \rho(x, y, z) dz}, \qquad (7) $$

where $\rho$ is the gas density, $v = (v_x, v_y, v_z)$ is the three-dimensional velocity field, and $V = (V_x, V_y)$ is the reduced two-dimensional field. We choose $z_0 = 1$ kpc over the whole galactic plane. Then $\omega(x, y)$ is given by

$$ \omega(x, y) = \frac{\partial V_y(x, y)}{\partial x} - \frac{\partial V_x(x, y)}{\partial y}. \qquad (8) $$

Note that we are considering all the gas in $z \in [-1 \text{ kpc}, 1 \text{ kpc}]$ to compute the integrated velocity fields, which is about 20 times the scale height of our simulated galaxies. We are not using a density threshold to integrate the velocity field. In observations, different tracers do not necessarily trace all the gas and are biased toward high-density regions.

### 2.2.2. Smooth Component

Since the definition of $\omega_{\text{rot}}(R)$ involves radial derivatives, we choose to parameterize $\omega_{\text{rot}}$ by an analytic function. To get $\omega_{\text{rot}}$, we fit a rotation curve of the form

$$ V_{\text{rot}}(R) = v_0 \arctan(R/R_1) \exp(-R/R_2) \qquad (9) $$

to the circular velocity field, and we apply Equation (2) to obtain $\omega_{\text{rot}}$. The arctan($x$) function provides a good fit to observed rotation curves (Courteau 1997), while the $\exp(-x)$ function recovers the decay of the rotation curve for our simulated galaxies. To fit the function in Equation (9), we divide the disk in radial bins of width 500 pc. For each radial bin we have a pair velocity uncertainty ($v_i, \delta v_i$), where $v_i$ is the median of the circular velocity and $\delta v_i$ is half of the difference between the 84th and 16th percentiles of the circular velocity field. Finally, we perform a least-squares optimization to fit the rotation curve. In Appendix D.1, we show how our results change using a different model for the rotation curve, e.g., adding a more sophisticated measure of the large-scale motions via harmonic decomposition. We find that our results are not very sensitive to the choice of the model of galactic rotation.

### 2.2.3. Random Component

The final step is choosing a model for $v_{\text{nc}}$ to obtain its contribution to the vorticity field. Here we choose $v_{\text{nc}}$ to be defined by a generating function $V(k)$ in Fourier space. The relation between $V(k)$ and the field $v_{\text{nc}}$ is derived in
Appendix B. For a two-dimensional field, $V(k)$ is related to the energy power spectrum $E(k)$ by

$$ \frac{1}{2} \langle v_{nc}^2 \rangle = \pi \int V(k)^2 k dk = \int E(k) dk. \quad (10) $$

This relation is shown in Equation (37). The function $V(k)$ is not unique for a whole galaxy. Each region of a galaxy is subjected to different conditions of stability, feedback, and dynamics, which will give rise to different noncircular velocity fields (turbulence, collapse, and gas flows). Given the geometry of disk galaxies, we expect that the dynamics of gas change across galactocentric radius. The simplest way to approach these differences is to separate the galaxy into radial bins. Then, each radial bin will be described by its own function $V(k)$.

We choose $V(k)$ to be of the form

$$ V(k) \propto \begin{cases} k^{-n_1} & \text{if } k_{\min} < k < k_c \\ k^{-n_2} & \text{if } k_c \le k < k_{\max}, \\ 0 & \text{elsewhere} \end{cases} \quad (11) $$

where the wavenumber $k$ is related to the spatial scale as $k = \lambda^{-1}$ (no $2\pi$ factor). We fix $k_{\min}$ and $k_{\max}$ to $4/L$ and $1/(4\Delta x)$, respectively, where $L$ is the box size and $\Delta x$ is the spatial resolution of the two-dimensional field, 30 pc for the simulations described in Section 3.1. This limits the dynamic range of $V(k)$ between 120 pc and 10 kpc (choosing $L = 40$ kpc). If we add the constraint of continuity for $V(k)$ at $k_c$, we need a parameter to set the amplitude of $V(k)$, which translates into the amplitude of the velocity field. We choose this last parameter to be the characteristic velocity dispersion of the random velocity field, given by

$$ \sigma_0^2 = 2\pi \int_{k_{\min}}^{k_{\max}} V(k)^2 k dk. \quad (12) $$

This is a property of GRFs (Lang & Potthoff 2011). To resume, the parameters defining the velocity field are ($n_1, n_2, k_c, \sigma_0$). The vorticity of this velocity field can be calculated using Equation (1) or directly from the parameters using Equation (39). We denominate the velocity and vorticity fields obtained from the $V(k)$ as $v_{nc}$ and $\omega_{nc}$. In summary, our model for the vorticity field is

$$ \omega(x, y) = \omega_{\text{rot}}(R) + \omega_{\text{nc}}(n_1(R), n_2(R), k_c(R), \sigma_0(R)). \quad (13) $$

We test two additional models for $V(k)$ with $n_1 = n_2$ i.e., a single power law, shown in Appendix D.2. In one model $k_{\max}$ is fixed to $(4\Delta x)^{-1}$, while in the second we allow $k_{\max}$ to vary. A single power law with fixed $k_{\max}$ that best fits the noncircular component of $\gamma$ does not match the behavior of $\gamma_{nc}$ as a function of scale. This implies that it is not a good representation of $\gamma_{nc}$. On the other hand, a single power law with variable $k_{\max}$ gets similar results compared to the piecewise power law used in this work. This is consistent with the possible values that we find for $n_2$ in Section 3.3 (i.e., it mimics the adopted model with a high value of $n_2$).

### 2.2.4. Distribution of $\gamma$

To find the best values for $n_1, n_2, k_c$, and $\sigma_0$, we measure the distribution of the circulation $\gamma(x, y, \ell)$

$$ \gamma(x, y, \ell) = \frac{1}{\ell^2} \int \omega(x, y) dS \quad (14) $$

at each scale $\ell$, integrating over square regions of size $\ell \times \ell$ centered on points $(x, y)$. According to our model, $\gamma = \gamma_{\text{rot}} + \gamma_{\text{nc}}(n_1 n, n_2, \lambda_c, \sigma_0)$. Since $v_{nc}$ is a GRF, $\omega_{nc}$ and $\gamma_{nc}$ are GRFs too. In Equation (B9) we show that the variance of $\gamma_{nc}(x, y, \ell)$ is determined by $\ell$ and the generating function $V(k)$. Then, we only need to compute the variance of $\gamma_{nc}$ as a function of $\ell$ and use Equation (B9) to find the parameters of $V(k)$.

### 2.2.5. The Scale at Which Gas Noncircular Motions Start to Dominate

The last and most relevant quantity in our framework is $\lambda_{eq}$, the scale at which the contributions from $\gamma_{\text{rot}}$ and $\gamma_{nc}$ to the measured circulation of gas are roughly the same. For $\ell > \lambda_{eq}$ the circulation of gas is dominated by galactic rotation, while at scales $\ell < \lambda_{eq}$ it is dominated by noncircular motions.

We define the scale $\lambda_{eq}$ in the following way. At each scale $\ell$ we measure the ratio $f_{\gamma}(\ell) = \sum \gamma_{nc, \ell}^2 / \sum \gamma_{\text{rot}, \ell}^2$. Then, $\lambda_{eq}$ is given by the equation

$$ f_{\gamma}(\lambda_{\text{eq}}) = 1.0 \quad \text{with} \quad f_{\gamma}(\ell) = \frac{\sum \gamma_{nc, \ell}^2}{\sum \gamma_{\text{rot}, \ell}^2}. \quad (15) $$

We compare their squared values since $\gamma$ can have negative values. For a random variable $x$ with mean value $\mu_x$ and standard deviation $\sigma_x$, the expected value of $\sum x^2$ is $\mu_x^2 + \sigma_x^2$. We can rewrite $f_{\gamma}(\ell)$ as

$$ f_{\gamma}(\ell) = \frac{\mu_{\gamma_{nc}}^2(\ell) + \sigma_{\gamma_{nc}}^2(\ell)}{\mu_{\gamma_{\text{rot}}}^2(\ell) + \sigma_{\gamma_{\text{rot}}}^2(\ell)} = \frac{\sigma_{\gamma_{nc}}^2(\ell)}{\mu_{\gamma_{\text{rot}}}^2(\ell) + \sigma_{\gamma_{\text{rot}}}^2(\ell)}, \quad (16) $$

where $\mu_{\gamma_{nc}}, \mu_{\gamma_{\text{rot}}}, \sigma_{\gamma_{nc}}$, and $\sigma_{\gamma_{\text{rot}}}$ are the mean values and standard deviations of $\gamma_{nc}$ and $\gamma_{\text{rot}}$ as a function of $\ell$. For a region with constant $\omega_{\text{rot}}$, $f_{\gamma}(\lambda_{\text{eq}}) = 1.0$ is equivalent to the equation $\sigma_{\gamma_{nc}}(\ell = \lambda_{\text{eq}}) = \omega_{\text{rot}}$. Given the parameters that define $\gamma_{\text{rot}}$ and $\gamma_{nc}$, $\lambda_{\text{eq}}$ is a function of the form $\lambda_{\text{eq}}(\omega_{\text{rot}}, n_1, n_2, \lambda_c, \sigma_0)$.

### 2.2.6. Deriving $\gamma_{nc}$ and $\lambda_{eq}$

In practice, we are looking for the parameters $n_1, n_2, k_c$, and $\sigma_0$ that best represent the equation $\gamma = \gamma_{\text{rot}} + \gamma_{nc}$. However, for a given set of parameters, the field $\gamma_{nc}$ is a random realization from a parent GRF and is not single valued. This means that the expression $\gamma = \gamma_{\text{rot}} + \gamma_{nc}$ has to be considered as the sum of two distributions rather than the sum of two fields or images. Then, to look for parameters that can model our data, we need to compare $\gamma, \gamma_{\text{rot}}$, and $\gamma_{nc}$ as distributions. We choose to compare histograms of $\gamma$ and $\gamma_{\text{rot}} + \gamma_{nc}$ as the final ingredient of our technique.

Our method can be summarized as follows:

1. Measure the two-dimensional vorticity field $\omega(x, y)$.

2. Model the large-scale and axisymmetric component of the velocity field, $V_{\text{rot}}$, using Equation (9) and compute $\omega_{\text{rot}}$.
3. Divide the disk into different radial annuli.

4. From $\omega$ and $\omega_{\text{rot}}$, measure the distributions of $\gamma$ and $\gamma_{\text{rot}}$ at each scale $l$ within each radial annulus.

5. At each spatial scale we compare the distributions of $\gamma$ and $\gamma_{\text{rot}} + \mathcal{N}(0, \sigma_{\gamma})$, where $\mathcal{N}(0, \sigma_{\gamma})$ is a random field with dispersion $\sigma_{\gamma}$. We fit $\sigma_{\gamma}$ using the least-squares method. This step creates an array $\sigma_{\gamma}(l)$ as a function of $l$ with its respective uncertainty $\epsilon_{\sigma}$.

6. Explore the parameter space $\theta = (n_1, n_2, k_c, \sigma_0)$ of the function $\mathcal{V}(k)$. Each parameter vector $\theta$ defines a different curve $\sigma_{\gamma_{\text{nc}}}(l, \theta)$. To find the posterior distributions of $\theta$ given our data $\mathbf{D}$, we use Bayes's theorem:

$$P(\theta|\mathbf{D}) = \frac{P(\mathbf{D}|\theta)P(\theta)}{P(\mathbf{D})} \qquad (17)$$

In this equation, $P(\mathbf{D}|\theta)$ is the likelihood to obtain $\mathbf{D}$ that in our case corresponds to the array $\sigma_{\gamma}(l)$. Our likelihood is given by

$$P(\mathbf{D}|\theta) = \prod_{l} \frac{1}{\sqrt{2\pi} \epsilon_{\sigma}} \exp\left(\frac{-(\sigma_{\gamma}(l) - \sigma_{\gamma_{\text{nc}}}(l, \theta))^2}{2\epsilon_{\sigma}^2}\right). \quad (18)$$

$P(\theta)$ is our prior knowledge of the parameters. We assume uniform prior distributions for each parameter. $P(\mathbf{D})$ is the Bayesian evidence of the data that ensures proper normalization. To sample the posterior distributions, we use Markov Chain Monte Carlo methods. We use 72 random walkers that are updated using the Metropolis-Hastings algorithm. Each random walker creates a chain with 15,000 values of $\theta$. For our analysis we ignore the first 5000 steps. This step creates samples of $\theta = (n_1, n_2, k_c, \sigma_0)$, from which we reconstruct the pdf's of the model parameters. This sample also establishes the parent distribution of $\gamma_{\text{nc}}$.

7. From $\gamma_{\text{rot}}$ and the distributions of $\gamma_{\text{nc}}$ we derive the distribution of the scale $\lambda_{\text{eq}}$ at each radial annuli.

### 2.2.7. Choice of Spatial Bin Size and Estimation of Uncertainties

To build the histograms and to compare them, we must choose the width of the bins, $\Delta\gamma$, and the uncertainty, $\epsilon_{\gamma}$, for each measurement with their propagation in the histogram bins. For the bin width we choose a conservative criterion: 8 times the bin width set by the Freedman-Diaconis rule (Freedman & Diaconis 1981), $\Delta\gamma = 16 \text{ IQR}(\gamma) N^{-1/3}$, where IQR($\gamma$) is the interquartile range of $\gamma$ and $N$ is the total number of data points. The Freedman-Diaconis rule attempts to minimize the integrated mean squared difference of the histogram model and the true underlying density.

For the uncertainty in $\gamma$ we set an uncertainty of $\epsilon_{\gamma} = 1 \text{ km s}^{-1}$ in the measured velocity at the resolution of the simulation $\Delta x$ ($\approx 30 \text{ pc}$), which is comparable to the precision of recent gas velocity measurements on nearby galaxies (Druard et al. 2014; Caldú-Primo & Schruba 2016; Koch et al. 2018; Sun et al. 2018). To propagate the uncertainty, we use Equation (3). For a square region with area $\Delta x \times \Delta x$, i.e., at the maximum resolution, the circulation is the sum of the integral along the four faces of the square. Then, $\Gamma = v_1 \Delta x + v_2 \Delta x + v_3 \Delta x + v_4 \Delta x$ (with the respective signs due to the dot product) and its uncertainty is $\epsilon_{\Gamma} = \sqrt{4\epsilon_{\gamma}^2 \Delta x} = 2\epsilon_{\gamma} \Delta x$. For $\gamma$ the uncertainty is $\epsilon_{\gamma} = 2\epsilon_{\gamma}/\Delta x$. A square region of size $l = N \Delta x$ is delimited by $4N$ linear segments, $N$ at each side. Then, $\Gamma = \sum_{i=1}^{4N} v_i \Delta x$ and its

uncertainty is $2\sqrt{N} \epsilon_{\gamma} \Delta x$, while for $\gamma$ it is $2\sqrt{N} \epsilon_{\gamma} \Delta x / (N \Delta x)^2 = (2\epsilon_{\gamma}/\Delta x) \times N^{-3/2} = \epsilon_{\gamma} N^{-3/2}$. The scaling of the uncertainty of $\gamma$ with $N$ is not straightforward if we use Equation (4). We have to recall that $\omega$ is the difference between two terms. If we add the vorticity of two neighbor cells, we are also subtracting the line integral along the line that both regions share. For that reason, when we compute the uncertainty in $\gamma$ for a region with $N \times N$ elements, we are adding $4N$ terms instead of $N \times N$.

The histogram counts of our distributions have Poisson noise, i.e., an uncertainty of $\sqrt{N}$, where $N$ is the number of data points lying in a given bin. We also have to propagate the uncertainties in $\epsilon_{\gamma}$ into the histogram. Let us consider the $j$th bin with endpoints $[l_j, u_j]$ and a data point with value $\gamma_i$ and uncertainty $\epsilon_i$. Assuming that $\gamma_i$ is the mean value of a Gaussian random variable with standard deviation $\epsilon_i$, the probability $p_{ij}$ that this data point lies within $[l_j, u_j]$ is

$$p_{ij} = \frac{1}{2} \left[ \text{erf}\left(\frac{u_j - \gamma_i}{\sqrt{2} \epsilon_i}\right) - \text{erf}\left(\frac{l_j - \gamma_i}{\sqrt{2} \epsilon_i}\right) \right], \quad (19)$$

where erf is the error function. Each bin acts like a Bernoulli random variable: we add 1 if the measurement lies in the bin or zero otherwise, with a probability $p_{ij}$ and $1 - p_{ij}$, respectively. The variance for the Bernoulli distribution is $p_{ij}(1 - p_{ij})$. Then, the total variance in the $j$th bin due to the uncertainties in $\gamma$ is $\sum_{i}^{N} p_{ij}(1 - p_{ij})$.

### 2.3. Parameter Distributions

To sample the distributions for $n_1, n_2, k_c$, and $\sigma_0$, we use the Metropolis-Hastings algorithm to generate Markov chains using Equation (17). To compute the likelihood function in Equation (18), we need to calculate the integral in Equation (B9) each time we test a new set of parameters. Since this step is computationally expensive, we divide the subspace ($n_1, n_2, k_c$) into a $72^3$ grid over which we pre-tabulate the integral. The parameter $\sigma_0$ works as a normalization of $\mathcal{V}(k)$ and can be handled independently. The intervals chosen are $n_1 \in (0.8, 2.5)$, $n_2 \in (2.0, 20.0)$, and $k_c \in (2k_{\text{min}}, k_{\text{max}}/2)$.

## 3. Application to Simulated Galaxies

We start this section by introducing the simulations used as a test bed for our method.

### 3.1. Simulations

We test our technique on three hydrodynamical simulations of disk galaxies. In each simulation we apply our circulation-based method on nine radial annuli, 3 kpc wide, centered on 1.5, 3.0, 4.5, 6.0, 7.5, 9.0, 10.5, 12.0, and 13.5 kpc. We vary the scale $l$ from the maximum resolution of the simulations 30 pc–5 kpc. Notice that $l$ can be larger than the width of a radial annulus. We might argue that within a radial annulus there is no information of scales larger than the width of the radial bin. However, the function $\mathcal{V}(k)$ has information about the correlation between two points in the velocity field, and within each radial annulus we can find points separated by distances larger than 3 kpc. The center of each square region is inside the 3 kpc annulus, but it can cover cells outside the annulus. Information from neighbor regions will affect the values of the parameters within each annular region. This might smooth the resulting radial profiles of our model parameters.
To create this set of simulated galaxies, we use the adaptive mesh refinement code Enzo (Bryan et al. 2014). We run simulations of spiral galaxies for 700 Myr with a coarse resolution of 3.7 kpc. We use two criteria to refine a given gas cell, both of which have to be fulfilled: refinement by baryon mass, and Jeans length. The Jeans length is at least resolved by four cells to prevent artificial fragmentation (Truelove et al. 1997). During the first 500 Myr, we use a maximum resolution of 60 pc until a quasi-steady state is reached. Then, the resolution is increased to 30 pc for 200 Myr. Over 80% of the mass in gas cells is found at the highest resolution. To obtain images and the velocity fields of the simulated galaxies, we use the yt python package<fnref n="13" /> (Turk et al. 2011).

These simulations are modeled as a four-component system that includes gas, star particles, and time-independent stellar and dark matter potentials. The time-independent stellar potential is given by the Miyamoto–Nagai profile (Miyamoto & Nagai 1975), with a gravitational potential $\Phi_{\text{stellar}}$ given by

$$ \Phi_{\text{stellar}} = - \frac{GM_{\text{stellar}}}{\sqrt{R^2 + (a + \sqrt{z^2 + b^2})^2}}, \quad (20) $$

where $M_{\text{stellar}}$ is the total stellar mass of the field and a and b are characteristic length scales. We choose $a = 5$ kpc and $b = 200$ pc, which are similar to the fitted values for the Milky Way stellar disk (Kafle et al. 2014). For the DM potential we use the Navarro–Frenk–White profile (Navarro et al. 1997). The virial mass of the DM potential and the total stellar mass of the stellar potential are different for each simulation, which will be defined later. Note here that the external potential is axisymmetric, and we do not include spiral structure as done in other works (e.g., Dobbs et al. 2015). The gaseous disks start with an exponential radial profile with a radial scale length of 3.5 kpc.

Star particles are formed when (i) the local number density $n_{\text{cell}}$ is greater than $100 \text{ cm}^{-3}$, (ii) the velocity field is converging $\nabla \cdot \mathbf{v} < 0$, (iii) the cooling time is shorter than the local freefall time, and (iv) the gas mass in the cell $m_{\text{cell}}$ is greater than the Jeans mass, $m_{\text{cell}} > m_{\text{J}} = \frac{4\pi}{3}(\lambda_{\text{J}}/2)^3$. If these criteria are satisfied, a star particle is formed with a mass equal to $0.1m_{\text{cell}}$. The typical mass of a star particle is of the order of $1 \times 10^4 M_{\odot}$ and represents a population of stars. Our simulations match the Kennicutt–Schmidt relation (Kennicutt 1998).

We include stellar feedback from supernovae (SNe), H II regions, and momentum injection from massive stars. To compute the energy or momentum injection as a function of time, for each type of stellar feedback we use tabulated results from STARBURST99 (Leitherer et al. 1999) assuming a Kroupa initial mass function, solar metallicity, and instantaneous star formation. We model SN feedback by injecting $10^{51}$ erg of thermal energy per every 55 $M_{\odot}$ of stars formed. This is 1.6 times more energy than the simulations of the Agora Project (Kim et al. 2016). Since star particles represent a population of stars, the energy is deposited continuously at the cell where the star particle lies. To add the effects of H II regions, we follow the approach of Renaud et al. (2013) and Goldbaum et al. (2016), heating the gas up to $10^4$ K within the Strömgren radius. If the volume of the Strömgren sphere, $V_S$, is

<table><caption>**Table 1**
Simulation Parameters</caption><thead><tr><th>Run</th><th><i>M</i><sub>gas</sub><br>(<i>M</i><sub>☉</sub>)</th><th><i>M</i><sub>*</sub><br>(<i>M</i><sub>☉</sub>)</th><th><i>M</i><sub>200</sub><br>(<i>M</i><sub>☉</sub>)</th></tr></thead><tbody><tr><td>G2E1</td><td>2 × 10<sup>10</sup></td><td>1 × 10<sup>11</sup></td><td>8 × 10<sup>11</sup></td></tr><tr><td>G1E1</td><td>1 × 10<sup>10</sup></td><td>1 × 10<sup>11</sup></td><td>8 × 10<sup>11</sup></td></tr><tr><td>G1E0.5</td><td>1 × 10<sup>10</sup></td><td>5 × 10<sup>10</sup></td><td>4 × 10<sup>11</sup></td></tr></tbody></table>

smaller than the cell volume, $V$, only a fraction $V_S/V$ of the thermal energy is deposited in the corresponding cell. In our simulations, most of the time $V_S < V$, and on average the gas is heated up to 7000 K. To compute the momentum injection, we consider stellar winds and radiation pressure (Agertz et al. 2013). To account for the radiation pressure, we compute the bolometric luminosity $L_{\text{bol}}$ for each active stellar particle, and we distribute the momentum $L_{\text{bol}}/c$ evenly in the six nearest cells. For stellar particles separated by one cell this causes some cancellation of the injected momentum (Hopkins & Grudić 2019). We underestimate the effect of radiation since we do not consider the scattering of IR photons. For our simulations, and assuming an IR opacity if $\kappa_{\text{IR}} = 10 \text{ cm}^2 \text{ g}^{-1}$, the optical depth of IR radiation is usually $\tau_{\text{IR}} \simeq 0.2$. To model the energy lost by radiation, we use the cooling curves of Sarazin & White (1987) for temperatures $T \ge 10^4$ K and those of Rosen & Bregman (1995) for $300 \text{ K} < T < 10^4$ K. This imposes a minimum value for the Jeans scale $\lambda_{\text{J}}$. For a surface density of $10 M_{\odot} \text{ pc}^{-2}$ and a temperature $T = 300 \text{ K}$ the $\lambda_{\text{J}}$ is of the order of 100 pc. This means that overdensities in our simulations are more representative of H I clouds.

We run a second set of simulations using only SN feedback to explore the effects of changing the feedback prescription. This second set has a higher density threshold of $2800 \text{ cm}^{-3}$, necessary to match the Kennicutt–Schmidt relation (Kennicutt 1998). Except from the expected decrease in the magnitude of noncircular motions due to the injection of less energy on small scales, most of the conclusions from the previous set of simulations hold also for these simulations, so they are presented and discussed in Appendix E.

The three simulations discussed here are designed as follows: run G2E1 corresponds to our reference simulation, with $2 \times 10^{10} M_{\odot}$ of gas, a stellar potential of $10^{11} M_{\odot}$, and a halo mass $M_{200} = 8 \times 10^{11} M_{\odot}$. The masses for the stellar and DM component are similar to Milky Way values (Kafle et al. 2014). Run G1E1 is identical to the reference run in all parameters, except that the gas mass is reduced by half (i.e., it has a lower gas fraction). Run G1E0.5 has the same gas fraction and disk scale length as the reference run, but half the mass in all components (stars, gas, and DM), therefore being a low surface density version of the reference run. For G2E1 and G1E1 the concentration parameter of the DM potential is set to $c = 21$, and for G1E0.5 it is tuned to maintain the same shape of the rotation curve, although the normalization can be different. The relevant parameters for these simulations are shown in Table 1. In our nomenclature G stands for the amount of gas and E for the magnitude of the external potential. We have to point out that these models sit a factor of three above the stellar mass–halo mass relation. These add changes in the magnitude of vertical acceleration, but we do not expect to produce major changes in the two-dimensional velocity field.

In Figure 3 we show projections of the gas density field across the z-axis. These galaxies do not present grand-design

<footnote n="13">¹³ yt: http://yt-project.org.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Gas surface densities integrated along the z-axis at $t = 700$ Myr. Dashed white circles delimit the region inside 15 kpc, which corresponds to the maximum radius of the defined annuli. The dotted line shows the maximum radius of the cells that are included in the analysis of circulation, which corresponds to 18.54 kpc.</figcaption></figure>

spiral patterns. Near the galactic center we see clumps with different sizes for each simulation, particularly larger clumps for G2E1 that show signs of tidal interactions. Structures in G1E0.5 appear to be less affected by shear owing to the lower magnitude of its rotation curve. The dashed white circle of 15 kpc radius in Figure 3 shows the outer edge of the outermost radial annulus where we measure the circulation of gas. When we measure the circulation within square regions of size 5 kpc and centered on a point at radius 15 kpc, we include points that are up to a distance of 18.54 kpc from the galactic center. We show a 18.54 kpc radius dotted circle for illustrative purposes in Figure 3.

## 3.2. Galactic Rotation

We show the rotation curves $v_{\text{rot}}(R)$ of each run in the top panel of Figure 4. The shaded regions show the variations in the tangential velocity field. The analytic rotation curves, obtained by fitting Equation (9), are shown as dashed lines. In the bottom panel of Figure 4 we also show the resulting vorticity $\omega_{\text{rot}}(R)$. The vorticity coming from galactic rotation, $\omega_{\text{rot}}(R)$, decreases as we get far from the galactic center. This means that $\gamma_{\text{rot}}$ decreases with galactocentric radius. If the parameters of $V(k)$, which define the magnitude of $\gamma_{\text{nc}}$, were constant across galactocentric radius, the relative contribution from $\gamma_{\text{rot}}$ to the total circulation $\gamma$ would decrease with galactocentric radius, and the spatial scale $\lambda_{\text{eq}}$ should increase with radius.

## 3.3. Characterization of Random Motions

To analyze our simulations, we choose a box size of $L = 40$ kpc. This implies that $k_c$ is bounded between values $(2k_{\text{min}} = 40 \text{ kpc}/8)^{-1} = (5 \text{ kpc})^{-1}$ and $\frac{1}{2}k_{\text{max}} = (8 \times 30 \text{ pc})^{-1} = (240 \text{ pc})^{-1}$. In Figure 5 we show the probability distributions of $n_1, n_2, \lambda_c$, and $\sigma_0$, where $\lambda_c = 1/k_c$ corresponds to the characteristic scale of the model $V(k)$. Keep in mind that, given the shape of the chosen function $V(k)$, the parameters $n_1, n_2$, and $\lambda_c$ are correlated parameters. Runs G2E1 and G1E1 show variations of $n_1$ between 1.0 and 1.9. The distribution of $\lambda_c$ ranges from its minimum value of ~240 pc to 1.1 kpc. In the annulus of G1E1 centered at 4.5 kpc $\lambda_c$ is unresolved. The parameter $\sigma_0$ shows narrow distributions. G2E1, the most massive galaxy, shows values of $\sigma_0$ above 30 km s$^{-1}$ at all radii. The distribution of $n_2$ gets flat for values over 5. The model is not sensitive to variations of $n_2$ over $n_2 = 5$. Equation (38) shows the contribution from each scale and the amplitude of the vorticity field through an

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Profiles of rotational velocities and large-scale vorticity. Top: radial profiles of the circular velocity $v_{\text{rot}}$, which are used to compute the large-scale vorticity component $\omega_{\text{rot}}$. Solid lines show the median values of the rotation curve. Shaded regions represent $1\sigma$ uncertainties due to variations in the velocity field. The analytic models of the rotation curves are shown as the dashed lines. Bottom: radial profiles of vorticity $\omega_{\text{rot}}(R)$. The solid and dashed lines correspond to the median values and the analytic models, respectively. The gray dotted line shows a function proportional to $\exp(-R/5\text{kpc})$ for comparison.</figcaption></figure>

integral of $V(k)$. As the values of $n_2$ increase, the contribution from $k > k_c$ to the vorticity field starts to get smaller. This suggests that the function $V(k)$ can be approximated as a single power law with a cut at $k_c$, i.e., $V(k > k_c) = 0$ as discussed in Section 2.2.3. We show in Figure 6 detailed distributions for G2E1. The off-diagonal plots of Figure 6 show two-dimensional histograms of the model parameters that help to visualize the correlation between these parameters. We can see a high
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Posterior distributions of model parameters. From top to bottom we plot the distributions for each of the three simulations. From left to right we show normalized distributions of $n_1$, $n_2$, $\lambda_c = 1/k_c$, and $\sigma_0$. The distribution for each annulus is shown with a characteristic color, going from purple to yellow with increasing galactocentric radius. The vertical black dashed line in the distributions for $\lambda_c$ corresponds to eight times the resolution of the simulations, the minimum allowed value for $\lambda_c$.</figcaption></figure>

correlation between $\lambda_c$ and $n_1$ and also a correlation between $\lambda_c$ and $n_2$ for low values of $n_2$. The explored range of values for $n_2$ in Figure 6 is extended to $n_2 = 50$ to show that its distribution remains uniform beyond $n_2 = 20$. For the annular regions, 1.5–4.5 kpc and 3–6 kpc in G2E1, 0–3 kpc in G1E1, and 3–6 kpc in G1E0.5, the distributions of $n_2$ show peaks in their lower limits imposed by our prior.

Radial variations of the parameters are summarized in Figure 7. If we look at the first and second panels of Figure 7, we can see that $n_1$ and $\lambda_c$ are anticorrelated for runs G2E1 and G1E1. This is also true for G1E0.5 but is not noticeable in Figure 7. G1E1 and G1E0.5 show similar values of $\sigma_0$ at large galactocentric radius besides having different rotation curves. Their profiles of $\sigma_0$ have large magnitudes compared to the velocity dispersion profiles measured in nearby galaxies from CO emission lines (Sun et al. 2018). For velocity dispersions derived from H I in Mogotsi et al. (2016) our profiles are also higher. However, derived velocity dispersion profiles in Romeo & Mogotsi (2017) for some nearby galaxies show similar magnitudes, reaching up to $50 \text{ km s}^{-1}$ in H I and CO. In our model, $\sigma_0$ models the velocity dispersion of the whole annular region, with velocities measured at the maximum resolution and without a density cut. Then, $\sigma_0$ has not to be understood as the average velocity dispersion for clouds in an annular region.

### 3.4. Scales at Which Gas Dynamics Transitions from Galactic Rotation to Noncircular Motions

What is the role of galactic rotation on small scales? We want to know down to which scales galactic rotation still

dominates the dynamics of gas or even molecular clouds. In our framework this information is encapsulated in the scale $\lambda_{eq}$, the scale at which the contributions to the circulation field from galactic rotation and noncircular motions are roughly the same.

The bottom panel of Figure 7 shows the radial profiles of $\lambda_{eq}$. Within 8 kpc from the galactic center, $\lambda_{eq}$ increases with galactocentric radius and varies between the resolution limit of 30 pc and 3 kpc. This shows that as we get farther from the galactic center, gas dynamics at the scale of clouds is predominantly dominated by noncircular motions.

In the radial annulus centered at 7.5 kpc, G2E1 and G1E1 have similar values of $\omega_{rot}$, $n_1$, and $\lambda_c$, while G2E1 has a higher $\lambda_{eq}$ by about a factor of 2. This suggests that for the same rotation curve differences in $\lambda_{eq}$ are mainly driven by differences in $\sigma_0$. The fundamental change between these two simulations is their gas surface density. Run G2E1 has the largest values of $\sigma_0$, and likewise it has the largest values of $\lambda_{eq}$. On the other hand, G1E1 and G1E0.5 show similar profiles for $\sigma_0$, but $\lambda_{eq}$ is larger for G1E0.5. This illustrates the effect of the rotation curve, which for G1E0.5 has a lower magnitude. In the central region of G1E1, $\lambda_{eq}$ goes to zero, below the resolution of our simulations. This means that $\lambda_{eq}$ is not resolved in these regions and that galactic rotation is the dominant source of circulation down to the resolution limit.

### 3.5. Distribution of $\gamma$

Now we show how the measured distributions of $\gamma$ change with spatial scale $\ell$ and how our model compares with them. Figure 8 shows percentiles of the circulation $\gamma$ and the model
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Distributions of parameters $n_1$, $n_2$, $\lambda_c = 1/k_c$, and $\sigma_0$ for run G2E1, in the radial annulus [4.5 kpc – 7.5 kpc]. Purple histograms in the diagonal show the marginal posterior distributions for $n_1$, $n_2$, $\lambda_c = 1/k_c$, and $\sigma_0$ in descending order. Off-diagonal plots show two-dimensional histograms of the model parameters. White contours show the 68%, 95%, and 99.7% confidence intervals. In these plots the explored range of values for $n_2$ has been extended to $n_2 = 50$.</figcaption></figure>

$\gamma_{\text{rot}} + \gamma_{\text{nc}}$ as a function of scale for three of the nine annuli. Each percentile is shown with a different color. The percentiles of the measured circulation are shown as solid lines, while the models $\gamma_{\text{rot}} + \gamma_{\text{nc}}$ are shown by the dotted lines with their respective $1\sigma$ intervals.

Let us first discuss the general characteristics of these distributions. At the smallest spatial scale $\gamma$ is equal to the vorticity measured at the spatial resolution of the simulation. Near the galactic center the distributions of $\gamma$ are much broader at every scale $\ell$. This is also true for the large-scale component of circulation, $\gamma_{\text{rot}}$. As shown in Figure 4, the slope of $\omega_{\text{rot}}$ decreases with radius, which means that within an annulus variations in $\omega_{\text{rot}}$ also decrease with radius. We can see how the width of the distribution of $\gamma$ changes from large to small scales: toward small scales it gets broader as the influence of noncircular motions becomes more important, while above scales of hundreds of parsecs it starts to converge toward a constant level, set by the galactic rotation component. Figure C2 in Appendix C shows how the distribution width of $\gamma$ looks for coherent rotation and a random field as a function of scale.

Figure 4 shows that $\omega_{\text{rot}}$ and consequently the pdf of $\gamma_{\text{rot}}$ are always positive.<fnref n="14" /> At galactic scales, $\gamma$ is greater than zero since $\gamma \simeq \gamma_{\text{rot}}$. On the other hand, the distribution of $\gamma_{\text{nc}}$, the

GRF component, is half positive/half negative at any scale. At the scales where noncircular motions start to become important, $\gamma$ starts to show negative values. This departure to negative values is not the same for every region; it depends on the magnitude of $\gamma_{\text{rot}}$ and the dispersion of $\gamma$ at the smallest scales, which depends on $\gamma_{\text{nc}}$. By looking at the percentile curves, we see that the percentage of regions with retrograde rotation varies between 20% and 40% at the smallest scales, with the highest fractions in G2E1.

We show the scales $\lambda_c$ and $\lambda_{\text{eq}}$ in Figure 8. For scales smaller than $\lambda_c$ the rate at which the distribution broadens starts to decrease until it stops. This is also illustrated in the examples of GRFs in Figure C1 in Appendix C. With regard to $\lambda_{\text{eq}}$, we can see how $\lambda_{\text{eq}}$ shifts from left to right depending on the average value of $\gamma$ at large scales and its variance at the largest and smallest scales.

We see that our model reproduces the shape of the distributions of $\gamma$ as a function of $\ell$, with some discrepancies at both extremes of the distributions. Figure 8 also displays $1\sigma$ uncertainties around the median value for our model of $V(k)$. Best agreement is seen at large galactocentric radii, where the distributions are better sampled since the number of cells in each annulus increases with radius. Near the galactic center we expect to observe large variations due to low sampling. In addition, it is more difficult to set a well-defined center of large-scale rotation near the center of the galaxy, since the

<footnote n="14">¹⁴ According to the chosen orientation of the z-axis.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Radial profiles of the median values of the model parameters $n_1$, $\lambda_c = 1/k_c$, $\sigma_0$, and $\lambda_{eq}$. Each run is shown with a different color. Error bars represent the 16th–84th percentile interval for each parameter.</figcaption></figure>

interactions with small structures can be comparable to or greater than the large-scale gravitational influence of the galaxy. Regardless, Figure 8 shows that our model can reproduce the trends in the measured distribution of $\gamma$ ($\ell$). This supports our assumption that the velocity field can be separated as the sum of large-scale circular motions and a GRF representing noncircular motions.

Although we are able to capture the general behavior of $\gamma$, there are noticeable discrepancies. In Figure 8 we see a systematic discrepancy at the 90th percentile, with measured values greater than the model. In some regions we also see discrepancies at the 10th percentile. In the first panel of G1E1 we see less symmetry with respect to the median and a

distribution that is broader than our model. This shows that the extreme values of $\gamma$ coming from noncircular motions are not represented in our model. One possibility is that within an annulus the model parameters change quickly. In our model, we are assuming a unique Gaussian distribution for each radial bin instead of a superposition of Gaussian distributions for each radius. However, in this scenario the distributions of $\gamma$ should be more symmetric. In some regions the deviation from the model of the 90th percentile is higher compared to the 10th percentile.

Some of the discrepancies could be explained by regions under collapse: when high-density regions collapse, their vorticity $\omega$ increases in magnitude. Kruijssen et al. (2019) found that simulated molecular clouds with higher densities show higher velocity gradients. By checking the velocity divergence in the x-y plane, i.e., $\nabla_{xy} \cdot \mathbf{v} = \frac{dv_x}{dx} + \frac{dv_y}{dy}$, we find that between the 20th and 80th percentiles of $\omega$ the median value of $\nabla_{xy} \cdot \mathbf{v}$ is close to zero. However, at the extremes of the distribution of $\omega$, $\nabla_{xy} \cdot \mathbf{v}$ drifts to negative values. Once a cloud of gas starts to collapse, $\nabla \cdot \mathbf{v} < 0$, the magnitude of $\omega$ increases. It might be possible to address this discrepancy by considering the conservation of angular momentum or the conservation of circulation, but that is beyond the scope of this work.

# 4. Discussion

We begin the discussion by commenting on the posterior distributions of the parameters $n_1$, $n_2$, $k_c$, and $\sigma_0$. Then, we discuss and interpret the scale $\lambda_{eq}$. Before going into details, we have to remind that the properties of the velocity field derived from our parameters correspond to properties of the solenoidal component of the velocity field since $\omega$ and consequently $\gamma$ have no information of the irrotational component of the velocity field.

## 4.1. Model Parameters

We start by discussing the behavior of the exponents $n_1$ and $n_2$. The exponent $n_1$ shows how the circulation is distributed on larges scales down to the scale $\lambda_c = 1/k_c$. Beyond $k_c$, i.e., at scales smaller than $\lambda_c$, the function $V(k)$ quickly drops, showing values of $n_2$ with no apparent upper limits. As we show in Appendix D.2, we get a similar results if $V(k) = 0$ for scales larger than $\lambda_c$. We also find in Appendix D.2 that this break is necessary to reproduce the distribution of $\gamma$ at smaller scales.

We have to point out that spatial correlations in the noncircular field are given by $V(k)$ and GRFs have coherent substructures unlike white-noise fields. The scale $\lambda_c$ could be showing the size of the coherent structures in the velocity field (Musacchio & Boffetta 2017). Coherent structures are long-lived structures that can be identified in the vorticity field (Ruppert-Felsot et al. 2005), which transport mass and energy across different scales. For scales smaller than the size of these structures, i.e., $k > k_c$, the function $V(k)$ decays quickly, meaning that there is little information about the circulation field. One interpretation is that circulation or rotation is transferred from large scales down to the scale $\lambda_c$. Below the scale $\lambda_c$ the redistribution of circulation from large scales stops and the distribution of $\gamma$ starts to converge. At these scales the solenoidal component of the velocity field starts to show
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Percentiles of circulation $\gamma$ at different annuli, as a function of scale $\ell$. Solid lines represent the percentiles of $\gamma$ in the simulation. Each color corresponds to a different percentile. Dotted lines show the percentiles of $\gamma_{\text{rot}} + \gamma_{\text{pc}}$ for the median values of $n_1$, $\lambda_c = 1/k_c$, and $\sigma_0$. Shaded regions represent $1\sigma$ uncertainty intervals from the posterior distributions of the model parameters. Vertical purple and orange lines show the spatial scales $\lambda_c = 1/k_c$ and $\lambda_{\text{eq}}$, respectively, with their corresponding uncertainty illustrated by the shaded regions. Black dashed horizontal lines represent $\gamma = 0$.</figcaption></figure>

coherent structures that are decoupled from the random behavior of the noncircular component. An observational example of this scale might be found in Rosolowsky et al. (2003), where the velocity gradients of massive clouds within regions of 500 pc are preferentially aligned. Rosolowsky et al. (2003) show an observational correlation in the velocity field for scales smaller than 500 pc, which is similar to the values we find for $\lambda_c$.

The next step is to link the structures in the velocity field with structures in the density field. One of the first scales that should affect the behavior of the velocity field is the scale height. However, the scale height of the gas density field for the three runs is of the order of 100 pc, which is usually smaller than $\lambda_c$ by a factor of 4–5. We continue analyzing different structures that exist in the plane of galactic disks.

Let us assume that the scale $\lambda_c$ is related to the fragmentation of gas with a characteristic scale $\lambda_f$. If $\Sigma_{\text{gas}}(R)$ is the gas surface density profile, we can estimate a characteristic mass for the fragments or clumps as $M_c \approx \pi\Sigma_{\text{gas}}(\lambda_f/2)^2$, and then the scale of fragmentation is approximated by

$$ \lambda_f \approx \left( \frac{M_c}{\Sigma_{\text{gas}}} \right)^{1/2} . \qquad (21) $$

The scale $\lambda_c$ could also be related to interactions between clumps (Dobbs & Pringle 2013); however, this interpretation is not straightforward since we do not find an energy cascade at smaller scales. If cloud–cloud interactions play a role, it is worth studying the typical distance between clumps of gas at each radius.

To compute the two aforementioned scales, we need to define a clump of gas. We use the following criteria to define a gas clump:

1. Each cell in the clump has a number density above $100 \text{ cm}^{-3}$.

2. Clumps are gravitationally bound.

3. Clumps have a minimum of 20 cells (a $3 \times 3 \times 3$ cube without the corners has 19 cells).

The clumps identified by these criteria have typical sizes of 100 pc, while in nearby galaxies sizes range between 10 and 100 pc (Rosolowsky 2007; Heyer et al. 2009; Colombo et al. 2014). These sizes are of the order of the scale height across the disk for the three runs. We measure $M_c$ as the total gas content in the clump. To compute the typical distance, we do the following: for each clump we average the distance to the three closest clumps, and we average that quantity at each annulus. We refer to this scale as $\lambda_{\text{int}}$.

We show $\lambda_f$ and $\lambda_{\text{int}}$ in Figure 9. For most regions, $\lambda_{\text{int}}$ and $\lambda_f$ are of the same order of magnitude, and they appear to be correlated. The spatial scale $\lambda_c$ seems to lie closer to $\lambda_f$. However, the scatter of these scales is too large to derive any strong conclusion. Each annulus has a characteristic $\Sigma_{\text{gas}}$, $\Omega$, and $\sigma_v$, quantities that are dynamically correlated. In that regard it is not surprising that spatial scales defined by these quantities show similar behaviors.

Another scale related to a change of the behavior of the velocity field is the epicyclic scale $\sigma_v/\kappa$ (Meidt et al. 2018),
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Comparison between $\lambda_c$ and characteristic spatial scales of gas clumps measured in the simulations. The orange line shows the epicyclic scale as a function of galactocentric radius. The purple line represents the length scale of fragmentation $\lambda_f = (M_c/\Sigma_{\text{gas}})^{1/2}$. The blue line shows the characteristic separation between clumps, $\lambda_{\text{int}}$. Error bars represent the 16th–84th percentile interval for each parameter.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Spatial scales as a function of galactocentric radius: solid pink and yellow lines correspond to $\lambda_{\text{eq}}$ and $\lambda_c$, respectively. The shaded regions correspond to $1\sigma$ uncertainties. The classical instability scales for two-dimensional disks, $\lambda_{\text{rot}}$ and $\lambda_J$, are shown as solid light-blue lines and black dashed lines, respectively. The dotted-dashed line shows the effective values $\lambda_{\text{rot}}$ after adding the effects of resolution.</figcaption></figure>

where $\sigma_v$ is the velocity dispersion of a gas cloud and $\kappa$ is the epicyclic frequency. Epicyclic motions correspond to the evolution to small perturbations of circular orbits under the gravitational potential of the galaxy. Structures larger than their corresponding epicyclic scale are ensured to be affected by the galactic potential. Figure 9 shows $\sigma_z/\kappa$ as the orange solid line, where $\sigma_z$ corresponds to the dispersion velocity in the z-axis in a radial bin. The choice of $\sigma_v = \sigma_z$ assumes that once a structure has formed, its velocity dispersion is nearly isotropic. Like the scales $\lambda_f$ and $\lambda_{\text{int}}$, the epicyclic scale lies close to $\lambda_c$.

The physical correlations between all these scales and their level of uncertainty make it difficult to compare them with $\lambda_c$. Therefore, we are not able to elucidate the fundamental physical origin of $\lambda_c$. We can only conclude that $\lambda_c$ is related to the formation of structure in our simulations and that the details of gas dynamics below such structures do not significantly affect the overall circulation of gas.

## 4.2. Circulation Scales: $\lambda_c$ and $\lambda_{\text{eq}}$

We start this section by discussing the role that gravitational instabilities can play in setting the distribution of circulation, particularly their effect in $\lambda_c$ and $\lambda_{\text{eq}}$. First, we recall the Toomre parameter, $Q = \kappa\sigma_v/\pi G\Sigma$, which for marginal stable

systems ($Q \approx 1$) and with constant $\kappa$ requires $\sigma_v \propto \Sigma$. At the scales of clouds, for virial parameters $\alpha \approx 2$, $\sigma_v$ also grows with $\Sigma$ ($\sigma_v \propto \sqrt{\Sigma_{\text{gas}}}$; Sun et al. 2018). For any of these two pictures, we expect that galaxies with more gas have more randomness in their velocity fields. This is illustrated by the run G2E1 in Figure 7, which shows higher values of $\sigma_0$ and $\lambda_{\text{eq}}$ compared to the other runs.

Since we are dealing with galactic disks, the first step to visualize relations with gravitational instabilities is to compute $\lambda_{\text{rot}}$ and the two-dimensional thermal Jeans scale $\lambda_J$, given by $c_s^2/G\Sigma$, where $c_s$ is the gas sound speed and $\Sigma$ is the gas surface density. The thermal Jeans scale sets the size of the smallest structures that can be formed. Both length scales, $\lambda_{\text{rot}}$ and $\lambda_J$, are shown in Figure 10.

Since the two-dimensional stability is affected by the disk thickness (or the resolution in the case of simulations), we have to consider the dispersion relation $\omega_p^2(k) = \kappa^2 - 2\pi G\Sigma|k|e^{-k\epsilon}$, where $\omega_p$ is the frequency of perturbations and $H$ is the disk thickness with a minimum value set by the numerical resolution (Binney & Tremaine 2008, p. 552). Perturbations where $\omega_p^2(k) < 0$ correspond to instabilities. To obtain the correct values if $\lambda_{\text{rot}}$, which we call $\lambda_{\text{rot,c}}$, we solve the equation $\omega_p^2(k) = 0$. We plot $\lambda_{\text{rot,c}}$ in Figure 10 as the dotted-dashed blue line. Once the disk thickness is considered, near the galactic
centers $\omega_p$ is always real and any radial perturbation is stable. This means that in the inner regions of these galaxies gravitational instabilities are not resolved. At these scales we do not expect to see an important injection of energy due to gravitational instabilities. The effect of the disk thickness is more noticeable in our additional simulations in Figure E2: if $\lambda_{rot,c}$ is not resolved, $\lambda_{eq}$ falls below the resolution of the simulations.

For this set of simulations, $\lambda_c$ lies between $\lambda_j$ and $\lambda_{rot,c}$, that is, within the range of scales in which gravitational instabilities can exist, consistent with this scale being associated with scales of structure formation. For all runs we see that $\lambda_{eq}$ is above $\lambda_j$, increases with radius, and in some regions is higher than $\lambda_{rot,c}$. Since $\lambda_{rot,c}$ is the maximum size of unstable perturbations, clouds formed in regions where $\lambda_{eq} > \lambda_{rot,c}$ will be predominantly dominated by noncircular motions.

For all runs we see that $\lambda_{eq}$ can show values up to kiloparsec scales. For G2E1, $\lambda_{eq}$ is higher than $\lambda_{rot}$ in most regions. This is probably caused by the high star formation rate, due to its higher gas content and associated increase in feedback-induced noncircular motions. In Figure E2 in Appendix E, our simulations with only SN feedback show values of $\lambda_{eq}$ lower than $\lambda_{rot}$ and lower than $\lambda_c$. This last point shows that only considering the difference in gas content between simulations without taking into account stellar feedback is insufficient to explain differences in $\lambda_{eq}$ across the different runs. These additional simulations also show an apparent correlation between $\lambda_{eq}$ and $\lambda_{rot}$. However, once we use a more energetic type of feedback, $\lambda_{eq}$ grows, and this apparent correlation disappears. This suggests that there might be two different regimes where the distribution of circulation is set by gravitational instabilities or by stellar feedback. This idea goes in line with results from numerical simulations and analytical models showing that turbulence can be powered by gravity or stellar feedback, and that the dominant driver of turbulence changes across the evolution of the universe (Krumholz & Burkert 2010; Goldbaum et al. 2015, 2016; Krumholz et al. 2018).

Since we are discussing the turbulent behavior of gas and its dynamical stability, we can also discuss the relevance of the turbulent Jeans scale $\lambda_{turb} = \sigma_v^2/G\Sigma$, where $\sigma_v$ is the velocity dispersion of gas considering nonthermal motions. However, $\sigma_v$ is a function of scale $\ell$ (Elmegreen & Scalo 2004; Romeo et al. 2010), and to properly take into account the effects of turbulence, we need to know how $\sigma_v$ changes with $\ell$. Here we take a first-order approach considering $\sigma_v = (c_s^2 + \sigma_z^2)^{1/2}$, where $\sigma_z$ is the mass-weighted vertical dispersion velocity in the disk of the galaxy. In this approximation we are assuming that the velocity dispersion at the scale of the disk scale height is a representative value of the velocity dispersion for bound structures in the presence of turbulence. To compute $\lambda_{turb}$, we use a temperature cut of $T < 5000$ K to avoid considering gas that is currently affected by stellar feedback. We show $\lambda_{turb}$ in Figure 10. For runs G2E1 and G1E1 $\lambda_{turb}$ is of the order of $\lambda_c$, which suggests that $\lambda_c$ could be tracing the scales at which turbulent structures are affected by their self-gravity. This is not shown by run G1E0.5, but we need to keep in mind that we are assuming that $\sigma_z^2$ is a good proxy for the turbulent velocity in the plane of the galaxy for self-gravitating structures.

In the field of fluid dynamics, it is known that in turbulent fluids coherent structures naturally arise, and that these structures are fundamental for the transport of angular momentum across different scales (Kraichnan 1967; Ruppert-Felsot et al. 2005).

This motivates us to look for structures that can be defined by kinematics only. One alternative is to look for structures whose behavior is defined by the ratio between the galactic angular velocity, which traces the galactic potential, and the local noncircular motions at cloud scales. A spatial scale that goes in that direction and compares the magnitudes of the velocity dispersion of gas and the galactic potential is the epicyclic scale $\sigma_v/\kappa$ (Meidt et al. 2018). However, as shown in Figure 9, if we use the velocity dispersion across the z-axis, the epicyclic scale is of the order of $\lambda_c$, which in most regions is smaller than $\lambda_{eq}$. Another approach is to consider the velocity dispersion in the x-y plane, $\sigma_{xy}$, within a radial bin. To compute $\sigma_{xy}$, we subtract the circular velocity model from the velocity field. The scale $\sigma_{xy}/\kappa$ compares the energy in the noncircular velocity field with respect to epicyclic motions given by the galactic potential. We show the spatial scale $\sigma_{xy}/\kappa$ in Figure 10 as the orange line.

Figure 10 shows that $\sigma_{xy}/\kappa$ is similar to $\lambda_{eq}$. This result would suggest that the ratio $\sigma_{xy}/\kappa$ is a good proxy for $\lambda_{eq}$. However, this might be valid only for our feedback prescription and for galaxies with an average star formation rate according to the Kennicutt–Schmidt relation (Kennicutt 1998; Daddi et al. 2010). For our second set of simulations described in Appendix E with only SN feedback, the values of $\lambda_{eq}$ are usually lower than $\sigma_{xy}/\kappa$ and $\lambda_c$. This again shows the effect of using different feedback prescriptions. Momentum feedback prescriptions change the velocity field more aggressively; $\lambda_{eq}$ shows large values and is similar to $\sigma_{xy}/\kappa$.

We can interpret this difference as differences in the sources of noncircular motions. Simulations with only thermal SN feedback can lose this source of energy quickly owing to our resolution of 30 pc, producing a lower effect in the velocity field. In this scenario, gravitational instabilities might become a relevant source of noncircular motions or turbulence. On the other hand, the stellar feedback prescription used in our main simulations changes explicitly the velocity field at the smallest scales and increases directly the magnitude of $V(k)$. Feedback might erase the correlation between gravitational instabilities and the noncircular motions at small scales.

This point has implications for the analysis of star formation in numerical simulations. In simulations with mechanical stellar feedback, i.e., momentum injection, the velocity field is explicitly changed and the amount of kinetic energy at small scales increases. This reduces the coupling between the dynamics of clouds and galactic rotation, as well as the coupling between the efficiency of star formation and the galactic environment. The magnitude of correlations between star formation and galactic properties found in simulations might depend on the specific stellar feedback prescription used.

Up to this point we have shown that the distribution of circulation is affected by stellar feedback and gravitational instabilities. It is interesting to discuss what other studies show with respect to the distribution of circulation or rotation. Here we mention the works of Tasker & Tan (2009) and Ward et al. (2016), which use simulations to measure how the rotation of molecular clouds aligns with respect to the rotation of their galaxies. These simulations have similar surface gas densities and the same shape of the velocity curve, but with different magnitudes. Both simulations have weak forms of feedback; Tasker & Tan (2009) did not include stellar feedback in their simulations, and the simulations in Ward et al. (2016) had a reduced feedback efficiency of 10%. Hence, we expect that their distribution of circulation is set by gravitational instabilities. For
comparison, at a galactocentric radius of 8 kpc, $\lambda_{\text{rot}} \simeq 2100$ pc in Tasker & Tan (2009) and $\lambda_{\text{rot}} \simeq 1500$ pc in Ward et al. (2016). In addition, the simulation of Ward et al. (2016) shows spiral structures, while the density field in Tasker & Tan (2009) is more random. The simulation of Tasker & Tan (2009) is more unstable than the one from Ward et al. (2016), and consequently the former should have higher values of $\lambda_{\text{eq}}$, or a higher fraction of molecular clouds with retrograde rotation with respect to their galaxy. These simulations effectively find different fractions of retrograde clouds, 30% in Tasker & Tan (2009) and 13% in Ward et al. (2016). This shows that more unstable systems are more dominated by noncircular motions and have higher values of $\lambda_{\text{eq}}$. Tasker & Tan (2009) also analyzed the effects of resolution, which directly influences the size of molecular clouds and the stability of gas dynamics as shown by $\lambda_{\text{rot,c}}$. Tasker & Tan (2009) show that as the resolution is increased, more molecular clouds present retrograde rotation. In summary, these studies show that in the absence of strong feedback, gravitational instabilities play a role in setting how circulation is distributed at smaller scales and the relevance of the spatial resolution used in numerical simulations.

### 4.3. Turbulence and the Power Spectrum

Since we are studying two-dimensional velocity fields, we discuss how our results are compared with known properties of two-dimensional turbulence. A turbulent velocity field $\delta v$ is characterized by its kinetic energy spectrum $E(k)$ such that the mean turbulent kinetic energy per unit mass is $\frac{1}{2}\langle\delta v^2\rangle = \int_0^\infty E(k)dk$. The power spectrum $E(k)$ of two-dimensional turbulence is characterized by the existence of two inertial regimes: (i) an inverse energy cascade $E(k) \propto k^{-5/3}$ for $k < k_f$ and (ii) a direct enstrophy cascade with $E(k) \propto k^{-3}$ for $k > k_f$, where $k_f$ is the wavenumber of the forcing scale (Kraichnan 1967; Wada et al. 2002; Bournaud et al. 2010; Musacchio & Boffetta 2017). In our model, the random component of the velocity $v_{nc}$ field is characterized by $V(k)$. From Appendix B it can be shown that $E(k) \propto V(k)^2k$ for a two-dimensional field, which is the case of $v_{nc}$. For two- and three-dimensional turbulence, the relation between $n$ and the energy spectrum $E(k)$ is given by

$$ (2D) \quad E(k) \propto v(k)^2 k \propto k^{-2n+1}. \quad (22) $$

In this formalism, the inverse energy and the direct enstrophy cascade are represented by $n = 4/3 \approx 1.33$ and $n = 2$. As shown in Figure 7, the distribution for the exponent $n_1$ ranges between 1.1 and 1.8. This exponent lies close to the expected values of both regimes. On the other hand, $n_2$, which should be associated with the enstrophy cascade, has no upper limit in our model, and the scale where $V(k)$ breaks, $\lambda_c$, is of the order of 240 pc to 1 kpc. If we look at the middle panel of Figure C1 in Appendix C, we see that for high values of $n$ the distribution of $\gamma$ starts to be less sensitive to changes in $n$. It is likely that most of the information of the distribution of $\gamma$ is given by $n_1$ and $k_c$. If this is the case, $k_c$ is related to the turbulent forcing scale $k_f$. Experiments of thin layer fluids show strong long-lived vortices at the scale $\lambda_f = 1/k_f$ (Musacchio & Boffetta 2017). Then, the turbulent picture also suggests that $\lambda_c$ might be related to the formation of structure in the turbulent velocity field.

Numerical and observational studies suggest that information about the turbulent velocity field can be extracted from the spectrum of the gas surface density (Elmegreen et al. 2001;

Bournaud et al. 2010; Combes et al. 2012). The numerical work of Bournaud et al. (2010) shows that the spectrum of the surface density field $\Sigma(k)$ may be described by a broken power law with a critical scale that is interpreted as the disk scale height. The slope of the power law in Bournaud et al. (2010) simulations changes from -2 at large scales to -3 small scales. Measuring the density spectrum $\Sigma(k)$ for our simulations, we find that the slope changes from around zero for scales larger than 500 pc (about the same order of $\lambda_c$) to a continuous decaying function at small scales with no clear critical scale. In our simulations we are not able to link the density and the power spectrum of gas.

## 5. Limitations and Caveats

The work presented here is largely exploratory and aimed at establishing the basic concepts associated with modeling the spatially dependent distribution of gas circulation in disk galaxies. Here we discuss some of the limitations associated with this modeling. In a future work we expect to address several of these limitations in order to apply the presented methods to extract information from observed galaxy velocity fields.

### 5.1. Velocity Model

The main assumption in our model is that the velocity field in the plane of the disk can be approximated as the contributions of two different fields: $V = V_{\text{rot}} + V_{\text{nc}}$, where $V_{\text{rot}}$ corresponds to the galactic velocity curve and $V_{\text{nc}}$ corresponds to a Gaussian random velocity field. In real galaxies, we find other types of coherent motions that are different from galactic rotation and pure random motions. Among these, we find induced motions by galactic bars and spirals, and epicyclic motions. At the scale of epicycles, gas is still affected by the tidal forces exerted by the galactic potential (Meidt et al. 2018). According to Meidt et al. (2018), depending on the strength of self-gravity, the dynamical structure of clouds shows preferred orientations in radial or azimuthal coordinates, which does not occur in our model of $V_{\text{nc}}$. In addition, galactic bars and spirals would also produce deviations from global galactic rotation, which we are implicitly including in $V_{\text{nc}}$.

Figure 8 shows that our simple model can successfully fit the distribution of circulation across different spatial scales in general terms, but there are some clear deviations in particular regimes (e.g., small-scale, prograde rotating regions with high values of $\gamma$ at intermediate galactocentric radii), which probably signal more complex types of motions not recovered by the model. Moreover, due to conservation of angular momentum, the vorticity in high-density regions is enhanced. It is unclear how to statistically model these types of motions.

### 5.2. Full Velocity Field

In this work we have made use of isolated galaxy simulations whose rotation axis is aligned with the z-axis of the simulation box by default. To compute the circulation, we used the two-dimensional velocity field $V(x, y)$, the density-weighted projection across the z-axis of the three-dimensional velocity field. We can separate the total circulation into two
terms, $\Gamma_x$ and $\Gamma_y$:

$$ \Gamma = \oint \mathbf{v} \cdot d\mathbf{r} = \oint v_x \cdot dx + v_y \cdot dy = \Gamma_x + \Gamma_y, \quad (23) $$

where x and y are coordinates on the plane of the disk. In observations we only have access to the velocity along the line of sight, which will be the sum of one of the velocities in the plane of the galaxy, $v_x$ or $v_y$ and $v_z$, motions vertical to the disk midplane. Consider a disk with inclination i such that the line of sight lies in the x-z plane. The coordinates in the plane of the sky are $z' = z \sin i + x \cos i$ and $y' = y$. The axis of the line of sight is $x' = x \sin i + z \cos i$, and the velocity is $v_{LOS} = v'_x = v_x \sin i + v_z \cos i$. On the midplane $z=0$ and the projected position $z' = x \cos i$. From the observed quantities we can compute

$$ \Gamma' = \oint v_{LOS} dz' = \sin i \cos i \oint v_x dx + \cos^2 i \oint v_z dx \\ \Gamma' = \sin i \cos i \Gamma_x + \cos^2 i K_{xz}, \quad (24) $$

where $\Gamma_x$ is the component of $\Gamma$ in the x-axis and $K_{xz}$ is the sum of vertical motions along the x-axis. The term $K_{xz}$ should be of the order of $\langle v_z \rangle l$. This component has to be treated as an additional term in the assumed decomposition of the velocity field. In this work we have assumed that the velocity field in the plane is the sum of galactic rotation plus a random component. For $v_x$ this means $v_x = -\Omega y + \delta v_x$, where $\Omega$ is the galactic angular velocity and $\delta v_x$ is the random velocity term. For $v_z$ we can assume that $v_z = \delta v_z$. Then, $\Gamma' \approx -\cos i [y \sin i \oint \Omega dx + l(\langle v_x \rangle \sin i - \langle v_z \rangle \cos i)]$. At large scales $\langle v_x \rangle$ and $\langle v_z \rangle$ are approximately zero, while at small scales both terms behave like random variables. The sum of these two terms would be the observed random component. Since we want to compare them with galactic rotation, the best inclination has to maximize the contribution of $\Omega$ to $\Gamma'$ that occurs at $i = 45^\circ$.

### 5.3. Surface Brightness Limits and Recovery of Velocity Information

A major limitation comes from the observational detection limits for different transition lines, which lead to an incomplete sampling of the velocity field. For example, the CO (1–0) transition has a critical density $\simeq 10^3$ cm$^{-3}$, tracing the distribution of molecular gas in galaxies. This implies that we can only observe a small fraction of the velocity field at the scales of molecular clouds. Proper ways of dealing with noise and censored data in faint regions will also need to be implemented.

### 5.4. Simulations

In this paper we use hydrodynamical simulations of disk galaxies to test our method. The results presented here are valid to our set of simulations, with their defined prescriptions for star formation and stellar feedback. However, caution should be taken before directly extrapolating our findings to the environments of real galaxies. Here we list what we consider are the most important aspects in which our simulations and observed galaxies differ:

1. Resolution: The maximum spatial resolution corresponds to 30 pc. In practice, this means that we are able to

resolve structures and instabilities of the order of 100 pc, corresponding to approximately four times our resolution. In nearby galaxies, the size of molecular clouds typically ranges from tens to hundreds of parsecs. Although we see formation of structures, this resolution is not enough to resolve the inner turbulence of molecular clouds, their gravitational collapse, and the interactions of clouds smaller than 100 pc. A higher resolution would imply more interactions and a higher velocity dispersion at the smallest scales studied here. Then, we might expect a change in the values of $n_2$ that sets the behavior of $V(k)$ at large wavenumber $k$, i.e., at smaller spatial scales. Despite this caveat, $\lambda_{eq}$ and $\lambda_c$ are well resolved almost everywhere.

2. Temperature: Gas is allowed to cool owing to radiation down to a temperature of 300 K. This means that the smallest structures in our simulations are more similar to H I clouds. Also, this temperature floor sets a minimum Jeans scale as a function of gas surface density

$$ \lambda_{J,\min} = \frac{5 k_B T_{\min}}{3 \mu m_p} \frac{1}{G \Sigma_{\text{gas}}} = 133 \text{pc} \left( \frac{\Sigma_{\text{gas}}}{10 M_{\odot} \text{ pc}^{-2}} \right)^{-1} \quad (25) $$

assuming a mean molecular weight $\mu = 0.6$. In local galaxies, the Jeans length is of the order of a few parsecs, about two orders of magnitude below the average Jeans scales found in our simulations. The temperature floor leads to an overestimation of the relevance of the Jeans scale. It is important to mention that to compute the radial profiles of $\lambda_J$ shown here, we are considering all the gas in an annulus and its respective average temperature instead of the average $\lambda_J$ for cold and dense gas. This makes sense for our analysis since we are computing the circulation for all the gas within the physical volume described in the paper. However, for regions with densities below $10 M_{\odot} \text{ pc}^{-2}$ the minimum value for $\lambda_J$ is larger than four resolution elements. Then, even for our resolution the values of $\lambda_J$ are likely overestimated and should be considered as upper limits.

3. Stellar feedback: In our recipe of stellar feedback, we include the direct injection of momentum from radiation pressure and stellar winds to the six nearest cells. Although the amount of added momentum $p_*$ does not explicitly change with spatial resolution, the typical masses, $m_{cell}$, of cells around star particles do change with different spatial resolution. This translates in different magnitudes for the change of the velocity field around star particles, since the velocity $v_* = p_*/m_{cell}$. We have not tested how sensitive to resolution is this feedback prescription.

4. Spiral arms and bars: A relevant difference between observations, other simulations, and our runs is that our simulated galaxies lack grand-design spiral arms. The main difference is that the old stellar population in this work is represented by an external axisymmetric potential, whereas other studies use particles (Renaud et al. 2013) or spiral potentials (Dobbs et al. 2015). Only new stars are particles; hence, only this stellar component can respond to perturbations making the stellar disk more stable. The impact of spiral arms and bars in our analysis can be separated by their effect on large and small scales. At large scales, the bulk motion $v_{rot}$ must be a function of radius and the azimuthal angle. At small scales, spiral and
bars can induce vortex motions by Kelvin–Helmholtz instabilities, Rayleigh–Taylor instabilities, or tidal fields (Dobbs & Bonnell 2006; Renaud et al. 2013). These structures create new sources of turbulence; therefore, the velocity field at the scale of molecular clouds has different properties. In this work, we have divided disk galaxies into annular regions and measured radial profiles of the parameters that define the small-scale velocity field $v_{nc}$. To test the effect of spirals and arms, we also need to separate regions according to their azimuthal distance to these structures.

## 6. Summary and Conclusions

In this study we characterize the rotation of gas in galaxies at different scales by measuring the circulation $\Gamma$, a macroscopic measure of fluid rotation. We develop a method to measure the contributions of large-scale motions, i.e., galactic rotation, and noncircular motions in the observed distribution of circulation at different spatial scales. Noncircular motions are modeled as random Gaussian velocities described by a generating function in Fourier space $V(k)$. We apply this method on three hydrodynamical simulations of galactic disks, performed with the AMR code Enzo, which includes star formation, SN feedback, and momentum feedback from stellar winds, and with a spatial resolution of 30 pc.

We summarize the major points of this work:

1. We model the velocity field of galaxies with two components: a galactic component given by the circular velocity profile, and a Gaussian random component. The random component is obtained from a function $V(k)$, whose functional form corresponds to a broken power law with exponents $n_1$ and $n_2$, transitioning at the wavenumber $k = k_c$. The amplitude of $V(k)$ is defined by the characteristic velocity dispersion of the random field $\sigma_0$. We apply the model to hydrodynamical simulations and confirm that motions can be well modeled by two components with different circulation, as hypothesized. The model successfully reproduces the distribution of circulation as a function of scale, except when regions are under gravitational collapse.

2. We find that a sharp transition in the behavior of gas dynamics at the scale $\lambda_c = 1/k_c$ is necessary to fit the circulation distribution. This may correspond to the scale at which kinematics transition from being coupled to the galaxy to more disordered motion, associated with feedback-driven turbulence or gravity-driven turbulence. However, the resolution of the current simulations limits our ability to probe this in greater detail.

3. The scale $\lambda_c$ is similar to the scale at which gas fragments and to the epicyclic scale $\sigma_z/\kappa$ that defines the scale at which self-gravity and the potential of the galaxy are equally important to determine the internal dynamics of clouds. The scale $\lambda_c$ is also similar to the scale of fragmentation and the distance between clumps, suggesting that $\lambda_c$ shows the formation of structure in the density field.

4. We introduce a dynamical spatial scale $\lambda_{eq}$. At spatial scales similar to $\lambda_{eq}$ the contributions of galactic circular motions and noncircular motions to the observed circulation or local rotation of gas are roughly the same. For regions larger than $\lambda_{eq}$ galactic rotation dominates the

circulation of gas and consequently the measured rotation. At these scales the distribution of circulation shows largely positive values, which means that gas rotates in the same orientation of the galaxy. For patches of gas smaller than $\lambda_{eq}$, noncircular and random motions start to dominate the observed circulation and retrograde rotating regions can be found.

5. We find that $\lambda_{eq}$ depends on the local properties of gas. From the center of the galaxies, $\lambda_{eq}$ increases with galactocentric radius. This shows that the spatial scale at which gas dynamics is dominated by noncircular motions depends on the position in the galactic disk. We see different behaviors in the central regions and outskirts of galaxies. Galactic rotation appears to be more important or dominant toward the center of galactic disks.

6. The scale $\lambda_{eq}$ is similar to the ratio $\sigma_{xy}/\kappa$, as predicted by models about balance of rotation and turbulence, where $\sigma_{xy}$ is the in-plane velocity dispersion of the two-dimensional velocity field within a radial annulus and $\kappa$ is the epicyclic frequency. However, when suppressing momentum feedback in simulations, $\lambda_{eq}$ can be lower than $\sigma_{xy}/\kappa$, and its radial profile might be correlated with the spatial scales of gravitational instabilities since self-gravity becomes a relevant source of turbulence.

7. In some regions $\lambda_{eq}$ is greater than $\lambda_{rot}$. The formation of structures in such regions will be dominated by noncircular motions. Depending on the sources of feedback, turbulence, and local dynamics, $\lambda_{eq}$ can be larger or smaller than $\lambda_{rot}$. Particularly, for strong modes of stellar feedback $\lambda_{eq}$ can be larger than $\lambda_{rot}$.

8. Stellar feedback changes $\lambda_{eq}$ by injecting momentum at smaller scales, increasing motions that do not follow galactic rotation. Different prescriptions of stellar feedback will produce different velocity fields that show different coupling between the dynamics of gas at small scales and the large-scale galactic rotation. This can also produce changes in the coupling between star formation and galactic rotation.

This work shows that the characterization of the ISM circulation, from the modeling of velocity fields of galaxies, opens the possibility of directly measuring scales associated with gravitational collapse and structure formation and studying how it changes with galaxy properties and local conditions. It also shows that rotation is dynamically important in some environments like the centers of galactic disks. In the future, the analysis of circulation can be extended to real galaxies to study the relevance of galactic rotation at the scale of molecular clouds.

Some of the computing for this project was performed on the Memex cluster. We would like to thank Carnegie Institution for Science and the Carnegie Sci-Comp Committee for providing computational resources and support that contributed to these research results. Powered@NLHPC: This research was partially supported by the supercomputing infrastructure of the NLHPC (ECM-02). The Geryon/Geryon2 cluster housed at the Centro de Astro-Ingenieria UC was used for part of the calculations performed in this paper. The BASAL PFB-06 CATA, Anillo ACT-86, FONDEQUIP AIC-57, and QUIMAL 130008 provided funding for several improvements to the Geryon2 cluster. J.U. acknowledges support from Programa Nacional de Becas de Postgrado, CONICYT (grant
D-21140839). A.E. acknowledges support from CONICYT project Basal AFB-170002 and Proyecto FONDECYT Regular grant 1181663. F.B. acknowledges funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement No. 726384). J.M.D.K. gratefully acknowledges funding from the German Research Foundation (DFG) in the form of an Emmy Noether Research Group (grant No. KR4801/1-1) and from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program via the ERC Starting Grant MUSTANG (grant agreement No. 714907). A.H. is a fellow of the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD).

Software: yt (Turk et al. 2011), Enzo (Bryan et al. 2014).

## Appendix A
### Vorticity and Angular Momentum

Let us approximate the velocity field around a point $\mathbf{r}_0 = x_0\hat{x} + y_0\hat{y}$ in Cartesian coordinates. We choose to do this analysis in Cartesian coordinates since cylindrical or spherical coordinate systems are noninertial reference frames. The unit vector $\hat{x}$ is instantly aligned with the radial unit vector $\hat{R}$, i.e., $\hat{x} = \hat{R}$ and $\hat{y} = \hat{\phi}$, where $\hat{\phi}$ is the azimuthal cylindrical coordinate. Then, $x_0 = R_0$ and $y_0 = 0$. For a velocity field $\mathbf{v} = R\Omega\hat{\phi}$,

$$v_x = -y_0\Omega(R_0) = 0 \qquad (A1)$$

$$v_y = x_0\Omega(R_0) = R_0\Omega(R_0). \qquad (A2)$$

We expand $\Omega(R)$ in terms of $x$ and $y$, around the point $(R_0 + x, 0 + y)$. Since the Cartesian coordinates are aligned with the cylindrical coordinates $R$ and $\phi$, we can use the gradient in cylindrical coordinates:

$$\Omega(R_0 + x, y) = \Omega_0 + \left. \frac{\partial \Omega}{\partial R} \right|_{R=R_0} x + \left. \frac{\partial \Omega}{R \partial \phi} \right|_{R=R_0} y \quad (A3)$$

$$= \Omega_0 + \left. \frac{\partial \Omega}{\partial R} \right|_{R=R_0} x. \qquad (A4)$$

We expand the velocity at first order in $x$ and $y$,

$$v_x = -(y_0 + y)\Omega \approx -y\Omega_0 \qquad (A5)$$

$$v_y = (R_0 + x)\Omega \approx R_0\Omega_0 + \left( \Omega_0 + R_0 \left. \frac{\partial \Omega}{\partial R} \right|_{R_0} \right) x. \qquad (A6)$$

Using this, we can compute the local angular momentum in the z-axis, $L_z = \hat{z}L$, where $L$ is the angular momentum vector and $\hat{z}$ is the unit vector in the z-axis:

$$L_z = \hat{z} \int \rho \mathbf{r} \times \mathbf{v} \, d^3\mathbf{x} \qquad (A7)$$

$$= \int \left[ R\Omega_0 x + x^2 \left( \Omega_0 + R \left. \frac{\partial \Omega}{\partial R} \right|_{R_0} \right) + y^2 \Omega_0 \right] \rho d^3\mathbf{x} \qquad (A8)$$

$$= MR\Omega_0 x_{CM} + I_{11}\Omega_0 \left( 1 + \left. \frac{\partial \ln \Omega}{\partial \ln R} \right|_{R_0} \right) + I_{22}\Omega_0,$$

$$I_{11} = \int \rho x^2 d^3 \mathbf{x}, \quad I_{22} = \int \rho y^2 d^3 \mathbf{x}, \qquad (A9)$$

where $x_{CM}$ is the position of the center of mass in the x-axis and $I_{11}$ and $I_{22}$ are the components of the moment of inertia tensor. If the coordinates $x$ and $y$ are centered on the center of mass and $I_{11} = I_{22} = \frac{1}{2}I$, the angular momentum is reduced to

$$L_z = \frac{1}{2} I \Omega_0 \left( 1 + \left. \frac{\partial \ln \Omega}{\partial \ln R} \right|_{R_0} \right) = \frac{1}{2} I \omega. \qquad (A10)$$

## Appendix B
### Random Fields: Velocity, Vorticity, and Normalized Circulation

In this section we show how to create GRFs for the noncircular velocity component $v_{nc}$, using a function $\mathcal{V}(k)$ in Fourier space. Let $\mathbb{W}(\mathbf{r})$ be a Gaussian white-noise field with $\mu = 0$ and $\sigma = 1$. From $\mathbb{W}(\mathbf{r})$ we create a Gaussian random velocity field $v(\mathbf{r})$ by

$$\begin{aligned} v(\mathbf{r}) &= \mathcal{F}^{-1}(\mathcal{V}(\mathbf{k})\mathcal{F}(\mathbb{W}(\mathbf{r}))) \\ &= \int e^{2\pi i \mathbf{k} \cdot \mathbf{r}} \mathcal{V}(\mathbf{k}) \, d\mathbf{k} \int e^{-2\pi i \mathbf{k} \cdot \mathbf{r}'} \mathbb{W}(\mathbf{r}') \, d\mathbf{r}' \end{aligned} \qquad (B1)$$

where $\mathcal{V}(k)$ is an even and positive function in Fourier space that generates the field $v(\mathbf{r})$, and $\mathcal{F}$ is the Fourier transform (Lang & Potthoff 2011). The velocity dispersion of this field is given by

$$\sigma_v^2 = \int_{k_{\min}}^{k_{\max}} \mathcal{V}(\mathbf{k})^2 d\mathbf{k} = 2\pi \int_{k_{\min}}^{k_{\max}} \mathcal{V}(k)^2 k dk. \qquad (B2)$$

This velocity field generates a vorticity field $\omega = \nabla \times \mathbf{v}$. The Fourier transform $\mathcal{F}$ of a field $f(\mathbf{r})$ is given by $\hat{f}(\mathbf{k}) = \mathcal{F}[f(\mathbf{r})] = \int e^{-2\pi i \mathbf{k} \cdot \mathbf{r}} f(\mathbf{r}) d\mathbf{r}$. Then, in Fourier space $\hat{\omega}(\mathbf{k}) = -2\pi i \mathbf{k} \times \hat{\mathbf{v}}(\mathbf{k})$. More directly, $\omega$ will be generated by the function $\mathbb{W}(k) = 2\pi\mathcal{V}(k)k$. At the smallest scales the dispersion of both fields $\omega$ and $v$ will be related through

$$\frac{\sigma_{\omega}^2}{\sigma_v^2} = \frac{\int_{k_{\min}}^{k_{\max}} \mathcal{W}(k)^2 dk}{\int_{k_{\min}}^{k_{\max}} \mathcal{V}(k)^2 k dk} = 4\pi^2 \frac{\int_{k_{\min}}^{k_{\max}} \mathcal{V}(k)^2 k^3 dk}{\int_{k_{\min}}^{k_{\max}} \mathcal{V}(k)^2 k dk}. \qquad (B3)$$

Now that we have a way to relate $\omega(x, y)$ with the function $\mathcal{V}(k)$ that generates the random velocity field, let us compute the normalized circulation $\gamma(x, y, l)$,

$$\begin{aligned} \omega(\mathbf{r}) &= \mathcal{F}^{-1}(\mathcal{W}(\mathbf{k})\mathcal{F}(\mathbb{W}(\mathbf{r}))) \\ &= \int e^{2\pi i \mathbf{k} \cdot \mathbf{r}} \mathcal{W}(\mathbf{k}) e^{-2\pi i \mathbf{k} \cdot \mathbf{r}'} \mathbb{W}(\mathbf{r}') \, dk \, d\mathbf{r}' \end{aligned} \qquad (B4)$$

$$\begin{aligned} \gamma(\mathbf{r}, l) &= \frac{1}{l^2} \int_{x-\frac{l}{2}}^{x+\frac{l}{2}} \int_{y-\frac{l}{2}}^{y+\frac{l}{2}} \omega(\mathbf{r}'') d\mathbf{r}'' \\ &= \frac{1}{l^2} \int_{x-\frac{l}{2}}^{x+\frac{l}{2}} \int_{y-\frac{l}{2}}^{y+\frac{l}{2}} d\mathbf{r}'' \int e^{2\pi i \mathbf{k} \cdot \mathbf{r}''} \mathcal{W}(\mathbf{k}) e^{-2\pi i \mathbf{k} \cdot \mathbf{r}''} \mathbb{W}(\mathbf{r}'') d\mathbf{k} d\mathbf{r}'' \end{aligned} \qquad (B5)$$

Changing $\mathbf{r}'' \rightarrow \mathbf{r} + \mathbf{r}''$, we get

$$\begin{aligned} \gamma(\mathbf{r}, l) &= \frac{1}{l^2} \int_{-\frac{l}{2}}^{\frac{l}{2}} \int_{-\frac{l}{2}}^{\frac{l}{2}} d\mathbf{r}'' \\ & \times \int e^{2\pi i \mathbf{k} \cdot \mathbf{r}} e^{2\pi i \mathbf{k} \cdot \mathbf{r}''} \mathcal{W}(\mathbf{k}) e^{-2\pi i \mathbf{k} \cdot \mathbf{r}''} \mathbb{W}(\mathbf{r}'') d\mathbf{k} d\mathbf{r}'' \end{aligned} \qquad (B6)$$
We can integrate the term $e^{2\pi i k \cdot r''} dr''$ over the rectangular square

$$ \int_{-\frac{\ell}{2}}^{\frac{\ell}{2}} \int_{-\frac{\ell}{2}}^{\frac{\ell}{2}} e^{2\pi i k \cdot r''} dr'' = \int_{-\frac{\ell}{2}}^{\frac{\ell}{2}} e^{2\pi i k_x x''} dx'' \\ \times \int_{-\frac{\ell}{2}}^{\frac{\ell}{2}} e^{2\pi i k_y y''} dy'' = \frac{\sin(\pi k_x \ell) \sin(\pi k_y \ell)}{\pi k_x \pi k_y}. \quad (B7) $$

Now $\gamma(\mathbf{r}, \ell)$ has the form $F^{-1}(M(\mathbf{k})F(\mathbf{W}(\mathbf{r})))$,

$$ M(\mathbf{k}) = \frac{1}{\ell^2} W(\mathbf{k}) \frac{\sin(\pi k_x \ell) \sin(\pi k_y \ell)}{\pi^2 k_x k_y}. \quad (B8) $$

Then, $\gamma(\mathbf{r}, \ell)$ is a random field in two dimensions, and its variance is given by

$$ \sigma_{\gamma, \ell}^2 = \int M(\mathbf{k})^2 d\mathbf{k} = \frac{1}{\ell^4} \int W(\mathbf{k})^2 \\ \times \frac{\sin(\pi k_x \ell)^2 \sin(\pi k_y \ell)^2}{\pi^4 k_x^2 k_y^2} dk_x dk_y. \quad (B9) $$

This formulation has been tested by creating two-dimensional arrays from $V(k)$ representing the velocity field. From these simulated fields we calculate empirically $\sigma_{\gamma, \ell}^2$, which is in agreement with the equation derived in this section.

## Appendix C

### Examples and Toy Models

In this section we use toy models to show how the parameters of the $V(k)$ change the distribution of $\gamma$ at different scales $\ell$. For these models, we set $\sigma_\omega = 1$ (Equation (B3)), the box size $L=1$, and a resolution $L/N$ with $N=10^4$. We use a function $V(k)$ of the form

$$ V(k) \propto \begin{cases} k^{-n_1} & \text{if } k_{\min} < k < k_c \\ k^{-n_2} & \text{if } k_c < k < k_{\max} \\ 0 & \text{elsewhere.} \end{cases} \quad (C1) $$

We show the dispersion in $\gamma_\ell$ as a function of $\ell$ for different choices of $V(k)$ in Figure C1. The dispersion is calculated by the integral (B9). In the left panel we vary $k_{\max}$ while $n_1 = n_2 = 0$, which represents white noise. The parameter $k_{\max}$ displaces the curve as a function of scale. In the middle panel we vary $n_1 = n_2 = n$ from -0.5 to 3.0 while $k_{\min} = 1$ and $k_{\max} = 500$. Fields with exponents between -0.5 and 1 show

similar profiles. There is a notorious degeneracy between spectra with low exponents. In the right panel $k_{\min} = 1$, $k_{\max} = 500$, $n_1 = 1.0$, $n_2 = 2.5$, and we vary $k_c$ from $k_{\min}$ to $k_{\max}$.

### C.1. Toy Models

To illustrate the behavior of the normalized circulation distribution, we create toy models with different velocity fields. First, we define the circular velocity fields, which are given by

$$ v_{\text{rot}} = R\Omega(R)\hat{\phi} \quad \Omega(R) \propto (R + r_0)^{-\beta}. \quad (C2) $$

The toy models are computed over a 2000 × 2000 grid with a box size $L=1$, and we choose $r_0 = 1/2000$. The left panel of Figure C2 shows the 16th, 50th, and 84th percentiles of $\gamma_{\text{rot}}$, the circulation of the circular velocity field, as a function of scale $\ell$. We normalize the values of $\gamma_{\text{rot}}$ such that the mean value of $\gamma_{\text{rot}}$ equals 1 at the highest resolution. We shows percentiles for $\beta \in [0.0, 0.5, 1.0, 1.5]$. For solid rotation $\beta = 0.0$ the distribution of $\gamma_{\text{rot}}$ is single valued. For other values of $\beta$ each percentile converges to an specific value at small $\ell$. Since $\Omega(R)$ is a decreasing function of radius, each percentile corresponds to $\gamma_{\text{rot}}$ measured along a unique radius in the field. It is important to notice that for these models $\gamma_{\text{rot}}$ is always positive. Variations at large scale are due to low sampling.

In the middle panel we plot the percentiles for five models of Gaussian random velocity fields defined by their spectra $V(k)$ according to Equation (C1). The parameters of each model are shown in Table C1.

From these random velocity fields we obtain $\gamma_{\text{nc}}$ for each model. We normalize $\gamma_{\text{nc}}$ fields such that their variance $\sigma_\gamma^2 = 1$ at the highest resolution. Models 1 and 2 show similar distributions as a function of $\ell$. As seen before in Figure C1, at low values of the exponent $n$ the distribution of $\gamma_{\text{nc}}$ does not depend strongly on $n$. Model 3, a single power law from $k \in [4, 64]$, lies close to model 5, which corresponds to a broken power law with $k_c = 64$. The only difference between these two models is the behavior of $V(k)$ for values of $k > 64$: for model 3 $V(k) = 0$, which is equivalent to $n_2 = \infty$, while for model 5 $V(k) \propto k^{-3}$. These show that functions $V(k)$ with high values of $n_2$ are similar to single power laws with $k_{\max} = k_c$. The right panel of Figure C2 shows the distribution of four different composite models using a velocity field with

<figure><img src="image_12.png" /><figcaption>Figure C1. Dispersion of the distribution of $\gamma(\mathbf{r}, \ell)$ as a function of spatial scale. Left: the parameters $n_1$ and $n_2$ are fixed to zero, and $k_{\min} = 1$. The parameter $k_{\max}$ is varied with values $k_{\max} \in [10, 10^2, 10^3, 10^4]$. Middle: we fix the parameters $k_{\min} = 1$ and $k_{\max} = 500$. The parameters $n_1$ and $n_2$ are equal to $n_1 = n_2 = n$, which corresponds to a function $V(k)$ with a single power law. The parameter $n$ takes values from the set [-0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0]. Right: the fixed parameters are $k_{\min} = 1$, $k_{\max} = 500$, $n_1 = 1.0$, and $n_2 = 2.5$. We vary the parameter $k_c \in [1, 2, 10, 50, 100, 200, 300, 500].$</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure C2.** Percentiles of $\gamma$ for coherent galactic rotation, $\gamma_{\text{rot}}$, and for random velocity fields, $\gamma_{\text{nc}}$, and different toy models. Each field has been computed in a 2000 × 2000 grid. Left: solid lines show median values for $\gamma_{\text{rot}}$ as a function of the scale $\ell$, while dashed lines correspond to the 16th and 84th percentiles. Middle: percentiles of $\gamma_{\text{nc}}$ as a function of $\ell$, for different models of Gaussian random fields defined in Table C1. Dashed lines show the 16th and 84th percentiles, which correspond to $1\sigma$ uncertainties for $\gamma_{\text{nc}}$. Solid lines show the median values of $\gamma_{\text{nc}}$. Right: percentiles of toy models for $\gamma_{\text{rot}} + \gamma_{\text{nc}}$. The model of $\gamma_{\text{rot}}$ has $\beta = 1$ for the four lines. The red and green lines correspond to $\gamma_{\text{rot}}$ plus model 1 and model 4, respectively. Yellow and blue lines have the same models for $\gamma_{\text{nc}}$ but with twice the magnitude.</figcaption></figure>

<table><caption>**Table C1**
Random Field Parameters</caption><thead><tr><th>Model</th><th>n</th><th>m</th><th>k<sub>min</sub></th><th>k<sub>max</sub></th><th>k<sub>c</sub></th></tr></thead><tbody><tr><td>1</td><td>0.0</td><td>...</td><td>4</td><td>500</td><td>...</td></tr><tr><td>2</td><td>1.0</td><td>...</td><td>64</td><td>500</td><td>...</td></tr><tr><td>3</td><td>1.0</td><td>...</td><td>4</td><td>64</td><td>...</td></tr><tr><td>4</td><td>1.0</td><td>3.0</td><td>4</td><td>500</td><td>16</td></tr><tr><td>5</td><td>1.0</td><td>3.0</td><td>4</td><td>500</td><td>64</td></tr></tbody></table>

$\beta = 1.0$ plus models 1 and 4 times a factor of 1 or 2. Since model 1 changes mostly at small scales, for most of the spatial scales the circulation is given by $\gamma_{\text{rot}}$ until the width of both distributions is comparable. Then, if we increase the magnitude of the random field, the transition where $\Delta\gamma_{\text{rot}} \sim \Delta\gamma_{\text{nc}}$ moves to larger scales as shown by the yellow dashed lines in Figure C2. It is important to mention that as we increase the magnitude of the random field, the number of regions with negative circulation, i.e., with retrograde rotation with respect to the galaxy, also increases.

# Appendix D
## Results with Different Models

### D.1. Rotation Curve

Our method to compute $\omega_{\text{rot}}$, and consequently $\gamma_{\text{rot}}$, is to measure a radial profile for the circular velocity field. The resulting radial profiles depend on the chosen size of the radial bins. In the main text we use the analytic function in Equation (9) to compute $\gamma_{\text{rot}}$. To test the sensitivity with respect to the chosen rotation curves, we calculate $\lambda_{\text{eq}}$ for different velocity models. The simplest model corresponds to measuring the median values of the rotation curve at intervals of 500 pc. To compute the derivatives in Equation (2), we fit a fourth-order polynomial function to $\ln(\nu(R))$ as a function of $\ln(R)$. We show the resulting $\lambda_{\text{eq}}$ as the blue line in Figure D1.

We also test the effects of adding information in the azimuthal component by means of a Fourier series expansion and include the radial component of the velocity field. The

velocity field and its vorticity are given by

$$v_{\theta}(R, \theta, m) = A_0 + \sum_{j=1}^{m=4} A_j \cos(j\theta) + B_j \sin(j\theta) \\ v_R(R, \theta, m) = C_0 + \sum_{j=1}^{m=4} C_j \cos(j\theta) + D_j \sin(j\theta) \quad (\text{D1})$$

$$\omega_z = \frac{1}{R} \left( \frac{\partial(Rv_{\theta})}{\partial R} - \frac{\partial v_R}{\partial \theta} \right), \quad (\text{D2})$$

where the coefficients $A_j$, $B_j$, $C_j$, and $D_j$ are functions of radius and $m$ is the order of the expansion. To obtain $A_j(R)$, $B_j(R)$, $C_j(R)$, and $D_j(R)$ as functions of $R$, we separate the disk in radial bins of constant width. We tested bin widths of 200, 500, and 800 pc, but they led to almost the same results in $\lambda_{\text{eq}}$. We show the effect of this complex velocity field on $\lambda_{\text{eq}}$ in Figure D1 using a bin width of 200 pc. Since $\gamma_{\text{rot}}$ has more information about the circulation of the fluid, the effect of $\gamma_{\text{nc}}$ is being noticed at smaller scales.

### D.2. Different Models for the Random Velocity Field

For our model, the spatial scale $\lambda_{\text{eq}}$, at which large-scale motions and noncircular motions contribute equally to the circulation of gas, depends on the function $V(k)$. We test different choices for $V(k)$. The first model of $V(k)$ consists of a power law ($n_1 = n_2$), and the wavenumber $k$ is bounded between the values $k_{\text{min}} = 4/L$ and $k_{\text{max}} = 1/4\Delta x$. The second model also consists of a single power law, but allowing $k_{\text{max}}$ to vary. We show these models in the right panel of Figure D1. We can see that the model used for $V(k)$ in the main text shows similar results for $\lambda_{\text{eq}}$ if we use a single power law with a variable $k_{\text{max}}$. However, if we fix $k_{\text{max}} = 1/4\Delta x$, the values of $\lambda_{\text{eq}}$ are higher. We show in Figure D2 an example of the percentiles obtained by using a function $V(k)$ with a single power law from $k_{\text{min}} = 4/L$ to $k_{\text{max}} = 1/4\Delta x$. We see that without the inclusion of the break at $k_c$ or a variable $k_{\text{max}}$ the model is unable to reproduce the plateau of the percentiles at the smallest scales.
<figure><img src="image_14.png" /></figure>

**Figure D1.** Radial profiles of $\lambda_{eq}$. Left: different choices of the large-scale velocity field. The magenta line shows $\lambda_{eq}$ obtained in the main text, with its respective $1\sigma$ uncertainties as a shaded region. The yellow line shows the resulting profile of $\lambda_{eq}$ using a Fourier series expansion for the large-scale velocity field with a radial bin of 200 pc. The blue line shows the results for the median rotation curve. Right: different models for $V(k)$. The magenta line corresponds to the model of $V(k)$ used in the main text. The model with a unique power law and fixed $k_{max}$ is shown in yellow. The blue line corresponds to a single power law and variable $k_{max}$.

<figure><img src="image_15.png" /></figure>

**Figure D2.** Percentiles of $\gamma = \gamma_{rot} + \gamma_{nc}$ as a function of the scale $\ell$, for a single power law $V(k)$ from $k_{min} = 4/L$ to $k_{max} = 1/4\Delta x$. Solid lines show the percentiles of $\gamma$ measured in the simulations. Dotted lines show the percentiles of $\gamma_{nc}$, and the shaded regions correspond to $1\sigma$ uncertainties.

## Appendix E
### Low Feedback Simulations

In this section we show the results for the simulations without early stellar feedback. Since these simulations have a lower stellar feedback, we add the suffix “-low” to distinguish them from the simulations presented in the main text.

Runs G2E1, G1E1, and G1E0.5 have 2.9, 1.3, and 1.4 times more gas than G2E1-low, G1E1-low, and G1E0.5-low,

respectively. Particularly, G2E1 has $1.8 \times 10^{10} M_{\odot}$ of gas in the disk, while all the other runs have below $1 \times 10^{10} M_{\odot}$ in gas. This makes G2E1 the most unstable disk at this point of time. This might explain why G2E1 shows the largest values of $\sigma_0$, which translates into higher values of $\lambda_{eq}$.

If we compare the distribution of G2E1 and G2E1-low, as we show in Figure E1, the former shows a broader distribution of $\gamma$ and a slightly higher fraction of regions with retrograde rotation. However, as pointed out before, the comparison is not straightforward since G2E1-low has about half the mass in gas compared to G2E1.

In Figure E2 we show the scales $\lambda_{eq}$, $\lambda_c$, and the scales of gravitational instability for G2E1-low, G1E1-low, and G1E0.5-low.

Each galaxy shows a similar behavior; at some particular radii $R^*$, $\lambda_{eq}$ decreases beyond the resolution of the simulations. We see from Figures 4 and 7 that $\omega_{rot}(R)$ decays exponentially while $\sigma_0$ decays somewhat linearly, which also holds for these simulations. At the center, where $\omega_{rot}$ peaks, $\lambda_{eq} \ll \Delta x$. Since $\omega_{rot}$ decays faster, at some particular radius $\lambda_{eq} \simeq \Delta x$, and $\lambda_{eq}$ starts to be resolved.

The regions where $\lambda_{eq}$ is resolved overlap with the regions where $\lambda_{rot}$ is resolved, once we consider the effect of the disk thickness. This shows that $\lambda_{rot}$ must be resolved in order to study $\lambda_{eq}$ in simulations.
<figure><img src="image_16.png" /><figcaption>**Figure E1.** Percentiles of circulation $\gamma$ within 6–9 kpc for G2E1 and G2E1-low, as a function of scale $\ell$. Solid lines represent the percentiles of $\gamma$ in the simulation, while dashed regions represent $1\sigma$ uncertainty intervals for the model $\mathcal{V}(k)$ around the median values, shown as dotted lines. Vertical lines show the spatial scale $\lambda_c = 1/k_c$ and its corresponding uncertainty illustrated by the dashed region. Black dashed horizontal lines show $\gamma = 0$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure E2.** Spatial scales as a function of galactocentric radius for simulations without early stellar feedback. Solid pink and yellow lines correspond to $\lambda_{eq}$ and $\lambda_c$, respectively. The shaded regions correspond to $1\sigma$ uncertainties. The classical instability scales, $\lambda_{rot}$ and $\lambda_j$, are shown as a solid light-blue line and a black dashed line, respectively. The dotted-dashed line shows the effects of the spatial resolution of the simulation on $\lambda_{rot}$.</figcaption></figure>

## ORCID iDs

José Utreras https://orcid.org/0000-0002-9226-0736

Sharon Meidt https://orcid.org/0000-0002-6118-4048

Eric Emsellem https://orcid.org/0000-0002-6155-7166

Frank Bigiel https://orcid.org/0000-0003-0166-9745

Erik Rosolowsky https://orcid.org/0000-0002-5204-2259

Eva Schinnerer https://orcid.org/0000-0002-3933-7677

## References

Agertz, O., Kravtsov, A. V., Leitner, S. N., & Gnedin, N. Y. 2013, ApJ, 770, 25

Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton, NJ: Princeton Univ. Press)

Bournaud, F., Elmegreen, B. G., Teyssier, R., Block, D. L., & Puerari, I. 2010, MNRAS, 409, 1088

Braine, J., Rosolowsky, E., Gratier, P., Corbelli, E., & Schuster, K. F. 2018, A&A, 612, A51

Bryan, G. L., Norman, M. L., O’Shea, B. W., et al. 2014, ApJS, 211, 19

Caldú-Primo, A., & Schruba, A. 2016, AJ, 151, 34

Colling, C., Hennebelle, P., Geen, S., Iffrig, O., & Bournaud, F. 2018, A&A, 620, A21

Colombo, D., Hughes, A., Schinnerer, E., et al. 2014, ApJ, 784, 3

Combes, F., Boquien, M., Kramer, C., et al. 2012, A&A, 539, A67

Courteau, S. 1997, AJ, 114, 2402

Daddi, E., Elbaz, D., Walter, F., et al. 2010, ApJL, 714, L118

Dobbs, C. L., & Bonnell, I. A. 2006, MNRAS, 367, 873

Dobbs, C. L., & Pringle, J. E. 2013, MNRAS, 432, 653

Dobbs, C. L., Pringle, J. E., & Duarte-Cabral, A. 2015, MNRAS, 446, 3608

Druard, C., Braine, J., Schuster, K. F., et al. 2014, A&A, 567, A118

Elmegreen, B. G. 1997, RMXAC, 6, 165

Elmegreen, B. G., Kim, S., & Staveley-Smith, L. 2001, ApJ, 548, 749

Elmegreen, B. G., & Scalo, J. 2004, ARA&A, 42, 211

Escala, A., & Larson, R. B. 2008, ApJL, 685, L31

Federrath, C., & Klessen, R. S. 2012, ApJ, 761, 156

Freedman, D., & Diaconis, P. 1981, Wahrscheinlichkeitstheorie verw Gebiete, 57, 453

Goldbaum, N. J., Krumholz, M. R., & Forbes, J. C. 2015, ApJ, 814, 131

Goldbaum, N. J., Krumholz, M. R., & Forbes, J. C. 2016, ApJ, 827, 28

Heyer, M., Krawczyk, C., Duval, J., & Jackson, J. M. 2009, ApJ, 699, 1092

Hopkins, P. F., & Grudić, M. Y. 2019, MNRAS, 483, 4187

Jeffreson, S. M. R., & Kruijssen, J. M. D. 2018, MNRAS, 476, 3688

Kafle, P. R., Sharma, S., Lewis, G. F., & Bland-Hawthorn, J. 2014, ApJ, 794, 59

Kennicutt, R. C., Jr. 1998, ApJ, 498, 541

Kim, J.-h., Agertz, O., Teyssier, R., et al. 2016, ApJ, 833, 202
Koch, E. W., Rosolowsky, E. W., Lockman, F. J., et al. 2018, MNRAS, 479, 2505

Koda, J., Sawada, T., Hasegawa, T., & Scoville, N. Z. 2006, ApJ, 638, 191

Kraichnan, R. H. 1967, PhFl, 10, 1417

Kruijssen, J. M. D., Dale, J. E., Longmore, S. N., et al. 2019, MNRAS, 484, 5734

Krumholz, M., & Burkert, A. 2010, ApJ, 724, 895

Krumholz, M. R., Burkhart, B., Forbes, J. C., & Crocker, R. M. 2018, MNRAS, 477, 2716

Krumholz, M. R., Dekel, A., & McKee, C. F. 2012, ApJ, 745, 69

Lang, A., & Potthoff, J. 2011, arXiv:1105.2737

Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3

Leroy, A. K., Schinnerer, E., Hughes, A., et al. 2017, ApJ, 846, 71

Matsumoto, Y., & Seki, K. 2010, JGRA, 115, A10231

Meidt, S. E., Glover, S. C. O., Kruijssen, J. M. D., et al. 2020, arXiv:2001.07459

Meidt, S. E., Leroy, A. K., Rosolowsky, E., et al. 2018, ApJ, 854, 100

Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533

Mogotsi, K. M., de Blok, W. J. G., Caldú-Primo, A., et al. 2016, AJ, 151, 15

Musacchio, S., & Boffetta, G. 2017, PhFl, 29, 111106

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Padoan, P., Haugbølle, T., & Nordlund, Å. 2012, ApJL, 759, L27

Padoan, P., Haugbølle, T., Nordlund, Å, & Frimann, S. 2017, ApJ, 840, 48

Pedlosky, J. 1992, Geophysical Fluid Dynamics, Springer Study Edition (New York: Springer)

Pranav, P., van de Weygaert, R., Vegter, G., et al. 2019, MNRAS, 485, 4167

Renaud, F., Bournaud, F., Emsellem, E., et al. 2013, MNRAS, 436, 1836

Romeo, A. B., Burkert, A., & Agertz, O. 2010, MNRAS, 407, 1223

Romeo, A. B., & Mogotsi, K. M. 2017, MNRAS, 469, 286

Rosen, A., & Bregman, J. N. 1995, ApJ, 440, 634

Rosolowsky, E. 2007, ApJ, 654, 240

Rosolowsky, E., Engargiola, G., Plambeck, R., & Blitz, L. 2003, ApJ, 599, 258

Ruppert-Felsot, J. E., Praud, O., Sharon, E., & Swinney, H. L. 2005, PhRvE, 72, 016311

Sarazin, C. L., & White, R. E., III 1987, ApJ, 320, 32

Schruba, A., Kruijssen, J. M. D., & Leroy, A. K. 2019, ApJ, 883, 2

Semenov, V. A., Kravtsov, A. V., & Gnedin, N. Y. 2016, ApJ, 826, 200

Silk, J. 1997, ApJ, 481, 703

Sormani, M. C., Treß, R. G., Ridley, M., et al. 2018, MNRAS, 475, 2383

Sun, J., Leroy, A. K., Schruba, A., et al. 2018, ApJ, 860, 172

Tan, J. C. 2000, ApJ, 536, 173

Tasker, E. J., & Tan, J. C. 2009, ApJ, 700, 358

Toomre, A. 1963, ApJ, 138, 385

Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, ApJL, 489, L179

Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9

Utreras, J., Becerra, F., & Escala, A. 2016, ApJ, 833, 13

Wada, K., Meurer, G., & Norman, C. A. 2002, ApJ, 577, 197

Ward, R. L., Benincasa, S. M., Wadsley, J., Sills, A., & Couchman, H. M. P. 2016, MNRAS, 455, 920