Dynamical Properties of Molecular-forming Gas Clumps in Galaxies at the Epoch of Reionization

TL;DR

This study uses cosmological simulations to analyze the properties and stability of molecular-forming gas clumps in a high-redshift galaxy, revealing their turbulence support, stability, and implications for future observations.

Contribution

First detailed simulation-based analysis of molecular gas clumps at reionization epoch, comparing their properties to local galaxies and assessing their stability and observability.

Findings

MGCs have masses ~10^6.5 Msun and sizes 45-100 pc, similar to local galaxies.

MGCs are highly turbulent with velocity dispersions 20-100 km/s.

Densest star-forming regions have lower virial and Toomre-Q parameters.

Abstract

We study the properties of molecular-forming gas clumps (MGCs) at the epoch of reionization using cosmological zoom-in simulations. We identify MGCs in a z=6 prototypical galaxy ("Althaea") using an H2 density-based clump finder. We compare their mass, size, velocity dispersion, gas surface density, and virial parameter (alpha_vir) to observations. In Althaea, the typical MGC mass and size are Mgas=10^6.5 Msun and R=45-100 pc, which are comparable to those found in nearby spirals and starburst galaxies. MGCs are highly supersonic and supported by turbulence, with rms velocity dispersions of sigma_gas=20-100 km s^-1 and pressure of P/k_B=10^7.6 K cm^-3 (i.e., >1000x with respect to the Milky Way), similar to those found in nearby and z~2 gas-rich starburst galaxies. In addition, we perform stability analysis to understand the origin and dynamical properties of MGCs. We find that MGCs are globally stable in the main disk of Althaea. Densest regions where star formation is expected to take place in clouds and cores on even smaller scales instead have lower alpha_vir and Toomre-Q values. Detailed studies of the star-forming gas dynamics at the epoch of reionization thus require a spatial resolution of <40 pc (=0.01"), which is within reach with the Atacama Large (sub-)Millimeter Array and the Next Generation Very Large Array.

Equations

$$\frac{1}{2}\ddot{\mathcal{I}}=2(\mathcal{T}-{\mathcal{T}}_{\mathrm{ext}})+\mathcal{W},$$

$$\frac{1}{2}\ddot{\mathcal{I}}=3M{\sigma }^2-4\pi P_{\mathrm{ext}}{R}^3-\Gamma \frac{G{M}^2}{R},$$

$$\frac{{\sigma }_{\mathrm{gas}}^2}{R}=\frac{1}{3}\left(\frac{4P_{\mathrm{ext}}}{\Sigma }+\frac{3}{5}\pi G\Sigma \right),$$

$$\frac{{\sigma }_{\mathrm{gas}}^2}{R}=\frac{\pi }{5}G\Sigma ,$$

$${\alpha }_{\mathrm{vir}}\equiv \frac{5{\sigma }_{\mathrm{gas}}^{2}R}{G{M}_{\mathrm{gas}}}=\frac{5{\sigma }_{\mathrm{gas}}^2}{\pi G\Sigma R},$$

$${\alpha }_{\mathrm{vir},\mathrm{tot}}\equiv \frac{5R}{G(M_{\mathrm{gas}}+M_{\star })}\frac{M_{\mathrm{gas}}{\sigma }_{\mathrm{gas}}^2+M_{\star }{\sigma }_{\star }^2}{M_{\mathrm{gas}}+M_{\star }},$$

$${\omega }^2={\kappa }^2-\frac{2\pi G\Sigma |k|}{1+|k|h}+{\sigma }_{\mathrm{disk}}^2{k}^2,$$

$${\kappa }^2\equiv \frac{2\Omega }{\varpi }\frac{d}{d\varpi }\left({\varpi }^2\Omega \right).$$

$$Q_{\star }\equiv \frac{{\sigma }_{\star }\kappa }{3.36G{\Sigma }_{\star }}.$$

$$Q_{\mathrm{gas}}\equiv \frac{{\sigma }_{\mathrm{gas}}\kappa }{\pi G{\Sigma }_{\mathrm{gas}}}.$$

$$T_{x}=\left\{\begin{array}[]{lccr}{\displaystyle 0.8+0.7\left(\frac{\sigma_{z}}{\sigma_{r}}\right)}&&&\mbox{if\ }\sigma_{z}\gtrsim 0.5\times\sigma_{r}\\[12.50002pt] {\displaystyle 1+0.6\left(\frac{\sigma_{z}}{\sigma_{r}}\right)^{2}}&&&\mbox{if\ }\sigma_{z}\lesssim 0.5\times\sigma_{r}\\ \end{array}\right.$$

$$Q_x^{\mathrm{thick}}=T_{x}Q,$$

$$Q^{-1}_{\rm eff}=\left\{\begin{array}[]{lccr}{\displaystyle\frac{w}{Q^{\rm thick}_{\star}}+\frac{1}{Q^{\rm thick}_{\rm gas}}}&&&\mbox{if\ }Q^{\rm thick}_{\star}\geq Q^{\rm thick}_{\rm gas}\\[7.5pt] {\displaystyle\frac{1}{Q^{\rm thick}_{\star}}+\frac{w}{Q^{\rm thick}_{\rm gas}}}&&&\mbox{if\ }Q^{\rm thick}_{\star}\leq Q^{\rm thick}_{\rm gas}\\ \end{array}\right.$$

$$w\equiv \frac{2{\sigma }_{\star }{\sigma }_{\mathrm{gas}}}{{\sigma }_{\star }^2+{\sigma }_{\mathrm{gas}}^2},$$

$${\sigma }_{\mathrm{gas}}^2={\sigma }_{\mathrm{bulk}}^2+c_s^2+{\sigma }_{\mathrm{NT}}^2.$$

$$⟨x⟩\equiv \frac{{\sum }_i{\rho }_i{x}_i}{{\sum }_i{\rho }_i},$$

$${\sigma }_{\mathrm{bulk}}^2=\frac{1}{3}⟨{\left|\mathbf{v}-⟨\mathbf{v}⟩\right|}^2⟩,$$

$$c_s^2=⟨\frac{{\mathrm{k}}_{\mathrm{B}}{P}_{\mathrm{TH}}}{{\mathrm{m}}_{\mathrm{p}}n}⟩,$$

$${\sigma }_{\mathrm{NT}}^2=⟨\frac{{\mathrm{k}}_{\mathrm{B}}{P}_{\mathrm{NT}}}{{\mathrm{m}}_{\mathrm{p}}n}⟩.$$

$$\mathcal{M}=⟨\sqrt{1+\frac{P_{\mathrm{NT}}}{P_{\mathrm{TH}}}}⟩.$$

$$f_{\mathrm{gas}}=\frac{M_{\mathrm{gas}}}{\left(M_{\mathrm{gas}}+M_{\star }\right)}.$$

Full text

††slugcomment: Accepted to the ApJ.††software: Python (Van Rossum & de Boer 1991), Astropy (Astropy Collaboration et al. 2013), Cython (Behnel et al. 2011), Matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011), PyMSES (Labadens et al. 2012), SciPy (Jones et al. 2001–), and yt (Smith et al. 2009; Turk et al. 2011)

1. Introduction

The growth of galaxies and their subsequent evolution are governed by the baryon cycle—galaxies accrete gas from the intergalactic medium (IGM) either directly from the cosmic web, or through mergers with other galaxies. This gas fuels star formation and feeds central supermassive black holes. Subsequent feedback replenishes and enriches the circumgalactic medium by expelling some part of this material. Existing studies indicate that early galaxies are more gas-rich with molecular fractions higher than those of present-day galaxies (e.g., van de Voort et al. 2011; Decarli et al. 2016, 2017). Massive gas inflows from the IGM trigger gravitational instabilities that lead to the formation of dense gas clumps in which molecules form quickly (hereafter molecular-forming gas clumps111Note that by clumps, we refer to kiloparsec-scale gas concentrations in galaxies, rather than molecular gas concentrations small than molecular clouds and containing cores.; MGCs) that are typically more massive ($M_{\rm gas}\simeq 10^{9}$ $M_{\odot}$) and extended ($\simeq$ sub-kpc) than those observed in nearby galaxies (e.g., Gabor & Bournaud 2013; Hopkins et al. 2014; Inoue et al. 2016). Some theoretical works argue that the migration of such giant massive clumps through disks is largely responsible for the buildup of the bulges of massive galaxies at redshift $z$ $\sim$ 0 (e.g., Ceverino et al. 2010).

Early galaxies have higher star formation rates (SFR; Behroozi et al. 2013; Sparre et al. 2015; Maiolino et al. 2015; Dunlop et al. 2017) and smaller sizes (e.g., Bouwens et al. 2011; Ono et al. 2013) compared to the local population (see also a review by Stark 2016). As a consequence, we expect them to be significantly more ionized, and have intense and hard interstellar radiation fields (ISRF). Since their metallicity and dust content are also expected to be lower in these early evolutionary stages, shielding of UV photons—responsible for the photoheating of the gas— is strongly reduced. Such differences in turn affect the regulation of the thermal and chemical state of the multi-phase interstellar medium (ISM).

Studying ISM properties of early galaxies is essential for understanding how star formation proceeds under more extreme conditions. Even in the local Universe, where detailed observations can be made, variations in molecular cloud properties have been observed between different galaxy populations (see e.g., Hughes et al. 2010, 2013). Given that high-$z$ galaxies statistically represent the early evolutionary stages of present-day galaxies, it is thus reasonable to pose the question: what are the physical properties of MGCs in early galaxies, and how do they differ from those found in local galaxy populations?

FIR fine-structure lines (e.g., [CII], [NII], and [OIII]), and CO and [C i] lines are key diagnostics for constraining the ISM conditions of galaxies. They also provide highly complementary information on different ISM phases (ionized, atomic, molecular; e.g., Scoville & Solomon 1974; Rubin 1985; Malhotra et al. 2001). Global measurements of these diagnostics in high-$z$ galaxies have provided preliminary information on their global properties (e.g., gas masses, gas temperature, and radiation field intensity). However, spatially resolving their ISM is necessary to fully understand many aspects of galaxy evolution and the physics behind their intense star formation (SFR $\sim$ 100$-$3000 $M_{\odot}$  yr*-1*; see e.g., the review by Carilli & Walter 2013). To date, spatially resolved ISM properties of high-$z$ galaxies have only been mapped in a handful of galaxies at high redshift using tracers such as dust continuum, CO, and [CII] lines (e.g., Swinbank et al. 2011; Hodge et al. 2015; Ferkinhoff et al. 2015; Hodge et al. 2016; Leung et al. 2019). These studies find that galaxies close to the peak of cosmic star formation ($z$ $\sim$ 2) are more molecular gas-rich, turbulent, and clumpy than nearby galaxies, although luminous compact blue galaxies may represent local analogs of these gas-rich, clumpy systems (Garland et al. 2015).

Earlier epochs still represent an essentially uncharted territory for ISM investigations. At present, it remains unclear how star formation proceeds in the (sub-)$L^{*}$ galaxy population at $z$ $\gtrsim$ 6 which is responsible for producing the bulk of the UV photons that reionized the Universe. High resolution hydrodynamics simulations have been carried out to investigate the global properties, structures, and morphologies of galaxies out to the epoch of reionization (EoR) and their importance for providing the ionizing photons (see e.g., Ceverino et al. 2017; Katz et al. 2017; Ma et al. 2018; Trebitsch et al. 2018; Rosdahl et al. 2018). Compared to the zoom-in galaxies presented in these works, the simulation used here focuses on a Lyman-break galaxy (Pallottini et al. 2017b) from the serra simulation suite (see Pallottini et al. 2019 for a comparison between these zoom-in simulations). Most notably, we include a non-equilibrium thermo-chemical network to follow the formation of H2, which is crucial to examining the properties of MGCs.

While ALMA has enabled the detection of [CII] 158 $\micron$ and CO line emission in normal (SFR$<$ 100 $M_{\odot}$  yr*-1*) galaxies at $z$ $>$ 6 over the past few years (e.g., Carniani et al. 2018; D’Odorico et al. 2018), the first spatially resolved observations are just starting to become available (e.g., Jones et al. 2017). On the theoretical side, a schematic investigation by Behrendt et al. (2016a) revealed that kiloparsec-scale clumps likely consist of many smaller clouds formed by gravitational instability, a hierarchical structure that persists in an idealized disk model by Behrendt et al. (2019).

To understand the physical properties of MGCs in early galaxies, we have undertaken a detailed numerical study whose aim is to characterize the dynamical properties of the star-forming MGCs in prototypical (i.e., $L^{*}$) galaxies in the EoR.

The paper is structured as follows222Throughout this paper, we adopt a concordance cosmology, with total matter, vacuum and baryonic densities in units of the critical density $\Omega_{\Lambda}$ = 0.692, $\Omega_{m}$ = 0.308, $\Omega_{b}$ = 0.0481, Hubble constant $H_{0}$ = 100 $h$ km s*-1*  Mpc*-1* with $h$ = 0.678, spectral index $n$ = 0.967 and $\sigma_{8}$ = 0.826 (Planck Collaboration et al. 2014).. We start by providing some physical background in §2. In §3, we describe the setup of our simulation and properties of our main galaxy (Althæa). In §4, we describe the method used to identify MGCs, and present the formalism within which we interpret the results. In §5, we present the results and characteristic properties of the MGCs. We then interpret the results and discuss the implications of our findings in §6, and give our conclusions in §7.

2. Physical Background

We begin by introducing some empirical relations commonly mentioned in observational studies of molecular clouds in the literature, as well as describing the gravitational instability of galactic disks, which might be driving MGC formation. These concepts will be used in the subsequent analysis of our simulations.

2.1. Larson’s Relations

Larson (1981) discussed a number of relations among Galactic molecular cloud properties, namely the linewidth-size, density-size, and mass-size relations. Larson relations are routinely used for comparing properties of molecular structures in different galactic environments. They also represent a useful framework to analyze our results as they have been argued to arise from the interplay between gravity and velocity dispersion given by virial theorem (note, however, the alternative interpretation involving gravitational collapse by e.g., Ballesteros-Paredes et al. 2011a).

The virial theorem for a distribution of unmagnetized gas can be written as (McKee & Zweibel 1992)

[TABLE]

where $\ddot{\mathcal{I}}$ is the second time derivative of the Lagrangian moment of inertia, $\mathcal{T}$ is the internal energy of the gas (including thermal, turbulent, and bulk motions), $\mathcal{T}_{\rm ext}$ is external pressure support, and $\mathcal{W}$ is the gravitational energy. Let us specialize to the case of a spherical self-gravitating molecular cloud of mass $M_{\rm gas}$, radius $R$, and root-mean-square velocity dispersion $\sigma$, accounting for both thermal and turbulent contributions. Defining $P_{\rm ext}$ as the external pressure, Equation 1 can be written as

[TABLE]

where $\Gamma$ is a geometrical factor that is equal to 3/5 for a uniform sphere; in Equation 2 the terms on the right-hand side represent the kinetic energy, external pressure, and gravitational potential energy terms.

Motivated by Larson’s linewidth-size relation (Larson 1981) and the work by Heyer et al. (2009), we assume equilibrium (i.e., ${\ddot{\mathcal{I}}}=0$), define the cloud surface density as $\Sigma$ = $M/\pi R^{2}$, and rewrite the previous equation as

[TABLE]

which further reduces to

[TABLE]

if the external pressure $P_{\rm ext}=0$. For this case (often referred to as simple virial equilibrium) from the balance between kinetic and gravity terms we can define the virial parameter as

[TABLE]

Based on Equation 4, a one-to-one mapping between $\sigma_{\rm gas}^{2}/R$ and $\Sigma$ is therefore expected for a virialized cloud, since $\sigma_{\rm gas}^{2}/R\propto\Sigma$. (Heyer et al. 2009) pointed out that the original size-linewidth relation implies constant surface density in this interpretation. The limited dynamic range in column density of CO observations appears to account for this (see also Ballesteros-Paredes et al. 2011a). Deviations from this relation are often attributed to a significant contribution from external pressure as per Equation 3 (see e.g., Heyer et al. 2009; Hughes et al. 2010, 2013; Meidt et al. 2013). Similar conclusions are also reached from the analysis of clouds forming in a Milky-Way like galaxy simulated with a $\simeq 4\,{\rm pc}$ resolution (i.e., Grisdale et al. 2018).

Summarizing, the virial parameter can be used to quantify the stability/boundedness of a molecular cloud. Accounting for the external pressure, a virial parameter of $\alpha_{\rm vir}\lesssim 2$ would be unstable (Bertoldi & McKee 1992). Such a criterion is often used in observations (see e.g., Kauffmann et al. 2017b).

In cases where a stellar component plays an important role in the dynamics, the virial parameter becomes

[TABLE]

where $M_{\star}$ is the stellar mass enclosed within the MGC volume.

2.2. Toomre Stability

The onset of gravitational instability is tightly connected to star formation (e.g., Kennicutt 1989; Martin & Kennicutt 2001; Wang & Silk 1994; Li et al. 2005, 2006). For axisymmetric modes, the dispersion relation for the growth of density perturbations in a rotating, turbulent disk of finite thickness $h$ is described by

[TABLE]

where $k$ is the wavenumber and $\kappa$ is the epicyclic frequency, defined as:

[TABLE]

(Romeo 1992), where $\varpi$ is the galactic radius. In Equation 7 the terms on the right hand side are related to rotation, self-gravity and internal pressure, respectively. Heuristically, the instability can be understood by considering the scale at which gravitational potential overcomes the internal energy. Gravity dominates at scales $L>L_{J}$, where $L_{J}$ is the Jeans length. However, differential rotation in disk galaxies can stabilize perturbations that might otherwise collapse for $L>L_{\rm rot}$, where $L_{\rm rot}$ is set by $\kappa$. As a result, disks are unstable to gravitational collapse on scales between $L_{J}<L<L_{\rm rot}$.

From the dispersion relation, a parameter $Q$ can be derived such that $Q<1$ when instability occurs, that reproduces this inequality to order unity. For a collisionless fluid—such as an ensemble of stars— this parameter is (Toomre 1964)

[TABLE]

The equivalent parameter for a collisional gas was derived by Goldreich & Lynden-Bell (1965)

[TABLE]

In the thin disk approximation ($kh\ll 1$), instability occurs on scales $k$ such that $Q<Q_{\rm crit}\simeq 1$ (or equivalently $\omega^{2}<0$ in Equation 7). A frequently used observable proxy for $Q$ is the ratio of disk circular velocity to root-mean-square velocity dispersion $v_{\rm circ}/\sigma_{\rm disk}$ (e.g., García-Burillo et al. 2003; Genzel et al. 2011; Kassin et al. 2012; Leung et al. 2019).

In our stability analysis, we account for the combined effect of gas and stars (derived exactly by Rafikov 2001), and for the non-negligible disk thickness. This is done by adopting an approximation for an effective two-component $Q_{\rm eff}$ parameter (i.e., Romeo & Wiegert 2011; Romeo & Falstad 2013, see also Inoue et al. 2016). The effect of disk thickness modifies the $Q$ parameter for gas and stars by accounting for the vertical velocity dispersion

[TABLE]

and

[TABLE]

with $x$ indicating either gas or stars. The combined effect of gas and stars can then be accounted for by writing

[TABLE]

where the relative weight $w$ is defined as

[TABLE]

Conceptually, the finite disk thickness reduces the gravity in the vertical direction, thereby making it easier for a system to maintain stability, and thus lowering the critical Toomre $Q_{\rm crit}$ from $\simeq$ 1 to 0.67 (Goldreich & Lynden-Bell 1965). On the other hand, including the contribution of the stellar component promotes gravitational instability, and thus increases $Q_{\rm crit}$, more-so if the stars have low velocity dispersion. As a rule of thumb, $Q_{\rm crit}\,=\,1.34$ for $Q_{\rm gas}$ = $Q_{\star}$.

3. Numerical Simulations

The simulations used in this work are described by Pallottini et al. (2017b, a) and are briefly summarized here. Serra333Greenhouse in Italian. is a suite of cosmological zoom-in simulations performed using Eulerian hydrodynamics and adaptive mesh refinement (AMR) techniques to achieve high spatial resolution in regions of interest (e.g., regions of high density). In particular, it uses a modified version of ramses (Teyssier 2002) as the AMR backend. The simulation used here covers a comoving box of 20 Mpc $h$-1 in size. The simulation zooms in on a target halo of mass M_{\rm DM}\simeq 10^{11}\mbox{M_{\odot}} at $z$ = 6. The Lagrangian region of the halo (2.1 Mpc $h$-1) has a dark matter mass resolution of \simeq 6\times 10^{4}\mbox{M_{\odot}}, equivalent at initial density to a baryonic mass resolution of 1.2\times 10^{4}\mbox{M_{\odot}}. This region is spatially refined with a quasi-Lagrangian criterion based on a mass threshold, so that a cell is refined if its total (dark+baryonic) mass exceeds the mass resolution by a factor of eight. The finest refined cell allowed in the zoom in region has size $l_{\rm cell}\simeq$ 30 pc (at $z=6$), i.e., sizes are comparable to the sizes of local giant molecular clouds (e.g., Sanders et al. 1985; Federrath & Klessen 2013; Goodman et al. 2014),

The models include a non-equilibrium chemical network (Krome) following e*-, H+, H-, He, He+, He++*, H2, and H${}_{2}^{+}$ (Grassi et al. 2014; Bovino et al. 2016). Of particular importance here is that the abundances are calculated using an on-the-fly non-equilibrium formation of molecular hydrogen scheme described by Pallottini et al. (2017b). In the simulation, the UV radiation field is assumed to be uniform spatially. Self-shielding of H2 from photo-dissociation is accounted for using the Richings et al. (2014) prescription, while the formation rate on dust grain is computed following Jura (1975). Photo-dissociation via Lyman-Werner photons is included as part of the Krome network. We do not assume a clumping factor in the chemistry solver (i.e., $C_{\rho}$ = 1, cf. Lupi & Bovino 2020). In the simulation, star formation follows an H2-based Schmidt-Kennicutt relation.

Star formation is modeled using an H2-based prescription of the Schmidt-Kennicutt relation (Krumholz et al. 2009). We adopt stellar tracks from starburst99 and include stellar feedback from supernovae (SNe) as well as winds from OB and asymptotic giant branch stars. Coupling to the gas is implemented via a sub-grid model for blastwaves, which accounts for radiative energy losses inside the cell. The remaining energy is injected into the ISM in both thermal and non-thermal form (Agertz et al. 2013); the latter is a term that mimics unresolved turbulence. The non-thermal energy is not affected by radiative cooling but rather is dissipated on the eddy turn-over time scale (Mac Low 1999). See Pallottini et al. (2017a) for further details. Radiation pressure on dust and gas is also included. Photoionization from local sources is neglected in this simulation (see Pallottini et al. 2019; Decataldo et al. 2019, for its effect).

The simulation zooms in on a galaxy named after the flower Althæa, whose properties are given in Pallottini et al. (2017b) and are briefly summarized in the following. Althæa is a Lyman-break galaxy that at $z\simeq 6$ is hosted by a dark matter halo of mass $M_{\rm DM}\simeq$ = 3.5 $\times$ 1011 $M_{\odot}$ at the center of a cosmic web knot, and accretes mass from the IGM mainly via three filaments of length $\simeq$ 100 kpc. Althæa has a stellar mass of $M_{\star}\simeq$ 2.6$\times 10^{10}$ $M_{\odot}$, a metallicity of $Z\simeq$ 0.5 $Z_{\odot}$, a molecular gas mass of $M_{\rm H2}\simeq$ 5$\times 10^{7}$ $M_{\odot}$, and a SFR of 30 – 80 $M_{\odot}$  yr*-1*. The specific SFR of Althæa is ${\rm sSFR}\equiv{\rm SFR}/M_{\star}\simeq$ 4–40  Gyr*-1*, compatible with the sSFR vs. $M_{\star}$ relation observed at high redshift (Jiang et al. 2016). The effective stellar radius of Althæa is $\simeq$ 0.5 kpc and the dark matter virial radius is $r_{\rm 200}\simeq$ 15 kpc. The stellar-to-halo mass ratio of Althæa is $M_{\star}/M_{\rm DM}\,\simeq\,$0.07. This is a factor of about three higher than the typical value inferred by abundance matching models, which typically give values of $\sim 0.02$. But studies show that this ratio varies as the halo mass becomes much larger or smaller than from $M_{h}$ $\sim$ 1012 $M_{\odot}$ (see Moster et al. 2013; Behroozi et al. 2013; Behroozi & Silk 2015; Moster et al. 2018)444 This ratio from abundance matching can vary; for instance, Behroozi et al. (2013) (main method) finds a ratio of 0.015 – 0.025 for a halo mass between $M_{\rm halo}$ $\sim$ 1010.5–$10^{12}$ $M_{\odot}$, whereas in the same paper the authors find 0.01–0.05 in the same mass range from direct abundance matching (see the Appendix therein). As another example, from an earlier abundance matching study, Moster et al. (2013) find a ratio 0.005, also differing by a factor three from Behroozi et al. (2013), though in opposite direction.. The ratio of Althæa is higher than e.g., Katz et al. (2017) and Ceverino et al. (2017) which could be a result of different feedback implementation.

3.1. Star Formation History

One of the main advantages of studying galaxies in simulations is that we can examine how their dynamical properties evolve with time, which is interesting in order to understand the physical processes that determine the morphology and dynamical properties of galaxies, which in turn affect their star formation. This is advantageous especially at early cosmic epochs, when the densest structures are beginning to form; gas is constantly being accreted onto the central galaxy from the cosmic web and satellite galaxies, thereby leading to bursts of star formation. Meanwhile, tidal forces resulting from interactions with these surrounding galaxies can disrupt the main disk and arms, likely leading to different dynamical states for the molecular structures compared to more evolved galaxies found at a later cosmic time (e.g., some molecular structures may disperse while others may agglomerate into more massive ones). The star formation history555Note that here the SFR is calculated based on the stellar mass formed in the last 4 Myr within 3.5 kpc from the galaxy center of mass. The SFR plotted in Figure 2 of Pallottini et al. (2017a) is a factor of two higher since there the SFR accounts for the contribution from massive satellite galaxies within the virial radius ($\approx$15 kpc). of Althæa is shown in Figure 1. The SFR of Althæa varies between $\sim$30–80 $M_{\odot}$  yr*-1* as it evolves from an actively accreting phase to a starburst phase after a merger, and then back to a relatively quiescent phase, over the simulated $\approx 700$ Myr.

In §5.1, we show the importance of rotation support from large-scale motions in the MGC dynamics by comparing their properties in the disturbed phase of Althæa and in the disk-like ordered phase. Given the stochastic nature of Althæa in its star formation history, we mostly focus on few of its most extreme evolutionary stages in this work (see §5.2). These phases correspond to (a) an intensely accreting phase and (b) a starburst phase (Figure 1). We are interested in determining whether MGC properties are sensitive to these different dynamical conditions. For completeness, we also show relationships of MGC properties examined for other evolutionary stages traced in the simulation (see §5.3).

4. Molecular-forming Gas Clumps

4.1. Identification

To identify the molecular complexes, we use a customized version of the clump-finding algorithm available in the python package yt (Turk et al. 2011), which was initially described in Smith et al. (2009) and modified since. The latest version of the default yt clump finder decomposes the zones of the simulation into non-overlapping tiles, which are stored in a three-dimensional tree that can be processed using $k$-dimensional tree algorithms. It then identifies the contours of a variable field (here, the density field) within a tile and connects them across the tiles. In the customized version used for this study, we enhanced the stability of the code. Due to the nature of our AMR simulation, we regrid the simulation data into uniform grids. The grid size is defined based on the highest resolution of the simulation data, i.e., the less refined regions are supersampled in the resulting uniform grids.

In the clump-finding process (in position-position-position space; PPP), we employ a set of different density thresholds defined based on the molecular hydrogen density of Althæa at different evolutionary stages ($z=6.0$–7.2). We note that this process is the three-dimensional analog to identifying molecular structures based on the noise levels in position-position-velocity (PPV) maps that observers obtain with telescopes, using molecular line tracers such as CO, CS, and HCN. This is commonly done by identifying clumps based on/after applying S/N-clipping, using tools such as the aips tasks serch, clumpfind, and cprops; (e.g., Williams et al. 1994; Oka et al. 2001; Rosolowsky & Leroy 2006; Rosolowsky et al. 2008; Donovan Meyer et al. 2013)). Existing studies find a good correspondence in the dynamical properties extracted in PPV- versus PPP-space for well-isolated structures (Ballesteros-Paredes & Mac Low 2002a; Heitsch et al. 2009; Shetty et al. 2010; Beaumont et al. 2013; Pan et al. 2015, but see Ballesteros-Paredes & Mac Low 2002a and Shetty et al. 2010 for a discussion on caveats and limitations).

In Figure 2, we show the distribution of H2 number density $n_{\rm H2}$ for Althæa during its accreting phase, including the contribution from the gas within 3.5 kpc from the galaxy center. We note that the distribution is almost flat for $n_{\rm H2}\gtrsim 1$ ${\rm cm}^{-3}$ and it samples the range of densities where clumps are found based on morphological analysis666For $n_{\rm H2}\gtrsim 1$ ${\rm cm}^{-3}$, the fourth Minkowski functional of the H2 density field is significantly larger than zero. This implies that the field is made of isolated components. See Fig. 6 in Pallottini et al. (2017a). (Pallottini et al. 2017a). In each evolutionary stage, we identify MGCs by applying density cuts to the H2 density distribution, i.e., MGCs are selected at $n_{\rm cut}\leq n_{\rm H2}$. We select 10 equally-spaced cuts in H2 density in log scale: ${(n_{\rm cut}/{1{\rm cm}^{-3})}}\,=\,[0.32,0.53,0.88,1.45,2.45,4.08,6.81,11.36,18.96,31.62]$. Note that with these choices we are including MGCs that are not fully molecular ($n=0.5n_{\rm H2}$).

To visually display the clump finding procedure, we overplot the molecular structures identified using a subset of the H2 density cuts ($n_{\rm cut}$ = 0.53, 5.81, and 18.96 ${\rm cm}^{-3}$) on the H2 density maps (Figure 3). Since the molecular structures are identified in the 3D H2 density field, they can appear as overlapping structures depending on the viewing angle; thus we also plot them in different projections so that the identification can be more easily appreciated. We repeat this identification process for 14 evolutionary stages between redshift $z$ $\in$[6.0, 7.2], spaced by $\Delta t$ = 15 Myr.

We impose the additional constraint that an identified structure must be composed of at least 10 cells. We caution that an important caveat of such a constraint is that we can only examine the parameter space of cloud complexes of radius $R\gtrsim 40$ pc, because of the resolution limit of the simulation.

4.2. MGC Properties

Upon identifying the molecular structures, we extract properties such as the gas mass $M_{\rm gas}$, effective size $R$, Mach number $\mathcal{M}$, velocity dispersion $\sigma_{\rm gas}$, and gas surface density $\Sigma_{\rm gas}$ to examine their dynamics.

The mass of an MGC is calculated from the uniformly-gridded 3D density field, integrating over the MGC volume $V$. The effective size is defined assuming spherical geometry, i.e., $R\equiv(3V/4\pi)^{1/3}$. The full velocity dispersion of MGCs is calculated from the bulk velocity field ($\sigma_{\rm bulk}$), thermal sound speed ($c_{s}$), and the non-thermal velocity dispersion ($\sigma_{\rm NT}$)

[TABLE]

In observations of MGCs, the linewidth contribution of dense gas exceeds that from the diffuse gas. Therefore, when calculating the global quantities of MGCs, we always perform a mass-weighting. Since we operate on data on a uni-grid, this is equivalent to density averaged quantities. In general, for the quantity $x$ in a MGC we write

[TABLE]

where the sum is done for the cells indexed by $i$ composing the MGC. We use the definition Equation 16 to write each term of the right-hand side of Equation 15 as follows. The bulk velocity dispersion is

[TABLE]

where Equation 16 is applied to each velocity component. The thermal sound speed is calculated from the thermal pressure ($P_{\rm TH}$) through

[TABLE]

where the pressure is in units of K ${\rm cm}^{-3}$, $\rm k_{B}$ is the Boltzmann constant, $m_{p}$ is the mass of a proton. Similarly, contribution from non-thermal energy is calculated from non-thermal pressure ($P_{\rm NT}$) as follows

[TABLE]

Finally, the Mach number is related to the pressure terms as follows

[TABLE]

5. Results

5.1. MGC Basic Properties

We start by considering the distributions of MGC total gas mass $M_{\rm gas}$, radius $R$, and gas mass fraction

[TABLE]

We show in Figure 4 the distributions from the combination of the three considered evolution stages. We show results for a low density cut of $n_{\rm cut}$ = 0.32 ${\rm cm}^{-3}$ (top panels) and a high cut777The highest density threshold we used of $n_{\rm cut}$ = 31.62 ${\rm cm}^{-3}$ corresponds to a minimum MGC mass of the order of 105.5 $M_{\odot}$ for the densest structure. However, not all evolutionary stages considered have at least one MGC at this highest density threshold, so we consider MGCs identified with the second highest density threshold instead. of $n_{\rm cut}$ = 18.96 ${\rm cm}^{-3}$ (lower panels). Overall, the mass of all MGCs identified ranges from $M_{\rm gas}\simeq$105.5-8.5 $M_{\odot}$, whereas the gas fraction ranges from $f_{\rm gas}\simeq$ 0.1 – 1.

The most massive structures ($M_{\rm gas}>$107 $M_{\odot}$) correspond to the molecular disk of the galaxy (see e.g., Figure 3), which is identified as a single component at low $n_{\rm cut}$. This main disk component occupies the top left corner of Figure 5, with a low gas mass fraction of $f_{\rm gas}<5\%$. The velocity dispersion of such a structure is dominated by the bulk motion term888The bulk velocity dispersion term is likely associated with the rotational velocity, see Kohandel et al. (2019) for kinematical analysis of Althæa., as shown in Figure 6, where we compare the velocity dispersions of MGCs resulting from bulk ($\sigma_{\rm bulk}$) versus non-thermal turbulent motions ($\sigma_{\rm NT}$) in the accreting and quiescent phases of Althæa. The accreting phase of Althæa displays a disturbed morphology, whereas the quiescent phase displays a disk-like morphology (Figure 3). Excluding the structure corresponding to the molecular disk, the velocity dispersions of most MGCs in the accreting phase of Althæa are dominated by turbulent motions, while in the quiescent disk-like phase, $\sigma_{\rm bulk}$ is at least comparable to $\sigma_{\rm NT}$ for almost half the MGCs. For these reasons, the most massive MGCs are excluded in the discussion of the dynamics (§6), as they are dominated by large-scale shear; however, we include them in the scaling relation plots (Figure 7 and Figure 8), as would be done for an observational study.

Considering only the MGCs identified at the highest density threshold, the mass of MGCs is $\simeq 10^{6.5}$ $M_{\odot}$. We note that similarly massive molecular structures (few times 107 $M_{\odot}$ 999See §3.1 of Behrendt et al. (2016b).) have been reported in idealized closed-box isolated galaxy simulations done at higher resolution (e.g., a maximum resolution of 3 pc in the 48 kpc box studied by Behrendt et al. 2016b). However, the stellar feedback, IGM, merger and accretion histories are not properly modeled in such simulations, e.g. see Grisdale et al. (2017, 2019) for the importance of feedback in shaping the properties of clouds forming in local simulated galaxies. That said, the similar mass range found in MGCs of Althæa is reassuring — our results are not far off in spite of the limited resolution ($l_{\rm cell}\simeq$ 30 pc).

The local sound speed of all MGCs identified is typically much smaller than their non-thermal (turbulent) velocities, as non-thermal pressure dominates thermal pressure for dense gas (Pallottini et al. 2017a), i.e., $c_{s}^{2}\ll\sigma_{\rm NT}^{2}$. In particular, the average Mach number for MGCs identified at the highest density threshold is $\bar{\mathcal{M}}\,\simeq 6$. This is consistent with the analysis done in Vallini et al. (2018), that finds a global Mach number of $\mathcal{M}\sim 10$ for Althæa.

5.2. Single Evolutionary Stage

Here, we focus on the MGC properties in the two most extreme evolutionary stages of Althæa — the accreting and starburst phases (see §3.1). These properties are plotted in Figure 7 together with observational results of nearby and $z$ $\simeq$2 molecular structures for comparison.

During the accreting phase, MGCs in Althæa are characterized by large velocity dispersions ($\sigma\simeq 100\mbox{ km s}^{-1}$) and sizes ($R\simeq 100$ pc). These values are comparable to those found in starburst galaxies, such as the nearby gas-rich galaxy M64 and the $z$ $\sim$ 2 starbursting disk galaxy SMM J2135-0102 (Rosolowsky & Blitz 2005; Swinbank et al. 2011), but are higher than those found in the Milky Way and extragalactic GMCs by an order of magnitude (Heyer & Brunt 2004; Bolatto et al. 2008).

While the virial parameter defined as Equation 5 is typically used in observational studies, as they are designed to probe structures that are composed of mainly molecular gas only, we calculate the virial parameter of MGCs including the influence from the stellar component (i.e., $\alpha_{\rm vir,tot}$ via Equation 6) given the high stellar-to-gas mass ratios of some molecular structures identified (see Figure 5). The biggest MGC identified at low $n_{\rm cut}$ and with the highest stellar-to-gas mass ratio ($\gtrsim$ 1; Figure 5) corresponds to the central main disk of Althæa (see top panels of Figure 3 and §4.2). The large amount of stellar mass already assembled in Althæa explains the low total virial parameter seen in the middle left panel of Figure 7, where the stellar gravitational potential influences the overall stability. On the other hand, the MGCs in this phase are mostly stable101010We check that the gas-only virial parameter of MGCs is comparable to the total virial parameter., with virial parameter within the range observed in the molecular structures of the Milky Way but higher than its average value, though the MGCs are more massive. The higher than typical MW value is linked to the higher turbulence of these systems.

Motivated by observational studies, we also plot the $\sigma_{\rm gas}^{2}/R$ ratio and gas surface density of MGCs in the accreting phase of Althæa in the bottom left panels of Figure 7. In the same plot, we show those observed in the Galactic Ring Survey (GRS) of the Milky Way (Heyer et al. 2009) for comparison. Dashed lines in the figure show the loci along which the annotated external pressures are needed for any molecular clouds in equilibrium to have certain linewidths for a given set of surface densities (see §2.1). At a given gas surface density, MGCs display a range of $\sigma_{\rm gas}^{2}/R$. Such variation in $\sigma_{\rm gas}^{2}/R$ is observed in GMCs in the central and outer part of the Milky Way (Oka et al. 2001; Heyer et al. 2009).

If constant column density (Larson’s third relation) were truly a fundamental property of molecular clouds, and if virial equilibrium were a universal property of molecular clouds with negligible surface pressure, one would expect a single point in the $\sigma_{\rm gas}^{2}/R$–$\Sigma$ relation at that mass surface density $\Sigma$. The variation in $\sigma_{\rm gas}^{2}/R$ with $\Sigma$ seen here suggests that column density is not a fundamental property of molecular clouds in agreement with Heyer et al. (2009). Previous studies suggest that the observed mass-size relation of molecular clouds is solely a result of $\rho\propto MR^{-3}$ and may be an artifact of the limited range of column densities a specific molecular line tracer is sensitive to (see e.g., Ballesteros-Paredes & Mac Low 2002b; Ballesteros-Paredes et al. 2011b).

In the starburst phase, MGCs have velocity dispersion and sizes similar to those in the accreting phase; however, MGCs in the starburst phase span a wider range in gas surface density (bottom panels of Figure 7) and have lower virial parameters. That is, MGCs in the starburst phase are more susceptible to collapse, as expected. The gas surface density of MGCs in both accreting and starburst phases is higher than that in the solar neighborhood of the Milky Way and the Large Magellanic Cloud (LMC), but comparable to those observed in starburst galaxies and the nearby ultra-luminous IR galaxies (ULIRGS; Boulares & Cox 1990; Scoville et al. 1991; Weiß et al. 2001; Hughes et al. 2010; Leroy et al. 2015).

5.3. Adopted Density Threshold Dependence

As seen in §5.1, the cuts used to identify the molecular structures are based on the H2 number density, which ranges from nH2 $\sim$ 0.1 – 30 ${\rm cm}^{-3}$. We do not explicitly use a minimum H2 fraction to select the structures. Among the various cuts, two physical regimes are identified, corresponding to materials in the disk/interarm material and MGCs (see Figure 5). The former have n${}_{\rm gas}<100$ ${\rm cm}^{-3}$ and n${}_{\rm H2}<5$ ${\rm cm}^{-3}$, whereas the latter have ngas $\sim$ $100$ ${\rm cm}^{-3}$ and n${}_{\rm H2}\gtrsim 10$ ${\rm cm}^{-3}$ (see also Figs. 5 and 6 of Pallottini et al. (2017a) for details111111In particular, their analysis of the Minkowski function of the molecular density field (see Fig. 6 therein)). Considering that $f_{\rm H2}$ = $\rho_{\rm H2}/\rho_{\rm H,total}$ = $\rho_{\rm H2}/(\mu\rho_{\rm gas})$ = $m_{\rm H2}n_{\rm H2}/(m_{p}\,\mu\,n_{\rm gas})$ $\sim$ 2 $n_{\rm H2}/(0.7~{}n_{\rm gas})$, for the disk/interarm region $f_{\rm H2}$ $\sim$ 14%, while for the MGCs $f_{\rm H2}$ $\sim$ 28%. Note that for the lowest value of the H2 cuts, these fractions might be lower; however, removing regions selected with such cuts yields no qualitative change to the physical interpretation resulting from the analysis.

We investigate possible variations in the dynamics of the molecular structures of Althæa by adopting different values of $n_{\rm cut}$ to test the robustness of our results against the choice of density threshold. That is, how sensitively are the structure properties, and thus, the results presented in §5.2, dependent on the choice of density thresholds? Observationally, this effect would be mimicked by e.g., adopting different molecular gas tracers since they have different critical densities.

The sizes of MGCs are dependent on the choice of $n_{\rm cut}$ in the following ways. As mentioned in §5.1, the most massive cloud identified at low H2 gas density thresholds corresponds to the molecular disk of Althæa, which is broken down into multiple smaller MGCs at higher $n_{\rm cut}$121212That said, there are fewer MGCs at highest $n_{\rm cut}$ than the lowest $n_{\rm cut}$ since there are fewer cells in the simulations with correspondingly high H2 gas densities.. Second, excluding such structure, some MGCs are further broken down into smaller MGCs at higher $n_{\rm cut}$ (see top panels of Figure 8). At the resolution limit of our simulation, most MGCs studied in this work have $R\simeq$ 50 pc (see also Figure 4). The gas velocity dispersion $\sigma_{\rm gas}$ of MGCs, on the other hand, is rather insensitive to the actual value of $n_{\rm cut}$. This lack of variation is reassuring: our inference on the velocity dispersion of $z$ $\sim$ 6 MGCs in relation to those observed in nearby and $z$ $\sim$2 galaxies is not biased by our choice of $n_{\rm cut}$ in this work.

At the highest density cut, only the densest gas structures are identified131313In the simulation, the densest molecular regions reaches gas density of $n_{\rm gas}$ $\sim$ 103 ${\rm cm}^{-3}$, which is comparable to the critical density of low-$J$ CO lines (see e.g., Vallini et al. 2018); however, we emphasize that we do not resolve the internal structure of molecular clouds at the resolution of the simulation.. The virial parameter of these MGCs are, on average, lower than those identified at low $n_{\rm cut}$, and are thus, more unstable against collapse (see middle right panel of Figure 8).

6. Discussion

6.1. Synthetic vs Observed MGC Properties

In this section, we compare the observable quantities of MGCs obtained from the simulations (e.g., $\sigma_{\rm gas}$, $R$, $M_{\rm gas}$, $\Sigma_{\rm gas}$) with those of molecular gas structures from observations at lower redshift.

The velocity dispersion of MGCs in Althæa is similar to those observed in $z$ $\sim$ 2 spatially resolved studies of gas-rich star-forming galaxies, spanning a range of $\sigma_{\rm gas}\simeq$ 20$-$80 km s*-1* (Figure 7; see e.g., Swinbank et al. 2011). However, the sizes of MGCs are almost a factor of two smaller. This suggests that the sizes of high-$z$ molecular gas structures reported in the literature are probably limited by the spatial resolution. The velocity dispersions of some MGCs are also comparable to the most turbulent gas clouds observed in the inner Milky way and nearby gas-rich galaxies (e.g., M64; Oka et al. 2001; Rosolowsky & Blitz 2005; Heyer et al. 2009; Leroy et al. 2015), which lie along the locus of $\sigma_{\rm gas}\propto R^{0.56}$. Such high velocity dispersion in Althæa is dominated by non-thermal energy (Figure 6) and results from the injection of kinetic energy from recent star formation as Althæa is assembling its stellar mass (given the high Mach numbers of MGCs; §4.2). In fact, by the accreting stage (see Figure 1a) at $z\simeq$ 7.2, Althæa has assembled a stellar mass of $M_{\star}$ = 7.5$\times 10^{9}$ $M_{\odot}$. This contrasts with Krumholz et al. (2018), who show that numerical models and analytic arguments point to feedback generally not being able to produce velocity dispersions of more than 10–20 km s*-1*, with higher velocities likely being produced by Toomre instability driven radial inflow.

The total pressure141414Non-thermal pressure dominates the total pressure in these high density regions. of MGCs identified at the highest $n_{\rm cut}$ ranges from $P/k_{B}\simeq 10^{7-9}$ K cm*-3*, with a median of $10^{7.6}$  K cm*-3*. Such high pressures are comparable to those observed in local ULIRGs; however, the molecular clouds in these local galaxies are concentrated within their central regions (Downes & Solomon 1998; Sakamoto et al. 2008). In our simulated galaxy, the high pressure MGCs are found throughout the disk. This difference likely stems from the different physical mechanisms giving rise to the highly turbulent nature of these molecular structures. While MGCs in both local ULIRGs and Althæa may form via gravitational instability (see e.g., McKee & Ostriker 2007), MGCs in the local merger-driven ULIRGs are likely formed by shock compression and cloud-cloud collisions that funnel large amounts of gas from the progenitor galaxies towards the central region (Tan 2000; Wu et al. 2018; see §6.2 for Toomre-$Q$ analysis on Althæa). As shown by Kohandel et al. (2019) and Gallerani et al. (2018), the highly turbulent nature seen in the MGCs of Althæa result from extra-planar flows and higher velocity accretion/SN-driven outflows. The presence of extra-planar flows may also be the dominant mode for forming the highly supersonic massive MGCs observed in gas-rich star-forming galaxies at $z$ $\sim$ 2. For instance, Swinbank et al. (2011) report ISM pressure of $P/k_{B}\sim$107 K ${\rm cm}^{-3}$ and molecular gas mass of $M_{\rm gas}\simeq 10^{8-9}$ $M_{\odot}$ in the star-forming regions of a gas-rich starburst galaxy at $z$  = 2.3 based on spatially resolved CO line observations.

We do not find any major quantitative differences (see Figure 7) in the MGCs of Althæa with respect to those observed in the nearby Universe between its accreting and starburst phase, which are separated by $\simeq$300 Myr. However, MGCs in the starburst phase have lower $\alpha_{\rm vir}$ compared to the accreting phase. This is expected for star formation to proceed.

6.2. Toomre Parameter and Stability of MGCs

We define Toomre $Q$ parameters for the gas, stars, and the effective two-component $Q_{\rm eff}$ in §2.2. Maps of the corresponding Toomre $Q$ are shown in Figure 10 (see also Figure 2 for example of MGCs identified in this evolutionary stage of Althæa using different $n_{\rm cut}$).

Close resemblance of the $Q_{\rm eff}$ and $Q_{\star}$ maps indicates that contributions from the stellar component play an important role in governing the stability of the MGCs against perturbations. This can be understood since stars in Althæa dominate the central part of the galaxy in mass, so their gravitational potential provides a non-negligible contribution to the instability. Similarly, contribution from the thickness of the disk is important in Althæa since its disk is warped and has a scale height-to-radius ratio of $h/r_{\rm gal}$ $\sim$ 150 pc/1 kpc $\simeq$ 0.15. That is, some MGCs are found in regions of $Q_{\rm gas}$ $\sim$ 1, which is consistent with the expectation that they correspond to regions of high surface densities that are gravitationally unstable. Note, however, that when including the stabilizing effects due to the stellar potential of Althæa and the thickness of its disk (via $\sigma_{z}$, i.e., $Q_{\rm eff}$; see Equation 13), some of these MGCs are consistent with $Q_{\rm eff}>1$. This demonstrates the importance of accounting for stellar contribution and disk thickness when examining the stability of molecular gas structures. This consideration is especially relevant for the relatively evolved and enriched systems at high redshift that are preferentially being imaged at high resolution with ALMA. In the outer regions of Althæa, on the other hand, $Q_{\rm eff}$ resembles $Q_{\rm gas}$, with both $\log{Q_{\rm gas}}$ and $\log{Q_{\rm eff}}<-1$. Notably, these are the more gas-rich regions (see Figure 9).

The large virial parameters ($\alpha_{\rm vir}>$ 2) seen in some MGCs of Althæa can be understood by first noting that fragmentation can happen in regions of low $Q$ (if we only consider instability against axisymmetric perturbations), but further evolution and gas collapse depends on the equation of state of the gas. Such fragmentation is expected to take place at the critical scale length $\lambda_{\rm crit}<2\pi^{2}G\Sigma/\kappa^{2}$. Second, this fragmentation scale is greater than the typical size of GMCs. This could be interpreted as a result of instability setting the scales for fragmentation, but the truly star-forming regions correspond to the collapsing, denser, and cooler molecular structures that are on smaller scales. Thus, the high virial parameter found for some MGCs in Althæa indicates that they are not the collapsing structures and are found in regions of $Q_{\rm eff}>1$. On the other hand, MGCs in (the denser) regions with a lower Toomre $Q_{\rm eff}$ parameter are more unstable, and star formation may take place within its star-forming clumps and cores on smaller scales, where energy quickly dissipates.

7. Summary and Conclusions

Properties of the star-forming ISM of galaxies near the Epoch of Reionization are now within reach with ALMA. While it is possible to obtain sensitive and high fidelity imaging that reveals their gas and dust morphology on sub-kpc scales and even smaller, such observations remain challenging; logistically, data from multiple cycles pushing to increasing resolution and sensitivity are needed. As such, observational studies of MGCs on such scales are still missing in the literature. In this work, we aim to understand the origin and dynamical properties of MGCs in prototypical galaxies at the EoR in numerical simulations to provide a framework within which upcoming observations can be compared against to aid in the interpretation.

We study the dynamics of MGCs and their temporal evolution in Althæa, a $z$ $\sim$ 6 prototypical galaxy obtained from the state-of-the-art cosmological zoom-in simulations Serra, which include a chemical network to determine the formation of molecular hydrogen, heating and cooling of the ISM by UV radiation and metal lines, and detailed stellar feedback. We use a three-dimensional clump-finding algorithm to identify MGCs. We decompose the molecular structures into non-overlapping objects by using a set of H2 density contours ($n_{\rm cut}$) at multiple evolutionary stages, in particular focusing on the accreting phase, starburst phase, and quiescent phase. We extract properties such as mass, size, Mach number, velocity dispersion, gas surface density, and virial parameter ($M_{\rm gas},R,\mathcal{M},\sigma_{\rm gas},\Sigma_{\rm gas},\alpha_{\rm vir}$) for each MGC and perform a Toomre-$Q$ stability analysis on Althæa.

Excluding the main structure (disk of Althæa) identified, the typical mass and size of MGC is $M_{\rm gas}\simeq\,10^{6.5}$ $M_{\odot}$ and $R\simeq\,50$ pc, respectively. Similarly massive molecular structures have been observed in nearby star-forming and starburst galaxies (e.g., Keto et al. 2005; Donovan Meyer et al. 2013; Colombo et al. 2014; Leroy et al. 2015), and reported in idealized (no stellar feedback) isolated (no cosmological initial conditions) galaxy simulations done at higher resolution (e.g., $l_{\rm cell}\simeq 3$ pc in the 48 kpc box studied by Behrendt et al. 2016b). That said, the similar mass range found in MGCs of Althæa is reassuring — our results are not far off despite the resolution limit ($l_{\rm cell}\simeq$ 30 pc). On the other hand, the gas velocity dispersion of $\sigma_{\rm gas}\simeq$ 20-10 km s*-1* is rather insensitive to the adopted density threshold.

Velocity dispersion and gas surface density of MGCs of Althæa are systematically higher than Milky Way clouds regardless of the density threshold $n_{\rm cut}$ adopted. These MGCs are, in fact, highly supersonic, with high Mach number of $\bar{\mathcal{M}}\simeq 6$. Their velocity dispersions are comparable to those observed in $z$ $\sim$ 2 starburst galaxies. A comparison between the bulk and non-thermal velocity dispersions of MGCs indicates that MGCs are supported by turbulence motions.

High pressure MGCs are found throughout the disk of Althæa, with a median pressure of $\bar{P}\simeq 10^{7.6}$ K ${\rm cm}^{-3}$. This is in contrast to the local ULIRGs, where such high pressure has only been observed in molecular clouds concentrated in their central regions (Downes & Solomon 1998; Sakamoto et al. 2008). In ULIRGs, MGCs are formed by shock compression and cloud-cloud collisions that funnel large amounts of gas from progenitor galaxies towards the central region (Tan 2000; Wu et al. 2018). On the other hand, the highly turbulent MGCs of Althæa results from extra-planar flows and high velocity accretion/SN-driven outflows. The presence of such extra-planar flows may also be the dominant mode for forming the highly supersonic massive MGCs observed in gas-rich star-forming galaxies at $z$  $\sim$ 2, for instance, those reported in Swinbank et al. (2011). The variation in $\sigma_{\rm gas}^{2}/R$ with $\Sigma$ seen in this work is in line with previous studies, suggesting that “Larson’s third relation” is an artifact arising from observational bias artifact of the limited range of column densities a specific molecular line tracer is sensitive to (Ballesteros-Paredes & Mac Low 2002b; Ballesteros-Paredes et al. 2011b). This would imply column density is not a fundamental properties of molecular clouds and not all molecular clouds are virialized.

We perform virial analysis, as motivated by observations, to assess the stability of MGCs. On average, the virial parameter of MGCs in the starburst phase of Althæa is lower than in the accreting phase, as expected for star formation. Similarly, the virial parameter of MGCs identified at the highest density thresholds are lower than those identified at lower density thresholds. Close resemblance of $Q_{\rm eff}$ and $Q_{\rm star}$ maps indicates that contribution from the stellar component plays an important role in governing the stability of the MGCs against axisymmetric perturbations, especially in the central part of Althæa. Similarly, stabilizing effect due to the thickness of its disk is also non-negligible. This illustrates the importance of accounting for stellar contribution and disk thickness when examining the stability of molecular gas structures, especially in relatively evolved and enriched systems at high redshift that are preferentially being observed now. Star formation is expected to take place within its star-forming clumps and cores on smaller scales, where energy quickly dissipates. This is consistent with the notion that collapsing structures result from gravitational instability occurring within globally stable structures, which are supported by turbulence and rotation on large scale.

This also implies that observations with spatial resolution better than $\simeq$40 pc are needed to examine the truly star-forming structures (cores), and thus, star formation in the first galaxies. Such resolution is in principle within reach using ALMA, for which 0$\farcs$01 resolution images can be attained using the highest frequency bands. With the Next Generation VLA151515http://library.nrao.edu/public/memos/ngvla/NGVLA_21.pdf, which is approximately ten times more sensitive than ALMA, direct observations of molecular gas clouds in galaxies at the EoR on cloud scales will be possible via the higher-$J$ rotational transitions of CO in the 2030s. In the meantime, cosmological zoom-in simulations, such as Serra, while inherently limited in galaxy statistics and depend on the sub-grid models adopted, serve as a useful tool for examining and making predictions on the morphology and dynamics of the molecular ISM of the first galaxies.

Acknowledgements

We thank Jens Kauffmann, Thushara Pillai, and Mark Swinbank for sharing their data, Caitlin Casey for useful discussions, and the referee for providing constructive comments on the manuscript. TKDL gratefully acknowledges support by the NSF through award SOSPA4-009 from the NRAO and support from the Simons Foundation. AF acknowledges support from the ERC Advanced Grant INTERSTELLAR H2020/740120. M-MML acknowledges partial support by the NSF through award AST18-15461, and hospitality from the Insitute for Theoretical Astrophysics at the University of Heidleberg. This work is based on a project developed at the Kavli Summer Program in Astrophysics (KSPA) held at the Center for Computational Astrophysics of the Flatiron Institute in 2018. The program was co-funded by the Kavli Foundation and the Simons Foundation. We thank the KSPA Scientific and Local Organizing Committees, and the program founder, Pascale Garaud, for supporting the genesis of this work. We also thank the Center for Cosmology and Particle Physics at the New York University for their hospitality in hosting us after the steam pipe explosion in NYC during the KSPA. The Flatiron Institute is funded by the Simons Foundation. This research has made use of NASA’s Astrophysics Data System Bibliographic Services.