The formation of low metallicity globular clusters in dwarf galaxy mergers

TL;DR

This study uses detailed hydro-dynamical simulations to explore how low-metallicity globular clusters form rapidly in dwarf galaxy mergers, highlighting the physical conditions and processes involved.

Contribution

It presents a novel high-resolution simulation model that captures the formation, feedback, and early evolution of star clusters in dwarf galaxy mergers, linking to observed globular clusters.

Findings

Clusters form rapidly within 6-8 Myr in dense gas flows.

Star formation rates exceed 10 M_sun/yr per kpc^2 during formation.

Gas expulsion by supernovae terminates cluster formation.

Abstract

We present a hydro-dynamical simulation at sub-parsec and few solar mass resolution of a merger between two gas-rich dwarf galaxies. Our simulation includes a detailed model for the multi-phase interstellar medium (ISM) and is able to follow the entire formation history of spatially resolved star clusters, including feedback from individual massive stars. Shortly after the merger we find a population of $\sim 900$ stellar clusters with masses above $10^{2.5}\; \rm{M_\odot}$ and a cluster mass function (CMF), which is well fitted with a power-law with a slope of $\alpha=-1.70\pm0.08$. We describe here in detail the formation of the three most massive clusters (M$_{*} \gtrsim 10^5$ M$_\odot$), which populate the high-mass end of the CMF. The simulated clusters form rapidly on a timescale of $6$-$8$ Myr in converging flows of dense gas. The embedded merger phase has extremely high star formation rate surface densities of $\Sigma_\mathrm{SFR}>10\; \mathrm{M}_\odot\; \mathrm{yr}^{-1}\; \mathrm{kpc}^{-2}$ and thermal gas pressures in excess of $P_{\rm th}\sim10^7 \; \mathrm{k}_{\rm B}\;(\rm K\;\mathrm{cm}^{-3})^{-1}$. The formation process is terminated by rapid gas expulsion driven by the first generation of supernovae, after which the cluster centers relax and both their structure and kinematics become indistinguishable from observed local globular clusters. The simulation presented here provides a general model for the formation of metal-poor globular clusters in chemically unevolved starbursting environments of low-mass dwarf galaxies, which are common at high redshifts.

Figures

Equations

$$M_J=\frac{{\pi }^{5/2}{c}_s^3}{6G^{3/2}{\rho }^{1/2}}$$

$$\Sigma (r)={\Sigma }_0{\left(1+\frac{r^2}{a^2}\right)}^{-\gamma /2}$$

Full text

galaxies: star clusters: general — galaxies: dwarf — galaxies: interactions — galaxies: ISM — methods: numerical

††journal: ApJL

1 Introduction

Globular clusters (GCs) are the densest gravitationally bound stellar systems in the Universe, and they are found in all types of galaxies, even low-mass dwarf galaxies (Harris, 1991). GCs have old ages ($\gtrsim 100$ Myr), with some of the oldest clusters being nearly as old as the Universe (Marín-Franch et al., 2009), and they have typical masses of several times $10^{5}$ M⊙, with typical effective radii of several parsecs. The elemental abundances of GCs range from $0.3\%$ solar to about solar metallicity and they show diverse chemical compositions indicative of multiple stellar populations (e.g. Bastian & Lardo 2018).

It is very challenging to resolve the formation sites of GCs at high redshifts with present-day telescopes, making their star formation environments largely unknown (e.g. Vanzella et al. 2017). However, the present-day structure of GCs suggests that extreme gas pressures, about $\sim 2-3$ orders of magnitude higher than in the present-day Milky Way, were required to make them at their high stellar densities and gravitational binding energy densities (Elmegreen & Efremov, 1997; Kruijssen, 2012).

Contrary to the high pressure requirement, present-day dwarf galaxies typically have the lowest interstellar pressures (e.g. Elmegreen & Hunter 2000). This contradiction has led to the idea that at least some metal-poor GCs formed during mergers of dwarf galaxies, where fast shock fronts compressed the gas and caused large-scale gravitational collapse into dense clouds (Brodie & Huchra, 1991).

Numerically the formation of massive stellar clusters have been studied using both isolated simulations of collapsing giant molecular clouds (e.g. Dale & Bonnell 2011; Howard et al. 2018), in mergers of both spiral (e.g. Bekki et al. 2002; Bournaud et al. 2008) and dwarf galaxies (e.g. Saitoh et al. 2010; Hopkins et al. 2013) and recently also in cosmological simulations (e.g. Kimm et al. 2016; Kim et al. 2018).

The pioneering simulations of Bekki et al. (2002) and Bournaud et al. (2008) resolved the forming stellar clusters with particle masses of $\sim 7000\;\mathrm{M}_{\odot}$ at a typical spatial resolution of a few hundred parsecs. More recently, improved numerical merger simulations by Hopkins et al. (2013) and Renaud et al. (2015) resolved masses down to $\sim 100\;\mathrm{M}_{\odot}$ at $\sim 1$ parsec scales, whereas current galaxy-scale simulations of isolated dwarf galaxies are even able to resolve the ISM down to sub-parsec scales and follow the feedback from individual stars (e.g. Hu et al. 2017; Su et al. 2018; Emerick et al. 2019).

The formation of a GCs is fundamentally a two-step process (Kruijssen, 2014). Firstly, enough mass must be accumulated to form a massive stellar cluster and secondly this cluster must survive for at least $\sim 100$ Myr and up to a Hubble time in order to be observable in the present-day Universe. In this Letter we concentrate on the formation phase of massive stellar clusters, and present a hydro-dynamical simulation of a dwarf galaxy merger resolved at stellar-mass and sub-parsec spatial resolution.

2 Simulation setup

2.1 Simulation code

The simulations were run using a modern smoothed particle hydrodynamics (SPH) solver SPHGal (Hu et al., 2014, 2016) based on the Gadget -code (Springel, 2005). The code employs a pressure-energy formulation of the hydrodynamic equations (Hopkins, 2013) using the Wendland C4 smoothing kernel (Dehnen & Aly, 2012) smoothed over 100 neighbors, and also includes well-tested implementations of artificial viscosity and thermal conduction in converging gas flows (Price, 2008).

In the temperature range $10\ \mathrm{K}<T<3\times 10^{4}\ \rm K$ gas is cooled using the model described in Hu et al. (2016), which is based on a non-equilibrium chemical network (Glover & Mac Low, 2007) that follows six chemical species (H2, H*+, H, CO, C+*, O), in addition to the free electron density. For the dust-to-gas mass ratio we assume a constant value of $0.1\%$. At temperatures above $T>3\times 10^{4}\ \rm K$, we adopt an equilibrium cooling model using a metallicity-dependent cooling function (Wiersma et al., 2009). Metallicity in the ISM is followed by tracing the mass of $12$ elements (H, He, N, C, O, Si, Mg, Fe, S, Ca, Ne, Zn) in both gas and star particles (Aumer et al., 2013), which evolve through stellar feedback. In this temperature regime the ISM is assumed to be optically thin and in ionisation equilibrium with a cosmic UV background (Haardt & Madau, 2001).

We follow the mass, metal (Chieffi & Limongi, 2004) and radiation (Hu et al., 2016, 2017) output into the ISM along the lifetime of each massive star. The stellar radiation includes the far-ultraviolet (FUV) radiation with photon energy below $13.6$ eV and the hydrogen-ionising radiation from massive stars. Shielding against the FUV radiation is accounted for by using the TreeCol algorithm (Clark et al., 2012).

The star formation model samples stellar masses from a Kroupa (Kroupa, 2001) initial mass function and statistically realises each individual star more massive than $4$ M⊙ as a stellar particle. If the sampled mass exceeds the available mass in the gas particle, the mass deficit is stochastically taken from nearby gas particles in order to satisfy mass conservation. Stars with masses below $4$ M⊙ are grouped into stellar population particles.

Star formation (SF) is taking place at a $2\%$ efficiency on a dynamical time-scale for gas particles which are present in a converging flow and have a Jeans mass

[TABLE]

below eight SPH kernel masses, where $G$ is the gravitational constant, $\rho$ is gas density, and $c_{s}$ is the local sound speed of the gas. In addition, we enforce instantaneous star formation when the Jeans mass falls below half a kernel mass, i.e. $\sim 200$ M⊙ (see Fig. 1).

At the end of their lifetimes, stars more massive than $8$ M⊙ undergo type II supernova (SNII) explosions with mass-dependent delay times (Hu et al., 2017), with the resulting mass and metal yields distributed within the weighted SPH kernel. At this low metallicity and dust content, feedback from supernovae and ionising radiation will dominate and the additional contribution from radiation pressure on dust (e.g. Howard et al. 2018) and from line driven winds (e.g. Puls et al. 2008) are expected to be weak and are therefore not included in our present model.

2.2 Initial conditions

The initial conditions consist of two identical dwarf galaxies setup at z=0 with virial radii of $44$ kpc including identical dark matter halos, stellar disks and extended gas disks. The dark matter halo is set up with a Hernquist profile with virial mass of $M_{\rm vir}=2\times 10^{10}$ M⊙, an NFW-equivalent (Springel et al., 2005) concentration parameter of $c=10$ and a Bullock spin-parameter of $\lambda=0.03$.

The rotationally supported stellar disk with mass $M_{*}=2\times 10^{7}\;\mathrm{M}_{\odot}$ has an exponential surface density distribution with a scale radius of $r_{\rm D,star}=0.73\ \rm kpc$ and a constant scale height of $h_{\rm star}=0.35\ \rm kpc$. The exponential gas disk with a mass of $M_{\rm gas}=4\times 10^{7}\;\mathrm{M}_{\odot}$ is set up with a scale radius of $r_{\rm D,gas}=1.46\ \rm kpc$ and a gas scale height initially set by the condition of hydrostatic equilibrium. The initial metallicity of the gas disk is set to a constant value of $0.1$ Z⊙.

The galaxies are oriented off the orbital plane with the inclination and argument of pericenter chosen as $\{i_{1},i_{2}\}=\{60^{\circ},60^{\circ}\}$ and $\{\omega_{1},\omega_{2}\}=\{30^{\circ},60^{\circ}\}$. The parabolic orbit of the dwarfs, with an initial distance of $5$ kpc and a pericentric separation of $1.46$ kpc, results in a relatively rapid coalescence time for the two gas disks. The first passage occurs $\sim 50$ Myr from the start of the simulation, followed by the second encounter $\sim 100$ Myr later, with the final coalescence of the gas disks occurring $\sim 20$ Myr after the second passage. The simulation was run for more than $100$ Myr past the coalescence for a total simulation time of $\sim 300$ Myr.

The gravitational forces are softened with a spline kernel with a softening length of $\epsilon=0.1\ \rm pc$ for stars and gas, and $\epsilon=62\ \rm pc$ for dark matter, respectively. The time evolution is resolved down to $\sim 10$ yr. The mass resolution for the stellar and gaseous particles is typically $4$ M⊙ per particle. The dark matter component has a mass resolution of $10^{4}$ M⊙ per particle.

Our simulation is limited by our gravitational softening length of $0.1$ pc and mass resolution of $\gtrsim 4$ M⊙. In this study we thus concentrate on the formation and early evolution of GCs. However, we stress that unsoftened simulations with direct N-body codes have shown that massive clusters similar to ours evolve into present-day GCs after $10$ Gyrs of evolution, even in the tidal field of the Milky Way (Wang et al., 2016).

3 Results

3.1 The formation of stellar clusters

When evolved in isolation, the model dwarf galaxies have a low SFR of $10^{-3}\;\mathrm{M}_{\odot}\;\mathrm{yr}^{-1}$, in agreement with the long gas depletion times observed locally in these types of systems (Hu et al., 2017). However, during the galaxy merger the global SFR increases by two orders of magnitude to SFR $=0.3$ M⊙ yr*-1*, with a concurrent increase in the fraction of stars forming in bound star clusters (see e.g. Kruijssen 2012).

In Fig. 1 we show the phase diagrams for all gas in the merging dwarf system at the time of peak SFR. The instantaneous SF threshold (solid line) corresponds to densities of $n_{\rm H}=10^{3.5}-10^{4.5}\ \rm cm^{-3}$, at temperatures in the range of $T=10-100\ \rm K$, with the most extreme values only reached during the merger.

We use the Subfind algorithm (Springel et al., 2001) to identify stellar clusters in the simulation. In this study we define clusters to be objects with more than $50$ bound stellar particles, which results in a minimum cluster mass of $M_{*}\sim 200$ M⊙.

After the coalescence of the galactic disks the merger remnant harbors a considerable population of more than 1000 young stellar clusters. The distribution of the $\sim 900$ clusters with masses above a typical observed cut-off of $M_{*}>10^{2.5}\;\mathrm{M}_{\odot}$ follows a mass function of the form $dN\propto M^{\alpha}dM$ with a best-fit power-law index of $\alpha=-1.70\pm 0.08$, as also seen in observed dwarf galaxies (e.g. Chandar et al. 2010). At the peak epoch of star formation we identify three compact regions with extremely high SFR surface densities exceeding $\Sigma_{\rm SFR}>10$ M⊙ yr*-1* kpc*-2*, comparable to values in massive high-redshift starburst galaxies (Tacconi et al., 2013). Each of these regions produces a massive dense star cluster on a timescale of only a few Myr, after which the clusters remain on orbits with apocenters of $\sim 500$ pc from the center of mass of the merger remnant, similarly to the globular clusters in local dwarf galaxies (Mackey & Gilmore, 2003).

We show in Fig. 2 the environment of the most massive cluster, where the average gas surface density exceeds $\Sigma_{\rm gas}\sim 10^{2}$ M⊙ pc*-2* with thermal pressures of $P_{\rm th}\sim 10^{5}$–$10^{7}$ k${}_{\mathrm{B}}\;(\rm K\;\mathrm{cm}^{-3})^{-1}$. The top panels depict the region when the cluster has assembled half of its final stellar mass, $5$ Myr before the formation is completed at time t0. The bottom panels show the time when the maximum thermal pressure is reached, approximately $2$ Myr later. At these times the entire formation region is still embedded in a dense gas filament. Over a total time span of $\sim 6$–$8$ Myr, nearby newly formed stellar clusters with a wide range of masses along the CMF merge together while embedded in a converging gas flow that funnels gas at a typical infall rate of $0.2\;\mathrm{M}_{\odot}\;\mathrm{yr}^{-1}$ into a region $50$ pc in radius.

The slow rotation of the galaxies allows for such a prolonged collapse phase in the absence of destructive shear forces. Additional compression is provided by a super-shell generated by the first massive cluster that formed. The central density of gas within the cluster formation region reaches $10^{3}$–$10^{4}$ M⊙ pc*-3* ( $n_{\rm H}>10^{4}$ cm*-3*) where the gas is efficiently transforming into stars on the local dynamical timescale of $\tau_{\mathrm{dyn}}\sim 10^{6}$ yr. The mean SFR in this region is $0.1$ M⊙ yr*-1* at a mean efficiency of $\epsilon_{\rm{SF}}\sim\dot{M}_{\rm*}/\dot{M}_{\rm infall}\sim 50\%$ with respect to the infalling gas, with peak values of up to $\epsilon_{\rm{SF}}\sim 90\%$.

3.2 The mass assembly of stellar clusters

Fig. 3 shows the rapid build-up of stellar density in the most massive cluster, where the newly formed stars concentrate at stellar densities of $\rho_{*}>10^{6}$ M⊙ pc*-3*. This maximum stellar density exceeds the density of the gas by more than two orders of magnitude and has a steeper radial density profile, in good agreement with what is observed in young massive clusters (Walker et al., 2016).

The central 10 half-mass radii (see Section 3.3) around the center of mass of each of the three most massive clusters retains a fairly constant gas infall rate of $0.05$–$0.2\;\mathrm{M}_{\odot}\;\mathrm{yr}^{-1}$ during their assembly, while on average $25\%$–$50\%$ of the infalling gas mass is converted into stars which mostly end up in the central cluster in each region.

Following the increasing SFR, the SNII rate increases as the cluster mass builds up. Typically, not more than $10\%$ of the supernova energy couples to the surrounding dense ISM even if the medium is ionised (Walch & Naab, 2015). However, at high enough supernova rates ($\sim 500$ Myr*-1*), even compensating for the short cooling times of $\tau_{\rm cool}\sim 2\times 10^{3}$ yr in the dense gas, the supernovae end up depositing as much as the gas binding energy of $4\times 10^{52}$ erg. As a result, all of the remaining gas is blown out of the cluster about $7$ Myr into its formation, soon after the epoch shown in the bottom panel of Fig. 2.

At this time the most massive cluster has reached a mass of $\sim 8\times 10^{5}$ M⊙ and has consumed $\sim 30\%$ of the available gas within the nearby $50$ pc region; the remaining gas has either been expelled ($63\%$) or converted into stars ($7\%$) mostly bound in smaller mass clusters ($5\%$). By the time the gas is expelled, about $80\%$ of the cluster stars have formed in-situ within $30$ pc, with the remaining $20\%$ of the stellar mass coming from smaller clusters and unbound stars that accreted from the surrounding $100$ pc neighborhood.

3.3 Comparison of cluster properties with observations

Once star formation is terminated, the cluster evolution is governed by its tidal environment and by collisional processes such as two-body relaxation. Even in a softened system, the number of particles is the main driver of the two-body relaxation timescale (Binney & Tremaine, 2008). At our gravitational softening of $\epsilon=0.1\ \rm pc$ we find that a constant density stellar core forms roughly on the local relaxation timescale of $t_{\rm relax}\sim 40$ Myr, which is within a factor of a few of the true relaxation time given our stellar-mass particle resolution.

In Fig. 4 we show the radial and 2D stellar surface density profiles of the most massive cluster $\sim 100$ Myr after assembly, together with the radial surface densities of the two next-most massive clusters, which have similar formation histories. Evolving the simulated clusters several Gyrs would further lower the central stellar density (Wang et al., 2016), bringing especially the lower-mass clusters closer to their observed counterparts.

A common model for fitting the surface brightness and mass profiles of especially young and intermediate-age stellar clusters, is the Elson-Fall-Freeman power-law (EFF, Elson et al. 1987; McLaughlin & van der Marel 2005) which is defined as

[TABLE]

and parametrised by the central surface density $\Sigma_{0}$, the scale radius $a$ and the power-law index $\gamma$. All three simulated clusters in Fig. 4 resemble the EFF profiles observed for similar mass clusters with ages in excess of $100$ Myr in the Local Group (McLaughlin & van der Marel, 2005), and they are well fitted by similar power-law exponents of $\gamma\sim-2.2$ to $\gamma\sim-2.5$ above the limit of $\Sigma_{*}\sim 10$ M⊙ pc*-2* set by the ambient stellar background.

In Fig. 5 we show how the velocity dispersions, half-mass radii, and central surface densities of the three model GCs are indistinguishable from massive stellar clusters in the dwarf galaxies Fornax and the Small and Large Magellanic Clouds (SMC and LMC), and the Milky Way (McLaughlin & van der Marel, 2005). Out of the cluster sample shown in Fig. 5, only the SMC (1 cluster) and LMC (16 clusters) harbor stellar clusters younger than 100 Myr, respectively. However, the properties of these younger clusters are inseparable from the other older clusters in the sample shown in Fig. 5. The average stellar metallicities of the simulated clusters are close to the initial value of $0.1\;Z_{\odot}$, indicating that the clusters have not been significantly enriched by star formation during their short formation times. The scatter in stellar ages for the simulated GCs is $4$ Myr or less, also in good agreement with the observations.

4 Conclusions

We have presented results from a detailed numerical simulation in which three dense massive stellar clusters form on timescales of several Myr in a dwarf galaxy merger with low metallicity. After the gas is expelled by supernovae, the clusters relax and appear similar to present-day GCs.

We propose that such a formation scenario results in GCs in chemically unevolved, low-mass galaxies (M${}_{*}<10^{9}$ M⊙) at high redshift ($z>3.5$) (e.g. Elmegreen et al. 2012). The peak SFR of $0.3$ M⊙ yr*-1* corresponds to a relatively low H$\alpha$ luminosity (Hunter et al., 2010) of L${}_{\rm{H}_{\alpha}}=4.3\times 10^{40}$ erg s*-1*. The most massive clusters can either survive to the present-day in these dwarfs, where tidal forces are weak or, alternatively, they can be accreted onto more massive systems such as the Milky Way (Muratov & Gnedin, 2010). Finally, this formation scenario is also consistent with observations of young massive clusters in the local Universe, which are thought to be the counterparts of old globular clusters (Portegies Zwart et al., 2010).

The computations were carried out at CSC – IT Center for Science Ltd. in Finland and at Max-Planck Institute for Astrophysics in Germany. N.L. acknowledges the financial support by the Jenny and Antti Wihuri Foundation. N.L. and P.H.J. acknowledges support by the Academy of Finland grant 274931 and the European Research Council via ERC Consolidator Grant KETJU (no. 818930). C.-Y.H. acknowledges The Center for Computational Astrophysics, supported by the Simons Foundation. S.W. acknowledges support by the European Research Council via ERC Starting Grant RADFEEDBACK (no. 679852) and by the German Science Foundation via CRC956, Project C5.