#  The Supersonic Project: The Early Evolutionary Path of SIGOs

William Lake Department of Physics and Astronomy, UCLA, Los Angeles, CA 90095
Mani L. Bhaumik Institute for Theoretical Physics, Department of Physics and
Astronomy, UCLA, Los Angeles, CA 90095, USA
Smadar Naoz Department of Physics and Astronomy, UCLA, Los Angeles, CA 90095
Mani L. Bhaumik Institute for Theoretical Physics, Department of Physics and
Astronomy, UCLA, Los Angeles, CA 90095, USA
Blakesley Burkhart Department of Physics and Astronomy, Rutgers, The State
University of New Jersey, 136 Frelinghuysen Rd, Piscataway, NJ 08854, USA
Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue,
New York, NY 10010, USA
Federico Marinacci Department of Physics & Astronomy “Augusto Righi”,
University of Bologna, via Gobetti 93/2, 40129 Bologna, Italy
Mark Vogelsberger Department of Physics and Kavli Institute for Astrophysics
and Space Research, Massachusetts Institute of Technology, Cambridge, MA
02139, USA
Gen Chiaki Astronomical Institute, Tohoku University, 6-3, Aramaki, Aoba-ku,
Sendai, Miyagi 980-8578, Japan  Yeou S. Chiou Department of Physics and
Astronomy, UCLA, Los Angeles, CA 90095  Mani L. Bhaumik Institute for
Theoretical Physics, Department of Physics and Astronomy, UCLA, Los Angeles,
CA 90095, USA
Naoki Yoshida Kavli Institute for the Physics and Mathematics of the Universe
(WPI), UT Institute for Advanced Study, The University of Tokyo, Kashiwa,
Chiba 277-8583, Japan  Department of Physics, The University of Tokyo, 7-3-1
Hongo, Bunkyo, Tokyo 113-0033, Japan  Yurina Nakazato Department of Physics,
The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan  Claire E.
Williams Department of Physics and Astronomy, UCLA, Los Angeles, CA 90095
Mani L. Bhaumik Institute for Theoretical Physics, Department of Physics and
Astronomy, UCLA, Los Angeles, CA 90095, USA

(Received ApJ)

###### Abstract

Supersonically Induced Gas Objects (SIGOs) are a class of early Universe
objects that have gained attention as a potential formation route for globular
clusters. SIGOs have recently begun to be studied in the context of molecular
hydrogen cooling, which is key to characterizing their structure and
evolution. Studying the population-level properties of SIGOs with molecular
cooling is important for understanding their potential for collapse and star
formation, and for addressing whether SIGOs can survive to the present epoch.
Here, we investigate the evolution of SIGOs before they form stars, using a
combination of numerical and analytical analysis. We study timescales
important to the evolution of SIGOs at a population level in the presence of
molecular cooling. Revising the previous formulation for the critical density
of collapse for SIGOs allows us to show that their prolateness tends to act as
an inhibiting factor to collapse. We find that simulated SIGOs are limited by
artificial two-body relaxation effects that tend to disperse them. We expect
that SIGOs in nature will be longer-lived compared to our simulations.
Further, the fall-back timescale on which SIGOs fall into nearby dark matter
halos, potentially producing a globular-cluster-like system, is frequently
longer than their cooling timescale and the collapse timescale on which they
shrink through gravity. Therefore, some SIGOs have time to cool and collapse
outside of halos despite initially failing to exceed the critical density.
From this analysis we conclude that SIGOs should form stars outside of halos
in non-negligible stream velocity patches in the Universe.

Globular star clusters — High-redshift galaxies — Giant molecular clouds —
Star formation — Galactic and extragalactic astronomy

##  1 Introduction

The standard model of structure formation, invoking dark energy and dark
matter that dominate the Universe’s energy density budget, has had great
success in explaining a wide variety of observations. The $\Lambda$CDM model
successfully explains the anisotropies in the Universe’s radiation field and
the large-scale distribution of galaxies, as well as more generally explaining
properties of the Universe’s structure on scales larger than $100$ Mpc
(Springel, 2005; Vogelsberger et al., 2014a, b; Schaye et al., 2015). This
model predicts that the Universe’s structure formed hierarchically, with very
early primordial baryon overdensities at $z\lessapprox 30$ collapsing to form
larger objects, which eventually form galaxies and other structures.

These early baryon overdensities evolved in turn from interactions with
existing DM overdensities following recombination. Because baryon objects had
been inhibited from collapse prior to recombination by the photon field, DM
overdensities at the time of recombination were about $10^{5}$ times larger
than baryon overdensities (e.g., Naoz & Barkana, 2005a). Because of this
imbalance, the formation and growth of baryon overdensities at this time was
driven by these large existing DM overdensities.

However, this process was complicated by the fact that there was a significant
relative velocity between baryons and DM following recombination
(Tseliakhovich & Hirata, 2010). This relative velocity had an rms value at
recombination of $30$ km s-1 (which was about $5$ times the average speed of
sound in the Universe at this time) and was coherent on few-Mpc scales.
Because of its coherence on these scales, it is also known as the stream
velocity.

Because the stream velocity was so highly supersonic, it induced significant
spatial offsets between baryon overdensities and their parent DM
overdensities, forming collapsed baryon objects outside the virial radii of
their parent DM halos (Naoz & Narayan, 2014). These collapsed baryon objects,
known as Supersonically Induced Gas Objects or SIGOs, have been proposed as a
progenitor of some modern-day globular clusters (GCs), and have been shown to
have similar properties to GCs (e.g., Chiou et al., 2018, 2019, 2021; Lake et
al., 2021; Nakazato et al., 2022). However, SIGOs remain poorly understood.
Early simulations of SIGOs, aimed at confirming their existence and putting
basic constraints on their properties (e.g., Popa et al., 2016; Chiou et al.,
2018, 2019, 2021; Lake et al., 2021), included only adiabatic and sometimes
atomic hydrogen cooling. However, molecular hydrogen cooling has important
effects on the structure of SIGOs (Nakazato et al., 2022), and on the
abundance of gas objects in general in the early Universe (e.g. Glover, 2013;
Schauer et al., 2021; Nakazato et al., 2022; Conaboy et al., 2022). For
example, SIGOs in molecular cooling simulations tend to be far more
filamentary than those in atomic cooling simulations (Nakazato et al., 2022),
potentially for reasons analogous to the Zel’dovich approximation, which
states that collapse will occur along successive axes, starting with the
shortest (Zel’dovich, 1970). H2 cooling also lowers the temperature in these
SIGOs to $\sim 200$ K, which lowers the Jeans mass to about $10^{3}$ M⊙,
potentially leading to star formation in SIGOs (Nakazato et al., 2022, Lake et
al. in prep).

In order to understand and contextualize the properties of SIGOs with
molecular hydrogen cooling, it is useful to have an analogous class of objects
to compare to, which share many of the properties of SIGOs. A natural class of
objects to compare to here are giant molecular clouds (GMCs). Like SIGOs, GMCs
form in environments shaped by supersonic turbulence (e.g., Larson, 1981;
Padoan & Nordlund, 1999; Mac Low & Klessen, 2004; Krumholz & McKee, 2005;
McKee & Ostriker, 2007; Krumholz & Tan, 2007; Bergin & Tafalla, 2007; Padoan &
Nordlund, 2011; Burkhart & Lazarian, 2012; Semenov et al., 2016; Mocz et al.,
2017; Burkhart, 2018, 2021; Appel et al., 2022). The energy from this
turbulence cascades from the driving scale of the turbulence ${\rm L}_{\rm
Drive}$ down to the sonic scale $\lambda_{s}$ which marks the boundary between
supersonic and subsonic turbulence:

$\lambda_{s}=\frac{\rm{L}_{\rm Drive}}{\mathcal{M}^{2}}\ ,$ (1)

where $\mathcal{M}$ is the mach number of the turbulent flow on the driving
scale, and where we have assumed Larson’s law (Larson, 1981). This energy
cascade is important in both cases, because it creates high density peaks in
the turbulent medium on scales comparable to the sonic scale (e.g., Krumholz &
McKee, 2005), which is key to understanding both objects’ structure.
Ultimately, this process leads to a critical density for star formation that
can be expressed in the form

$\rho_{\rm crit,GMC}=\frac{\pi c_{s}^{2}\mathcal{M}^{4}}{GL_{\rm Drive}^{2}}\
,$ (2)

where $c_{s}$ is the speed of sound and $G$ is the gravitational constant.

GMCs are also interesting as a point of comparison for SIGOs because they may
be a major formation channel for globular clusters at early times (e.g.,
$z\sim 6$ Harris & Pudritz, 1994; Elmegreen & Efremov, 1997; Kravtsov &
Gnedin, 2005; Shapiro et al., 2010; Grudić et al., 2022). GMCs form with a
wide range of masses and densities, and it has been shown that the high-
density, high-mass end of this formation channel may be capable of efficiently
forming globular clusters. This formation process is supported by observations
of the merging Antennae system and its population of massive young clusters
(e.g., Whitmore & Schweizer, 1995; Whitmore et al., 1999). Therefore, a better
understanding of the comparison between SIGOs and GMCs may also be important
to understanding the link between SIGOs and globular clusters.

Motivated by these factors, in the present paper we present a suite of
simulations using the cosmological simulation code AREPO, including primordial
chemistry and accounting for the effects of molecular hydrogen cooling. We
constrain SIGOs’ structural and kinematic properties, including the effects of
the stream velocity and supersonic turbulence on the size and density of
SIGOs, and compare these properties to those of GMCs to better contextualize
them.

This paper is organized as follows: in Section 2 we discuss the setups of the
simulations used, as well as modifications made to definitions of SIGOs used
in prior works. In Section 3 we expand upon the importance of molecular
hydrogen cooling to SIGOs’ properties. In Section 4, we discuss the
similarities and differences between SIGOs and GMCs, and revisit previous
density-driven analyses of SIGO evolution. In Section 5, we overview the
evolution of SIGOs from a timescale perspective, primarily discussing the
potential for SIGOs to collapse through the lens of molecular hydrogen
cooling. Lastly, in Section 6 we overview future avenues of exploration
towards understanding star formation in SIGOs and summarize our results.

For this work, we have assumed a $\Lambda$CDM cosmology with
$\Omega_{\rm\Lambda}=0.73$, $\Omega_{\rm M}=0.27$, $\Omega_{\rm B}=0.044$,
$\sigma_{8}=1.7$, and $h=0.71$.

##  2 Numerical Setup and Object Classification

Using the cosmological and hydrodynamic simulation code AREPO, we present a
suite of $4$ simulations demonstrating the effect of molecular hydrogen
cooling on the properties of SIGOs. The main parameters of these simulations
are listed in Table 1. Here, we indicate runs with molecular hydrogen cooling
using “H$2$”, and indicate runs without molecular hydrogen cooling (simply
using adiabatic and atomic hydrogen cooling) with “H”. Runs without the stream
velocity are labelled with “$0$v,” whereas runs with a $2\sigma$ stream
velocity (a 2 $v_{\rm rms}$ relative velocity between baryons and dark matter)
are labelled as $2$v. The stream velocity in the latter runs is implemented as
a uniform boost to baryon velocities in the x direction in the initial
conditions–at the initial redshift of z$=200$, this is a boost of $11.8$ km
s-1.

Our initial conditions were generated using transfer functions calculated
using a modified version of CMBFAST (Seljak & Zaldarriaga, 1996) taking into
account the first-order correction of scale dependent temperature fluctuations
(e.g., Naoz & Barkana, 2005b; Naoz & Narayan, 2013) and second-order
corrections to the equations presented in Tseliakhovich & Hirata (2010),
describing the evolution of the stream velocity. We use two transfer
functions, one each for the baryons and DM, as the evolution of the gas
fraction strongly depends on the initial conditions of the baryons (e.g., Naoz
et al., 2009, 2011, 2012, 2013; Park et al., 2020).

Our simulations use $512^{3}$ DM particles with a mass of $1.9\times 10^{3}$
M⊙ and $512^{3}$ Voronoi mesh cells corresponding to a gas mass of $360$ M⊙,
in a box $2$ comoving Mpc on a side. Following Chiou et al. (2019, 2021); Lake
et al. (2021); Nakazato et al. (2022), we choose $\sigma_{8}=1.7$,
representing a rare, overdense region where structure forms early in a large
volume, in order to increase the number of gas objects in our simulation,
allowing increased statistical power. The simulations begin at z$=200$, and
they run to z$=20$.

In runs $0$vH$2$ and $2$vH$2$, using the chemistry and cooling library GRACKLE
(Smith et al., 2017; Chiaki & Wise, 2019), we track nonequilibrium chemical
reactions and their associated radiative cooling explicitly in the gas. This
run includes H2 and HD molecular cooling, as well as chemistry for 15
primordial species: e-, H, H+, He, He+, He++, H-, H2, H${}_{2}^{+}$, D, D+,
HD, HeH+, D-, and HD+. The radiative cooling rate of H2 follows both
rotational and vibrational transitions (Chiaki & Wise, 2019).

Figure 1: Comparison of star formation in SIGOs with and without molecular
cooling at $z=20$. We show the gas density field in our simulation box for
simulation $2$vH$2$ (molecular cooling, right panel), compared to simulation
$2$vH (atomic cooling, left panel). SIGOs are marked with Xs. Note here that
SIGOs trace gas and halo abundances on these scales. Note also that molecular
hydrogen cooling dramatically increases the abundance of SIGOs.

For comparison, in runs $0$vH and $2$vH, we consider only atomic hydrogen
cooling, following the species e-, H, H+, He, He+, and He++ in equilibrium
with a spatially uniform, redshift-dependent photoionizing background, as
described in Vogelsberger et al. (2013).

We use the object classifications laid out in Chiou et al. (2018). We identify
DM halos with an FOF (friends-of-friends) algorithm that uses a linking length
of $20\%$ of the mean particle separation of the DM component, which is about
$780$ comoving pc. Assuming sphericity for simplicity, we use this to
calculate positions and radii of all DM halos in the simulation output (though
note that DM halos at this time resemble triaxial ellipsoids, e.g. Sheth et
al., 2001; Lithwick & Dalal, 2011; Vogelsberger & White, 2011; Schneider et
al., 2012; Vogelsberger et al., 2020).

We then run the FOF finder on the gas component, allowing us to identify gas-
primary objects that contain over 100 gas cells. Especially when molecular
hydrogen is included, SIGOs at these redshifts are distinctly filamentary,
similar to in the Zel’dovich approximation (Zel’dovich, 1970). We cannot,
therefore, assume sphericity when determining these objects’ properties, as we
did with DM halos. To address this issue, we follow the procedure introduced
in Popa et al. (2016). In particular, we fit these gas-primary objects to
ellipsoids determined as the smallest ellipsoidal surface that surrounds all
of the constituent gas cells. We then reduce the major axis of this ellipsoid
by $5$% until either the ratio of the axes lengths of the tightened ellipsoid
to that of the original ellipsoid is greater than the fraction of gas cells
remaining inside the tightened ellipsoid, or until 20% of their particles have
been removed. Finally, in order to distinguish SIGOs from other classes of gas
objects, we require that the center of mass of the gas objects be located
outside the virial radius of all nearby halos, and that the SIGOs have a
baryon fraction above $60$%111This condition was introduced as $40$% in
initial work by Chiou et al. (2018), and modified by Nakazato et al. (2022) to
account for false identifications of SIGOs arising from molecular cooling, as
discussed in Section 3.

Run | $v_{\mathrm{bc}}$ | H2 Cooling
---|---|---
0vH2 | 0 | Yes
0vH | 0 | No
2vH2 | 2$\sigma$ | Yes
2vH | 2$\sigma$ | No
Table 1: Simulation Parameters

##  3 The Importance of Molecular Cooling to SIGOs’ Properties

Molecular hydrogen cooling plays a vital role in the formation of the first
stars (e.g., Saslaw & Zipoy, 1967; Haiman et al., 1996; Abel et al., 1998;
Yoshida et al., 2003; Stacy et al., 2011; Greif et al., 2011; Yoshida et al.,
2008; Bromm, 2013; Glover, 2013; Hummel et al., 2016; Nakazato et al., 2022).
Previous studies of SIGOs have argued that adiabatic and atomic cooling may be
a sufficient approximation to understand SIGOs’ formation and morphology
(e.g., Popa et al., 2016; Chiou et al., 2018, 2021; Lake et al., 2021). This
is because the main factor in the formation of SIGOs is the phase shift of
gravitational fluctuations between DM and baryons (Naoz & Narayan, 2014).
However, it was recently shown that molecular cooling may be able to increase
the efficiency of SIGO formation, as well as playing an important role in
their density (Schauer et al., 2021; Nakazato et al., 2022).

Quantifying the population-level differences in SIGOs’ properties with
molecular cooling compared to atomic cooling can help to illuminate the
physical processes important to their later evolution. Figure 1 shows a
comparison of SIGO abundances with and without molecular hydrogen cooling
(runs $2$vH and $2$vH$2$). SIGOs in this figure are marked with white Xs,
displayed against a backdrop of the gas density in the full simulation box.
With molecular hydrogen cooling, SIGOs are dramatically more efficient at
condensing from overdensities in the gas. Thus, we find $84$ SIGOs at z$=20$
in run $2$vH$2$, compared to $27$ in run $2$vH. With molecular cooling, SIGOs
also reach higher densities, due to their increased cooling efficiencies. In
run $2$vH$2$, the most dense SIGO at z$=20$ had an overall gas density of
$39.1$ cm-3, compared to a maximum gas density in any SIGO of $2.5$ cm-3 in
run $2$vH. Therefore, SIGOs are more abundant and denser by redshift $20$ with
molecular cooling.

Figure 2: Changes in SIGO abundance with redshift. Plotted here is the
evolution of SIGO abundances with redshift in Run $2$vH$2$. For comparison,
SIGO abundances from Run $2$vH are plotted in red. Note that the time
resolution of Run $\mathbf{2}$vH$\mathbf{2}$ is higher than that of Run
$\mathbf{2}$vH, for the purposes of this study.

The increase in SIGO abundances through molecular cooling seen in Figure 1 can
be seen as a function of redshift in Figure 2. This figure shows the time
evolution of SIGO abundances with molecular cooling, as well as showing SIGO
abundances at integer redshifts in Run $2$vH for comparison. As one can see,
the abundance of SIGOs increases with time, and consistent with Nakazato et
al. (2022), abundances are generally enhanced through molecular cooling,
through the process outlined above. This process is most important at later
redshifts: as one can see in the figure, the abundance of SIGOs in Run
$2$vH$2$ begins to significantly diverge from the abundance without molecular
cooling after $z\sim 25$ in our simulations.

Molecular cooling affects the structure around SIGOs as well: SIGOs form
embedded in gas filaments, and molecular hydrogen cooling permits these
filaments to condense much more efficiently than does atomic hydrogen cooling.
As mentioned in Section 2, the gas fraction cutoff for identification of SIGOs
with our FOF algorithm had to be revised upwards when working with molecular
hydrogen cooling compared to prior work, because gas filaments condense so
efficiently that they sometimes appear to structure-finding algorithms (even
in Run $0$vH$2$ where no SIGOs should form) as collapsed baryon objects
outside of halos. Because SIGOs are embedded in these gas filaments, this
process is inextricably linked to their formation: the mass and particularly
the density of cooled gas around SIGOs in these filaments increases in Run
$2$vH$2$ compared to Run $2$vH, potentially allowing SIGOs to accrete external
gas and grow further with molecular cooling.

Figure 3: The velocity field of a typical SIGO (in red) with molecular
cooling. This SIGO is associated with a nearby parent DM halo, and there is a
larger DM halo at the bottom right of the field. Here, we see that the
velocity dispersion inside the SIGO is quite small compared to that outside of
it. The sonic scale here is larger than the scale of the SIGO, so turbulent
flow plays a small role in the SIGO’s potential further collapse. As can be
seen the SIGO is embedded in a shock front, which has formed on a scale where
the Mach number is about unity. See Figure 4.

##  4 SIGOs as GMC analogues

As mentioned above (§1), at face value, GMCs and SIGOs seem to share several
key properties. At the most fundamental level, the impact of molecular
hydrogen is critical for understanding the evolution of both classes of
objects, which also have similar masses ($\sim 10^{6}$ M⊙) and scales ($\sim
100$ pc). Additionally, these classes have similar substructures, with SIGOs
often exhibiting local density fluctuations such as cores and filaments in
simulations, similar to those in GMCs (Mac Low & Klessen, 2004; Krumholz &
McKee, 2005; Mocz & Burkhart, 2018). In addition, the formation of both
classes of object is dominated by supersonic turbulence from their environment
(Burkhart et al., 2009; Padoan & Nordlund, 2011; Burkhart & Lazarian, 2012;
Semenov et al., 2016; Mocz et al., 2017; Burkhart, 2018; Chiou et al., 2019;
Burkhart, 2021; Appel et al., 2022). In considering star formation within both
classes of objects, it is often impossible to consider only the Jeans mass–
instead, one must use a critical density that incorporates the impact of
turbulence.

To illustrate the nature of the turbulence around SIGOs, in Figure 3, we see
the velocity field around one SIGO from Run $2$vH$2$, including that of its
shock front. Mach numbers are labelled based on the local temperature and the
velocity relative to the center of mass of the SIGO. The SIGO itself here has
a maximum radius of $118$ pc, which is slightly smaller than the sonic scale.
The volume-averaged Mach number within the SIGO is $0.97$. One way to think
about this SIGO is as a density fluctuation induced by this supersonic
turbulence: the strongest density perturbations from the supersonic motions
take place at the sonic scale, appearing as gas overdensities whose positions
are correlated with those of nearby dark matter overdensities. This creates
gas objects (SIGOs) offset from their parent halos, with sizes that are by
necessity comparable to the sonic scale. Once this turbulence-induced
structure formation occurs, molecular cooling helps to cool the gas within the
SIGO to temperatures of order $200$ K (Nakazato et al., 2022), which enables
further collapse. Collapse of the SIGO reduces its axis length from the sonic
scale, so a typical SIGO will have a maximum axis length that is comparable to
but smaller than the sonic scale, unless it either accretes additional gas
from its surroundings (permitting it to grow while accreting mass) or
collapses.

Figure 4 shows a more general comparison of the sonic scales and the maximum
axis lengths of the SIGOs in Run $2$vH$2$. In the top panel, each faded line
shows the velocity dispersion within a SIGO as a function of its radius, and
in the bottom panel each faded line shows the mach number dispersion as a
function of radius. The black line is the average of all SIGOs’ dispersions in
log space. The spatial distributions of each are similar, because SIGOs are
close to isothermal: their sound speeds don’t change much on the edges of the
SIGOs compared to the centers. Typical temperatures at the edge of a SIGO are
of order $10\%$ higher than those in the central regions. Therefore, their
mach dispersions roughly follow their velocity dispersions.

Figure 4: SIGOs’ linewidth-size relation:  Panel (a) shows the linewidth-size
relation for SIGOs, computed over a sphere centered on the center of mass of
each SIGO, with increasing radius (scaled for each SIGO by the length of the
longest axis of that SIGO). The translucent lines represent individual SIGOs,
and the black opaque line represents an average of all SIGOs in log space.
Panel (b) shows the mach number dispersion, showing that most SIGOs are
somewhat smaller than the sonic scale. The vertical red line in each panel
displays the length of the longest axis of the SIGO, and the horizontal red
line in panel b indicates Mach 1.

A typical SIGO has a volume-averaged velocity dispersion of about Mach 0.8, or
about $3$ km/s (with fairly high variance) on the scale of its largest axis.
Put another way, a typical SIGO has a sonic scale that is about $1.6$ times
larger than its longest axis. This assumes a driving scale for the turbulence
on the order of the maximum axis length of the SIGO, as discussed by Chiou et
al. (2019). Because the sonic scale is comparable to the size of a SIGO (and
in some cases is smaller than the SIGO), it is important to account for this
turbulence when considering the ability of a SIGO to undergo gravitational
collapse (see also Hirano et al. (2017) for similar effect in non-SIGOs).
Chiou et al. (2019) accounted for this effect by adapting a critical density
(previously used in GMCs by e.g., Krumholz & McKee, 2005), defined by Equation
(2). This critical density defines the scale on which turbulence-induced
fluctuations collapse through Jeans collapse, in effect equating the Jeans
scale to the sonic scale. However, the analysis in the paper treated SIGOs as
spherical objects with uniform density $\rho_{\rm SIGO}$ within a radius
$R_{\rm max}$. SIGOs are not spherical (particularly with molecular hydrogen
cooling accounted for), and have a typical shortest axis length that is of
order $10$ times smaller than their longest axis in Run $2$vH$2$; therefore,
this overestimates SIGOs’ masses. An updated treatment for the critical
density of collapse for SIGOs, then, equates the mass of the SIGO to the Jeans
mass, giving a new definition of the critical density of a SIGO for collapse
$\rho_{\rm crit}$:

$\rho_{\rm crit}\approx\frac{\pi c_{s}^{2}\mathcal{M}^{4}}{4GR_{\rm min}R_{\rm
int}},$ (3)

where Rint is the second-largest principal axis of the SIGO and where we have
assumed $L_{\rm drive}\approx 2\times R_{\rm max}$.

At z$=20$, no SIGOs in run $2$vH exceed this critical density. This result
contrasts with that of Chiou et al. (2019), which found that SIGOs can
collapse through only atomic hydrogen cooling, because of the additional
factor of the prolateness of the SIGOs considered in this analysis. With
atomic cooling alone, the cooling timescale for these SIGOs is long, so it is
unlikely that SIGOs could form stars outside of halos without considering
molecular hydrogen cooling. However, when considering molecular hydrogen
cooling, looking only at these SIGOs’ densities at early redshifts does not
paint the whole picture of SIGOs’ ability to form stars outside of halos.

While the critical density and the Jeans density are both important to SIGOs,
their location outside nearby halos allows their evolution to be
characteristically slow. At z=$20$, about $50$ Myr or less after SIGOs begin
to form in meaningful quantities, SIGOs tend to be fairly isolated objects.
They are outside of the immediate vicinity of the DM halos that birthed them,
and have characteristically long fall-back times to their nearest halo.
Because of this, even though their cooling times are long (Schauer et al.,
2021), SIGOs can have shorter cooling times than fall-back times, allowing
them to collapse in spite of their early Jeans stability (though typically
later than z=$20$, because these timescales tend to be $\gtrapprox 50$ Myr).
Subsequently, they may be able to accrete onto halos on the fall-back
timescale, and potentially survive as identifiable clusters due to their
compact nature and existing population of stars. In order for this to happen,
they must have cooling and collapse timescales that are shorter than their
fall-back timescale to nearby halos, as discussed in Section 5. At any given
time, then, SIGOs may be destined for collapse through cooling despite not
exceeding the critical density for collapse. Therefore, to study the evolution
of SIGOs, we argue that a timescale analysis is more appropriate than a
density analysis, as follows in Section 5.

Another of the primary characteristics of GMC populations is their power law
mass spectrum. GMCs in the inner disk of the Milky Way follow a power law mass
spectrum given by dN$\propto$MγdM with $\gamma\sim-1.5$. In contrast, the mass
spectra of GMCs in the outer Milky Way ($\gamma\sim-2.1\pm 0.2$) and M$33$
($\gamma\sim-2.9\pm 0.4$) are notably steeper (e.g. Rosolowsky, 2005).
Characterizing this mass spectrum is vital to understanding both how the
clouds formed, as well as their overall contribution to large-scale star
formation. Establishing how the mass spectrum of SIGOs compares to these power
laws is useful not only for furthering theoretical work with SIGOs, but also
for contextualizing similarities between SIGOs and these different populations
of GMCs.

Figure 5 shows the mass spectrum of SIGOs in Run $2$vH$2$. The blue histogram
bins show the probability density of finding a given SIGO in the labelled mass
bin. The green line shows a gaussian kernel density estimate of this spectrum,
with a bandwidth of $.15$, in log space. However, immediately translating this
estimate of the mass spectrum to a power law is complicated by the artificial
$2$-body relaxation present in the simulation (see Section 5.4 for more
details). At low masses, this mass spectrum is distorted, or truncated, by the
limited resolution of our simulation. Therefore, we only use mass points above
MSIGO = $10^{5}$ M⊙ to construct a best-fit power spectrum. This cutoff
reflects SIGOs whose artificial $2$-body relaxation timescales are longer than
the age of the oldest SIGO in the simulation, which are therefore less
affected by relaxation. We introduce an upper bound mass cutoff of $10^{6}$
M⊙, to avoid over-fitting to the SIGO with a mass of $\sim 10^{7}$ M⊙, because
we do not have enough data points to determine whether or not the spectral
index varies over the range $10^{6}$ M${}_{\odot}\leq$ MSIGO $\leq 10^{7}$ M⊙.
The best-fit power law is expected to diverge from the data in this range due
to Poisson fluctuations.

Figure 5: The Mass Spectrum of SIGOs:  Here we show the mass spectrum of SIGOs
from run $2$vH$2$ (the blue histogram, plotted as a probability density in log
space), as well as showing (in orange) a best-fit power law mass spectrum for
high-mass SIGOs. Overlaid in green is a gaussian kernel density estimate to
the mass spectrum. The mass spectrum of SIGOs in this run is consistent with a
power law index of $-2.3$. Figure 6: Timescales of SIGOs: here we show a
number of important timescales to the evolution of SIGOs at z$=20$ in Run
$2$vH$2$. Black dots mark the age and baryon mass of SIGOs found at this
redshift. The green region and lines in these figures show the range of fall-
back timescales to the nearest halo, as a function of mass. The dashed line
indicates the same for a SIGO with median properties. The blue region shows
the same, but for the range of cooling timescales, and the gray region shows
the timescale for gravitational collapse of a SIGO. The black line shows the
relaxation timescale of SIGOs in the simulation at the resolution of the
simulation (much shorter than in the real Universe). SIGOs with shorter
cooling and collapse timescales than their fall-back timescales are marked in
red as having the potential to form stars outside of a DM halo. As depicted,
the main limitation of SIGOs’ growth is the numerical evaporation process
(i.e., 2-body relaxation, see text). Therefore, in our adopted resolution high
mass SIGOs can collapse and form stars, while in the Universe we expect that
more SIGOs will form stars.

With these conditions in place, the orange line in Figure 5 shows a best-fit
power law to our kernel density estimate with $10^{5}$ M${}_{\odot}\leq$ MSIGO
$\leq 10^{6}$ M⊙. We find a best-fit mass spectrum index for SIGOs in this run
of $\gamma\sim-2.3$, which is, interestingly, consistent with the mass
spectrum of GMCs in the outer Milky Way, and is broadly consistent with the
range of mass spectra seen in GMC populations.

##  5 Timescale Analysis of the Evolution of SIGOs

The similarity between SIGOs and GMCs is limited when considering the main
challenge for SIGOs: their location outside of a DM halo, and eventual fall
back into their parent halo. In considering this challenge, and as a starting
point for our timescale analysis, we identify the following evolutionary
channels for SIGOs in a simulation (see Figure 9 for the relevant processes).

  * •

Gravitational Collapse: SIGOs can undergo gravitational collapse in their
overdense cores, or on the scale of the SIGO. We expand on this physical
process in §5.1.

  * •

Cooling: As mentioned above, molecular hydrogen cooling lowers the temperature
of the gas, reducing the gas pressure and assisting with collapse. See §5.2
for more details.

  * •

Fall-back Into Halos: SIGOs are inherently found near DM halos (Naoz &
Narayan, 2014). Therefore, SIGOs are likely to eventually fall into a DM halo.
Accretion is the most likely final state of even potential star-forming SIGOs,
and could lead to the formation of accreted clusters. See §5.3 for an estimate
of this timescale.

  * •

$2$-body Relaxation: Gas in AREPO numerical simulations has an associated mass
that depends on the simulation resolution. This gas interacts gravitationally,
causing it to undergo artificial two-body relaxation processes. These
artificial interactions result in the evaporation of gas structures on a
characteristic timescale. We estimate the timescale of this process in §5.4.

  * •

Growth: SIGOs can grow by accreting gas and dark matter from their
environment. Over time, this changes their mass and gas fraction, and can
impact their ability to collapse. We overview this process in §5.5

###  5.1 Gravitational Collapse

The timescale for gravitational collapse is important to determining whether a
SIGO can achieve the overdensities needed to form stars. In order for this to
occur, a SIGO needs to be able to collapse within a shorter time span than the
fall-back timescale to nearby DM halos. If this condition is not satisfied,
the SIGO will fall into a nearby DM halo before it forms stars, becoming
disrupted by tides or external pressure. As a starting point for this
timescale comparison, we calculate a free collapse timescale:

$t_{\rm collapse}\approx\frac{1}{\sqrt{G\rho_{\rm peak}}}\ ,$ (4)

depicted in Figure 6 for a representative SIGO, and in Appendix A, Figure 9
for individual SIGOs. For the collapse timescale for a typical SIGO in Figure
6, we assume a representative range of peak densities. The upper bound density
is taken as

${\rm log}_{10}(\rho_{\rm max})=0.7\times{\rm log}_{10}\left(\frac{{\rm
M}_{\rm SIGO}}{{\rm M}_{\odot}}\right)-26.3.$ (5)

The lower bound density is taken to be

${\rm log}_{10}(\rho_{\rm min})=0.4\times{\rm log}_{10}\left(\frac{{\rm
M}_{\rm SIGO}}{{\rm M}_{\odot}}\right)-25.6,$ (6)

where $\rho$ has units of g/cm3. Here we use peak rather than mean densities
to reflect that the SIGO’s star formation is likely primarily occurring in its
most dense regions: the process of star formation in a SIGO is not $100\%$
efficient.

As seen in the figure, the collapse timescale has a negative dependence on
mass, driven by larger central overdensities in more massive SIGOs. This
contributes to the enhanced ability of massive SIGOs to form stars.

###  5.2 Cooling

However, this timescale for free gravitational collapse does not paint the
whole picture of the collapse of SIGOs. In order for a SIGO to collapse, it
also must be able to efficiently cool, allowing its Jeans mass to decline and
permitting gravity to overcome gas pressure. There are 2 phases to this
cooling in a SIGO that has an initial Jeans mass above its actual mass: in the
first phase, the SIGO isochorically cools, lowering its Jeans mass until it
drops below the SIGO’s mass. Secondly, the SIGO begins to collapse, with the
dynamics of its collapse determined by the cooling and collapse timescales. To
account for this, we can define a cooling timescale ${\rm t}_{\rm cool}$,
given as the time it would take for a gas clump to cool isochorically at a
rate $\Lambda(T,n_{\rm H})$ to its Jeans temperature, or, for SIGOs with
supersonic velocity dispersions, to a temperature at which its density exceeds
the critical density. We then add the timescale for the second phase of this
collapse: the cooling time as defined in Schauer et al. (2021), taken from the
temperature at which the SIGO initially collapses.

$t_{\rm cool}=\frac{nk_{\rm B}T_{\rm collapse}}{(\gamma-1)\Lambda(T_{\rm
collapse})}+t_{\rm iso},$ (7)

where tiso is the additional cooling time a SIGO that initially does not
exceed its Jeans mass takes to isochorically cool to lower its Jeans mass to
MSIGO and begin to collapse:

${\rm t}_{\rm iso}=\left\\{\begin{array}[]{ c l }0&\quad{\scriptstyle{\rm
M}_{\rm SIGO}\geq{\rm M}_{\rm J}({\rm T}_{\rm init})}\\\ &\\\ \int_{{\rm
T}_{\rm Collapse}}^{{\rm T}_{\rm init}}\frac{nk_{\rm
B}}{(\gamma-1)\Lambda(T)}dT\ .&\quad{\scriptstyle{\rm M}_{\rm SIGO}<{\rm
M}_{\rm J}({\rm T}_{\rm init})}\end{array}\right.$ (8)

If the SIGO initially exceeds the Jeans mass, its collapse temperature is
taken as its volume-averaged temperature. Otherwise, its collapse temperature
(after the SIGO has isochorically cooled), is taken as the Jeans temperature
that the SIGO has cooled to:

${\rm T}_{\rm J}=\frac{3{\rm G}\mu{\rm m}_{\rm H}}{5{\rm k}_{\rm
B}}\sqrt[3]{\frac{36{\rm M}_{\rm SIGO}^{2}\rho_{\rm SIGO}}{\pi^{2}}}$ (9)

In SIGOs, $\Lambda$(T) is dominated by and approximated as the molecular
cooling rate defined in Galli & Palla (1998):

$\Lambda({\rm T})_{\rm H_{2}}=\frac{\Lambda({\rm LTE})}{1+[n^{\rm cr}/n({\rm
H})]},$ (10)

consistent with the molecular cooling prescription used in GRACKLE. Note here
that because the SIGO may accrete gas from its surroundings and gain mass as
it cools (raising its Jeans temperature), and because this does not account
for overdense regions within the SIGO which enhance cooling, this cooling
timescale assuming an isochoric process is an upper bound on the true cooling
time relevant to collapse.

We show the representative timescale based on Equation 8 in Figure 6 (see also
Figure 9, for this timescale for individual SIGOs ). To depict the typical
value of this timescale we take the median initial gas number density of a
SIGO, $1.08$ cm-3. To represent the maximum cooling time, we take a lower
bound gas density of $0.7$ cm-3. For the lower bound cooling time, we take
n${}_{\rm H}\sim 2.0$ cm-3 (see Figure 10 from Appendix A). We assume an
initial temperature of $500~{}{\rm K}$ or the Jeans temperature (Eq 9),
whichever is higher, for this representative line. As depicted in Figure 6,
most SIGOs at $z=20$ have not yet cooled substantially.

###  5.3 Fall-back Into Halos

Figure 7: Fall-back of a SIGO (the creation of a GC-like system at high
redshift): here we show the time evolution of a single SIGO in Run $2$vH$2$.
This SIGO, the left-most branch of this plot, is accreted by a nearby halo
(the gas component of which is shown in purple, the right branch of the plot,
with a very low baryon fraction). This accretion produces a more massive
merged object, containing the SIGO in its substructure. This SIGO (labelled
“GC-like” in the figure, for globular cluster-like candidate), will be subject
to future evolution within the halo, and is identifiable as a distinct
component, with the potential to continue to evolve into a star cluster. For a
movie of this evolution see here.

Eventually, the majority of SIGOs will be accreted onto a halo (most often
their parent halo), forming globular-cluster-like systems. Like GCs, these are
mostly expected to reside in the halos of their host galaxies (e.g., Chaboyer,
1999; Benedict et al., 2002; Chen et al., 2018). Naoz & Narayan (2014)
calculated the timescale on which a SIGO free-falls to its parent halo:

$t_{\rm ff}=0.27{\rm Gyr}\left(\frac{\Delta r_{\rm phys}}{0.59{\rm
kpc}}\right)^{\frac{3}{2}}\left(\frac{\Omega_{m}}{\Omega_{b}}\frac{{\rm
M}_{\rm b}}{10^{6}{\rm M}_{\odot}}\right),$ (11)

where $\Delta r_{\rm phys}$ is the physical separation between the SIGO and
the halo, and ${\rm M}_{\rm b}$ is the baryon mass of the SIGO. In order for a
SIGO to form stars outside of a DM halo, this timescale must be longer than
both the cooling and collapse timescales for that SIGO, permitting the object
to cool and reach high central densities prior to being affected by tides or
ram pressure stripping near the larger halo.

We show the representative timescale from Eq. (11) in Figure 6, where we
assumed a median physical separation of $0.4$ kpc between a SIGO and its
parent halo, with a range of $0.1-0.9$ kpc for a characteristic SIGO. As seen
in this figure, some SIGOs “fall” into a DM halo. Figure 7 shows the
evolutionary path of a single such SIGO (left of the figure with high baryon
fractions) that falls back into a dark matter halo (the gas component of which
forms the right branch of the figure, with low gas fractions). Notably, this
particular SIGO was accreted by a larger nearby DM halo rather than its parent
halo, allowing its lifetime outside of a halo to be unusually short for this
process, only $15$ Myr. We label the final outcome as “GC-like,” indicating
the potential for this object to form a star cluster at the outskirts of the
halo, similar to a present-day GC.

###  5.4 2-body Relaxation

In a simulation with limited resolution, each Voronoi cell has an associated
gas mass. Therefore, gravitational interactions between these cells can lead
to changing the object’s energy by an order of itself. This processes, called
$2$-body relaxation, has been shown to be a possible limiting factor on
simulations of SIGOs (Naoz & Narayan, 2014; Popa et al., 2016). This process
essentially “evaporates” the object as a function of time. The $2$-body
relaxation timescale for $N$ particles can be written as

$t_{\rm rlx}\sim\frac{{0.1\rm N}}{{\rm ln}({\rm N})}\times\frac{{\rm
r}}{\sigma},$ (12)

(Binney & Tremaine, 2008) where $\sigma$ is the $1$D velocity dispersion
within the SIGO and $r$ is taken to be the minimum principle axis length.
SIGOs can have as few as $\sim 100$ gas particles in these simulations, which
combined with their small sizes and fairly high velocity dispersions leads to
artificial relaxation times on the order of tens of Myr. This underscores the
need for caution when following the evolution of SIGOs to low redshifts in
these cosmological simulations: at our resolution, low-mass SIGOs with around
$10^{4}$ M⊙ of material are already beginning to be destroyed by artificial
relaxation at redshift $20$.

We depict the timescale from Eq. (12) in Figure 6 (see also 9), as a solid
black line. As seen in the figure, the main effect shortening the lifetime of
our SIGOs is the relaxation timescale generated by our limited resolution.
Thus, we expect SIGOs in the Universe to have increased longevity not captured
by our simulation.

###  5.5 Growth

Figure 8: Growth of a SIGO: here we show the time evolution of a single SIGO
in Run $2$vH$2$. This SIGO increases its mass by nearly a factor of $10$
between its formation and the final redshift of our simulation, accreting gas
from its surroundings with a similar baryon fraction to the SIGO itself and
thereby maintaining its baryon fraction over time. For a movie of this
evolution see here.

While not itself associated with a generally-defined timescale, it is
important to also note that SIGOs can grow through gas and dark matter
accretion. As can be seen in Figure 3, SIGOs form embedded in gas
streams–typical gas densities surrounding the forming SIGOs can be as high as
$\sim 20\%$ that of the SIGO itself. In addition, the material around the SIGO
may be itself enriched in baryons (though its gas fraction will not
necessarily be as high as that of the SIGO). This gas, as well as some dark
matter, may be accreted from a SIGO’s environment, increasing its mass over
time and allowing its gas fraction to change.

Figure 8 shows an example of this growth for a single SIGO in Run $2$vH$2$.
This SIGO formed with a mass of about $6\times 10^{4}$ M⊙ and a gas fraction
near $60\%$. The SIGO was able to accrete gas from its surroundings with a
similar baryon fraction to the SIGO itself, growing by a factor of nearly 10
while maintaining a roughly constant gas fraction. This has the potential to
impact the SIGO’s future evolution, as the more massive SIGO at z$=20$ is much
more likely to quickly cool and collapse to form stars through the process
outlined above in §5.2.

It is important to note, however, that this process of growth does not have to
maintain a SIGO’s gas fraction. Over time, accretion of surrounding material
can cause the gas fraction in SIGOs to drop below the $60\%$ threshold set for
SIGO identification, even though the SIGOs still exist (and are commonly still
enriched in gas). In Run $2$vH$2$, about $50\%$ of formed SIGOs at all
redshifts were found to eventually drop below a gas fraction of $60\%$ through
this process before z$=20$, while generally maintaining a gas fraction above
$40\%$ and maintaining or increasing their density. In future studies
modelling star formation in SIGOs to lower redshifts, it will be important to
follow these evolved SIGOs, which have higher masses than newborn SIGOs and
may be promising grounds for star formation in SIGOs at lower redshifts.

###  5.6 Timescale Comparison

Figure 6 shows these characteristic timescales for the evolutionary channels
for SIGOs mentioned above. In grey, one can see the characteristic timescale
for gravitational collapse of SIGOs identified at z$=20$ and at least one
other snapshot222Here we use multiple snapshots in order to limit false
detections of SIGOs and better constrain their dynamical properties. of Run
$2$vH$2$. In green, we show the range of fall-back timescales to the nearest
halo in the simulation. In blue, we show the range of cooling times for the
SIGOs to cool through molecular hydrogen cooling to Jeans collapse (for SIGOs
with a Mach dispersion less than 1) or collapse through the critical density
(for SIGOs with a Mach dispersion greater than 1). For all of these regions,
we also add a dotted line, showing the timescale for a SIGO of a given mass
with median properties (principle axis lengths, distance from the nearest
halo, temperature, gas fraction, and velocity dispersion). We also add a black
line showing a typical timescale for 2-body relaxation in the simulation at
simulation resolution. The figure also shows the age and mass of each SIGO in
Run $2$vH$2$ at z=$20$. Note that SIGOs with an age of $0.45$ Myr are those
that have just formed in the simulation.

As one can see from the plot, the cooling timescale has a strong (roughly
M-2.6) dependence on mass. Because of this, molecular cooling alone is not
sufficient to efficiently cool SIGOs below about $\sim 10^{5}$ M⊙ to collapse.
In essentially all cases with such SIGOs, the cooling timescale will be
dramatically longer than the collapse and fall-back timescales, so the SIGO
will fall into a halo before forming stars.

Above this mass limit, one must also compare the timescales for a SIGO to fall
back into a nearby halo or collapse gravitationally to determine SIGOs’ fates.
In isolation, high-mass SIGOs are frequently capable of collapsing to form
stars; however, because of the wide variance of these timescales, it is
possible for the fall-back timescale for a SIGO to be shorter than the
collapse timescale. Because of this, the fall-back and collapse timescales
must be individually compared, to determine which SIGOs will form stars before
accretion onto nearby DM objects.

One final consideration in studying the evolution of SIGOs apparent from this
plot is the simulation’s relaxation timescale. At the resolution of this
study, even $10^{5}$ M⊙ SIGOs dissipate due to $2$-body relaxation within $80$
Myr–a similar length of time to that between the first SIGO’s formation and
our z$=20$ snapshot. This implies that even here at redshift $20$, we are
undercounting the lowest-mass SIGOs, as well as potentially underestimating
the densities of those that exist. In order to eventually follow SIGOs through
the process of forming stars and merging with halos at later times, future
studies will need to work at higher resolutions ($\sim 200$ ${\rm M}_{\odot}$
or better), increasing the relaxation timescale.

##  6 Discussion

Understanding the behavior of supersonic turbulence in and around SIGOs is key
to understanding their evolution and potential for star formation (Chiou et
al., 2019). Here we show how supersonic turbulence acts together with
molecular cooling in SIGOs to form density peaks, with high enough densities
to become sites of star formation. This process is similar to that in GMCs,
where turbulence has a scale-dependent effect on star formation, boosting
densities on the sonic scale (e.g., Krumholz & McKee, 2005; Burkhart, 2018).

The structure and population-level properties of SIGOs are also similar to
those of GMCs. For example, we find that SIGOs have similar substructures and
masses to GMCs. In particular, we find that SIGOs have core-like structures
and filaments (Figure 10, or see the evolved SIGO in Figure 8), similar to
those in GMCs (Mac Low & Klessen, 2004; Krumholz & McKee, 2005; Mocz &
Burkhart, 2018). Further, we find that the mass spectrum of SIGOs is
consistent with the mass spectrum of GMCs in the outer Milky Way (see Figure 5
in §4).

In particular, we show for the first time in a simulation that SIGOs form on
scales comparable to the sonic scale (Figure 4). Supersonic turbulence aid in
the formation and collapse of SIGOs, with SIGOs forming as high density peaks.
Molecular cooling lowers these objects’ Jeans scale. When this process is
sufficiently efficient, this lowers the Jeans scale below the sonic scale and
permits collapse of these newly formed SIGOs.

Chiou et al. (2019) found that most SIGOs should form stars through atomic
cooling alone. However, this result was called into question by Schauer et al.
(2021), which included atomic and molecular cooling and did not find star
formation sites outside of DM halos in a simulation. Nakazato et al. (2022),
on the other hand, followed a SIGO in a zoom-in simulation and confirmed that
the SIGO experienced Jeans collapse. The modifications in this paper to the
prescriptions in Chiou et al. (2019) are a step towards reconciling these
disparate results. Correcting the prescription for collapse in SIGOs based on
their asphericity, we come to the conclusion that SIGOs should not form stars
when considering only atomic hydrogen cooling. With molecular hydrogen cooling
accounted for, SIGOs should be able to form stars. However, not all SIGOs
undergo global collapse, exceeding the Jeans scale on the physical scale of
the entire SIGO rather than in substructure (neglecting metal line cooling,
which has the potential to enhance star formation in SIGOs but has not yet
been studied in a simulation; see e.g. Schauer et al., 2021).

At $z=20$, when most SIGOs have not yet cooled, about $5$ SIGOs will be
capable of undergoing Jeans collapse. Thus, at lower redshifts, we expect many
SIGOs to collapse. We estimate a star-forming SIGO abundance of $0.63$ Mpc-3
at z$=20$. This abundance is comparable to the present-day local density of
low-metallicity globular clusters (i.e., $0.44$ Mpc-3, Rodriguez et al. 2015).
Because of the high Poisson uncertainty inherent to this figure, this suggests
an explanation for the lack of star-forming SIGOs in Schauer et al. (2021). A
smaller simulation box size and $\sigma_{8}$, combined with lower star-forming
SIGO abundances than could be inferred from the literature at the time, may
allow a $<1.5\sigma$ Poisson fluctuation to yield a simulation box with no
star-forming SIGOs outside of halos, neglecting metal line cooling. However,
these results indicate that the abundance of star-forming SIGOs in the
Universe may be quite high in regions with significant stream velocities.

The effect of $\sigma_{8}$ on the redshift of SIGO formation must also be
considered in contextualizing these results. Here, we have assumed an elevated
$\sigma_{8}$, representative of a $\sim 2\sigma$ density fluctuation. This
elevated value of $\sigma_{8}$ yields more power; i.e., it helps structure
develop earlier (e.g., Greif et al., 2011; Park et al., 2020). As shown in
Figure 10 (see Appendix A), the comoving density of SIGOs before they collapse
is not significantly dependent on redshift, so the pre-collapse physical
density and therefore the cooling rate of SIGOs is dependent on their redshift
of formation. In overdense regions such as those that form galaxy clusters,
then, SIGOs can more efficiently cool and collapse than in underdense regions.
This factor explains why SIGOs are not seen in regions of the Universe outside
of galaxies: SIGOs formed too late to efficiently cool in the underdense
regions of the early Universe that gave rise to these volumes. On the other
hand, SIGOs form earlier and are more prone to form star clusters in overdense
regions.

As highlighted in Figure 6, the resolution limitation yields an artificial age
limit on SIGOs, which should not exist in the case of gas objects(e.g., Naoz &
Narayan, 2014). Thus, we use a semi-analytical formulation to show that the
evolution of SIGOs generally proceeds to one of two states: either
gravitational collapse to form stars outside of halos (followed by accretion
onto halos as a cluster), or fall back into a nearby DM object, where the SIGO
can evolve as substructure within the halo. This can be expressed as a
timescale problem: SIGOs with a shorter collapse timescale and cooling
timescale than fall-back timescale (to the nearest DM object) should form
stars outside of halos. Otherwise, we suggest that the SIGO can form a GC-like
object, since GCs are often found near the edges of DM halos. Star formation
may take place in this overdense gas at the outskirts of the halo. An example
evolutionary path of a SIGO that has fallen back into a nearby halo is shown
in Figure 7. This SIGO had a low mass and an unusually low fall-back timescale
(owing to a massive nearby halo), and therefore was unable to collapse outside
of a DM halo. We also show, however, that a low initial mass may not prevent a
SIGO from reaching the masses needed from star formation outside of halos. In
Figure 8, we show an example of a SIGO that grows in mass by a factor of
nearly $10$ over the course of the simulation, lowering its cooling time by
$2$ orders of magnitude. This presents the possibility that even SIGOs that
form outside of halos at lower masses may be capable of cooling to form stars
outside of halos through growth, and suggests the need in future work for an
analytic model for understanding SIGO growth rates.

Future studies aiming to investigate the collapse of SIGOs should carefully
consider the cooling timescale (Figure 9) in determining a final redshift and
resolution for their simulations. While individual SIGOs with unusually short
cooling timescales may collapse prior to z$=20$, a typical early SIGO forming
at z$=25$ with a typical cooling timescale of $100$ Myr will not fully
collapse until around z$=17$ (or potentially later with lower values of
$\sigma_{8}$). For similar reasons, it is important to account for the two-
body relaxation timescale in simulations beyond z$=20$, which may be
comparable to $100$ Myr at similar resolutions to those in this paper.

Explicit treatments of star formation with metal line cooling are needed to
study the further evolution of SIGOs, including the size of their stellar
populations. Exciting observational predictions may become possible with this
next step, such as determining a half-light radius for SIGOs, and comparing it
not only to present day GCs, but also to expected higher-redshift GC
observations and potentially GC progenitor observations from JWST. These
simulations will also begin to establish the efficiency of star formation in
SIGOs, allowing us to understand characteristic maximum stellar population
masses from the process. As these simulations develop, zoom-in simulations
including radiative feedback should also be run to the present day, in order
to give us a picture of the entire evolutionary history of a SIGO, allowing us
to much more firmly connect these high-redshift objects to the structures we
see in today’s Universe.

W.L., S.N., Y.S.C, B.B., F.M., and M.V. thank the support of NASA grant No.
80NSSC20K0500 and the XSEDE AST180056 allocation, as well as the Simons
Foundation Center for Computational Astrophysics and the UCLA cluster Hoffman2
for computational resources. F.M. acknowledges support from the Program “Rita
Levi Montalcini” of the Italian MUR. S.N thanks Howard and Astrid Preston for
their generous support. Y.S.C thanks the partial support from UCLA
dissertation year fellowship. B.B. also thanks the the Alfred P. Sloan
Foundation and the Packard Foundation for support. MV acknowledges support
through NASA ATP grants 16-ATP16-0167, 19-ATP19-0019, 19-ATP19-0020,
19-ATP19-0167, and NSF grants AST-1814053, AST-1814259, AST-1909831 and
AST-2007355. NY acknowledges financial support from JST AIP Acceleration
Research JP20317829.

##  Appendix A Comparison of Important Timescales

Figure 9: Timescale Comparisons: here we show a number of important timescales
to the evolution of SIGOs at z$=20$ in Run $2$vH$2$. In the first $3$ panels
(a-c), a black dot indicates a SIGO from the simulation. A red star indicates
a SIGO that is likely to form stars outside of a DM halo, with cooling and
collapse timescales shorter than its fallback timescale. Note that $2$ star-
forming SIGOs could not be pictured in this figure, due to their very short
cooling and collapse timescales indicative that collapse has already occurred.
The red lines in these figures indicate equality between the two timescales
being compared. The final panel (d) shows the ages of z=$20$ SIGOs at the
resolution of the simulation plotted against their masses. The red line shows
the relaxation timescale at simulation resolution of a SIGO with typical
properties as a function of mass.

For another comparison of all of these timescales, the top left panel of
Figure 9 shows the collapse timescale of SIGOs in Run $2$vH$2$. This is
plotted against the SIGOs’ median fall-back timescale (across time) to their
nearest halo. This median corrects for errors in the FOF finder for DM halos,
ensuring that all nearby halos are accounted for, even if they are erroneously
not found or their centers are mislocated at z$=20$. As one can see, most
SIGOs have a shorter collapse timescale than fall-back timescale (and are
below the red line in the panel reflecting equal timescales). This indicates
that the cooling timescale is key to determining whether or not SIGOs will
collapse before falling into nearby DM halos.

Panels (b) and (c) of Figure 9 show two other comparisons of the timescales
involved in this evolution. Panel (b) shows a comparison of the cooling
timescale and the fall-back timescale for all z$=20$ SIGOs. This panel further
demonstrates that while the collapse timescale is important for the evolution
of SIGOs, that timescale alone is usually insufficient to determine whether or
not a SIGO will collapse independently of a DM halo– the cooling timescale
tends to be a bigger driver of such SIGOs’ evolution, especially at lower
masses. Panel (c) compares the cooling timescale for SIGOs to their collapse
timescale, showing that the cooling timescale for SIGOs is, in fact, always
longer than the collapse timescale at formation. Note that this cooling
timescale includes an isochoric cooling term for the SIGO to begin to collapse
(Eq. 8): this does not reflect the balance between cooling times and collapse
times once SIGOs achieve the Jeans mass and begin to collapse.

Figure 10: Densities: here we show the mean (left panel) and mean (right
panel) comoving densities of SIGOs at various redshifts in Run $2$vH$2$. In
each panel, a mark indicates a SIGO from the simulation. The peak density is
calculated as the volume-weighted average density of the most dense Voronoi
gas cell in the SIGO and its $9$ nearest neighbours.

SIGOs generally form in a similar range of densities regardless of mass. Their
early evolution occurs primarily through cooling, and sometimes mass growth
through gas accretion, rather than through immediate global collapse. However,
SIGOs’ collapse timescales depend on their mass, as their peak densities in
their cores have a mass dependence. In the right panel of Figure 10, we show
this effect. As one can see, even SIGOs that have evolved for some time at
z=$20$ do not have mean gas densities that exhibit a strong dependence on
mass. This is expected: most SIGOs at this time are around $10$ Myr old, and
nearly all of these SIGOs have not yet had a chance to globally cool and
collapse on $\gtrapprox 50$ Myr timescales. However, there are hints of
substructure formation in high-mass SIGOs. In the left panel, we show the
trend in peak density with mass, excepting $4$ z$=20$ SIGOs that have already
partially or fully collapsed and reached much higher densities. A linear
regression yields an approximate formula for the relation between peak density
and mass at z=$20$: log10($\rho_{\rm peak}$ (g/cm3)) =
$0.6\times$log10(MSIGO(M⊙)) - $26.2$. As one can see, there is a clear
positive correlation between peak density and mass–one that becomes stronger
at later redshifts.

## References

  * Abel et al. (1998) Abel, T., Anninos, P., Norman, M. L., & Zhang, Y. 1998, ApJ, 508, 518, doi: 10.1086/306410
  * Appel et al. (2022) Appel, S. M., Burkhart, B., Semenov, V. A., Federrath, C., & Rosen, A. L. 2022, ApJ, 927, 75, doi: 10.3847/1538-4357/ac4be3
  * Benedict et al. (2002) Benedict, G. F., McArthur, B. E., Fredrick, L. W., et al. 2002, AJ, 123, 473, doi: 10.1086/338087
  * Bergin & Tafalla (2007) Bergin, E. A., & Tafalla, M. 2007, ARA&A, 45, 339, doi: 10.1146/annurev.astro.45.071206.100404
  * Binney & Tremaine (2008) Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press)
  * Bromm (2013) Bromm, V. 2013, Reports on Progress in Physics, 76, 112901, doi: 10.1088/0034-4885/76/11/112901
  * Burkhart (2018) Burkhart, B. 2018, ApJ, 863, 118, doi: 10.3847/1538-4357/aad002
  * Burkhart (2021) —. 2021, PASP, 133, 102001, doi: 10.1088/1538-3873/ac25cf
  * Burkhart et al. (2009) Burkhart, B., Falceta-Gonçalves, D., Kowal, G., & Lazarian, A. 2009, ApJ, 693, 250, doi: 10.1088/0004-637X/693/1/250
  * Burkhart & Lazarian (2012) Burkhart, B., & Lazarian, A. 2012, ApJ, 755, L19, doi: 10.1088/2041-8205/755/1/L19
  * Chaboyer (1999) Chaboyer, B. 1999, in Astrophysics and Space Science Library, Vol. 237, Post-Hipparcos Cosmic Candles, ed. A. Heck & F. Caputo, 111, doi: 10.1007/978-94-011-4734-7_7
  * Chen et al. (2018) Chen, S., Richer, H., Caiazzo, I., & Heyl, J. 2018, ApJ, 867, 132, doi: 10.3847/1538-4357/aae089
  * Chiaki & Wise (2019) Chiaki, G., & Wise, J. H. 2019, MNRAS, 482, 3933, doi: 10.1093/mnras/sty2984
  * Chiou et al. (2019) Chiou, Y. S., Naoz, S., Burkhart, B., Marinacci, F., & Vogelsberger, M. 2019, ApJ, 878, L23, doi: 10.3847/2041-8213/ab263a
  * Chiou et al. (2021) —. 2021, ApJ, 906, 25, doi: 10.3847/1538-4357/abc88f
  * Chiou et al. (2018) Chiou, Y. S., Naoz, S., Marinacci, F., & Vogelsberger, M. 2018, MNRAS, 481, 3108, doi: 10.1093/mnras/sty2480
  * Conaboy et al. (2022) Conaboy, L., Iliev, I. T., Fialkov, A., Dixon, K. L., & Sullivan, D. 2022, arXiv e-prints, arXiv:2207.11614.  https://arxiv.org/abs/2207.11614
  * Elmegreen & Efremov (1997) Elmegreen, B. G., & Efremov, Y. N. 1997, ApJ, 480, 235, doi: 10.1086/303966
  * Galli & Palla (1998) Galli, D., & Palla, F. 1998, A&A, 335, 403.  https://arxiv.org/abs/astro-ph/9803315
  * Glover (2013) Glover, S. 2013, in Astrophysics and Space Science Library, Vol. 396, The First Galaxies, ed. T. Wiklind, B. Mobasher, & V. Bromm, 103, doi: 10.1007/978-3-642-32362-1_3
  * Greif et al. (2011) Greif, T. H., Springel, V., White, S. D. M., et al. 2011, ApJ, 737, 75, doi: 10.1088/0004-637X/737/2/75
  * Grudić et al. (2022) Grudić, M. Y., Hafen, Z., Rodriguez, C. L., et al. 2022, arXiv e-prints, arXiv:2203.05732.  https://arxiv.org/abs/2203.05732
  * Haiman et al. (1996) Haiman, Z., Rees, M. J., & Loeb, A. 1996, ApJ, 467, 522, doi: 10.1086/177628
  * Harris & Pudritz (1994) Harris, W. E., & Pudritz, R. E. 1994, ApJ, 429, 177, doi: 10.1086/174310
  * Hirano et al. (2017) Hirano, S., Hosokawa, T., Yoshida, N., & Kuiper, R. 2017, Science, 357, 1375, doi: 10.1126/science.aai9119
  * Hummel et al. (2016) Hummel, J. A., Stacy, A., & Bromm, V. 2016, MNRAS, 460, 2432, doi: 10.1093/mnras/stw1127
  * Kravtsov & Gnedin (2005) Kravtsov, A. V., & Gnedin, O. Y. 2005, ApJ, 623, 650, doi: 10.1086/428636
  * Krumholz & McKee (2005) Krumholz, M. R., & McKee, C. F. 2005, ApJ, 630, 250, doi: 10.1086/431734
  * Krumholz & Tan (2007) Krumholz, M. R., & Tan, J. C. 2007, ApJ, 654, 304, doi: 10.1086/509101
  * Lake et al. (2021) Lake, W., Naoz, S., Chiou, Y. S., et al. 2021, ApJ, 922, 86, doi: 10.3847/1538-4357/ac20d0
  * Larson (1981) Larson, R. B. 1981, MNRAS, 194, 809, doi: 10.1093/mnras/194.4.809
  * Lithwick & Dalal (2011) Lithwick, Y., & Dalal, N. 2011, ApJ, 734, 100, doi: 10.1088/0004-637X/734/2/100
  * Mac Low & Klessen (2004) Mac Low, M.-M., & Klessen, R. S. 2004, Reviews of Modern Physics, 76, 125, doi: 10.1103/RevModPhys.76.125
  * McKee & Ostriker (2007) McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565, doi: 10.1146/annurev.astro.45.051806.110602
  * Mocz & Burkhart (2018) Mocz, P., & Burkhart, B. 2018, MNRAS, 480, 3916, doi: 10.1093/mnras/sty1976
  * Mocz et al. (2017) Mocz, P., Burkhart, B., Hernquist, L., McKee, C. F., & Springel, V. 2017, ApJ, 838, 40, doi: 10.3847/1538-4357/aa6475
  * Nakazato et al. (2022) Nakazato, Y., Chiaki, G., Yoshida, N., et al. 2022, ApJ, 927, L12, doi: 10.3847/2041-8213/ac573e
  * Naoz & Barkana (2005a) Naoz, S., & Barkana, R. 2005a, MNRAS, 362, 1047, doi: 10.1111/j.1365-2966.2005.09385.x
  * Naoz & Barkana (2005b) —. 2005b, MNRAS, 362, 1047, doi: 10.1111/j.1365-2966.2005.09385.x
  * Naoz et al. (2009) Naoz, S., Barkana, R., & Mesinger, A. 2009, MNRAS, 399, 369, doi: 10.1111/j.1365-2966.2009.15282.x
  * Naoz & Narayan (2013) Naoz, S., & Narayan, R. 2013, Physical Review Letters, 111, 051303, doi: 10.1103/PhysRevLett.111.051303
  * Naoz & Narayan (2014) —. 2014, ApJ, 791, L8, doi: 10.1088/2041-8205/791/1/L8
  * Naoz et al. (2011) Naoz, S., Yoshida, N., & Barkana, R. 2011, MNRAS, 416, 232, doi: 10.1111/j.1365-2966.2011.19025.x
  * Naoz et al. (2012) Naoz, S., Yoshida, N., & Gnedin, N. Y. 2012, ApJ, 747, 128, doi: 10.1088/0004-637X/747/2/128
  * Naoz et al. (2013) —. 2013, ApJ, 763, 27, doi: 10.1088/0004-637X/763/1/27
  * Padoan & Nordlund (1999) Padoan, P., & Nordlund, Å. 1999, ApJ, 526, 279, doi: 10.1086/307956
  * Padoan & Nordlund (2011) —. 2011, ApJ, 730, 40, doi: 10.1088/0004-637X/730/1/40
  * Park et al. (2020) Park, H., Ahn, K., Yoshida, N., & Hirano, S. 2020, ApJ, 900, 30, doi: 10.3847/1538-4357/aba26e
  * Popa et al. (2016) Popa, C., Naoz, S., Marinacci, F., & Vogelsberger, M. 2016, MNRAS, 460, 1625, doi: 10.1093/mnras/stw1045
  * Rodriguez et al. (2015) Rodriguez, C. L., Morscher, M., Pattabiraman, B., et al. 2015, Phys. Rev. Lett., 115, 051101, doi: 10.1103/PhysRevLett.115.051101
  * Rosolowsky (2005) Rosolowsky, E. 2005, PASP, 117, 1403, doi: 10.1086/497582
  * Saslaw & Zipoy (1967) Saslaw, W. C., & Zipoy, D. 1967, Nature, 216, 976, doi: 10.1038/216976a0
  * Schauer et al. (2021) Schauer, A. T. P., Bromm, V., Boylan-Kolchin, M., Glover, S. C. O., & Klessen, R. S. 2021, ApJ, 922, 193, doi: 10.3847/1538-4357/ac27aa
  * Schaye et al. (2015) Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521, doi: 10.1093/mnras/stu2058
  * Schneider et al. (2012) Schneider, A., Smith, R. E., Macciò, A. V., & Moore, B. 2012, MNRAS, 424, 684, doi: 10.1111/j.1365-2966.2012.21252.x
  * Seljak & Zaldarriaga (1996) Seljak, U., & Zaldarriaga, M. 1996, ApJ, 469, 437, doi: 10.1086/177793
  * Semenov et al. (2016) Semenov, V. A., Kravtsov, A. V., & Gnedin, N. Y. 2016, ApJ, 826, 200, doi: 10.3847/0004-637X/826/2/200
  * Shapiro et al. (2010) Shapiro, K. L., Genzel, R., & Förster Schreiber, N. M. 2010, MNRAS, 403, L36, doi: 10.1111/j.1745-3933.2010.00810.x
  * Sheth et al. (2001) Sheth, R. K., Mo, H. J., & Tormen, G. 2001, MNRAS, 323, 1, doi: 10.1046/j.1365-8711.2001.04006.x
  * Smith et al. (2017) Smith, B. D., Bryan, G. L., Glover, S. C. O., et al. 2017, MNRAS, 466, 2217, doi: 10.1093/mnras/stw3291
  * Springel (2005) Springel, V. 2005, MNRAS, 364, 1105, doi: 10.1111/j.1365-2966.2005.09655.x
  * Stacy et al. (2011) Stacy, A., Bromm, V., & Loeb, A. 2011, ApJ, 730, L1, doi: 10.1088/2041-8205/730/1/L1
  * Tseliakhovich & Hirata (2010) Tseliakhovich, D., & Hirata, C. 2010, Phys. Rev. D, 82, 083520, doi: 10.1103/PhysRevD.82.083520
  * Vogelsberger et al. (2013) Vogelsberger, M., Genel, S., Sijacki, D., et al. 2013, MNRAS, 436, 3031, doi: 10.1093/mnras/stt1789
  * Vogelsberger et al. (2020) Vogelsberger, M., Marinacci, F., Torrey, P., & Puchwein, E. 2020, Nature Reviews Physics, 2, 42, doi: 10.1038/s42254-019-0127-2
  * Vogelsberger & White (2011) Vogelsberger, M., & White, S. D. M. 2011, MNRAS, 413, 1419, doi: 10.1111/j.1365-2966.2011.18224.x
  * Vogelsberger et al. (2014a) Vogelsberger, M., Genel, S., Springel, V., et al. 2014a, MNRAS, 444, 1518, doi: 10.1093/mnras/stu1536
  * Vogelsberger et al. (2014b) —. 2014b, Nature, 509, 177, doi: 10.1038/nature13316
  * Whitmore & Schweizer (1995) Whitmore, B. C., & Schweizer, F. 1995, AJ, 109, 960, doi: 10.1086/117334
  * Whitmore et al. (1999) Whitmore, B. C., Zhang, Q., Leitherer, C., et al. 1999, AJ, 118, 1551, doi: 10.1086/301041
  * Yoshida et al. (2003) Yoshida, N., Abel, T., Hernquist, L., & Sugiyama, N. 2003, ApJ, 592, 645, doi: 10.1086/375810
  * Yoshida et al. (2008) Yoshida, N., Omukai, K., & Hernquist, L. 2008, Science, 321, 669, doi: 10.1126/science.1160259
  * Zel’dovich (1970) Zel’dovich, Y. B. 1970, A&A, 5, 84


# Non-linear frequency warping using constant-Q transformation for speech
emotion recognition

Premjeet Singh1, Goutam Saha2 Dept of Electronics and ECE
Indian Institute of Technology Kharagpur, Kharagpur, India
<EMAIL_ADDRESS><EMAIL_ADDRESS>Md Sahidullah3 Université de
Lorraine
CNRS, Inria, LORIA, F-54000, Nancy, France
<EMAIL_ADDRESS>
###### Abstract

In this work, we explore the constant-Q transform (CQT) for speech emotion
recognition (SER). The CQT-based time-frequency analysis provides variable
spectro-temporal resolution with higher frequency resolution at lower
frequencies. Since lower-frequency regions of speech signal contain more
emotion-related information than higher-frequency regions, the increased low-
frequency resolution of CQT makes it more promising for SER than standard
short-time Fourier transform (STFT). We present a comparative analysis of
short-term acoustic features based on STFT and CQT for SER with deep neural
network (DNN) as a back-end classifier. We optimize different parameters for
both features. The CQT-based features outperform the STFT-based spectral
features for SER experiments. Further experiments with cross-corpora
evaluation demonstrate that the CQT-based systems provide better
generalization with out-of-domain training data.

###### Index Terms:

Speech emotion recognition (SER), Constant-Q transform (CQT), Mel frequency
analysis, Cross-corpora evaluation.

##  I Introduction

The _speech emotion recognition_ (SER) is the task for recognizing emotion
from human speech. The potential applications of SER include _human-computer
interaction_ , _sentiment analysis_ and _health-care_ [1, 2, 3, 4]. Humans
naturally sense the emotions in speech while machines find it difficult to
characterize them [5, 6]. Techniques proposed till date have significantly
increased the machine’s ability to recognize speech emotions. However, the
task is still challenging mainly due to the presence of large _interpersonal
and intrapersonal variability_ and the _differences in speech quality_ used to
train and evaluate the system. The goal of this work is to develop an improved
SER system by considering emotion-specific acoustic parameters from speech
that are assumed to be more robust to unwanted variabilities.

Previous studies in SER research have shown that spectral and prosodic
characteristics of speech contain emotion-related information. Spectral
features include _lower formants frequencies_ (F1 and F2), _speech amplitude
and energy_ , _zero crossing rate_ (ZCR) and spectral parameters, e.g., like
_spectral flux_ and _spectral roll-off_ [7, 8, 9]. Prosodic features include
_pitch_ , _pitch harmonics_ , _intonation_ , and _speaking rate_ [8, 9]. The
acoustic front-ends are used with _Gaussian mixture model_ (GMM) or _support
vector machines_ (SVMs) as back-end classifiers for SER tasks [10]. Studies
with prosody reveal that high arousal emotions, such as _Angry_ , _Happy_ and
_Fear_ , have higher average pitch values with abrupt pitch variations whereas
low arousal emotions like _Sadness_ and _Neutral_ have lower pitch values with
consistent pitch contours [1, 11, 12, 13, 14]. The authors in [15] have
reported that recognition accuracy of Anger is higher near F2 (1250-1750 Hz)
and that of Neutral is higher near F1 (around 200-1000 Hz). Authors in [16]
report that center frequencies of F2 and F3 are reduced in depressed
individuals. In [17], the authors report that high arousal emotions have
higher mean F1 and lower F2 and high (positive) valence emotions have high
mean F2. In [18], authors report discrimination between idle and negative
emotions using temporal patterns in formants. In [19], the authors have
demonstrated that non-linear frequency scales, such as _logarithmic_ , mel and
_equivalent rectangular bandwidth_ (ERB), have considerable impact in SER
performance over linear frequency scale.

Recent works with deep learning methods such as _convolutional neural
networks_ (CNNs) or CNN with _recurrent neural networks_ (CNN-RNNs), use
spectrogram or raw waveform as input and have shown impressive results [20,
21, 22]. These data-driven methods automatically learn the emotion-related
representation, however, the role of individual speech attributes in the
decision making process is not clear due to the lack of _explainability_. On
the other hand, the generalization of these methods remains an open problem,
especially when the audio-data for train and test are substantially different
in terms of language and speech quality [23].

Figure 1: F-ratios of spectrograms based on CQT (top), mel-filter (middle),
and standard STFT (bottom) corresponding to the frequency bins. We use the
speech sentences with fixed text ‘a02’ from EmoDB database (discussed in
Section III-A). We select the same text assuming spectra characteristics of
emotions to be text-dependent. First column shows the values over the entire
frequency range while the second column focuses only on the lower-frequency
regions.

We address this generalization issue by capturing emotion-related information
from speech before processing with a neural network back-end. Given the fact
that the low and mid frequency regions of speech spectrum contain pitch
harmonics and lower formants that are relevant for emotion recognition, we
propose to use a more appropriate approach for time-frequency analysis that
produces emotion-oriented speech representation in the first place. Even
though the processing with _mel frequency warping_ introduces non-linearity in
some sense, the power spectrum from the speech is essentially computed with a
uniform frequency resolution. We propose to use a time-frequency analysis
method called _constant-Q transform_ (CQT). This transformation offers higher
frequency resolution at low-frequency regions and higher time resolution at
high-frequency regions. As the pitch harmonics and lower formants reside in
the low-frequency regions of speech spectrum, we hypothesize that keeping high
resolution in this region may efficiently capture emotion-related information.

The CQT was initially proposed for music signal processing [24]. Then it was
also applied in different speech processing tasks, e.g., _anti-spoofing_ [25,
26], _speaker verification_ [27] and _acoustic scene classification_ [28].
Recently, the CQT has also been studied for SER [29], but without success.
This is possibly due to the lack of optimization of CQT parameters and/or the
applied end-to-end model fails to exploit the advantages of CQT. Recent
studies show that CNN-based models are suitable for SER including cross-
corpora evaluation [23, 30]. In this work, we also adopt a CNN-based approach
for modeling SER systems. Our main contributions in this work are summarized
as follows: (i) We propose a new framework for CQT-based SER by optimizing CQT
extraction parameters for reduced redundancy and improved performance, (ii) We
investigate CNN architecture known as _time-delay neural networks_ (TDNNs)
suitable for speech pattern classification tasks [31, 32] for SER, and (iii)
We perform cross-corpora evaluation with three different speech corpora to
assess the generalization ability of the proposed method. Our results
demonstrate that the optimized CQT features not only outperform short-time
Fourier transform (STFT) features but also provide better generalization.

##  II Methodology

In this section, we discuss the CQT-based feature extraction framework and the
TDNN architecture for emotion recognition.

###  II-A Constant-Q transform

The CQT of a time-domain signal $x[n]$ is defined as,

$X[k]=\frac{1}{N[k]}\sum_{n=0}^{N[k]-1}W[k,n]x[n]e^{-jw_{k}n},$ (1)

where $X[k]$ is the CQT coefficient for $k$-th frequency bin, $W[k,n]$ is the
time-domain window for $k$-th bin with duration $N[k]$, $x[n]$ denotes the
time samples and $w_{k}=\frac{2\pi Q}{N[k]}$, where $Q$ is the (constant) Q
factor of the filter banks [24]. In CQT computation, the window length $N[k]$
varies for different values of $k$. Hence, $x[n]$ is correlated with sinusoids
of different lengths with equal cycles of oscillation. This leads to
constant-Q filter bank representation with geometrically spaced center
frequencies over frequency octaves. Hence, we obtain a time-frequency
representation which has frequency resolution varying from high to low towards
increasing frequencies.

The CQT representation of an audio signal depends on the _number of octaves of
frequencies_ and the _number of frequency bins per each octave_. The number of
octaves depends upon the chosen minimum frequency ($F_{\mathrm{min}}$) and the
maximum frequency ($F_{\mathrm{max}}$) of operation, and this equals to
$\log_{2}{\frac{F_{\mathrm{max}}}{F_{\mathrm{min}}}}$ [25]. The CQT
representation with reduced number of total frequency bins over a fixed number
of octaves will provide detailed information for lower frequency region with
reduced redundancy. Conversely, due to linearly spaced frequency bins, _short-
time Fourier transform_ (STFT) does not offer this flexibility. Fixing
$F_{\mathrm{min}}$ to $32.7$ Hz and $F_{\mathrm{max}}$ to _Nyquist frequency_
gives us approximately $8$ octaves. _Hop length_ in CQT computation defines
the number of time samples by which CQT computation window moves. The CQT also
has resemblance with _continuous wavelet transform_ (CWT) which provides
variable time-frequency resolution and has been found helpful for SER [33].

During the CQT-based feature extraction process, the CQT coefficients are
uniformly resampled and then processed with _discrete cosine transform_ (DCT)
to compute speech features known as _constant-Q cepstral coefficients_
(CQCCs).

We perform _class separability_ analysis of the time-frequency representations
by computing the F-ratios [34]. The Fig. 1 shows the F-ratio obtained at
different frequency bins. The higher F-ratios at lower bins for CQT and STFT
show the presence of more discriminative information. The figure also
indicates that CQT-spectrogram has more number of discriminative coefficients
on an average over others due to higher resolution in low-frequency regions.

###  II-B CNN architecture

The time-frequency representation of speech-like signal is suitable to be used
with 1-D CNN, popularly known as TDNN in speech processing literature. Our
method is inspired by the TDNN-based _x-vector_ system [32] developed for
speaker verification task. This processes speech information at _frame_ and
_segment_ level. In frame level, the TDNN captures _contextual information_ by
applying _kernel_ over adjacent frames and by processing each speech frame in
an identical manner. This also applies _dilation_ in the temporal domain to
reduce redundancy and to make it computationally efficient. The frame-level
information is processed with several TDNN layers having different kernel
sizes and dilation parameters. Finally, _temporal pooling_ aggregates frame-
level information into segment-level and this is followed by processing with
_fully connected_ (FC) and _softmax_ layer for classification objective. The
standard x-vector system computes the segment-level intermediate
representation referred as _embeddings_ which are further processed with
another system for classification. In contrast, our proposed method trains the
network in an end-to-end fashion for which the emotion for a test speech is
obtained from the output of the trained network.

We empirically optimize the parameters for TDNN architecture. Finally, we use
four TDNN layers, followed by statistics pooling with mean and standard
deviation, and one FC layer before _softmax_. Table I describes the parameters
for different layers.

TABLE I: The parameters of CNN architecture for SER. Layer | Size | Kernel Size | Dilation
---|---|---|---
TDNN | 32 | 5 | 1
TDNN | 32 | 3 | 2
TDNN | 32 | 3 | 3
TDNN | 64 | 1 | 1
Statistics Pooling (Mean and SD) | 128 | - | -
Fully Connected | 64 | - | -
Softmax | #Classes | - | -

Figure 2: Performance comparison of CQT for different parameter values.
Optimized CQT shows the response of CQT with 3 bins per octave and hop length
of 64 samples. Baseline MFSC corresponds to MFSC extraction with standard
value of 128 mel-filters and 160 samples hop length, whereas, optimized MFSC
has 24 mel-filters with 64 hop. Optimized CQCC and MFCC shown are obtained
after applying DCT over optimized MFSC and CQT. The results shown in figure
are obtained over EmoDB database only.

##  III Experimental setup

###  III-A Speech corpora

In our experiments, we use three different speech corpora which are described
in Table II. We downsample speech files at sampling rate of 16 kHz when
required. The EmoDB is a German language corpora while RAVDESS and IEMOCAP are
in English. For IEMOCAP database, we select only four emotions (Angry, Happy,
Sad and Neutral) as some of the emotion class have inadequate data for
training neural network models [30]. We perform cross-corpora SER experiments
by selecting the same four emotions.

TABLE II: Summary of the speech corpora used in the experiments. (F=Female,
M=Male)

Databases Speakers Emotions Berlin Emotion Database (EmoDB) [35] 10 (5 F, 5 M)
7 (Anger, Sad, Boredom, Fear, Happy, Disgust and Neutral) Ryerson Audio-Visual
Database of Emotional Speech and Song (RAVDESS) [36] 24 (12 F, 12 M) 8 (Calm,
Happy, Sad, Angry, Neutral, Fearful, Surprise, and Disgust) Interactive
Emotional Dyadic Motion Capture Database (IEMOCAP) [37] 10 (5 F, 5 M) 4
(Happy, Angry, Sad and Neutral)

Figure 3: Confusion matrices for emotion classification experiment with
optimized CQT and MFSC features in EmoDB corpus. Given values are the ratio of
utterances identified in column class to the total number of utterances in
every corresponding row class.

###  III-B Experimental details & evaluation methodology

First, we optimize the parameters of the features on EmoDB. We perform
experiments on this corpus using _leave-one-speaker-out_ (LOSO) cross
validation by keeping one speaker in test. Out of the remaining speakers, we
use two of them in validation and seven in training. We also apply five-fold
data augmentation by corrupting training set with additive noises and room
reverberation effect following the _Kaldi_ recipe111https://github.com/kaldi-
asr/kaldi/tree/master/egs/voxceleb/v2 for x-vector training [32].

We extract features from each speech utterance and discard the non-speech
frames with a simple energy-based _speech activity detector_ (SAD). We apply
utterance-level _cepstral mean variance normalization_ (CMVN) before creating
the training and validation samples with chunks of $100$ consecutive frames.
We consider multiple non-overlapping chunks from the speech utterances
depending on the length. We use LibROSA222https://librosa.github.io/ python
library for feature extraction.

We do not apply chunking for testing and consider the full utterance for
computing the test accuracy. We report the final performances with accuracy as
well as _unweighted average recall_ (UAR). The accuracy is computed as the
ratio between the number of correctly classified sentences to the total number
of sentences in test. The UAR is given as [38],

$\mathrm{UAR}=\frac{1}{K}\sum_{i=1}^{K}\frac{A_{ii}}{\sum_{j=1}^{K}A_{ij}}$
(2)

where $A$ refers to the contingency matrix, $A_{ij}$ corresponds to number of
samples in class $i$ classified into class $j$ and $K$ is the total number of
classes. As accuracy is considered _unintuitive_ for databases with uneven
samples across different classes, we optimize the feature extraction
parameters based on the UAR metric.

In DNN, we use ReLU activation function and batch normalization for all the
hidden layers. For regularization, we apply dropout with probability $0.3$ on
the FC layer only. We use Adam optimizer with learning rate $0.001$. The mini-
batch size is $64$. We train the models for $50$ epochs and finally testing is
done with the model which achieves the highest UAR on the validation set. We
repeat each experiment multiple times and report the average performance.

##  IV Results

###  IV-A Experiments on EmoDB

First, we conduct experiments on EmoDB and optimize the CQT parameters. We
vary the number of bins per octave from $1$ to $96$. We also perform the
experiments with three different hop lengths: $64$, $128$, and $192$. The top
row of Fig. 2 shows the standard accuracy and UAR for CQT. We observe improved
performance for lower bins per octave and lower hop length. The performance
remains very similar for bins per octaves between $2$ and $5$. We select $3$
bins per octave as the optimum observing the consistency in different runs of
the experiment. We fix the hop size $64$ as the optimum since the performance
is consistently better with this hop size, especially, for lower bins per
octave. Since the optimized CQT features use $24$ filters and hop length as
$64$, we apply similar configuration for STFT-based _mel frequency cepstral
coefficients_ (MFCCs) as well as _mel frequency spectral coefficients_ (MFSCs)
(i.e., MFCCs without DCT). The SER performances with CQT and STFT-based
features are illustrated as a bar plot in Fig. 2. We observe that CQT
coefficients as well as CQCCs consistently outperform MFCCs and MFSCs. We also
notice that the optimized MFSC outperforms baseline MFSC. The DCT slightly
degrades performance in both CQT and STFT-based approaches. We chose the best
configuration for both features for the remaining experiments.

TABLE III: Cross-corpora results shown in accuracy / UAR. This uses optimized configuration of MFSC and CQT. To have equal number of classes, only four emotions (Happy, Angry, Sad and Neutral) are considered from every database here. All other parameter settings remain similar to other experiments. Train on | Test on | MFSC | CQT
---|---|---|---
EmoDB | RAVDESS | 0.41 / 0.44 | 0.44 / 0.46
IEMOCAP | 0.36 / 0.37 | 0.38 / 0.39
RAVDESS | EmoDB | 0.45 / 0.42 | 0.48 / 0.48
IEMOCAP | 0.30 / 0.32 | 0.32 / 0.34
IEMOCAP | EmoDB | 0.64 / 0.50 | 0.63 / 0.50
RAVDESS | 0.38 / 0.39 | 0.38 / 0.39

Figure 3 shows the confusion matrices obtained for CQT and MFSC in experiments
with EmoDB. We observe that CQT is better capable of discriminating emotions
such as Fear, Disgust, Sad, Anger and Neutral as compared to MFSC. The CQT-
based system yields improved accuracy for Sad, Neutral and Disgust because
those emotions are more prominent in low-frequency regions. Performance of
Boredom is slightly degraded. Among all the seven emotions, Fear shows the
highest gain in performance over MFSC. Happy shows the lowest classification
accuracy and a high confusion with Angry.

###  IV-B Cross-corpora evaluation

Table III shows the performance obtained after cross corpus testing. The
optimized CQT shows better performance than optimized MFSC for most cases
except when the train-test pair are IEMOCAP-EmoDB and IEMOCAP-RAVDESS. The
obtained results consolidate our hypothesis that CQT helps in better capturing
of emotion-dependent information leading to better generality across
databases.

##  V Discussion and Conclusion

We notice that increasing the frequency resolution at lower frequency regions
led to substantial improvement in SER performance. This also confirms that
low-frequency region containing pitch harmonics and lower formants convey
important emotion-specific information. At the same time, the CQT with lower
high-frequency resolution does not degrade the overall SER performance which
indicates that high-frequency regions are less important from emotion
perspective. Also, better performance with fewer frequency bins in both CQT
and MFSC indicates less redundant time-frequency representation is more
effective for emotion discrimination. Though STFT with optimized parameters
generates spectrograms with higher frequency resolution, the performance
degrades most likely due to increased redundancy caused by capturing details
of high-frequency region. Cross-corpora evaluation suggests that CQT-based
time-frequency representation provides better generalization for SER task with
different speech corpora in training and test.

We conclude that CQT is a better choice of time-frequency representation in
terms of both recognition performance and generalization ability. However, the
SER performance is still poor for real-world deployment. We also gain no
improvement over MFSC for all the seven emotions included in EmoDB corpus.
This indicates that the time-frequency representation needs further
investigation for SER. This work can also be extended by exploring CQT
representation with recurrent architecture and attention mechanisms which are
lacking within our TDNN framework but found useful for SER.

## References

  * [1] Mehmet Berkehan Akçay and Kaya Oğuz,  “Speech emotion recognition: Emotional models, databases, features, preprocessing methods, supporting modalities, and classifiers,”  Speech Communication, vol. 116, 2020.
  * [2] Sreenivasa Rao Krothapalli and Shashidhar G Koolagudi,  “Speech emotion recognition: a review,”  in Emotion recognition using speech features, pp. 15–34. Springer, 2013.
  * [3] Moataz El Ayadi, Mohamed S Kamel, and Fakhri Karray,  “Survey on speech emotion recognition: Features, classification schemes, and databases,”  Pattern Recognition, vol. 44, no. 3, pp. 572–587, 2011.
  * [4] Jesús B Alonso, Josué Cabrera, Manuel Medina, and Carlos M Travieso,  “New approach in quantification of emotional intensity from the speech signal: emotional temperature,”  Expert Systems with Applications, vol. 42, no. 24, pp. 9554–9564, 2015.
  * [5] Rosalind W Picard,  Affective computing,  MIT press, 2000.
  * [6] Rosalind W Picard,  “Affective computing: challenges,”  International Journal of Human-Computer Studies, vol. 59, no. 1-2, pp. 55–64, 2003.
  * [7] Lijiang Chen, Xia Mao, Yuli Xue, and Lee Lung Cheng,  “Speech emotion recognition: Features and classification models,”  Digital Signal Processing, vol. 22, no. 6, pp. 1154–1160, 2012\.
  * [8] Florian Eyben, Anton Batliner, and Bjoern Schuller,  “Towards a standard set of acoustic features for the processing of emotion in speech.,”  in Proceedings of Meetings on Acoustics 159ASA. Acoustical Society of America, 2010, vol. 9, p. 060006.
  * [9] Florian Eyben et al.,  “The Geneva minimalistic acoustic parameter set (GeMAPS) for voice research and affective computing,”  IEEE Transactions on Affective Computing, vol. 7, no. 2, pp. 190–202, 2015.
  * [10] M. Kockmann, L. Burget, and J.H. C̆ernocký,  “Application of speaker-and language identification state-of-the-art techniques for emotion recognition,”  Speech Communication, vol. 53, no. 9-10, pp. 1172–1185, 2011.
  * [11] Carl E Williams and Kenneth N Stevens,  “Emotions and speech: Some acoustical correlates,”  The Journal of the Acoustical Society of America, vol. 52, no. 4B, pp. 1238–1250, 1972.
  * [12] Rainer Banse and Klaus R Scherer,  “Acoustic profiles in vocal emotion expression.,”  Journal of Personality and Social Psychology, vol. 70, no. 3, pp. 614, 1996.
  * [13] Roddy Cowie and Ellen Douglas-Cowie,  “Automatic statistical analysis of the signal and prosodic signs of emotion in speech,”  in Proc. ICSLP. IEEE, 1996, vol. 3, pp. 1989–1992.
  * [14] Suman Deb and Samarendra Dandapat,  “Multiscale amplitude feature and significance of enhanced vocal tract information for emotion classification,”  IEEE Transactions on Cybernetics, vol. 49, no. 3, pp. 802–815, 2018\.
  * [15] Sahar E Bou-Ghazale and John HL Hansen,  “A comparative study of traditional and newly proposed features for recognition of speech under stress,”  IEEE Transactions on Speech and Audio Processing, vol. 8, no. 4, pp. 429–442, 2000.
  * [16] Daniel Joseph France, Richard G Shiavi, Stephen Silverman, Marilyn Silverman, and M Wilkes,  “Acoustical properties of speech as indicators of depression and suicidal risk,”  IEEE Transactions on Biomedical Engineering, vol. 47, no. 7, pp. 829–837, 2000.
  * [17] Martijn Goudbeek, Jean Philippe Goldman, and Klaus R Scherer,  “Emotion dimensions and formant position,”  in Proc. INTERPSEECH, 2009, pp. 1575–1578.
  * [18] Elif Bozkurt, Engin Erzin, Cigdem Eroglu Erdem, and A Tanju Erdem,  “Formant position based weighted spectral features for emotion recognition,”  Speech Communication, vol. 53, no. 9-10, pp. 1186–1197, 2011.
  * [19] Margaret Lech, Melissa Stolar, Robert Bolia, and Michael Skinner,  “Amplitude-frequency analysis of emotional speech using transfer learning and classification of spectrogram images,”  Advances in Science, Technology and Engineering Systems Journal, vol. 3, pp. 363–371, 2018.
  * [20] Shiqing Zhang, Shiliang Zhang, Tiejun Huang, and Wen Gao,  “Speech emotion recognition using deep convolutional neural network and discriminant temporal pyramid matching,”  IEEE Transactions on Multimedia, vol. 20, no. 6, pp. 1576–1590, 2017.
  * [21] Qirong Mao, Ming Dong, Zhengwei Huang, and Yongzhao Zhan,  “Learning salient features for speech emotion recognition using convolutional neural networks,”  IEEE Transactions on Multimedia, vol. 16, no. 8, pp. 2203–2213, 2014.
  * [22] Che-Wei Huang and Shrikanth Narayanan,  “Characterizing types of convolution in deep convolutional recurrent neural networks for robust speech emotion recognition,”  arXiv preprint arXiv:1706.02901, 2017.
  * [23] Jack Parry, Dimitri Palaz, Georgia Clarke, Pauline Lecomte, Rebecca Mead, Michael Berger, and Gregor Hofer,  “Analysis of deep learning architectures for cross-corpus speech emotion recognition,”  Proc. INTERSPEECH, pp. 1656–1660, 2019.
  * [24] Judith C Brown,  “Calculation of a constant Q spectral transform,”  The Journal of the Acoustical Society of America, vol. 89, no. 1, pp. 425–434, 1991.
  * [25] Massimiliano Todisco, Héctor Delgado, and Nicholas Evans,  “Constant Q cepstral coefficients: A spoofing countermeasure for automatic speaker verification,”  Computer Speech & Language, vol. 45, pp. 516–535, 2017.
  * [26] Monisankha Pal, Dipjyoti Paul, and Goutam Saha,  “Synthetic speech detection using fundamental frequency variation and spectral features,”  Computer Speech & Language, vol. 48, pp. 31–50, 2018.
  * [27] Héctor Delgado et al.,  “Further optimisations of constant Q cepstral processing for integrated utterance and text-dependent speaker verification,”  in Proc. IEEE SLT. IEEE, 2016, pp. 179–185.
  * [28] T. Lidy and A. Schindler,  “CQT-based convolutional neural networks for audio scene classification,”  in Proc. the Detection and Classification of Acoustic Scenes and Events 2016 Workshop (DCASE2016), 2016, vol. 90, pp. 1032–1048.
  * [29] Dengke Tang, Junlin Zeng, and Ming Li,  “An end-to-end deep learning framework for speech emotion recognition of atypical individuals,”  in Proc. INTERSPEECH, 2018, pp. 162–166.
  * [30] Dias Issa, M Fatih Demirci, and Adnan Yazici,  “Speech emotion recognition with deep convolutional neural networks,”  Biomedical Signal Processing and Control, vol. 59, pp. 101894, 2020\.
  * [31] A. Waibel, T. Hanazawa, G. Hinton, K. Shikano, and K. J. Lang,  “Phoneme recognition using time-delay neural networks,”  IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 37, no. 3, pp. 328–339, 1989.
  * [32] David Snyder, Daniel Garcia-Romero, Gregory Sell, Daniel Povey, and Sanjeev Khudanpur,  “X-vectors: Robust DNN embeddings for speaker recognition,”  in Proc. ICASSP. IEEE, 2018, pp. 5329–5333.
  * [33] P. Shegokar and P. Sircar,  “Continuous wavelet transform based speech emotion recognition,”  in Proc. ICSPCS, 2016, pp. 1–8.
  * [34] Simon Nicholson, Ben Milner, and Stephen Cox,  “Evaluating feature set performance using the F-ratio and J-measures,”  in Proc. EUROSPEECH, 1997, pp. 413–416.
  * [35] Felix Burkhardt, Astrid Paeschke, Miriam Rolfes, Walter F Sendlmeier, and Benjamin Weiss,  “A database of German emotional speech,”  in Proc. INTERSPEECH, 2005, pp. 1517–1520.
  * [36] Steven R Livingstone and Frank A Russo,  “The Ryerson audio-visual database of emotional speech and song (RAVDESS): A dynamic, multimodal set of facial and vocal expressions in North American English,”  PLOS One, vol. 13, no. 5, 2018.
  * [37] Carlos Busso et al.,  “IEMOCAP: Interactive emotional dyadic motion capture database,”  Language Resources and Evaluation, vol. 42, no. 4, pp. 335, 2008\.
  * [38] Andrew Rosenberg,  “Classifying skewed data: Importance weighting to optimize average recall,”  in Proc. INTERSPEECH, 2012, pp. 2242–2245.