Major mergers of galaxies are likely to trigger bursty star formation activities. Usually, the accumulation of dense gas and the boost of star formation efficiency (SFE) are considered to be the two main drivers of the starbursts. However, it is still unclear how each process operates on the scale of individual star-forming clouds. Here, we present a high-resolution (\(2\,\text{M}_{\astrosun}\)) radiation-hydrodynamic simulation of a gas-rich dwarf galaxy merger using the Realistic ISM modeling in Galaxy Evolution and Lifecycles (RIGEL) model to investigate how mergers affect the properties of the structure of dense star-forming gas and the cloud-scale SFE. With the unprecedented mass and temporal resolution of the simulations, we track the evolution of sub-virial dense clouds in the simulation by mapping them across successive snapshots spanning 200 Myr taken at intervals of 0.2 Myr. We found that the merger triggers a \(130\) times higher star formation rate (SFR) and shortens the galaxy-wide gas-depletion time by two orders of magnitude compared to those of two matched isolated galaxies. However, the depletion time of individual clouds and their lifetime distribution remained unchanged over the simulation period. The cloud life cycles and cloud-scale SFE are determined by the local stellar feedback rather than environmental factors like the tidal field regardless of the merger process, and the integrated SFE (\(\epsilon_\text{int}\)) of clouds in complex environments is still well described by an \(\epsilon_\text{int}\)–\(\Sigma_\text{tot}\) relation found in idealized isolated-cloud experiments. During the peak of the starburst, the median cloud-scale integrated SFE only changed by 0.17–0.33 dex lower compared to the value when the two galaxies are not interacting. The merger boosts the SFR primarily through the accumulation and compression of dense gas fueling star formation. Strong tidal torques assemble \(\gtrsim10^{5}\,\text{M}_{\astrosun}\) clouds that seed massive stellar clusters. The average separation between star-forming clouds decreases during the merger, which in turn decreases the cloud–cluster spatial decorrelation from \(\gtrsim1\) kpc to \(\sim0.1\) kpc depicted by tuning fork diagrams, a testable prediction for future observations of interacting low-mass galaxies.

1 Introduction↩︎

Star formation in the observable Universe is notably an inefficient process by any measure. On a galactic scale, only about \(2\%\) of the potentially star-forming gas is actually converted into stars during each dynamical time [1], [2]. This low galaxy-wide star formation efficiency (SFE) is generally linked to feedback mechanisms from massive stars as well as active galactic nuclei (AGNs). Specifically, stellar feedback shapes the evolution and life cycles of star-forming regions at the scale of molecular clouds (MCs) and results in a very short period of star formation for each MC of several megayears [3]–[6]. Consequently, this regulates the cloud-scale SFE to a low level in regular star-forming environments [7]–[11].

However, systems with significantly enhanced galaxy-wide SFE are ubiquitous throughout the Universe and are typically known as “starbursts” [12]. Major mergers of galaxies are commonly recognized as a triggering mechanism of starburst in observations [13]–[19] and numerical simulations ([20], [21]; [22], hereafter ; [23], [24]). The merger-driven starburst is an important stage for the galaxy growth across cosmic time and the formation and evolution of globular cluster populations [25]–[30].

Observationally, merger-driven bursts are evident by the decreased depletion time \(\tau_{\rm dep}\) (i.e. the inverse of the galaxy-wide SFE, see Sec. 3.2) on the Kennicutt–Schmidt plane and an enhanced amount of molecular gas traced by CO at the scale of the entire merger systems to kiloparsec [31]–[37]. Recently, the Atacama Large Millimeter/submillimeter Array (ALMA) observations have been able to resolve the molecular gas in nearby merger and starburst systems at a spatial resolution of \(<100\) pc [38]–[40]. The MCs in such galaxies like NGC 3256 and NGC 4038/9 are found to be large, dense, and highly turbulent. However, at the scale of individual molecular clouds and star-forming regions, it is still not clear whether the galaxy-wide starburst in highly turbulent, out-of-equilibrium mergers is primarily driven by an enhanced cloud-scale SFE, or just an increased gas supply, or a combination of both factors.

Previous numerical studies have demonstrated that the intense tidal forces during galaxy mergers can enhance the formation of massive star clusters [41], [42] while also quickly leading to their disruption [43], [44]. The enhancement of massive cluster formation is attributed to the tidal compression which assembles cold and dense gas, while the tidal disruption is a natural result of the stretching effect of the strong and variable external field. These dual effects are anticipated to influence the life cycle and ultimately the integrated SFE of the MCs within merger systems. While it is well-established that merger-driven enhancement of dense gas fuels star formation, grasping how these positive and negative effects of tides influence individual clouds is vital for physically comprehending merger-induced starburst events at the smallest scale. However, the complex hydrodynamical evolution of clouds and their fragility to stellar feedback render analytical modeling impractical.

In light of recent progress in observing and simulating the detailed multiphase interstellar medium (ISM) and star formation within various galaxies, we present an idealized simulation of a merger between dwarf galaxies. This simulation used the cutting-edge galaxy formation model, Realistic ISM modeling in Galaxy Evolution and Lifecycles (RIGEL, [45], henceforth referred to as ), to present a numerical view of how the dense clouds evolve and how the stars form in a merger-induced starburst system with both high mass and temporal resolutions of \(2\,\text{M}_{\astrosun}\) and \(0.2\,\)Myr, respectively. The aim of this paper is to quantify the driving mechanisms of merger-induced starburst, study the properties of dense cloud population and the relationship between cloud-scale SFE in different galactic environments.

The remainder of the paper is organized as follows. In Section 2, we briefly summarize the physical processes in the RIGEL model, describe the initial conditions (IC) of the dwarf–dwarf merger simulation, and illustrate the workflow to identify dense clouds and star clusters from the simulation snapshots. In Section 3, we give an overview of the global properties of the simulated galaxies. In Section 4, we zoom in to the cloud scale to study the properties of star-forming clouds and their integrated SFE. In Section 5, we use the tuning fork diagrams to provide sub-kpc scale observational predictions of the ISM structures in merger systems. Finally, we discuss the implications and caveats of our model and summarize the key results of the paper in Section 6.

2 Method↩︎

The simulations in this work are performed with arepo[46]–[48], a moving-mesh, finite-volume hydrodynamic code with a second-order Godunov scheme. Following , the idealized galaxy merger simulation used the RIGEL framework, which include self-gravity, hydrodynamics, radiative transfer, heating, cooling, star formation, and stellar feedback through radiation, stellar winds, and SNe. The initial condition (IC) for the merger is set the same as the one in for a direct comparison. We describe the RIGEL model and the ICs in detail in the following subsections.

2.1 The RIGEL model↩︎

RIGEL tracks the formation and evolution of individual massive stars drawn from the [49] initial mass function (IMF) in the resolved multiphase ISM. Stars are formed in cold (\(T < T_\text{th}\)), dense (\(n_\text{H} > n_\text{th}\)), contracting (\(\nabla\cdot {\boldsymbol{v}} < 0\)), self-gravitating, and marginally Jeans-resolved gas where the thermal Jeans length is smaller than 10 times the cell size (\(L_\text{J}=\sqrt{\pi c_s^2 / G \rho}<4\Delta x\)). Compared to , we used a new criterion to find the locally self-gravitating gas cells by \([||\nabla {\boldsymbol{v}}_i||^2+(c_{s,i}/\Delta x)^2)]/8\pi G\rho_i<f_\alpha\) [50], where \(||\nabla {\boldsymbol{v}}_i||\) is the Frobenius norm of velocity and \(f_\alpha\) is a free parameter. We adopted \(f_\alpha=0.5\) through numerical experiments and found it properly selected the self-gravitating cells in gas above the threshold density. In this work, given the mass resolution of \(2\,\text{M}_{\astrosun}\), the density and temperature thresholds for star formation were set as \(T_\text{th}=100\) K and \(n_\text{th}=3000\) cm\(^{-3}\), respectively.

Lifetimes, photon production rates, mass-loss rates, and wind velocities of massive stars (stellar mass \(M_\star > 8\,\text{M}_{\astrosun}\)) are determined by their initial masses and metallicities based on a library that incorporates a variety of stellar models. This stellar feedback model combines the stellar models covering the metallicity from \(10^{-8}\,\text{Z}_{\astrosun}\) to solar from [51], [52], [53], [54], and [55], and the details was described in Section 2.4 of . In our \(2\,\text{M}_{\astrosun}\) resolution simulation, the massive star mass drawn from IMF usually exceeds the mass of its host star particle. To ensure mass conservation at the star cluster scale, star particles continuously lose mass based on an IMF-averaged rate while still being able to release mass in discrete events like supernovae, with any inconsistency from large ejecta masses being counterbalanced by the collective mass loss from all other star particles.

The detailed cooling in star-forming ISM [56], ionization feedback from regions [57], and the Sedov–Taylor phase of SNe [58] are properly accounted for and resolved at a solar-mass-level resolution. The radiation field of seven spectral bands, including infrared (IR, \(0.1-1\) eV), optical (Opt., \(1-5.8\) eV), far-ultraviolet (FUV, \(5.8-11.2\) eV), Lyman–Warner (LW, \(11.2-13.6\) eV), hydrogen ionizing (EUV1, \(13.6-24.6\) eV), ionizing (EUV2, \(24.6-54.4\) eV), and ionizing (EUV3, \(54.4-\infty\) eV) bands [59], was explicitly modeled and coupled to gas hydrodynamics via the moment-based RHD solver arepo-rt [60], [61].

2.2 Initial condition↩︎

We created a merger IC of two identical disks with the parameters identical to the IC used in . Each dwarf galaxy consists of \(2\times10^{10}\,\text{M}_{\astrosun}\) dark matter (DM) particles, a \(4\times10^7\,\text{M}_{\astrosun}\) gas disk with a scale length of \(1.46\) kpc, and a \(2\times10^7\,\text{M}_{\astrosun}\) background stellar disk with a scale length of \(0.73\) kpc and a scale height of \(0.35\) kpc. The DM halo of each galaxy follows a Hernquist profile with an NFW-equivalent concentration parameter \(c=10\) and a spin parameter \(\lambda = 0.03\). The mass resolution was defined to \(2500\,\text{M}_{\astrosun}\) for DM and \(2\,\text{M}_{\astrosun}\) for baryonic components. The gravitational softening lengths were set to \(\epsilon_\text{DM}=39\) pc for DM and \(\epsilon_\star=0.05\) pc for star particles, while the softening of gas particles was determined adaptively with a minimum length of 0.01 pc. The two dwarf galaxies were set on parabolic orbits with a pericentric distance of \(d_\text{peri}=1.46\) kpc, an initial separation of \(d_\text{init}=5\) kpc, and the disc spin parameters of \((\theta_1,\phi_1,\theta_2,\phi_2)= (60^{\circ},30^{\circ},60^{\circ},240^{\circ})\). The initial metallicity of these two galaxies is \(0.1\,\text{Z}_{\astrosun}\).

We ran the simulation for \(200\) Myr of the simulation time. To accurately link the clouds between consecutive snapshots and build the cloud evolution graph (see Section 2.4), we output snapshots every 0.2 Myr.

2.3 Cloud and cluster identification↩︎

We used the CloudPhinder¹ package to identify the gas clouds in our simulations. CloudPhinder detects the local density peaks in galaxies and then identifies the largest self-gravitating gaseous structures by searching for neighboring cells.

To identify the potential and active star-forming clouds, we only considered dense gas, i.e. the gas cells having a density of \(n_\text{H}>100\) cm\(^{-3}\), similar to the density threshold in previous work [62]–[64]. Furthermore, the virial parameter \(\alpha\) for these detected gas structures must be less than 10. This set of selection criteria was relatively lenient, as our aim was to capture the evolution of clouds throughout their entire lifecycle. Similar to [64], our metal-poor dwarf galaxies are deficient in gas and the large scale reservoir for star formation is cold atomic gas. Thus, we did not use fraction to select star-forming clouds.

To identify young star clusters, we used the built-in subfind algorithm [65] of arepo to identify bound stellar systems. Bound stellar systems with at least 35 member stars were identified as star clusters. Clusters with a median age of member stars \(t_\text{50}<5\) Myr were defined as young star clusters.

2.4 Cloud and cluster evolution graph (merger tree) construction↩︎

After identifying the clouds and clusters in each simulation snapshot, we tracked their evolution as a function of time by linking their progenitors and descendants across successive snapshots. This procedure resembles constructing merger trees for halos in cosmological simulations; however, unlike halos, clouds and clusters frequently split and dissolve due to feedback and dynamical processes. In mathematics, such graphs are called directed acyclic graphs (DAGs). Here we refer to these graphs as an “evolution graph”.

Each cloud can spawn multiple children (cloud splits) or have multiple parents (cloud mergers). Following the method outlined in [63], we connected the clouds that have more than one particle ID in common and stored the parents and children of every cloud using the python package networkx [66]. To prevent the marginal connections, the parent and child cloud must have more than 1% of the mass in common.

To generate the evolution paths of the clouds within the graph, we started to walk through the graph from the clouds without any predecessors, referred to as “root nodes”. These root nodes could have one or more child clouds. Thus, each root node generated several possible evolution paths based on the options chosen from the multichild nodes. In general, we used the mass fraction of inherited gas in each child cloud to describe the probability that the child cloud \(k\) would become the next node in the path was \(P_k = m^\text{inh}_k/\sum_i m^\text{inh}_i\), where \(m^\text{inh}_k\) was the mass cloud \(k\) inherited and \(\sum_i m^\text{inh}_i\) is the total mass inherited by all children. We employed various techniques to examine the evolution graphs for particular queries. The specific methods will be presented when addressing relevant questions accordingly. The details of how to link clouds were described in Section 2.4 of [63].

3 Global properties of the dwarf galaxy merger↩︎

3.1 Global star formation rate↩︎

In this section, we present an overview of the star formation activities in our simulation. Fig. 1 provides an overview of the merger process through the gas surface density distribution. Fig. 2 shows the time evolution of the distance between the centers of the two galaxies (top panel) and the global star formation rate (bottom panel). The two galaxies have two pericenter passages during the 200 Myr simulation time, the first at \(\sim45\) Myr, and the second at \(\sim150\) Myr. During the two close encounters, intense galactic torques and tidal forces induce a considerable inflow of gas toward the central 2 kpc region of the merging system, resulting in significant star formation in these gas-rich regions, as demonstrated by the emergence of young stars in the central interacting regions. As expected, the close encounters of galaxies significantly boost the SFR of the galaxies and induce two starbursts onset at \(\sim60\) Myr and \(\sim160\) Myr, respectively. The SFR reaches \(\gtrsim0.02\,\text{M}_{\astrosun}\) yr\(^{-1 }\) during the first starburst and surges to a few \(0.1\,\text{M}_{\astrosun}\) yr\(^{-1}\) at the peak of the second starburst. The peak SFR during the starburst is 130 times the median SFR of two isolated dwarf galaxies. In our simulation, as the total gas mass decreases monotonically and is not significantly consumed by star formation, a 130 times increase in the SFR directly correlates to a 130 times enhancement in the galaxy-wide SFE.

As the SFR varies across two orders of magnitudes, we divide the 200 Myr of simulation time into four stages based on the SFR and merger progress to analyze the results. These phases are labeled as “Approach” for \([0,45)\) Myr, “\(1^\text{st}\) peak” for \([45,90)\) Myr, “Plateau” for \([90,135)\) Myr, and “\(2^\text{nd}\) peak” for \([135,200)\) Myr.

3.2 Galaxy-wide depletion time and star formation efficiency↩︎

To study the triggering mechanism of the starburst in our dwarf galaxy merger, here we check the dense gas fraction in the galaxy and the galaxy-wide SFE.

In the upper panel of Fig. 3, we plot the fraction of \(n_\text{H}>100\) cm\(^{-3}\) dense gas and gas in identified clouds. Similar to what [42] found in the merger of Milky Way (MW)-like galaxies, the merger significantly increases the amount of dense gas and enhances the formation of clouds. During the second SF peak, the total mass of dense gas is 56 times larger than that during the approach stage. The evolution of both gas fractions highly resembles the evolution of SFR, suggesting the SFE of dense gas is roughly constant.

Here, we use the galaxy-wide depletion time as an indicator of the galaxy-wide SFE. The depletion time \(\tau_\text{dep} = M_\text{gas}/\dot{M}_\star\) is the inverse of SFE per dynamical time. In Fig. 3, we present the depletion time for the total gas mass, dense gas mass, and total cloud mass, respectively. As expected, the total gas depletion time varies across two orders of magnitude from \(\gtrsim100\) Gyr to \(\sim 1\) Gyr. However, the depletion times of dense gas and, especially, clouds do not exhibit as significant variation as that of total gas. The depletion time for dense gas varies by up to a factor of 10. The depletion time of gas in clouds fluctuates around its median value of 7.4 Myr with a 16–84 percentile range of 4.7–12 Myr. At the moment when the total gas depletion time significantly decreases (i.e. the \(1^\text{st}\) and \(2^\text{nd}\) SF peak), the cloud depletion time even increases and decays rapidly to the median value.

This short but universal dense gas/cloud depletion time indicates that on the galaxy scale, the higher SFE in the sense of total gas is solely because of the increased amount of dense gas, while the efficiency of dense gas converting to stars remains unchanged. However, it is expected that the strong tides and compression during the merger should impact the properties of dense gas and clouds. To understand this phenomenon, we will zoom in to cloud scales to study how the merger affects the lifetimes and separations of clouds, and finally the cloud-scale integrated SFE.

4 Zoom in to cloud scale: life cycle of clouds and integrated SFE↩︎

We use the evolution graph of the identified clouds to study the star formation activities on the cloud scale. As an illustration, Fig. 4 shows the temporal evolution of the clouds in an evolution graph spanning 14 Myr from 150 Myr to 163 Myr. This period corresponds to the sharp surge of SFR before it reaches the summit of the second peak. As can be seen in 150–153 Myr (the \(3^{\text{rd}}\) to \(6^{\text{th}}\) panels), small, isolated clouds condense from the highly turbulent gas structure. These clouds keep accreting gas and move close to each other due to gravity. At 155 Myr (the \(8^{\text{th}}\) panel), a large cloud is assembled at the center and begin to spawn young stars. Star cluster appears at 156 Myr (the \(9^{\text{th}}\) panel) but the cloud is still accumulating mass from the upper-right, and forms a large single cloud at 159 Myr (the \(12^{\text{th}}\) panel). Another round of star formation happens in this large cloud. As star formation and stellar feedback are ongoing, this large cloud splits to several smaller clouds. Finally, stellar feedback then rapidly disperses the clouds and creates a bubble within 4 Myr.

This example demonstrates that cloud evolution dynamics is highly complex, potentially encompassing numerous merger and split events, as well as nurturing multiple cycles of star formation. These dynamic changes can only be captured owing to the high spatial and temporal resolutions provided by our simulation.

4.1 The cloud population in the merger system↩︎

We first give an overview of the physical properties of the clouds in our simulated merger system. In Fig. 5, we present the relation between the effective radius and velocity dispersion, a.k.a the Larson’s “law” [68]. Following [63], we define the effective radius \(R_\text{eff}\) and one-dimensional velocity dispersion \(\sigma_v\) of a cloud as \[R_\text{eff}=\sqrt{\frac{5}{3}\frac{\sum (m_i |\vec{x_i}-\vec{x_c}|^2)}{\sum m_i}},\] and \[\sigma_v=\sqrt{ \frac{\sum [m_i (|\vec{v_i}-\vec{v_c}|^2+c^2_{s,i})]}{3\sum m_i} },\] where \(\vec{x_c}\) and \(\vec{v_c}\) are the center of mass (COM) position and velocity vector of the cloud, \(m_i\), \(c_{s,i}\), \(\vec{v_i}\), and \(\vec{x_i}\) are the mass, sound speed, velocity, and position vector of the \(i\)-th member gas cell of each cloud, respectively.

The cloud sample for Fig. 5 is selected from simulation snapshots at different stages with a time interval of 10 Myr, as it is long enough compared to the lifetime of most clouds in our simulations to avoid duplicated counting. Each point is color-coded by the virial parameter \(\alpha = 2E_K/E_U\) (top panel) or the environmental tidal strength \(\lambda_\text{tid}\) (bottom panel) estimated at a scale of 50 pc following [69] and [42]. The tidal strength \(\lambda_\text{tid}\) is defined as the maximum of the absolute value of the eigenvalues of the tidal tensor \(\boldsymbol{T}_{ij}\) [70]. The tidal tensor of a cloud is defined at its COM \({\boldsymbol{x}}_\text{cl.}\) as \[\boldsymbol{T}_{ij} = \left.-\frac{\partial^2\Phi_\text{env.}(\boldsymbol{x})}{\partial x_i \partial x_j}\right|_{{\boldsymbol{x}} = {\boldsymbol{x}}_\text{cl.}}\,,\] where \(\Phi_\text{env.}(\boldsymbol{x})=\Phi(\boldsymbol{x}) - \Phi_\text{cl.}(\boldsymbol{x})\) is the environmental gravitational potential, \(\Phi(\boldsymbol{x})\) and \(\Phi_\text{cl.}(\boldsymbol{x})\) are the total potential and the potential generated by the cloud itself, respectively. We use the python gravity solver pytreegrav [71] to evaluate the potentials.

The overall trend of the clouds and the best-fitting result of our sample closely follow the MW observations by [67]. However, the clouds exhibit a slightly larger scatter of 0.31 dex on the Larson’s relation compared with the observed scatter of 0.22 dex. This is because we have a rather tolerant selection criterion so that we included dynamically very different clouds into our sample. The virial parameter is strongly correlated to the velocity dispersion as expected. The clouds with similar virial parameters roughly follow a similar relation with a similar slope.

Moreover, we find that clouds with large velocity dispersion usually present in strong tidal fields, while those with small velocity dispersion are in weak tidal fields. For a given size, clouds in stronger tidal environments have larger velocity dispersion and are highly unbounded, suggesting that the clouds are more turbulent. Though the clouds are unbound as entities, the overdense substructures can collapse and finally form stars. This indicates a highly turbulent star-forming environment during the merger, akin to the observed cloud properties in merger systems [40]. Large clouds usually appear in stronger tidal environments, suggesting that the formation of larger clouds is facilitated by the strong tidal fields due to the merger.

As the tidal force significantly affects the dynamics and boundness of the clouds [72], one may expect the lifetime of clouds in different merger stages to be different. To examine this, in Fig. 6, we show the lifetime PDFs of clouds in different merger stages. To statistically account for all identified clouds and their potential evolution paths, we use 300 Monte-Carlo walkers per root node to sample all possible evolution paths, resulting in 300 probability-weighted Monte-Carlo paths for each root node [73]. The lifetimes yielded by all the Monte-Carlo paths are considered as the sample for Fig. 6. Surprisingly, the normalized PDFs of different stages follow the same exponentially-decaying distribution. We fit the PDFs of all stages together and obtained a characteristic lifetime \(\tau_\text{life} = 1.48\) Myr.

To understand this, we first explain why the lifetime distribution is expected to be an exponentially-decaying distribution. The exponential distributions arise in systems where the probability of a disruptive event occurring is constant over time, regardless of how long the system has existed. For the clouds, its dispersal can be attributed to various physical mechanisms including dynamical shearing, tidal disruption, star formation, and stellar feedback (in the form of radiation, stellar winds, and supernovae). These mechanisms can be regarded as no memory and communication of each other in this context. The unified distribution of lifetimes across different stages suggests that the life cycles of our clouds barely depend on the galactic dynamics. On the contrary, the characteristic lifetime is dominated by the timescale of star formation and feedback. Given the free-fall time scale at our cloud identification criterion \(n_\text{H}=100\) cm\(^{-3}\) is 4 Myr, it is unlikely to be determined by gas depletion due to star formation. Therefore, the early stellar feedback through radiation and stellar winds is the key factor to determine the cloud lifetimes. This is also consistent with the observed feedback timescale of \(1.5\) Myr in NGC 300 [4].

4.2 Effects of compressive tidal fields on star-forming clouds↩︎

As the tidal field is unlikely to dominate the lifetime of individual clouds, now we examine its effect on other cloud properties like their masses and average separations.

Figure 7 shows the distribution of the median tidal strength \(\lambda_\text{tid}^\text{mid}\) and the maximum mass that clouds achieve throughout their lifetime, color-coded by the cloud lifetimes. We found the maximum cloud mass is positively correlated to the tidal strength, similar to the correlation for young cluster mass reported by [42]. Massive clouds preferentially reside in strong tidal fields, while the low-mass clouds have a spread distribution. This is because only the overdense regions in compressive tidal fields can keep accreting and assemble enough dense gas to collapse and form such massive bound structures.

For the low-mass clouds, those that reside in strong tidal fields exhibit shorter lifetimes. However, for the \(M_\text{cloud,max}>1000\,\text{M}_{\astrosun}\) clouds, this correlation is much less significant. Massive clouds located in strong tides can still be as long-lived as those in weak tides. For the most massive clouds, their lifetimes are also the longest.

This result is plausible because the clouds we selected are (quasi-)self-gravitationally bound systems, so that the more massive ones are less susceptible to the galactic tidal fields. Therefore, stellar feedback, which is much more rapid and powerful, still plays the key role in dispersing these clouds. Consequently, the cloud lifetimes exhibit a unified characteristic time scale of \(1.48\) Myr across all the merger stages.

To quantify the separation among individual clouds, we define the average separation \(l\) of a cloud as the mean distance to its 16 nearest neighbors. In the top panel of Fig. 8, we show the relation between the tidal strength \(\lambda_\text{tid}\) and average separation \(l\). The data points are color-coded by the nature of the tidal field (extensive or compressive) experienced by the clouds. The nature of the tidal field is characterized by \(\Lambda =(\lambda_2+\lambda_3 )/2\lambda_1\), where \(\lambda_{1,2,3}\) are the largest, middle, and smallest eigenvalues of the tidal tensor \({\boldsymbol{T}}_{ij}\). Basically, \(\Lambda>1\) indicates fully compressive tides, while \(\Lambda<0\) indicates extensive tides [26].

Apparently, the majority of clouds exist within completely compressive tidal fields. This is expected, as dense clouds predominantly form in central areas of gravitational potential where the gradient is relatively flat [74]. We observed a significant inverse relationship between \(\lambda_\text{tid}\) and \(l\), which can be described by the relationship \(\lambda_\text{tid}\propto l^{-0.43}\). This negative correlation indicates that intense tidal fields compress the gas accumulated and confine the effect of stellar feedback in the central regions, thereby shortening the distance between star-forming clouds.

The bottom panel of Fig. 8 shows the PDFs of average separation \(l\) in different merger stages. During the approach stage, the average separation \(l\) mainly distributes in the \(0.1\)–\(1\) kpc with a median of 288 pc. With the merger in progress and SFR increasing, \(l\) extends and shifts to smaller scales. The average separations in the first and second peaks are 46 pc and 32 pc, respectively. These results indicate the increase of volume density of individual clouds at the kiloparsec scale during the merger, which contributes to the reduced gas-star decorrelation found in Section 5.

From the results presented in Section 4.1 and 4.2 , we learned that the compressive tidal fields facilitate the formation of dense gas clouds by compressing and fueling gas in the central regions, so that we observed the reduced separation among clouds. However, once the self-gravitating clouds have formed, their lifetimes are still dominated by the feedback from the stars they nurture.

4.3 Star formation efficiency↩︎

We now study the cloud-scale SFE for our cloud sample. We use the integrated SFE, \(\epsilon_\text{int} = M_\star/M_0\), to characterize the efficiency of star formation for each cloud, where \(M_\star\) is the final mass of stars formed and \(M_0\) is the mass of the initial gas cloud.

Idealized simulations of isolated GMCs have revealed that the cloud-scale integrated SFE is positively related to the initial mass \(M_0\) or surface density \(\Sigma_0\) of the GMCs [75]–[77]. Analytically, this can be understood through the force balance between gravitational collapse and the feedback-driven momentum flux. Assuming a spherically symmetric geometry where \(M_\star\) stars at the center push a gas shell with a mass of \(M_\text{sh}\) through their feedback, the balance equation goes [77]: \[\frac{M_\star\dot{p}_\star}{4\pi R^2}= \frac{GM_\star\Sigma_\text{sh}}{R^2}+\frac{\beta GM_\text{sh}\Sigma_\text{sh}}{R^2}\,,\] where \(\Sigma_\text{sh}\) is the surface density of the gas shell, \(\dot{p}_\star\) is the specific feedback momentum flux and \(\beta\) is a geometric factor. By defining \(\Gamma = \pi G\Sigma_0/(4f_\text{boost}\dot{p}_\star)\), which describes the relative strength between gravity and feedback, we can solve \(\epsilon_\text{int} = M_\star/M_0\) as \[\label{equ:eps} \epsilon_\text{int} = \frac{\sqrt{\Gamma^2+(4\beta-2)\Gamma+1}-(2\beta-1)\Gamma-1}{2(1-\beta)\Gamma}\,,\tag{1}\] where \(f_\text{boost}\) is a free parameter reflecting the strength of feedback compared to the fiducial value. [77] provided the best-fit values for the free paremters of \(\beta = 1.83\pm0.89\) and \(\dot{p}_\star = (3.32\pm0.64)\times10^9\) cm s\(^{-2}\). This equation goes as \(\epsilon_\text{int}\rightarrow1\) when \(\Sigma_0\rightarrow\infty\), suggesting that stellar feedback will fail at very high gas surface density [76], [78]. However, these type of theoretical models have rarely been tested in a more realistic galactic environment like our merger system.

Here we use our high-resolution simulation to study the cloud-scale SFE. To do this, we must define the related quantities properly. As presented in Fig. 4, the dynamics of cloud evolution are highly complex. Clouds linked in a single graph can undergo multiple times of merging and splitting, and even multiple cycles of star formation. To obtain the integrated SFE, we first define the main evolution path of a cloud following [63]. Given a root node, its main evolution path is defined as the evolution path with the maximum probability of \[P(X_0 \rightarrow X_1 \rightarrow X_2 \rightarrow ... \rightarrow X_n)= \prod_{i=0}^{n-1} P(X_{i} \rightarrow X_{i+1})\] where \(P(X_{i} \rightarrow X_{i+1})\) is the probability of linking cloud \(X_{i}\) and \(X_{i+1}\).

The main evolution path will prefer the shorter paths, so that it avoids including multiple rounds of star formations. However, the main paths starting from different root nodes can merge together (i.e. small clouds assemble a large one), leading to duplicated counting.

To address this problem, we identify the node with the highest total baryon mass \(M_\text{tot} = M_\text{gas}+M_\star\) for every main path, and we group the main paths that achieve this maximum baryon mass at the same node. For example, the polygons with different colors in Fig. 4 belong to different main paths. These main paths finally merge together to one large cloud at \(158\) Myr, and this cloud is thus the shared node of these main paths with the maximum baryon mass. The total baryon mass \(M_\text{tot}\) and surface density for a cloud \(\Sigma_\text{tot}=M_\text{tot}/\pi R_\text{eff}^2\), are derived from the values of this shared node. Finally, the total new star mass \(M_{\star,\text{final}}\) is calculated by summing the masses of new stars formed from nodes belonging to this group of main paths. The integrated SFE is \(\epsilon_\text{int}=M_{\star,\text{final}}/M_\text{tot}\) as defined.

In Fig. 9, we present the distribution of the integrated SFE \(\epsilon_\text{int}\) as a function of the total baryon surface density \(\Sigma_\text{tot}\). The distributions of all cloud samples with star formation are shown in the top panels. We noticed that the theoretical models like equation (1 ) require at least one massive star to provide feedback. Thus, we selected a subset of clouds which form at least \(100\,\text{M}_{\astrosun}\) new stars, and the results are presented in the bottom panels of Fig. 9. We compare our results with the theoretical prediction given by equation (1 ) with the best-fitting parameters from [77] and \(f_\text{boost} = 1\), \(3\), and \(10\).

The highest integrated SFE within our simulated clouds is 84%, while this cloud only forms \(67\,\text{M}_{\astrosun}\) stars. Thus, its high SFE is a natural result as it is free from stellar feedback. The lowest integrated SFE is 0.1%, which also occurs in a low-mass cloud as a result of intense dynamical dispersal. Generally, the SFE of small clouds susceptible to dynamical effects strongly depends on the environments they reside in.

For all the clouds, the scatter of the simulated results is large, the median curve roughly follows the theoretical curve with \(f_\text{boost} = 3\), and the 16–84 confidence level lies between the \(f_\text{boost} = 1\) and \(10\) curves in all the stages. This large scatter is understandable as the SFE of numerous low-mass clouds is determined by environmental effects rather than stellar feedback.

On the contrary, once we confined our sample to the clouds with \(M_{\star,\text{final}}>100\,\text{M}_{\astrosun}\), the scatter becomes much smaller. By mass, 90% of the stars are formed from such clouds. For these clouds, the highest integrated SFE is 60% and the lowest is 1.7%. The median curves and 16–84 confidence levels now lie between the \(f_\text{boost} = 1\) and \(3\) curves for all four merger stages. Quantitatively, within the range \(30<\Sigma_\text{tot}/\text{M}_{\astrosun}\,\text{pc}^{2}<300\), the median SFE during the second peak is 0.17–0.33 dex (0.47 to 0.67 times) lower than that of the approach phase. This suggests that the strong tidal effects in merger only slightly influence the cloud-scale SFE, consistent with what we have found in the cloud lifetimes. As a result, the \(\epsilon_\text{int}\)–\(\Sigma_0\) relation found in idealized cases can still be valid in complicated galactic environments. Since 90% of the stars are formed in the clouds with \(M_{\star,\text{final}}>100\,\text{M}_{\astrosun}\), the majority of the star formation activity is regulated by stellar feedback. This feedback-determined SFE reflects on the galaxy scale and results in the constant galaxy-wide cloud depletion time we found in Fig. 3.

4.4 Tidal disruption of small clouds↩︎

We have seen in Fig. 9 that the stellar feedback determines the integrated SFE of large clouds with \(M_{\star,\text{final}}>100\,\text{M}_{\astrosun}\). Although the small or low SFE clouds with \(M_{\star,\text{final}}<100\,\text{M}_{\astrosun}\) only form 10% of the stars, it is interesting to examine if the galactic dynamics regulates their SFE.

In Fig. 10, we present the joint probability distribution of \(\epsilon_\text{int}\) and \(\Sigma_\text{tot}\) color-coded by the tidal strength for the \(M_{\star,\text{final}}<100\,\text{M}_{\astrosun}\) small clouds. As these clouds are free from stellar feedback, their \(\epsilon_\text{int}\) present a large scatter. Nonetheless, we can evidently see that for the \(\Sigma_\text{tot}<300\,\text{M}_{\astrosun}\) pc\(^{-2}\) clouds, the ones located in weak tidal fields present high integrated SFE reaching \(70\)%, while the ones located in strong fields present low integrated SFE down to \(0.1\)%. For these clouds, the strength of tidal disruption plays an essential role in determining their SFEs. The extensive distribution of \(\epsilon_\text{int}\) reflects the drastic variation of tidal strength in the merger system.

Interestingly, several outliers with high surface density of \(\Sigma_\text{tot}>300\,\text{M}_{\astrosun}\) pc\(^{-2}\), high SFE of \(\epsilon_\text{int}>10\%\) and high tidal strength of \(\lambda_\text{tid}>3.5\times10^4\) Gyr\(^{-2}\) can be found in Fig. 10. This may be due to the intense compressive tides causing these clouds to be sufficiently compact, allowing them to continue collapsing and forming stars without disturbance. We will present an order-of-magnitude analysis about this in Section 6.1.

4.5 Formation of massive star clusters from massive clouds↩︎

Lastly, we examine the mass function of clouds and clusters as results of the varying galactic tidal fields. In Fig. 11, we present the cloud mass functions (top panel) and young cluster (\(<5\) Myr) initial mass functions (CIMFs, bottom panel) of the merger system during different stages. To avoid duplicated counting, we select clouds in snapshots with a time interval of 10 Myr, which is much longer than the characteristic lifetime of our clouds (1.48 Myr, see Section 4.1). On the other hand, the young clusters here are defined as the clusters whose median stellar age is younger than 5 Myr, so the cadence to obtain CIMFs is 5 Myr.

We describe the cloud mass functions and CIMFs as a power-law distribution. Evidently, the cloud mass functions during the merger-induced starburst exhibit a flatter slope and extend to a mass larger than \(10^5\,\text{M}_{\astrosun}\). Quantitatively, we fit the mass functions with a power-law form \({\rm d}N\,/\,{\rm d}M\propto M^{\alpha}\). For this fitting, we consider only clouds with masses exceeding \(100\,\text{M}_{\astrosun}\) to ensure a complete mass sampling. The power-law indices \(\alpha\) and maximum cloud masses (captured by the selected snapshots) are \((-2.36,\,-1.84,\,-2.08,\,-1.79)\) and \((1.8,\,23,\,11,\,90)\times10^3\,\text{M}_{\astrosun}\) for the Approach, \(1^\text{st}\) peak, Plateau, and \(2^\text{nd}\) peak stages, respectively. As the merger proceeds, dense gas accumulated in a stronger tidal field tends to form larger clouds (see Fig. 7). Thus, the slope becomes shallower and the maximum cloud mass becomes larger. As a result, during the merger, there are more massive clouds that exist and provide the fuel for the intense star formation. Clouds dozens of times more massive than those in isolated galaxies appear during the merger and become the place to nurture massive star clusters.

As a result, the CIMF changes across stages with a similar trend, with the slopes of \((-2.23,\,-2.08,\,-2.32,\,-1.98)\) for each stage. This slope does not change as significantly as that of cloud mass functions, suggesting some self-regulation of the cluster formation. However, the maximum cluster mass does change significantly from \(926\,\text{M}_{\astrosun}\) within the approach stage to \(2.6\times10^4\,\text{M}_{\astrosun}\) within the \(2^\text{nd}\) peak. We notice that the massive star clusters are still deficient compared with the increase of massive clouds, especially compared with the results of . We will discuss this in Section 6.2.

An observation constraint of the cluster formation in different environments is the relation between the galaxy-wide SFR and maximum embedded cluster mass (\(M_\text{cluster,max}\)) in the galaxy [80], [81]. This relation can be interpreted as a consequence of how the shape of the CIMF is influenced by SFR [25], [82]. We examine this relation in our simulation in Fig. 12. Our simulated results roughly follow the theoretical relation derived by [79], who assumed a single slope power-law embedded cluster mass function with an SFR-dependent power-law index [82]. Although our simulation only covered the low SFR range with a handful of data, it matches both the slope and the magnitude of the observations. Moreover, the most intense SFR during the \(2^\text{nd}\) peak reaches the region with abundant observations, and the most massive young clusters formed in this period also present comparable masses with the observations. This result exhibits the ability and fidelity of the RIGEL model to capture and predict the self-regulated star formation in drastically variable environments. We noticed that the observations are mostly more massive galaxies and not specifically selected to be mergers. The future application of RIGEL across a broader range of systems will allow for a more equitable comparison with observations.

5 Spatial decorrelation between clouds and young clusters↩︎

In the previous section, we found that strong tides in merger decrease the average separation between dense clouds but do not significantly change the efficiency of stellar feedback. This finding can provide observational predictions for future observations.

In this section, we use the “tuning fork” diagram [83] to quantify the ISM structure and star-forming activities on scales of 0.1–1 kpc, where the internal structures of clouds are not yet resolved. This diagram illustrates the relative bias of molecular gas depletion time measured in the apertures focused either on dense gas or young stars, making it an observational tool to measure the spatial decorrelation between newly formed stellar regions and dense gas clouds.

We noticed that the observational lower limit for molecular gas surface density of \(\Sigma_\text{\ce{H2}}>13\,\text{M}_{\astrosun}\) pc\(^{-2}\) [4] is much higher than the typical \(\Sigma_\text{\ce{H2}}\) in our simulation. Thus, we do not apply realistic observational sensitivity cuts to our analysis. The results obtained here can be regarded as a theoretical exploration and more practical observational predictions require future simulations of larger systems.

5.1 Two-dimensional Tuning fork diagram↩︎

To generate the molecular gas map, we project the volume density of within an (8 kpc)\(^3\) box along the \(z\)-axis. In our \(0.1\,\text{Z}_{\astrosun}\) metal-poor galaxies, a significant fraction of the mass can be found in diffuse, -dominated gas where star formation does not happen [45], [84]. We apply a cut on the fraction of gas cells of \(2x_\text{\ce{H2}}>0.1\) to exclude such diffuse .

We follow the method outlined by [85] to generate the SFR map. We select young star particles with ages 2–5 Myr. We project the young stars to the same mesh as the molecular gas map and compute \(\dot{\Sigma}_\star\) in each pixel. The lower age cut mimics the typical timescale of the embedded phase of young clusters [86], [87] while the upper limit mimics the typical stellar age traced by H\(\alpha\) [88]–[90].

With the \(\Sigma_\text{\ce{H2}}\) and \(\dot{\Sigma}_\star\) maps at a native resolution of 5 pc, we smooth the maps using a 2D Gaussian filter with a width of 20 pc. We use the scikit-image package [91] to identify the local extrema of \(\Sigma_\text{\ce{H2}}\) and \(\dot{\Sigma}_\star\) maps. Following the method outlined by [83], we then smooth the maps with a top hat filter with the window size \(L\) and compute the average \(\Sigma_\text{\ce{H2}}\) and \(\dot{\Sigma}_\star\) at the locations of gas and SFR peaks identified above. By doing this, we can get a tuning fork diagram for a single snapshot.

We present the 2D tuning fork diagram from our simulation in the top panel of Fig. 13. In all four stages, our simulation results present the tuning-fork-like shape, while the scale at which the SFR branch and the gas branch converge differs with the stages. During the approach stage, the branches converge at the largest scale of \(L\gtrsim1\) kpc, which aligns with the observations of nearby disk galaxies [4], [5].

As the two galaxies merge, the convergence point shifts toward smaller scales. At the first SF peak, convergence scale decreases to \(\sim 300\) pc, and it increases again as the SFR decreases during the Plateau stage. Finally, the convergence point shifts to \(\sim100\) pc in the most intense starburst during the second SF peak. This corresponds to the average separation between the peaks of molecular gas and star formation shrinking as the SFR increases. The opening of the tuning fork at the smallest scale also decreases during the merger, reflecting that the contrast of spatial decorrelation between recent SFR and dense gas decreases.

When two galaxies collide, the projection effect could be substantial: several star-forming regions may appear to overlap along the line-of-sight so that it reduces both the convergence scale and opening. To justify our findings in the two-dimensional tuning fork, we extend the traditional 2D tuning fork analysis to three-dimensional space in the next subsection.

5.2 Three-dimensional tuning fork diagram↩︎

To generate the three-dimensional tuning fork diagram, we first use a cloud-in-cell (CIC) algorithm to assign the particle data to a Cartesian grid. Consistent with the 2D version, we only select the gas particles with \(2x_\text{\ce{H2}}>0.1\) for the molecular gas cube and 2–5 Myr age young stars for the SFR cube. We smooth the cubes using a 3D Gaussian fliter with a width of 20 pc and identify the local extrema of \(\rho_\text{\ce{H2}}\) and \(\dot{\rho}_\star\) in the smoothed cubes. We then smooth the cubes with a top hat filter and compute the average \(\rho_\text{\ce{H2}}\) and \(\dot{\rho}_\star\) at the locations of gas and SFR peaks. The resultant 3D tuning fork diagram is presented in the bottom panel of Fig. 13.

Similar to the 2D tuning fork diagram, the convergence point moves towards a smaller scale during the ongoing merger, reaching \(\sim100\) pc during the most intense starburst at the second SF peak. However, the opening of the tuning fork does not decrease as evident as in the 2D case. Therefore, this 3D tuning fork diagram confirms that the merger does, in fact, decrease the average separation between regions of star formation and dense clouds, whereas the reduction of the opening of the two branches is possibly a projection effect.

The reduced spatial decorrelation between young stars and dense clouds observed in the simulated merger can be explained by various physical causes. This might indicate that the continued accumulation of dense gas and high environmental pressure during a merger reduced the ability of stellar feedback to disperse clouds. Alternatively, stellar feedback may remain largely unchanged at the cloud scale, but the volume density of individual clouds at kiloparsec scale increases. Our findings in Section 4 favor the latter scenario, highlighting the importance of early stellar feedback. Future simulations of larger systems can further examine this finding and synergy with the ALMA observations of nearby mergers.

6 Discussions and conclusions↩︎

6.1 Energy of compressive tides↩︎

In the Section 4, we found that the strong tidal fields present during a merger boost the accumulation of cold clouds; however, they do not alter the SFE of these clouds that are forming massive clusters. Here, we present some order of magnitude estimates to illustrate these phenomena from an energy-based perspective.

Assuming an isotropic tidal field, the tidal energy of a spherical cloud with mass \(M\) and radius \(R\) can be estimated as \[\begin{align} E_\text{tid}&=\frac{1}{2}\lambda_\text{tid}MR^2\notag\\ &=10^{46}\,\text{erg}\left(\frac{\lambda_\text{tid}}{10^4\,\text{Gyr}^{-2}}\right)\left(\frac{M}{10^3\,\text{M}_{\astrosun}}\right)\left(\frac{R}{10\,\text{pc}}\right)^2\,. \end{align}\]

This is much smaller than the energy released by a single SN explosion (\(10^{51}\) erg), and even the early stellar wind and radiation feedback [92], [93]. Thus, even a single massive star can dominate the evolution of the clouds.

However, the binding energy of such a cloud is only \[\begin{align} U_\text{bind}&=-\frac{3}{5}\frac{GM^2}{R}\notag\\ &=-5\times10^{45}\,\text{erg}\left(\frac{M}{10^3\,\text{M}_{\astrosun}}\right)^2\left(\frac{R}{10\,\text{pc}}\right)^{-1}\,. \end{align}\] Therefore, the tidal energy is actually comparable to the binding energy of a cloud. If there is no massive star present in a cloud, the SFE will strongly depend on the tidal strength, as we have seen in Section 4.4.

In Section 4.4, we found several clouds in strong tidal fields present high integrated SFE, and we suggested that they are compact enough to be free from tidal disruption. To quantify this, we compare the tidal energy and binding energy by \[\left|\frac{U_\text{bind}}{E_\text{tid}} \right|\sim \frac{GM}{\lambda_\text{tid}R^3}\sim \frac{\beta G\bar{\rho}}{\lambda_\text{tid}}\,,\] where \(\bar{\rho}\) is the average density of the cloud and \(\beta\) is a geometric factor.

Given the maximum tidal strength in our simulation is \(\lambda_\text{tid}\sim10^5\) Gyr\(^{-2}\), we have the critical density of \(\bar{\rho}_\text{crit} = 22\beta\,\text{M}_{\astrosun}\,\text{pc}^{-3}\) for \({u_\text{bind}}/{u_\text{tid}}=1\). Clouds with lower average density should be susceptible to the tidal effects while much denser clouds are tightly bound. In our simulation, we found that the actual volume density of the outlier clouds in Fig.10 (i.e., clouds that have SFE within strong tidal fields) exceeds \(200\,\text{M}_{\astrosun}\) pc\(^{-3}\), consistent with the above estimate.

On the other hand, at a large scale, the tidal energy can far exceed the self-gravitational energy. In the warm neutral medium (WIM) which dominates the ISM in dwarf galaxies, the \(p{\rm d}V\) work done by tidal force can be expressed as \[{\rm d}U_\text{int}=4\pi nkTR^2{\rm d}R\,,\notag\] and the change of tidal energy is \[{\rm d} E_\text{tid}=\lambda_\text{tid}MR{\rm d}R\,.\notag\] Combining these two equations, we can get a characteristic scale analog to the Jeans length \[\begin{align} R_\text{tid}=\sqrt{\frac{3\pi kT}{\lambda_\text{tid}\mu m_\text{H}}}=280\,\text{pc}\left(\frac{\lambda_\text{tid}}{10^4\,\text{Gyr}^{-2}}\right)^{-0.5}\left(\frac{T}{10^4\,\text{K}}\right)^{0.5}\,. \end{align}\] As a comparison, the Jeans length in \(T=10^{4}\) K, \(n_\text{H}=1\) cm\(^{-3}\) WIM is \(1.6\) kpc, much larger that \(R_\text{tid}\). In the isolated dwarf galaxies or solar neighborhood tidal field with \(\lambda_\text{tid}=1600\) Gyr\(^{-2}\), \(R_\text{tid}\) is comparable to the Jeans length, while in the strong merger tidal field with \(\lambda_\text{tid}>10^5\) Gyr\(^{-2}\), \(R_\text{tid}\) decreases to \(<100\) pc. Therefore, the compressive tides can lead to the gravitational instability at a much smaller scale compared with the self-gravity and facilitate the formation of dense star-forming clouds. In the case of our simulation, we found an enhanced cloud formation during the merger with a mean separation of tens to hundreds of parsecs.

6.2 Comparison with other work↩︎

As the IC of our simulation is the same as , we can make a direct comparison with the series of GRIFFIN papers. The morphology and SFR of our merger system show good agreement with . The slope of their CIMFs varies between \(-1.7\) and \(-2\) during the merger, comparable but slighly flatter than our result of \(-1.98\). More importantly, the most massive cluster in reaches a mass of \(7\times10^5\,\text{M}_{\astrosun}\), while we only obtain \(\sim3\times10^4\,\text{M}_{\astrosun}\) clusters. Thus, our model deficits massive clusters compared to . As and we use the same algorithm to find young clusters (though they define young as \(<10\) Myr while our definition is \(<5\) Myr), the difference in CIMF and maximum cluster mass should be attributed to the star formation and feedback model. We notice that the present observation cannot constrain the models with such details [the only hint is that the maximum cluster mass of \(7\times10^5\,\text{M}_{\astrosun}\) in seems too large for their peak SFR of \(\sim0.2\,\text{M}_{\astrosun}\) yr\(^{-1}\) (see Fig.12)], yet it is valuable for a theoretical exploration.

The star formation model determines the site, timing, and dynamics of newly formed stars, thus determining the properties of clusters in terms of the mass, size, boundness of clusters and the CIMF [94]. As both and our model do not resolve the formation dynamics of stars, the size and boundness of clusters can be strongly affected by the underlying star formation density criterion and the gravitational softening. Since the gas mass resolution of our simulation is higher, we adopted a relatively high density criterion of \(3000\) cm\(^{-3}\) and a small softening for stars of \(0.05\) pc. As a comparison, the Jeans criterion adopted by typically allows stars to be formed at \(500\) cm\(^{-3}\) and their softening is \(0.1\) pc. Therefore, the clusters in our simulation are more compact, especially for the small cluster formed through a monolithic cloud collapse. Such compact clusters are thus hard to be accreted and assembled by large clusters. More sophisticated formation and dynamics models [28], [95], [96] are necessary to generate more realistic cluster populations.

Another major difference between our model and is that we use radiative transfer in arepo-rt to explicitly model the radiative feedback from massive stars, while adopted a Strömgren-type approximation method to model the effects of radiation. [97] noticed that their approximation method preferentially ionized mass concentration regions with high density. This can potentially reduce the momentum injection rate by region expansion [57], [98]. Thus, it can be interesting to test the dependence of CIMF on the intensity of radiative feedback. We run a low-resolution simulation using \(10\,\text{M}_{\astrosun}\) resolution and stellar radiation luminosities reduced by a factor of 0.1 across all bands, labeled as “low/wRT” (weak radiative feedback). As a controlled group, we also run a \(10\,\text{M}_{\astrosun}\) resolution simulation with normal luminosities. We show the CIMFs of these simulations in Fig. 14. The CIMF is not sensitive to the numerical resolution, but is very susceptible to the intensity of radiative feedback. The weaker radiative feedback leads to a flatter slope and the most massive cluster reaches a mass of \(>10^5\,\text{M}_{\astrosun}\). This result suggests that the process of cluster formation is influenced not merely by the presence of early feedback, but is also sensitive to the strength of this feedback, requiring an accurate modeling of all stellar feedback channels.

[64] studied the cold clouds in a similar dwarf–dwarf merger simulation. Their cloud selection criteria are slightly different from ours: they did not use the sub-virial criterion, and they adopted the same density criterion as ours but an additional temperature cut of \(T<100\) K. Opposite to [64], we found that the slope of the cloud mass function varies significantly across the merger stages. Although it can be related to the cloud selection criterion, it is also important to study whether the cloud mass function is sensitive to the star formation and feedback model in the future. Additionally, they found that the cloud lifetimes do not seem to be affected by the environment. Our results confirm this finding and explain it by the feedback-dominated cloud life cycle. They found that the clouds have an e-folding lifetime distribution with a characteristic timescale exceeding \(3\) Myr, potentially due to not restricting the virial status of the clouds, which allows them to encompass a longer period of accretion and collapse.

[42] studied the major merger of two MW-like galaxies with similar orbital configuration using the SMUGGLE model [50]. Similar to our results, they also found that both the mass function of GMCs and young clusters varies significantly across the merger stages. The slope of the GMC mass function varies from 2.47 in isolated galaxies to 1.99 in mergers, and the CIMF slope varies from 1.93 to 1.67. Since the MW-like galaxies are much more massive, the tidal strength in their simulation extends to \(\lambda_\text{tid}\gtrsim10^7\) Gyr\(^{-2}\), and they found a clear positive relation between tidal field strength and young cluster mass in this large \(\lambda_\text{tid}\) range. As we found the cloud-scale SFE is insensitive to the tidal environment, this cluster-tide correlation can be a direct result of the positive correlation between \(\lambda_\text{tid}\) and \(M_\text{cloud}\) we presented in Fig. 7. The SFR in the MW merger is enhanced by a factor of 10. This enhancement is an order of magnitude weaker compared to that observed in our dwarf-dwarf merger, whereas their galaxy-wide SFE, conversely, is comparable during the peak star formation phases and even presents a roughly factor of 2 higher than the maximum value. This suggests that major mergers can exert a more significant influence on star formation in low-mass galaxies. We can understand this by considering that star formation is less efficient in isolated dwarf galaxies compared to the MW disk. This inefficiency arises because the dwarf disk is predominantly composed of warm, diffuse HI gas. In contrast, during merger processes, this diffuse HI gas is accumulated into the central region, leading to the formation of star-forming dense clouds.

[63] studied the evolution of giant molecular clouds (GMCs) in a suite of Milky Way-mass galaxies. They found that the integrated SFE of the GMCs can also be described by a \(\epsilon_\text{int}\)–\(\Sigma_0\) relation. They also noticed that the galactic dynamics can introduce scatter in this relation. Our finding for the feedback-free clouds in Section 4.4 provides an extreme example for this: the SFE of clouds without massive stars is strongly correlated with the intensity of tidal fields.

Recently, [99] and [100] studied the cloud-scale SFE in the thesan-zoom cosmological simulation of high-redshift galaxies [101]. Although they did not track the evolution of individual GMCs, they found that GMCs exhibit universal statistics in terms of mass functions, size, surface density, and especially cloud-scale SFE, regardless of the environments. The feedback-regulated cloud-scale SFE identified in this work offers a direct explanation for their results, suggesting that our findings might be applicable to describe the star formation cycles throughout cosmic time.

Previous numerical experiments have demonstrated that the opening and convergence points of tuning fork diagrams are notably influenced by the feedback model, particularly the pre-SN feedback [85], [102]. [103] further proposed that they can be affected by galactic dynamics in MW-like galaxies. Our finding in Section 5 underscores the importance of galactic dynamics and predicts a significantly reduced spatial decorrelation in merger systems.

6.3 Limitations of this work↩︎

We note that there are a few limitations of this work. First, we use a highly idealized prograde major merger of two identical dwarf galaxies, which provides the maximum effects of tidal interaction and highest enhancement of gas accumulation. A comprehensive study requires fully exploring the properties of host galaxies and orbit parameters. Moreover, due to the low mass of the dwarf galaxies, the boost of the tidal field during the merger is within two orders of magnitudes, and the strongest tides only reach \(\lambda_\text{tid}\approx10^5\) Gyr\(^{-2}\). The major merger of MW-mass galaxies can further boost the tidal strength to \(\lambda_\text{tid}>10^7\) Gyr\(^{-2}\) [42], which potentially results in more drastic effects on the cloud evolution. The gas surface density and the star formation rate surface density of our system are also too low to directly compare with spatially-resolved observations of nearby massive galaxies. As discussed in Section 6.2, our model is not yet able to accurately capture the internal structure of star clusters and, therefore, cannot investigate their long-term dynamical evolution. The time and site of massive star formation do not have any preference in our model, which can be oversimplified. For example, if massive stars preferentially form in denser cores and at later times compared with low-mass stars [104], [105], the cloud-scale SFE can be higher with respect to our prediction. The compact cluster formation could also be attributed to this over-simplified star formation model. Additionally, the resolved collisional dynamics might considerably expand the cluster sizes, driven by stellar dynamical processes [96]. We are working on a new attempt to reproduce a more realistic star cluster population in the RIGEL model (Deng et al. in prep.).

6.4 Summary and conclusions↩︎

In this paper, we attempted to find the mechanisms that drive the starburst in the dwarf–dwarf major merger using a high-resolution numerical simulation.

Firstly, we found that the high-pressure, dense gas significantly increased due to the angular momentum cancellation and tidal compression. The total mass of \(n_\text{H}>100\) cm\(^{-3}\) dense gas increased by a factor of 56 during the merger.

However, the depletion time of dense star-forming gas, especially the gas in clouds, remains roughly constant across all merger stages. The depletion time of gas in clouds fluctuates around its median value of 7.4 Myr, with the 16–84 percentile range of 4.7–12 Myr. This suggests that the merger does not significantly change the cloud-scale SFE.

Indeed, by tracking the evolution of individual clouds we found that the lifetimes of clouds across all stages follow the same exponentially-decaying distribution with a characteristic timescale of \(\tau_\text{life}=1.48\) Myr. Moreover, the integrated SFE of clouds formed at least one massive star (\(M_{\star,\text{final}}>100\,\text{M}_{\astrosun}\)) can always be described by a \(\epsilon_\text{int}\)–\(\Sigma_0\) relation assuming the feedback-gravity balance. The median integrated SFE during the most intense starburst only changed by 0.17–0.33 dex lower compared to that in isolated galaxies. This suggests that the star formation activities in the merger system are regulated by stellar feedback rather than galactic dynamics. The integrated SFE in clouds without massive star present a strong correlation with the tidal strength, while these clouds only contribute 10% of the total stellar mass.

The compressive tidal fields in the merger system play an essential role in compressing and accumulating gas to form cold dense clouds. The cloud mass is positively correlated to the tidal strength, and massive clouds preferentially reside in strong tidal fields. As a result, the slope of the cloud mass function becomes flatter from \(-2.36\) in approach to \(-1.79\) in the \(2^\text{nd}\) peak. The presence of larger clouds serves as sites for the formation of massive clusters. Indeed, we have also noticed an alteration in the slope of CIMF from \(-2.23\) to \(-2.03\), although this variation is less pronounced compared to the change in the cloud mass function.

The average separation of clouds shows a strong negative correlation with the tidal strength: \(\lambda_\text{tid}\propto l^{-0.43}\). It is compressed from hundreds pc in isolated galaxies to \(<50\) pc in the merger. As a result, the spatial decorrelation between clouds and clusters observed on the tuning fork diagram also decreases from \(\gtrsim1\) kpc to \(\sim0.1\) kpc during the merger.

Based on these findings, we can conclude that the merger stimulates the formation of dense clouds fueling the star formation due to the strong compressive tidal fields. However, once massive stars form in the clouds, they determine the further evolution of their natal clouds and eventually the integrated SFE. Therefore, the cloud-scale SFE or depletion time remains unchanged during the merger.

References↩︎