The Formation and Evolution of Star Clusters in Interacting Galaxies

Gadget (Springel et al. 2001; Springel 2005) is a massively parallel N-body/SPH code. It handles the components of a galaxy in two distinct ways: it treats the motions and evolution of dark matter and stars as collisionless particles in an N-body problem, while the gas is dealt with using the SPH method (Gingold & Monaghan 1977; Hernquist & Katz 1989). The N-body particles are described by the collisionless Boltzmann and Poisson equations, and the hydrodynamics of the fluid is followed using properties of neighboring gas particles smoothed by a kernel function. The gravitational force of each particle is calculated with a tree algorithm in which particles are grouped together and their effect is taken as a single multipole force, which reduces the computation cost greatly to O(N log N) compared to the direct summation of each particle pair with complexity O(N²). In this code, an artificial viscosity term is introduced into the equation of motion of SPH to represent the viscosity that often arises in ideal gases due to shocks caused by microphysics. The Gadget-2 we use explicitly conserves energy and entropy in the SPH formulation (Springel & Hernquist 2002). This version and its variants have been widely used in a large number of applications, from large-scale cosmological simulations (e.g., Springel & Hernquist 2003b; Springel et al. 2005b; Feng et al. 2013; Schaye et al. 2015) to galaxy mergers (e.g., Springel 2000; Li et al. 2004; Hopkins et al. 2005; Springel & Hernquist 2005; Hopkins et al. 2006; Li et al. 2007; Cox & Loeb 2008; Hopkins et al. 2008; Hayward et al. 2014).

Gizmo (Hopkins 2015) is a new Lagrangian code developed to circumvent the many problems encountered by SPH methods (Agertz et al. 2007; Bauer & Springel 2012; Kereš et al. 2012; Sijacki et al. 2012; Vogelsberger et al. 2012; Nelson et al. 2013; Zhu et al. 2015; Zhu & Li 2016). It derives the hydrodynamic equations using a kernel function to partition the volume and a Riemann solver to evolve the equations at the Lagrangian face comoving with the mass. Gizmo implements strict conservation of mass, energy, and linear and angular momentum, and it does not require any artificial diffusion terms to deal with shocks, embodying the advantages of both SPH and grid-based methods. It captures the instabilities of fluid mixing well, greatly reduces numerical noise and artificial viscosity, and as a result calculates fluid physics at smaller Mach numbers more accurately. Gizmo treats the contact discontinuities and shocks more precisely and more efficiently, generally within one kernel length instead of two to three as in Gadget, and it does not have the zeroth-order and first-order errors that are present in SPH (Zhu et al. 2015), so it can attain higher accuracy with a much smaller number of neighbors, which results in a faster convergence. We have used the meshless finite-mass mode of GIZMO for our project. The mass of an individual gas element is conserved in this mode, which allows us to trivially trace the history of star particles to their progenitor gas particles (otherwise, one needs tracer particles to do so).

A detailed comparison between Gadget and Gizmo in galaxy simulations has been conducted by Zhu & Li (2016), who showed a general agreement between the two simulations, but there were notable differences in a number of galaxy properties, such as star formation history, gas fraction, and disk structures.

In this study, our motivation for using these two codes to perform the same merger simulation is to reduce the possibility of numerical artifacts affecting our results. As we will show in Section 3, although the detailed star formation history of the mergers is different between the two simulations, the overall distribution functions of the cluster mass and the precluster gas pressure agree well, which suggests that our results are physical and robust, because we can argue against the effects of numerical artifacts such as the number of neighbors and artificial viscosity on these results.

Our simulations consider major mergers of two equal-mass Milky Way–sized galaxies. The galaxy is constructed using the model of Mo et al. (1998), which has a dark matter halo with a Hernquist density profile (Hernquist 1990), a thin disk with gas and stars, and a black hole in the center. The galaxy properties are similar to those of the Milky Way: the total mass is ${10}^{12}\,{M}_{\odot }$, the gas fraction ${f}_{\mathrm{gas}}=0.2$, the radial scale length of the galaxy is 3 kpc, and the disk height is one-fifth of it. The disk contains 4% of its total mass, and the seed mass of the black hole is ${10}^{5}\,{M}_{\odot }$.

The ISM in these galaxies is modeled using a subgrid multiphase recipe, and the SFR follows the empirical Schmidt–Kennicutt law (Schmidt 1959; Kennicutt 1998) where the surface density of star formation is related to the surface density of gas (${{\rm{\Sigma }}}_{\mathrm{SFR}}\propto {{\rm{\Sigma }}}_{\mathrm{gas}}^{1.4}$). The radiative cooling and heating in the ISM are modeled with the assumption that the medium is in collisional equilibrium and there is an external UV background (Haardt & Madau 1996). We have also followed feedback processes from both supernovae (Springel & Hernquist 2003a) and active galactic nuclei (Springel et al. 2005a). Supernova feedback includes thermal energy and galactic winds. The wind energy efficiency is 5% of the supernova energy, and the wind direction is anisotropic: winds carry energy and matter perpendicular to the disk plane. The feedback from the black holes is in the form of thermal energy deposited isotropically into the surrounding gas.

We note that Renaud et al. (2015) have performed a simulation of an Antennae-like merger. In order to explore a more extreme merger environment, we simulate a head-on collision of two Milky Way-sized galaxies, in which the progenitors are initially placed on a parabolic orbit with the inclination of both with respect to the orbital plane as $\theta =0$ and $\phi =0$. In the simulations, each galaxy is initially started with 82,000 gas particles, 328,000 star particles, and 1,476,860 dark matter particles, which yield a mass resolution of $5\times {10}^{4}\,{M}_{\odot }$ for the gas and star particles and $6\times {10}^{5}\,{M}_{\odot }$ for each dark matter particle.

Finding groups or structures in a given set of data is a classic problem in data mining. There are many algorithms for group finding, which essentially differ in their notion of groups and in their methods. The two main classes are particle-based and density-based algorithms. The most widely used method of group finding in astronomy is the FOF algorithm (Davis et al. 1985). It is a particle-based algorithm where all particles within a given linking length are considered as a group. However, there are two significant downsides of this approach even with an adaptive linking length (Suginohara & Suto 1992): (1) if two groups have a linking bridge, they will be identified as one group; and (2) it cannot identify substructures within a structure.

The other class of group-finding algorithms is density-based, which identify overdensities in the field as groups (Bertschinger & Gelb 1991; Warren et al. 1992; Gill et al. 2004). These methods do not suffer from the problems of the FOF methods described above. For this reason, we adopt the density-based hierarchical group-finding algorithm AHF (Knollmann & Knebe 2009) to identify SCs in the simulations using the following procedures.

First, AHF divides the simulation box into grid regions. It determines the density inside each grid cell and compares the density to a threshold or background density value. If the computed density exceeds the threshold, it divides the grid into half of its initial size. It computes the densities in each of the refined grid cells and again compares them with the threshold. This process goes on recursively until all of the cells in the simulation box have densities less than the threshold value. Next, it starts from the finest grid and marks isolated overdense regions as possible clusters. It goes on to the next coarser level and again identifies possible regions as clusters. Importantly, it links the possible clusters in finer grids to their respective coarser parts (linking daughters to parents). This continues until it has reached the coarsest grid, and finally it builds a tree of clusters with subclusters.

Considering the observed physical properties of YMCs, we impose a few criteria on the groups identified by AHF to qualify them as SCs. Each group should be gravitationally bound and have no substructures in order to only include individual SCs. It should contain stars and have at least four baryonic particles, which means the minimum cluster mass is $2\times {10}^{5}\,{M}_{\odot }$. The upper limit of group mass is set at ${10}^{8}\,{M}_{\odot }$, and the baryonic fraction (mass ratio of baryons to dark matter) of each group should be larger than unity in order to distinguish the SCs from dwarf galaxies that may form in our simulations. Such an approach provides a holistic identification of SCs in our study.

In the simulations, the progenitor galaxies start moving toward each other at t = 0 on a parabolic orbit. They have their first close encounter at $t\sim 0.23\,\mathrm{Gyr}$ and the second one at $t\sim 0.9$ Gyr until they finally coalesce at $t\sim 1\,\mathrm{Gyr}$. During the close passages, vigorous star formation is triggered by the compression of gas by tidal forces and rising gas densities in the inner region of the galaxies due to gravitational torques (Barnes & Hernquist 1991, 1996). As shown in Figure 1, the first starburst occurs at $t\sim 0.2\mbox{--}0.5\,\mathrm{Gyr}$, when the SFR increases by nearly two orders of magnitude and reaches a peak of $\sim {10}^{3}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ at $t\sim 0.4\,\mathrm{Gyr}$. The second starburst takes place at $t\sim 0.9\mbox{--}1.1\,\mathrm{Gyr}$, and the SFR peaks at $\sim {10}^{2}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ at $t\sim 1\,\mathrm{Gyr}$.

Zoom In Zoom Out Reset image size

The star formation history of galaxy mergers depends strongly on the progenitor properties and orbital parameters. The SFR in our simulation is higher than typical mergers at the local universe. Renaud et al. (2015) estimated from their simulation that the SFR of the Antennae merger at its starburst phase is $\sim {10}^{2}{M}_{\odot }$. However, star formation in ultraluminous infrared galaxies (ULIRGs) in the nearby universe can have comparable intensity. For example, radio recombination line studies of merger-driven starburst galaxy Arp 220 (77 Mpc away) suggest a mean SFR of $\sim 240\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ or more plausibly short periods of intense starbursts with an SFR of $\sim {10}^{3}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ (Anantharamaiah et al. 2000; Thrall et al. 2008; Varenius et al. 2016). We note that the mass of Arp 220 is estimated as $\sim {10}^{10}\,{M}_{\odot }$ (Scoville et al. 1997), much lower than our modeled galaxies. The high SFR in our simulated galaxies may be a product of both their high-mass progenitors and their extreme orbital parameters with the head-on collision.

As demonstrated in Figure 1, there is a remarkable difference in the star formation histories between the two simulations in that Gizmo produces higher SFR peaks than Gadget by a factor of 3–5. This is due to the more accurate treatments of fluids and shocks in Gizmo. Similar differences have also been seen in the code comparison study of galaxy mergers using Gadget and the moving-mesh code Arepo by Hayward et al. (2014), who reported that Arepo produces higher SFRs than Gadget by up to a factor of 10 for mergers of Milky Way-sized galaxies.

The strong compression and shocks produced by the galaxy interaction fuel rapid formation of SCs during the starbursts. As demonstrated in Figure 2, most of the SCs form in the nuclear regions of the two merging galaxies, with a few spread in the tidal tails and the galactic bridges. Similar distributions of YMCs have also been observed in galaxy mergers, including nuclear region clusters by Whitmore & Schweizer (1995) and Miller et al. (1997), and tidal tail clusters by Barnes & Hernquist (1992), Knierman et al. (2003), Bastian et al. (2005), and Mullan et al. (2011).

Zoom In Zoom Out Reset image size

In order to investigate the triggering source of star formation during the merging process, we track the gas particles that form stars at the three star-formation peaks, minor peak 1a at 0.23 Gyr and major peak 1b at 0.4 Gyr during the first close passage, and major peak 2 at 1 Gyr during the final coalescence, as labeled in Figure 1. It was shown by Hernquist (1989) that the major star-formation episodes in galaxy mergers are marked by a rapid loss of the angular momentum of the star-forming gas driven by the gravitational torque. We follow the procedure of Barnes & Hernquist (1996) and calculate the gravitational torque, $\tau =r\times F$, exerted on these star-forming gas particles by the gas and stars in the same galaxy (internal torque), and by gas, stars, and halo particles of the other galaxy (external torque). As shown in Figure 3, the internal torque is higher than the external counterpart by orders of magnitude for all tracked star-forming gas particles during the galaxy interaction. Similar results have been reported by a number of theoretical studies of major mergers (e.g., Hernquist 1989; Mihos & Hernquist 1994; Barnes & Hernquist 1996; Hopkins et al. 2009), which show that the internal torque is the dominant source of torque that drives the loss of angular momentum in these interacting galaxies. The close encounter of the galaxies produces strong tidal forces that result in a nonaxisymmetric response in the galaxy disks. These forces deform the galaxy disks and form gaseous and stellar bars in the galaxies. These gas bars lead the stellar bars by a few degrees (Barnes & Hernquist 1991), which eventually produces a strong torque on the gas near the center that drives rapid gas inflow toward the nuclear region, resulting in a vigorous starburst. From Figure 2, the majority of clusters formed at the starburst phase are highly clustered in the center region of each galaxy, indicating their origin from the nuclear gas inflow, while the few clusters formed in the bridge and tidal tails may be triggered by tidal force, as suggested by Renaud et al. (2009) and Renaud et al. (2014).

Zoom In Zoom Out Reset image size

We note that the first major star-formation peak 1b (at 0.4 Gyr) takes place about 160 Myr after the first close pericentric passage at ∼0.24 Gyr. Similar time delays have been found in other simulations of galaxy mergers (e.g., Mihos & Hernquist 1994, 1996; Cox & Loeb 2008; Hayward et al. 2014), because of the timescale to build up the gas density driven by internal torque for star formation. In addition, we note that there is a minor star-formation peak of $\sim 15\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ at 1a (0.24 Gyr) preceding the major one of $\sim 800\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ at 1b (0.4 Gyr). We find that the ratio of external to internal torque peaks during the 1a phase, suggesting that external torque from the tidal force may contribute to the star formation as well. Studies by Renaud et al. (2009, 2014) have shown that tidal force during galaxy interaction may compress the gas and enhance the star formation.

The resulting mass functions of the SCs formed during the first close encounter are shown in Figure 4. Although the total number of clusters in the same snapshot differs between the Gizmo and Gadget simulations by a factor of ∼1.3, the range of $\sim 150\mbox{--}200$ is in good agreement with observations of galaxy mergers such as the Antennae (Whitmore et al. 1999; Larsen 2010). More interestingly, both simulations produce similar mass distributions, which resemble a peaked or a quasi-lognormal function, with a Gizmo peak around ${10}^{5.8-6}\,{M}_{\odot }$ and the Gadget peak at ${10}^{6-6.2}\,{M}_{\odot }$.

Zoom In Zoom Out Reset image size

Our mass functions do not show a purely declining power law, as suggested by many observations of YMCS (e.g., Zhang & Fall 1999, Bik et al. 2003; McCrady & Graham 2007; Fall & Chandar 2012). This could be due to the limited resolution in our simulation, as we cannot resolve clusters at mass lower than ${10}^{5}\,{M}_{\odot }$. However, Renaud et al. (2015) also reported lognormal-shape ICMFs from their Antennae simulation even though they have a much higher mass resolution ($\sim 70\,{M}_{\odot }$). It is also possible that the power-law phase is extremely short lived ($\lt 10$ Myr) in violent mergers and our snapshot interval (every 10 Myr) misses that phase. We note, however, some observations have suggested that the power-law ICMF for YMCs is not universal, but rather they have a turnover at low mass (Cresci et al. 2005; Anders et al. 2007).

The peak masses of these ICMFs are higher than those of the observed GCMFs, $\sim 1.5\mbox{--}3\times {10}^{5}\,{M}_{\odot }$, but it can be predicted that after an evolution of gigayears, stellar evolution, tidal stripping, and other disruptive processes will cause them to lose some of their mass (Kruijssen 2015; Webb & Leigh 2015). We discuss their evolution in more detail in Section 4. It is encouraging that both of our simulations, together with that by Renaud et al. (2015), produce similar quasi-lognormal ICMFs with similar peak positions, despite having vastly different hydrodynamic solvers, feedback processes, and numerical resolutions. This agreement suggests that the lognormal mass function is unlikely to be due to numerical artifacts but has a physical origin, which will be investigated in the next section.

In order to explore the physical origin of the quasi-lognormal ICMFs in Figure 4, we examine the physical conditions of cluster formation. As mentioned in Section 1, massive SCs form in molecular clouds with very high gas pressures (Elmegreen & Efremov 1997; Ashman & Zepf 2001). When a star-forming cloud is under high pressure, the efficiency of star formation increases and the gas velocity dispersion becomes higher, which in turn help to keep the cloud bound (Jog & Solomon 1992). These physical conditions create an ideal nursery to form massive, gravitationally bound clusters. Recently, Zubovas et al. (2014) performed an N-body simulation of a molecular cloud and concluded that high external pressures drive efficient star formation and can cause cloud fragmentation, leading to the formation of SCs. Quantitatively, the external pressure P on a nascent molecular cloud is given by Elmegreen (1989):

$$\begin{eqnarray}&&P=\displaystyle \frac{3{\rm{\Pi }}{M}_{\mathrm{cloud}}{\sigma }_{v}^{2}}{4\pi {r}^{3}},\end{eqnarray} \tag{ 1 }$$

where ${M}_{\mathrm{cloud}}$ is the mass of the cloud, ${\sigma }_{v}$ is the velocity dispersion, and r is the size of the cloud. The factor Π is given by the ratio of the density at the cloud edge and the average density, ${\rm{\Pi }}={n}_{e}/\langle {n}_{e}\rangle$. This Π ratio is dependent on the density profile of the parent molecular cloud. Locally the probability distribution function of ISM density due to turbulence can be approximated as lognormal, but at high-density regions ($\gt {10}^{3}\,{\mathrm{cm}}^{-3}$), such as at the center of a molecular cloud, the density profile can develop a power-law tail. At these dense places, the Π value can be $\gg 1$, and consequently it can increase the amount of gas above a density threshold that facilitates further star formation (Elmegreen 2011; Renaud et al. 2014).

Elmegreen & Efremov (1997) estimated that the pressure in the birth clouds of typical GC progenitors or SSCs is $\gtrsim {10}^{8}\,{\rm{K}}\,{\mathrm{cm}}^{-3}$, which is $\gtrsim {10}^{4}$ times higher than the ISM pressure in the Milky Way (Jenkins et al. 1983; Boulares & Cox 1990; Welty et al. 2016).

In the simulations, in order to determine the pressure of the clouds from which the SCs form, we track the cluster members back in time. We take the constituent star particles of a cluster and identify the gas particles from which they formed. We then measure the velocity dispersion of these gas particles and approximate the gas cloud radius as the average distance from gas particles to the center of mass of the cloud. Due to the limited spatial and mass resolutions, we cannot directly probe the cloud density profile for estimating Π, so we approximate ${\rm{\Pi }}=0.5$ following Johnson et al. (2015). We note that realistically Π can be higher, but it would increase all our pressure measurements similarly. With these parameters, we can then calculate the cloud pressure using Equation (1). In Figure 5, we show the resulting pressure distributions against the cluster mass distributions, in comparison with the pressure of a pre-SSC cloud observed in the Antennae by Johnson et al. (2015). We also calculate the pressures of observed SCs in the Antennae Galaxies using the velocity and radius data compiled in Portegies Zwart et al. (2010), Mengel et al. (2002), and Mengel et al. (2008) for comparison.

Zoom In Zoom Out Reset image size

As shown in Figure 5, the precluster cloud pressures from both the Gizmo and Gadget simulations fall in the range of $P/k\sim {10}^{8}\mbox{--}{10}^{12}\,{\rm{K}}\ {\mathrm{cm}}^{-3}$, in good agreement with observations of a proto-SSC cloud in Antennae (Johnson et al. 2015), but they are ${10}^{4}\mbox{--}{10}^{8}$ times higher than the typical pressure in the ISM (Jenkins et al. 1983; Boulares & Cox 1990; Welty et al. 2016). Our results support the theoretical expectations that massive SCs form in high-pressure clouds.

Moreover, the pressure distributions have near-lognormal-shape profiles in all panels. Such a quasi-lognormal-shape pressure distribution may be the cause of the quasi-lognormal-shape ICMF. If we assume a certain cluster-formation efficiency η (η varies with galactic environment, from 0.01 in quiescent galaxies to $\gt 0.4$ in interacting galaxies, as suggested by Goddard et al. 2010; Kruijssen 2015), then the mass of an SC ${M}_{\mathrm{cluster}}$ may be related to that of the birth clouds ${M}_{\mathrm{cloud}}$ as ${M}_{\mathrm{cluster}}\propto \eta {M}_{\mathrm{cloud}};$ then by inverting Equation (1) we get ${M}_{\mathrm{cluster}}\propto \eta {M}_{\mathrm{cloud}}\propto \eta P$. This qualitative relation, as can be inferred from Figure 5, suggests that the distribution of cluster mass depends on that of the cloud pressure, so a lognormal pressure distribution may lead to a lognormal mass distribution of the resulting clusters.

Furthermore, Figure 5 also shows that the pressure profiles have a peak at cluster mass around ${10}^{5.8-6}\,{M}_{\odot }$ in the Gizmo simulation and around ${10}^{6-6.2}\,{M}_{\odot }$ in the Gadget simulation, which are exactly the same peaks of their corresponding mass functions. This striking correlation may explain the preferred mass range around the characteristic peak ${10}^{6}\,{M}_{\odot }$, as determined by the cloud pressure, for the newly formed SCs in galaxy mergers.

In order to understand the possible physical processes behind the high pressure of the precluster clouds in our simulations, we compare the evolution of the total velocity dispersion of gas and stars in the central region of one galaxy (since the two merging galaxies are identical), that is, within 5 kpc from the galactic center, and that of the entire simulation box, as shown in Figure 6. We find that the peaks of the total velocity dispersion correspond to those of the star formation as in Figure 1, and that the velocity dispersion of the galactic central region is higher than that of the whole box during the major starburst phases at 0.4 and 1 Gyr. These results suggest that a high velocity dispersion around the galaxy center, which leads to high pressure and a strong circumnuclear starburst, may result from the same mechanism, compressive shocks driven by gravitational torques during galaxy merger.

Zoom In Zoom Out Reset image size

Our theoretical findings may provide an explanation for a number of observations. In addition to the recent observations of high pressure in a proto-SSC cloud in Antennae (Johnson et al. 2015), measurements of molecular clouds in the Antennae galaxies have revealed very high velocity dispersion in high star-forming regions (Zhang et al. 2010), which can be explained by compressive shocks (Wei et al. 2012). Herrera et al. (2011) carried out near-infrared imaging spectroscopy of the same region and found extended line widths in ${{\rm{H}}}_{2}$ emission, which indicates powerful shocks in the region. Similarly, measurements of the CO emission from the starbursting merger of M81/M82 by Keto et al. (2005) suggest that the molecular clouds undergoing star formation are driven by shock compression. Theoretical studies (Jog & Solomon 1992; Ashman & Zepf 2001) have also shown that, during the galaxy encounters, the giant molecular clouds undergo significant shock compression, which leads to an increase in the cloud pressure.

We also note that during the $1a$ phase at 0.23 Gyr, the velocity dispersion of the central region is similar to that of the whole box, which suggests a spatially extended star formation probably influenced by tidal forces, as discussed in Section 3.1. Simulations of Antennae galaxies by Renaud et al. (2014, 2015) have shown that compressive tides during the galactic encounter can cause high star formation over extended volumes.

We can see from Figure 5 that the pressures of our simulated clouds are somewhat higher than that of the observed systems. This can arise from the fact that in our simulation the two equal-mass, Milky Way-sized galaxies collide mainly head on and merge violently, whereas the Antennae Galaxies have a smaller mass and are on a milder pericentric passage (Renaud et al. 2015). The extreme conditions in the simulated galaxies produce more powerful shocks, which in turn help increase the gas pressure. Such an extreme high-pressure environment may preferentially form massive SCs in a narrow mass range, as shown in Figure 4, which may help to explain why we do not see a power-law ICMF in the simulations.

Our simulations bridge the observations and theories of cluster formation, and we confirm that massive SCs can form in gas-rich galaxy mergers due to the high gas pressure produced by strong gravitational interactions. Moreover, the special pressure range in such merger environments preferentially forms SCs in a narrow mass interval around the peak mass at $\sim {10}^{6}\,{M}_{\odot }$. Furthermore, the quasi-lognormal pressure distribution may lead to the quasi-lognormal ICMF of SCs formed in colliding galaxies. Our results, therefore, provide clues to the formation of GCs and their universal lognormal mass functions, which will be explored in the next section.