© The Authors 2025

1. Introduction

Star formation (SF) in the observable universe is notably an inefficient process by any measure. On a galactic scale, only about 2% of potential star-forming gas is actually converted into stars during each dynamical time (Kennicutt 1998; Kennicutt & Evans 2012). 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 SF for each MC of several megayears (Hollyhead et al. 2015; Kruijssen et al. 2019; Chevance et al. 2022, 2023). Consequently, this regulates the cloud-scale SFE at a low level in regular star-forming environments (Evans 2009; Murray 2011; Lee et al. 2016; Mattern et al. 2024; Leroy et al. 2025).

However, systems with significantly enhanced galaxy-wide SFE are ubiquitous throughout the universe and are typically known as “starbursts” (Weedman et al. 1981). Major mergers of galaxies are commonly recognized as a triggering mechanism for starbursts based on observations (e.g., Kennicutt et al. 1987; Sanders & Mirabel 1996; Barton et al. 2000; Ellison et al. 2008; Horstman et al. 2021; Gao et al. 2023; Kaviraj et al. 2025) and numerical simulations (Barnes & Hernquist 1991; Di Matteo et al. 2007; Lahén et al. 2020, hereafter L20; Renaud et al. 2019, 2022). The merger-driven starburst is an important stage for galaxy growth across cosmic time as well as the formation and evolution of globular cluster populations (Li et al. 2017; Renaud et al. 2017; Lahén et al. 2019, 2024; Keller et al. 2020; Pascale et al. 2025).

Observationally, merger-driven bursts are evident from the decreased depletion time τ_dep (i.e., the inverse of the galaxy-wide SFE; see Sect. 3.2) on the Kennicutt–Schmidt plane and an enhanced amount of molecular gas traced by CO on scales ranging from the entire merger systems down to kiloparsec (e.g., Daddi et al. 2010; Klaas et al. 2010; Martinez-Badenes et al. 2012; Pan et al. 2018; Bemis & Wilson 2019; Wilson et al. 2019; Kennicutt & De Los Reyes 2021). Recently, observations performed with the Atacama Large Millimeter/submillimeter Array (ALMA) have been able to resolve the molecular gas in nearby merger and starburst systems at a spatial resolution of < 100 pc (e.g., Brunetti et al. 2021; Brunetti & Wilson 2022; Brunetti et al. 2024). The MCs in such galaxies as NGC 3256 and NGC 4038/9 are found to be large, dense, and highly turbulent. However, at the scale of individual MCs and star-forming regions, it remains unclear whether the galaxy-wide starburst in highly turbulent, out-of-equilibrium mergers is primarily driven by an enhanced cloud-scale SFE, an increased gas supply, or both.

Previous numerical studies have demonstrated that intense tidal forces during galaxy mergers enhance the formation of massive star clusters (e.g., Renaud et al. 2008; Li et al. 2022) while also quickly leading to their disruption (e.g., Baumgardt & Makino 2003; Kruijssen et al. 2011). The enhancement of massive cluster formation is attributed to tidal compression, which assembles cold and dense gas, whereas 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 the merger-driven enhancement of dense gas fuels SF, 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 SF 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, Deng et al. 2024b, henceforth referred to as D24), to present a numerical view of how dense clouds evolve and how stars form in a merger-induced starburst system with both high mass and temporal resolutions of 2 M_⊙ and 0.2 Myr, respectively. The aim of this paper is to quantify the driving mechanisms of merger-induced starbursts and study both the properties of dense cloud populations 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 (ICs) 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 on 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-kiloparsec 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 (Springel 2010; Pakmor et al. 2016; Weinberger et al. 2020), a moving-mesh, finite-volume hydrodynamic code with a second-order Godunov scheme. Following D24, the idealized galaxy merger simulation used the RIGEL framework, which includes self-gravity, hydrodynamics, radiative transfer, heating, cooling, SF, and stellar feedback through radiation, stellar winds, and supernovae (SNe) explosion. The initial condition for the merger is set the same as the one in L20 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 Chabrier (2003) initial mass function (IMF) in the resolved multiphase ISM. Stars are formed in cold (T < T_th), dense (n_H > n_th), contracting (∇ ⋅ v < 0), self-gravitating, and marginally Jeans-resolved gas where the thermal Jeans length is smaller than four times the cell size ($L_{\mathrm{J}}=\sqrt{\pi c_s^2/G\rho }<4\mathrm{\Delta }x$). Compared to D24, we used a new criterion to find the locally self-gravitating gas cells by [||∇v_i||² + (c_s, i/Δx)²)]/8πGρ_i < f_α (Marinacci et al. 2019), where ||∇v_i|| is the Frobenius norm of velocity and f_α is a free parameter. We adopted f_α = 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 M_⊙, the density and temperature thresholds for SF were set as T_th = 100 K and n_th = 3000 cm⁻³, respectively.

Lifetimes, photon production rates, mass-loss rates, and wind velocities of massive stars (stellar mass M_★ > 8 M_⊙) 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⁻⁸ Z_⊙ to solar from Schaerer (2002), Lanz & Hubeny (2003), Emerick et al. (2019), Tanikawa et al. (2020), and Gessey-Jones et al. (2022), and the details are described in Section 2.4 of D24. In our 2 M_⊙ resolution simulation, the massive star mass drawn from IMF 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 such as 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 (Kim et al. 2023), ionization feedback from H II regions (Deng et al. 2024a), and the Sedov—Taylor phase of SNe (Kim & Ostriker 2015) were 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), He I ionizing (EUV2, 24.6 − 54.4 eV), and He II ionizing (EUV3, 54.4 − ∞ eV) bands (Kannan et al. 2020), was explicitly modeled and coupled to gas hydrodynamics via the moment-based RHD solver AREPO-RT (Kannan et al. 2019; Zier et al. 2024).

2.2. Initial condition

We created a merger IC of two identical disks with the parameters identical to the IC used in L20. Each dwarf galaxy consists of 2 × 10¹⁰ M_⊙ dark matter (DM) particles, a 4 × 10⁷ M_⊙ gas disk with a scale length of 1.46 kpc, and a 2 × 10⁷ M_⊙ 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 Navarro–Frenk–White (NFW)-equivalent concentration parameter c = 10 and a spin parameter λ = 0.03. The mass resolution was defined to 2500 M_⊙ for DM and 2 M_⊙ for baryonic components. The gravitational softening lengths were set to ϵ_DM = 39 pc for DM and ϵ_★ = 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_peri = 1.46 kpc, an initial separation of d_init = 5 kpc, and the disk spin parameters of (θ₁, ϕ₁, θ₂, ϕ₂) = (60° ,30° ,60° ,240° ). The initial metallicity of these two galaxies is 0.1 Z_⊙.

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, 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., gas cells with a density of n_H > 100 cm⁻³, similar to the density threshold in previous work (e.g., Grisdale et al. 2018; Ni et al. 2025; Fotopoulou et al. 2024). Furthermore, the virial parameter α for these detected gas structures had to 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 life cycle. Similar to Fotopoulou et al. (2024), our metal-poor dwarf galaxies are deficient in H₂ gas, and the large-scale reservoir for SF is cold atomic gas. Thus, we did not use H₂ fraction to select star-forming clouds.

To identify young star clusters, we used the built-in SUBFIND algorithm (Springel et al. 2001) 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₅₀ < 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 Ni et al. (2025), we connected the clouds that had more than one particle ID in common and stored the parents and children of every cloud using the PYTHON package NETWORKX (Hagberg et al. 2008). To prevent marginal connections, the parent and child cloud must have more than 1% of their mass in common.

To generate the evolution paths of the clouds within the graph, we began to walk through the graph from the clouds without any predecessors, referred to as “root nodes”. These root nodes can have one or more child clouds. Thus, each root node generated several possible evolution paths based on the options chosen from the multi-child 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_k^inh/∑_im_i^inh, where m_k^inh is the inherited mass cloud k and ∑_im_i^inh is the total mass inherited by all children. We employed various techniques to examine the evolution graphs for particular queries. The specific methods are presented when addressing the relevant questions raised in this paper. The details of how to link clouds are described in Section 2.4 of Ni et al. (2025).

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 (SFR; bottom panel). The two galaxies have two pericentric passages during the 200 Myr simulation time: the first at ∼45 Myr, and the second at ∼150 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 SF 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 ∼60 Myr and ∼160 Myr, respectively. The SFR reaches ≳0.02 M_⊙ yr⁻¹ during the first starburst and surges to a few 0.1 M_⊙ yr⁻¹ 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 SF, a 130 times increase in the SFR directly correlates to a 130 times enhancement in the galaxy-wide SFE.

Fig. 1. Gas surface density maps in different merger stages. The four panels show snapshots of the approaching stage (20 Myr, upper left), the first passage when the first SF peak appears (70 Myr, upper right), the first apocenter when the SFR exhibits a plateau (110 Myr, bottom left), and the second encounter when the second SF peak occurs (160 Myr, bottom right). The overset points represent young (< 40 Myr) and massive (> 8 M⊙) stars with ages color-coded from violet (youngest) to red (oldest). A movie of the simulation is available at https://www.bilibili.com/video/BV18vKRzzEnj.

Fig. 2. Evolution of the distance between two galaxies (upper panel) and the SF history of the merger system (bottom panel). The distance is measured between the centers of mass of the background stars in the galaxies. Bottom panel: SFR averaged over a timescale of 5 Myr. The blue curve represents the total SFR of the two dwarfs in the merger, while the orange curve represents twice the SFR of an isolated dwarf galaxy. The four stages of the merger–“Approach”, “1st peak”, “Plateau”, and “2nd peak” –are marked with shaded regions in various colors.

As the SFR varies across two orders of magnitudes, we divided 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^st peak” for [45, 90) Myr, “Plateau” for [90, 135) Myr, and “2^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_H > 100 cm⁻³ dense gas and gas in identified clouds. Similar to what Li et al. (2022) 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 the dense gas is 56 times higher than that during the approach stage. The evolution of these two gas fractions highly resembles the evolution of the SFR, suggesting that the SFE of the dense gas is roughly constant.

Fig. 3. Evolution of galaxy-wide ISM properties and depletion time. Upper panel: Evolution of the fraction of nH > 100 cm−3 dense gas (red curve) and gas in identified clouds (dense gas clumps with virial parameter α < 10, blue curve). Bottom panel: Galaxy-wide depletion times calculated by the total gas mass (black curve), dense gas mass (red curve), and total cloud mass (blue curve). The depletion time of the gas in the clouds fluctuates around its median value of 7.4 Myr (dashed blue line), with a 16–84 percentile range of 4.7–12 Myr.

Here, we use the galaxy-wide depletion time as an indicator of the galaxy-wide SFE. The depletion time τ_dep = M_gas/Ṁ_★ is the inverse of the 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 ≳100 Gyr to ∼1 Gyr. However, the depletion times of the dense gas and, especially, clouds do not exhibit as significant variation as that of the total gas. The depletion time for dense gas varies by up to a factor of 10. The depletion time of the 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^st and 2^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 zoomed-in on cloud scales to study how the merger affects the lifetimes and separations of clouds, as well as 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 SF 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 seen in 150–153 Myr (the 3^rd to 6^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^th panel), a large cloud is assembled at the center and begins to spawn young stars. Star clusters appear at 156 Myr (the 9^th panel), while the cloud continues to accumulate mass from the upper-right, forming a large single cloud at 159 Myr (the 12^th panel). Another round of SF occurs in this large cloud. As SF and stellar feedback occur, this large cloud splits into several smaller clouds. Subsequently, stellar feedback rapidly disperses the clouds and creates a bubble within 4 Myr.

Fig. 4. Evolution of clouds in a sequence spanning 14 Myr. The first two panels show the location of the zoom-in region in the merger system. The third to the last panels show the column density distribution inside a (200 pc)3 box. In these panels, the color maps refer to the gas surface density distribution of the zoom-in regions, and the dots represent the stars color-coded according to their age. The polygons in the zoom-in panels depict the 2D convex hulls of the clouds identified from a specific evolution graph. The evolution paths in this graph merge at 158 Myr and share a node (green polygon), where we evaluate the total baryon mass Mtot and surface density Σtot for the SFE calculation (see Section 4.3). The clouds form in the upper-right and lower-left regions at the beginning (see the panel of 153 Myr for an example) and accumulate at the center (155 Myr). Star clusters appear at 156 Myr, but the cloud continues to accumulate mass from the upper-right, forming a large single cloud at 159 Myr. The cloud continues to form stars while collapsing. Finally, the stellar feedback rapidly disperses the clouds, leaving bubbles within 4 Myr. The colorful polygons appearing between 155 Myr to 160 Myr indicate the main paths sharing the same maximum baryon mass node in this evolution graph.

This example demonstrates the high complexity of cloud evolution dynamics: it potentially encompasses numerous merger and split events, and exhibits multiple cycles of SF. 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” (Larson 1981). Following Ni et al. (2025), we define the effective radius R_eff and 1D velocity dispersion σ_v of a cloud as

$$\begin{array}{c}\hfill R_{\mathrm{eff}}=\sqrt{\frac{5}{3}\frac{\sum (m_i|{\mathit{x}}_{\mathit{i}}-{\mathit{x}}_{\mathit{c}}{|}^2)}{\sum m_i}},\end{array}$$(1)

Fig. 5. Cloud velocity dispersion as a function of effective radius. The clouds are selected from the snapshot spaced at 10 Myr intervals to avoid duplicate counts and eliminate bias toward the long-lived clouds. The points in the top panel are color-coded by their virial parameters, while those in the bottom panel are color-coded by tidal strength λtid. Clouds in stronger tidal environments typically have larger velocity dispersions and virial parameters. The red line represents the observations of MW clouds by Solomon et al. (1987), and the red-shaded region represents the observed 0.22 dex dispersion of log10σv. The gray line and shaded region represent our best-fitting relation and the log10σv dispersion of 0.31 dex. Our sample fitting closely follows the observations, although the scatter in the sample is slightly larger than that of the observed MW clouds, due to a rather tolerant selection criterion designed to include dynamically diverse clouds.

and

$$\begin{array}{c}\hfill {\sigma }_v=\sqrt{\frac{\sum [m_i{(|{\mathit{v}}_{\mathit{i}}-{\mathit{v}}_{\mathit{c}}|}^2+c_{s,i}^2)]}{3\sum m_i}},\end{array}$$(2)

where x_c and v_c are the center of mass (COM) position and velocity vector of the cloud, m_i, c_s, i, v_i, and 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 duplicate counts. Each point is color-coded by the virial parameter α = 2E_K/E_U (top panel) or the environmental tidal strength λ_tid (bottom panel) estimated at a scale of 50 pc following Li & Gnedin (2019) and Li et al. (2022). The tidal strength λ_tid is defined as the maximum of the absolute value of the eigenvalues of the tidal tensor T_ij (Renaud et al. 2011). The tidal tensor of a cloud is defined at its COM x_cl. as

$$\begin{array}{c}\hfill {\mathit{T}}_{\mathit{ij}}={-\frac{{\partial }^2{\mathrm{\Phi }}_{\mathrm{env}.}(\mathit{x})}{\partial x_i\partial x_j}|}_{\mathit{x}={\mathit{x}}_{\mathrm{cl}.}},\end{array}$$(3)

where Φ_env.(x) = Φ(x)−Φ_cl.(x) is the environmental gravitational potential, Φ(x) and Φ_cl.(x) are the total potential and the potential generated by the cloud itself, respectively. We used the PYTHON gravity solver PYTREEGRAV (Grudić & Gurvich 2021) to evaluate the potentials.

The overall trend of the clouds and the best-fitting result of our sample closely follow the MW observations by Solomon et al. (1987). However, the clouds exhibit a slightly larger scatter of 0.31 dex on the Larson’s relation compared to the observed scatter of 0.22 dex. This is because we had a rather tolerant selection criterion, which allowed us to include dynamically different clouds into our sample. The virial parameter strongly correlates with 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 higher velocity dispersions 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 (Brunetti et al. 2024). 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 (e.g., Lee et al. 2025), one may expect the lifetime of clouds in different merger stages to be different. To examine this, in Fig. 6, we show the lifetime probability distribution functions (PDFs) of clouds in different merger stages. To statistically account for all identified clouds and their potential evolution paths, we used 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 (see Section 5.1 in Jeffreson et al. 2021a, for details). The lifetimes yielded by all the Monte-Carlo paths are considered as the sample for Fig. 6. Surprisingly, the normalized PDFs of the different stages follow the same exponentially decaying distribution. We fit the PDFs of all stages together and obtained a characteristic lifetime τ_life = 1.48 Myr.

Fig. 6. Cloud lifetime distributions at different merger stages. The dashed gray line denotes the maximum likelihood estimate fitting of the lifetime for all stages. The clouds in all merger stages follow the same exponentially decaying relation with a characteristic time of 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, their dispersal can be attributed to various physical mechanisms including dynamical shearing, tidal disruption, SF, 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 SF and feedback. Given the free-fall timescale at our cloud identification criterion n_H = 100 cm⁻³ is 4 Myr, it is unlikely to be determined by gas depletion due to SF. 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 (see Kruijssen et al. 2019).

4.2. Effects of compressive tidal fields on star-forming clouds

As the tidal field is unlikely to dominate the lifetime of individual clouds, we now examine its effect on other cloud properties, such as their masses and average separations. Figure 7 shows the distribution of the median tidal strength ${\lambda }_{\mathrm{tid}}^{\mathrm{mid}}$ and the maximum mass that clouds achieve throughout their lifetime, color-coded by the cloud lifetimes. We find that the maximum cloud mass correlates positively with the tidal strength, as the young cluster mass reported by Li et al. (2022). Massive clouds preferentially reside in strong tidal fields, while low-mass clouds have a spread distribution. This is because only the overdense regions in compressive tidal fields can keep accreting and assembling enough dense gas to collapse and form such massive bound structures.

Fig. 7. Distribution of clouds in the ${\lambda }_{\mathrm{tid}}^{\mathrm{mid}}$–Mcloud, max plane, color-coded by their average lifetimes. ${\lambda }_{\mathrm{tid}}^{\mathrm{mid}}$ is the median tidal strength experienced by clouds throughout their lifetime, and Mcloud, max is the maximum mass the clouds ever reached during their lifetime. The solid black line and shaded region represent the median relation and the 16–84 percentile range, respectively.

The low-mass clouds that reside in strong tidal fields exhibit shorter lifetimes. However, for the M_cloud, max > 1000 M_⊙ 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 timescale of 1.48 Myr across all the merger stages.

To quantify the separation among individual clouds, we defined 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 λ_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 Λ = (λ₂ + λ₃)/2λ₁, where λ_1, 2, 3 are the largest, middle, and smallest eigenvalues of the tidal tensor T_ij. Basically, Λ > 1 indicates fully compressive tides, while Λ < 0 indicates extensive tides (see Renaud et al. 2017, Section 3.3 for details).

Fig. 8. Relation between average cloud separation l and tidal strength. The cloud sample is the same as that presented in Fig. 5. The top panel is color-coded by the nature of tides Λ. The dashed gray line represents the fitting of λtid ∝ l−0.43. The bottom panel shows the PDFs of the average separation l in different merger stages. Most of the clouds are in fully compressive tidal fields (Λ > 1). The average separation l extends and shifts to smaller scales as the merger progresses and the SFR increases.

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 (Renaud et al. 2009). We observed a significant inverse relationship between λ_tid and l, which can be described by the relationship λ_tid ∝ 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 an increase in volume density of individual clouds at the kiloparsec scale during the merger, which contributes to the reduced gas-star de-correlation found in Section 5.

From the results presented in Section 4.1 and 4.2, we learn that 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 remain 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, ϵ_int = M_★/M₀, to characterize the efficiency of SF for each cloud, where M_★ is the final mass of stars formed and M₀ 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₀ or surface density Σ₀ of the GMCs (e.g., Fall et al. 2010; Grudić et al. 2018; Li et al. 2019). 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_★ stars at the center push a gas shell with a mass of M_sh through their feedback, the balance equation goes (Li et al. 2019):

$$\begin{array}{c}\hfill \frac{M_{★}{\dot{p}}_{★}}{4\pi R^2}=\frac{G{M}_{★}{\mathrm{\Sigma }}_{\mathrm{sh}}}{R^2}+\frac{\beta G{M}_{\mathrm{sh}}{\mathrm{\Sigma }}_{\mathrm{sh}}}{R^2},\end{array}$$(4)

where Σ_sh is the surface density of the gas shell, ṗ_★ is the specific feedback momentum flux and β is a geometric factor. By defining $\mathrm{\Gamma }=\pi G{\mathrm{\Sigma }}_0/(4f_{\mathrm{boost}}{\dot{p}}_{★})$, which describes the relative strength between gravity and feedback, we can solve ϵ_int = M_★/M₀ as

$$\begin{array}{c}\hfill {ϵ}_{\mathrm{int}}=\frac{\sqrt{{\mathrm{\Gamma }}^2+(4\beta -2)\mathrm{\Gamma }+1}-(2\beta -1)\mathrm{\Gamma }-1}{2(1-\beta )\mathrm{\Gamma }},\end{array}$$(5)

where f_boost is a free parameter reflecting the strength of feedback compared to the fiducial value. Li et al. (2019) provided the best-fit values for the free parameters of β = 1.83 ± 0.89 and ṗ_★ = (3.32 ± 0.64) × 10⁹ cm s⁻². This equation goes as ϵ_int → 1 when Σ₀ → ∞, suggesting that stellar feedback will fail at very high gas surface density (Grudić et al. 2018; Dekel et al. 2023). However, these types of theoretical models have rarely been tested in a more realistic galactic environment such as 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 SF. To obtain the integrated SFE, we first define the main evolution path of a cloud following Ni et al. (2025). Given a root node, its main evolution path is defined as the evolution path with the maximum probability of

$$\begin{array}{c}\hfill P(X_0\to X_1\to X_2\to \dots \to X_n)={\displaystyle \prod _{i=0}^{n-1}}P(X_i\to X_{i+1}),\end{array}$$(6)

where P(X_i → X_i + 1) is the probability of linking cloud X_i and X_i + 1.

The main evolution path prefers the shorter paths, so that it avoids including multiple rounds of SF. 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_tot = M_gas + M_★ 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_tot and surface density for a cloud Σ_tot = M_tot/πR_eff², are derived from the values of this shared node. Finally, the total new star mass M_★, final is calculated by summing the masses of new stars formed from nodes belonging to this group of main paths. The integrated SFE is ϵ_int = M_★, final/M_tot as defined.

In Fig. 9, we present the distribution of the integrated SFE ϵ_int as a function of the total baryon surface density Σ_tot. The distributions of all cloud samples with SF are shown in the top panels. We noticed that the theoretical models such as equation (5) require at least one massive star to provide feedback. Thus, we selected a subset of clouds, which form at least 100 M_⊙ 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 (5) with the best-fitting parameters from Li et al. (2019) for f_boost = 1, 3, and 10.

Fig. 9. Distribution of integrated SFE ϵint as a function of the total baryon surface density Σtot. Top panels: Distribution of all cloud samples with SF. Bottom panel: Results for the clouds formed more than 100 M⊙ in newly formed stars. The data points are color-coded by the normalized density in log-space estimated by a Gaussian kernel. The green curves show the median values, and the shaded regions indicate the 16–84 confidence levels. The dashed curves denote Equation (5) with the best-fitting parameters from Li et al. (2019) for fboost = 1, 3, and 10. Although the scatter of the simulated results is large, the median curve roughly follows the theoretical curve with fboost = 3. The 16–84 confidence level lies between the fboost = 1 and 10 curves at all the stages.

The highest integrated SFE within our simulated clouds is 84%, while this cloud only forms 67 M_⊙ 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_boost = 3, and the 16–84 confidence level lies between the f_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_★, final > 100 M_⊙, 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_boost = 1 and 3 curves for all four merger stages. Quantitatively, within the range $30<{\mathrm{\Sigma }}_{\mathrm{tot}}/{\mathrm{M}}_{\odot }{\mathrm{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 mergers only slightly influence the cloud-scale SFE. This is consistent with what we find in cloud lifetimes. As a result, the ϵ_int–Σ₀ relation found in the idealized cases can still be valid in complicated galactic environments. Since 90% of the stars are formed in the clouds with M_★, final > 100 M_⊙, the majority of the SF 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 see in Fig. 9 that the stellar feedback determines the integrated SFE of large clouds with M_★, final > 100 M_⊙. Although the small or low SFE clouds with M_★, final < 100 M_⊙ only form 10% of the stars, it is interesting to examine whether the galactic dynamics regulates their SFE.

In Fig. 10, we present the joint probability distribution of ϵ_int and Σ_tot color-coded by the tidal strength for the M_★, final < 100 M_⊙ small clouds. As these clouds are free from stellar feedback, their ϵ_int present a large scatter. Nonetheless, we can evidently see that for the Σ_tot < 300 M_⊙ pc⁻² 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 ϵ_int reflects the drastic variation of tidal strength in the merger system.

Fig. 10. Joint probability distribution of ϵint and Σtot color-coded by the tidal strength. Here λtid is the value when the cloud reaches its maximum baryon mass.

Interestingly, several outliers with high surface density of Σ_tot > 300 M_⊙ pc⁻², high SFE of ϵ_int > 10%, and high tidal strength of λ_tid > 3.5 × 10⁴ Gyr⁻² 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 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 duplicate counts, we selected 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.

Fig. 11. Cloud mass functions (top panel) and young cluster (< 5 Myr) initial mass functions (CIMFs, bottom panel) of the merger system at different stages. The dashed gray line denotes the power-law relation with a slope of −1.5, −2, and −2.5, respectively. A starburst induced by the merger leads to the formation of massive clusters, especially during the second SFR peak. The maximum cloud masses are (1.8,  23,  11,  90)×103 M⊙ for the Approach, 1st peak, Plateau, and 2nd peak stages, respectively. The maximum cluster masses are (0.93,  3.4,  1.8,  26)×103 M⊙ for these four stages.

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⁵ M_⊙. Quantitatively, we fit the mass functions with a power-law form dN / dM ∝ M^α. For this fitting, we considered only clouds with masses exceeding 100 M_⊙ to ensure a complete mass sampling. The power-law indices α and maximum cloud masses (captured by the selected snapshots) are ( − 2.36,   − 1.84,   − 2.08,   − 1.79) and (1.8,  23,  11,  90)×10³ M_⊙ for the Approach, 1^st peak, Plateau, and 2^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 provide the fuel for the intense SF. 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 M_⊙ within the approach stage to 2.6 × 10⁴ M_⊙ within the 2^nd peak. We notice that the massive star clusters are still deficient compared with the increase in massive clouds, especially compared with the results of L20. We 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_{cluster, max}) in the galaxy (e.g., Weidner & Kroupa 2004; Schulz et al. 2015). This relation can be interpreted as a consequence of how the shape of the CIMF is influenced by SFR (Weidner et al. 2013; Li et al. 2017). We examine this relation in our simulation in Fig. 12. Our simulated results roughly follow the theoretical relation derived by Yan et al. (2017), who assumed a single slope power-law embedded cluster mass function with an SFR-dependent power-law index (e.g., Weidner et al. 2013). 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^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 SF in drastically variable environments. We notice that the observations are mostly from more massive galaxies and not specifically selected to be mergers. A future application of RIGEL across a broader range of systems will allow for a more equitable comparison with observations.

Fig. 12. Most massive young cluster mass vs. galaxy-wide SFR (Mcluster, max–Ṁ★) relation. The simulation data are color-coded by stage of the merger. The dashed red line represents the theoretical relation derived following Yan et al. (2017). The gray crosses represent observational data from Weidner & Kroupa (2004) and Schulz et al. (2015), which have a typical uncertainty of 0.3 dex.

5. Spatial de-correlation between clouds and young clusters

In the previous section, we found that strong tides in mergers 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 (Kruijssen et al. 2018) 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 on either dense gas or young stars, making it an observational tool for measuring the spatial de-correlation between newly formed stellar regions and dense gas clouds.

We noticed that the observational lower limit for molecular gas surface densities of Σ_H2 > 13 M_⊙ pc⁻² (Kruijssen et al. 2019) is much higher than the typical Σ_H2 in our simulation. Thus, we did 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. 2D tuning fork diagram

To generate the molecular gas map, we projected the volume density of H₂ within an (8 kpc)³ box along the z-axis. In our 0.1 Z_⊙ metal-poor galaxies, a significant fraction of the H₂ mass can be found in diffuse, H I-dominated gas where SF does not happen (e.g., Hu et al. 2021; Deng et al. 2024b). We applied a cut on the H₂ fraction of gas cells of 2x_H2 > 0.1 to exclude such diffuse H₂.

We followed the method outlined by Semenov et al. (2021) to generate the SFR map. We selected young star particles with ages 2–5 Myr. We projected the young stars to the same mesh as the molecular gas map and compute ${\dot{\mathrm{\Sigma }}}_{★}$ in each pixel. The lower age cut mimics the typical timescale of the embedded phase of young clusters (e.g., Kim et al. 2021; Deshmukh et al. 2024), while the upper limit mimics the typical stellar age traced by Hα (Flores 2021; Smith et al. 2022; Tacchella et al. 2022).

With the Σ_H2 and ${\dot{\mathrm{\Sigma }}}_{★}$ maps at a native resolution of 5 pc, we smoothed the maps using a 2D Gaussian filter with a width of 20 pc. We use the SCIKIT-IMAGE package (Van der Walt et al. 2014) to identify the local extrema of Σ_H2 and ${\dot{\mathrm{\Sigma }}}_{★}$ maps. Following the method outlined by Kruijssen et al. (2018), we then smoothed the maps with a top hat filter with the window size L and computed the average Σ_H2 and ${\dot{\mathrm{\Sigma }}}_{★}$ at the locations of gas and SFR peaks identified above. By doing this, we obtained 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 ≳ 1 kpc, which aligns with the observations of nearby disk galaxies (Kruijssen et al. 2019; Chevance et al. 2022).

Fig. 13. 2D (top panel) and 3D (bottom panel) “tuning fork” diagrams from our simulation at different phases. These diagrams reflect the deviation of the H2 depletion time centered on H2 (τH2, gas, upper branch) or recent SFR peaks (τH2, SFR, lower branch) from the galaxy-averaged value as a function of the spatial scale L. The curves of the 1st peak, Plateau, and 2nd peak are multiplied by constants of 10, 100, and 1000, respectively, to avoid overlapping. The vertical lines indicate the scale when τH2, gas exceeds τH2, SFR by 50%.

As the two galaxies merge, the convergence point shifts toward smaller scales. At the first SF peak, the convergence scale decreases to ∼300 pc, and it increases again as the SFR decreases during the Plateau stage. Finally, the convergence point shifts to ∼100 pc in the most intense starburst during the second SF peak. This corresponds to the average separation between the peaks of molecular gas and SF 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 de-correlation 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, reducing both the convergence scale and opening. To justify our findings in the 2D tuning fork, we extended the traditional 2D tuning fork analysis to 3D space in the next subsection.

5.2. 3D tuning fork diagram

To generate the 3D tuning fork diagram, we first used a cloud-in-cell (CIC) algorithm to assign the particle data to a Cartesian grid. Consistent with the 2D version, we only selected the gas particles with 2x_H2 > 0.1 for the molecular gas cube and 2–5 Myr age young stars for the SFR cube. We smoothed the cubes using a 3D Gaussian filter with a width of 20 pc and identify the local extrema of ρ_H2 and ${\dot{\rho }}_{★}$ in the smoothed cubes. We then smoothed the cubes with a top hat filter and compute the average ρ_H2 and ${\dot{\rho }}_{★}$ at the locations of gas and SFR peaks. The resulting 3D tuning fork diagram is presented in the bottom panel of Fig. 13.

Similar to the 2D tuning fork diagram, the convergence point moves toward a smaller scale during the ongoing merger, reaching ∼100 pc during the most intense starburst at the second SF peak. However, the opening of the tuning fork does not decrease as evidently 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 SF and dense clouds, whereas the reduction of the opening of the two branches is possibly a projection effect.

The reduced spatial de-correlation 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 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{array}{cc}\hfill E_{\mathrm{tid}}& =\frac{1}{2}{\lambda }_{\mathrm{tid}}M{R}^2\hfill \\ \hfill & ={10}^{46}\mathrm{erg}\left(\frac{{\lambda }_{\mathrm{tid}}}{{10}^4{\mathrm{Gyr}}^{-2}}\right)\left(\frac{M}{{10}^3{\mathrm{M}}_{\odot }}\right){\left(\frac{R}{10\mathrm{pc}}\right)}^2.\hfill \end{array}$$(7)

This is much smaller than the energy released by a single SN explosion (10⁵¹ erg), and even the early stellar wind and radiation feedback (e.g., Agertz et al. 2013; Jeffreson et al. 2021b). Thus, even a single massive star can dominate the evolution of the clouds.

However, the binding energy of such a cloud is only

$$\begin{array}{cc}\hfill U_{\mathrm{bind}}& =-\frac{3}{5}\frac{G{M}^2}{R}\hfill \\ \hfill & =-5\times {10}^{45}\mathrm{erg}{\left(\frac{M}{{10}^3{\mathrm{M}}_{\odot }}\right)}^2{\left(\frac{R}{10\mathrm{pc}}\right)}^{-1}.\hfill \end{array}$$(8)

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 seen in Section 4.4.

In Section 4.4, we found several clouds in strong tidal fields with 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

$$\begin{array}{c}\hfill \left|\frac{U_{\mathrm{bind}}}{E_{\mathrm{tid}}}\right|\sim \frac{\mathit{GM}}{{\lambda }_{\mathrm{tid}}{R}^3}\sim \frac{\beta G\overline{\rho }}{{\lambda }_{\mathrm{tid}}},\end{array}$$(9)

where $\overline{\rho }$ is the average density of the cloud and β is a geometric factor.

Given that the maximum tidal strength in our simulation is λ_tid ∼ 10⁵ Gyr⁻², we obtain the critical density of ${\overline{\rho }}_{\mathrm{crit}}=22\beta {\mathrm{M}}_{\odot }{\mathrm{pc}}^{-3}$ for u_bind/u_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 M_⊙ pc⁻³, 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 pdV work done by tidal force can be expressed as

$$\begin{array}{c}\hfill \mathrm{d}{U}_{\mathrm{int}}=4\pi nkT{R}^2\mathrm{d}R,\end{array}$$

and the change of tidal energy is

$$\begin{array}{c}\hfill \mathrm{d}{E}_{\mathrm{tid}}={\lambda }_{\mathrm{tid}}MR\mathrm{d}R.\end{array}$$

Combining these two equations, we can obtain a characteristic scale analogous to the Jeans length:

$$\begin{array}{c}\hfill R_{\mathrm{tid}}=\sqrt{\frac{3\pi kT}{{\lambda }_{\mathrm{tid}}\mu m_{\mathrm{H}}}}=280\mathrm{pc}{\left(\frac{{\lambda }_{\mathrm{tid}}}{{10}^4{\mathrm{Gyr}}^{-2}}\right)}^{-0.5}{\left(\frac{T}{{10}^4\mathrm{K}}\right)}^{0.5}.\end{array}$$(10)

As a comparison, the Jeans length in T = 10⁴ K, n_H = 1 cm⁻³ WIM is 1.6 kpc, much larger than R_tid. In the isolated dwarf galaxies or solar neighborhood tidal field with λ_tid = 1600 Gyr⁻², R_tid is comparable to the Jeans length, while in the strong merger tidal field with λ_tid > 10⁵ Gyr⁻², R_tid decreases to < 100 pc. Therefore, the compressive tides can lead to gravitational instability at a much smaller scale compared with 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 L20, we can make a direct comparison with the series of GRIFFIN papers. The morphology and SFR of our merger system show good agreement with L20. The slope of their CIMFs varies between −1.7 and −2 during the merger, comparable but slightly flatter than our result of −1.98. More importantly, the most massive cluster in L20 reaches a mass of 7 × 10⁵ M_⊙, while we only obtain ∼3 × 10⁴ M_⊙ clusters. Thus, our model lacks massive clusters compared to L20. As L20 we used 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 SF 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 × 10⁵ M_⊙ in L20 seems too large for their peak SFR of ∼0.2 M_⊙ yr⁻¹ (see Fig. 12)], yet it is valuable for a theoretical exploration.

The SF 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 (Hislop et al. 2022). As both L20 and our model do not resolve the formation dynamics of stars, the size and boundness of clusters can be strongly affected by the underlying SF 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⁻³ and a small softening for stars of 0.05 pc. As a comparison, the Jeans criterion adopted by L20 typically allows star formation at 500 cm⁻³, with a softening of 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 for large clusters to accrete and assemble. More sophisticated formation and dynamics models (e.g., Lahén et al. 2023, 2024, 2025) are necessary to generate more realistic cluster populations.

Another major difference between our model and L20 is that we used radiative transfer in AREPO-RT to explicitly model the radiative feedback from massive stars, while L20 adopted a Strömgren-type approximation method to model the effects of radiation. Hu et al. (2017) noticed that their approximation method preferentially ionized mass concentration regions with high density. This can potentially reduce the momentum injection rate by H II region expansion (e.g., Matzner 2002; Deng et al. 2024a). Thus, it can be interesting to test the dependence of CIMF on the intensity of radiative feedback. We ran a low-resolution simulation using 10 M_⊙ 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 ran a 10 M_⊙ 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⁵ M_⊙. 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.

Fig. 14. CIMFs for the simulation with weakened radiative feedback. The red curve shows the result of 10 M⊙ resolution stellar radiation luminosities reduced by a factor of 0.1 across all bands. Solid black and dashed gray curves represent the results from the simulation with normal radiative feedback at 2 M⊙ and 10 M⊙, respectively. The CIMF is not sensitive to the numerical resolution but is very susceptible to the intensity of radiative feedback.

Fotopoulou et al. (2024) 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. Contrary to Fotopoulou et al. (2024), we find 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 SF and feedback model in the future. Additionally, they find 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 find 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.

Li et al. (2022) studied the major merger of two MW-like galaxies with similar orbital configuration using the SMUGGLE model (Marinacci et al. 2019). Similar to our results, they also find that both the mass function of GMCs and young clusters vary 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 λ_tid ≳ 10⁷ Gyr⁻², and they find a clear positive relation between tidal field strength and young cluster mass in this large λ_tid range. Since we find the cloud-scale SFE to be insensitive to the tidal environment, this cluster-tide correlation can be a direct result of the positive correlation between λ_tid and M_cloud 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 SF 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 SF in low-mass galaxies. We can understand this by considering that SF 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.

Ni et al. (2025) studied the evolution of giant molecular clouds (GMCs) in a suite of MW-mass galaxies. They find that the integrated SFE of the GMCs can also be described by a ϵ_int–Σ₀ relation. They also notice 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, Shen et al. (2025) and Wang et al. (2025) studied the cloud-scale SFE in the THESAN-ZOOM cosmological simulation of high-redshift galaxies (Kannan et al. 2025). Although they did not track the evolution of individual GMCs, they find 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 SF 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 (Fujimoto et al. 2019; Semenov et al. 2021). Jeffreson et al. (2020) 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 de-correlation in merger systems.

6.3. Limitations of this work

We note that there are a few limitations in this work. First, we used 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 a thorough exploration of 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 λ_tid ≈ 10⁵ Gyr⁻². The major merger of MW-mass galaxies can further boost the tidal strength to λ_tid > 10⁷ Gyr⁻² (e.g., Li et al. 2022), which potentially results in more drastic effects on the cloud evolution. The gas surface density and the SFR 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 location of massive star formation do not have any preference in our model, which may be an oversimplification. For example, if massive stars preferentially form in denser cores and at later times compared with low-mass stars (Grudić et al. 2022; Farias et al. 2024), the cloud-scale SFE can be higher with respect to our prediction. The compact cluster formation could also be attributed to this oversimplified SF model. Additionally, the resolved collisional dynamics might considerably expand the cluster sizes, driven by stellar dynamical processes (Lahén et al. 2025). We are working on reproducing 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 starbursts in dwarf–dwarf major mergers using a high-resolution numerical simulation. Our key findings are summarized as follows:

First, we find that high-pressure, dense gas significantly increases due to angular momentum cancellation and tidal compression. The total mass of n_H > 100 cm⁻³ 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 find that the lifetimes of clouds across all stages follow the same exponentially decaying distribution with a characteristic timescale of τ_life = 1.48 Myr. Moreover, the integrated SFE of clouds that formed at least one massive star (M_★, final > 100 M_⊙) can always be described by a ϵ_int–Σ₀ relation assuming the feedback-gravity balance. The median integrated SFE during the most intense starburst was lower by only 0.17–0.33 dex compared to that in isolated galaxies. This suggests that the SF activities in the merger system are regulated by stellar feedback rather than galactic dynamics. The integrated SFE in clouds without massive stars presents a strong correlation with tidal strength, while these clouds only contribute 10% of their 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 upon reaching −1.79 in the 2^nd peak. Larger clouds serve as locations for the formation of massive clusters. Indeed, we also notice an alteration in the slope of the CIMF from −2.23 to −2.03, even though 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 tidal strength: λ_tid ∝ l^−0.43. It is compressed from hundreds of parsecs in isolated galaxies to < 50 pc during the merger. As a result, the spatial de-correlation between clouds and clusters observed in the tuning fork diagram also decreases from ≳1 kpc to ∼0.1 kpc during the merger.

Based on these findings, we can conclude that the merger stimulates the formation of dense clouds fueling the SF 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.

Data availability

Movie associated with Fig. 1 is available at https://www.aanda.org

Acknowledgments

We thank Volker Springel for giving us access to AREPO, and thank Natalia Lahén for helping us generate the initial condition. YD is grateful to Zhiyuan Li, Zhiqiang Yan, and Yulong Gao for useful discussions. HL is supported by the National Key R&D Program of China No. 2023YFB3002502, the National Natural Science Foundation of China under No. 12373006, and the China Manned Space Program with grant No. CMS-CSST-2025-A10. RK acknowledges support of the Natural Sciences and Engineering Research Council of Canada (NSERC) through a Discovery Grant and a Discovery Launch Supplement (funding reference numbers RGPIN-2024-06222 and DGECR-2024-00144) and York University’s Global Research Excellence Initiative. BL acknowledges the funding of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy EXC 2181/1 – 390900948 (the Heidelberg STRUCTURES Excellence Cluster). We use PYTHON packages NUMPY (Harris et al. 2020), SCIPY (Virtanen et al. 2020), ASTROPY (Astropy Collaboration 2013, 2018), MATPLOTLIB (Hunter 2007), and PAICOS (Berlok et al. 2024) to analyze and visualize the simulation data.

References

All Figures

Fig. 1. Gas surface density maps in different merger stages. The four panels show snapshots of the approaching stage (20 Myr, upper left), the first passage when the first SF peak appears (70 Myr, upper right), the first apocenter when the SFR exhibits a plateau (110 Myr, bottom left), and the second encounter when the second SF peak occurs (160 Myr, bottom right). The overset points represent young (< 40 Myr) and massive (> 8 M⊙) stars with ages color-coded from violet (youngest) to red (oldest). A movie of the simulation is available at https://www.bilibili.com/video/BV18vKRzzEnj.

In the text

Fig. 2. Evolution of the distance between two galaxies (upper panel) and the SF history of the merger system (bottom panel). The distance is measured between the centers of mass of the background stars in the galaxies. Bottom panel: SFR averaged over a timescale of 5 Myr. The blue curve represents the total SFR of the two dwarfs in the merger, while the orange curve represents twice the SFR of an isolated dwarf galaxy. The four stages of the merger–“Approach”, “1st peak”, “Plateau”, and “2nd peak” –are marked with shaded regions in various colors.

In the text

Fig. 3. Evolution of galaxy-wide ISM properties and depletion time. Upper panel: Evolution of the fraction of nH > 100 cm−3 dense gas (red curve) and gas in identified clouds (dense gas clumps with virial parameter α < 10, blue curve). Bottom panel: Galaxy-wide depletion times calculated by the total gas mass (black curve), dense gas mass (red curve), and total cloud mass (blue curve). The depletion time of the gas in the clouds fluctuates around its median value of 7.4 Myr (dashed blue line), with a 16–84 percentile range of 4.7–12 Myr.

In the text

Fig. 4. Evolution of clouds in a sequence spanning 14 Myr. The first two panels show the location of the zoom-in region in the merger system. The third to the last panels show the column density distribution inside a (200 pc)3 box. In these panels, the color maps refer to the gas surface density distribution of the zoom-in regions, and the dots represent the stars color-coded according to their age. The polygons in the zoom-in panels depict the 2D convex hulls of the clouds identified from a specific evolution graph. The evolution paths in this graph merge at 158 Myr and share a node (green polygon), where we evaluate the total baryon mass Mtot and surface density Σtot for the SFE calculation (see Section 4.3). The clouds form in the upper-right and lower-left regions at the beginning (see the panel of 153 Myr for an example) and accumulate at the center (155 Myr). Star clusters appear at 156 Myr, but the cloud continues to accumulate mass from the upper-right, forming a large single cloud at 159 Myr. The cloud continues to form stars while collapsing. Finally, the stellar feedback rapidly disperses the clouds, leaving bubbles within 4 Myr. The colorful polygons appearing between 155 Myr to 160 Myr indicate the main paths sharing the same maximum baryon mass node in this evolution graph.

In the text

Fig. 5. Cloud velocity dispersion as a function of effective radius. The clouds are selected from the snapshot spaced at 10 Myr intervals to avoid duplicate counts and eliminate bias toward the long-lived clouds. The points in the top panel are color-coded by their virial parameters, while those in the bottom panel are color-coded by tidal strength λtid. Clouds in stronger tidal environments typically have larger velocity dispersions and virial parameters. The red line represents the observations of MW clouds by Solomon et al. (1987), and the red-shaded region represents the observed 0.22 dex dispersion of log10σv. The gray line and shaded region represent our best-fitting relation and the log10σv dispersion of 0.31 dex. Our sample fitting closely follows the observations, although the scatter in the sample is slightly larger than that of the observed MW clouds, due to a rather tolerant selection criterion designed to include dynamically diverse clouds.

In the text

	Fig. 6. Cloud lifetime distributions at different merger stages. The dashed gray line denotes the maximum likelihood estimate fitting of the lifetime for all stages. The clouds in all merger stages follow the same exponentially decaying relation with a characteristic time of 1.48 Myr.
In the text

Fig. 7. Distribution of clouds in the ${\lambda }_{\mathrm{tid}}^{\mathrm{mid}}$–Mcloud, max plane, color-coded by their average lifetimes. ${\lambda }_{\mathrm{tid}}^{\mathrm{mid}}$ is the median tidal strength experienced by clouds throughout their lifetime, and Mcloud, max is the maximum mass the clouds ever reached during their lifetime. The solid black line and shaded region represent the median relation and the 16–84 percentile range, respectively.

In the text

Fig. 8. Relation between average cloud separation l and tidal strength. The cloud sample is the same as that presented in Fig. 5. The top panel is color-coded by the nature of tides Λ. The dashed gray line represents the fitting of λtid ∝ l−0.43. The bottom panel shows the PDFs of the average separation l in different merger stages. Most of the clouds are in fully compressive tidal fields (Λ > 1). The average separation l extends and shifts to smaller scales as the merger progresses and the SFR increases.

In the text

Fig. 9. Distribution of integrated SFE ϵint as a function of the total baryon surface density Σtot. Top panels: Distribution of all cloud samples with SF. Bottom panel: Results for the clouds formed more than 100 M⊙ in newly formed stars. The data points are color-coded by the normalized density in log-space estimated by a Gaussian kernel. The green curves show the median values, and the shaded regions indicate the 16–84 confidence levels. The dashed curves denote Equation (5) with the best-fitting parameters from Li et al. (2019) for fboost = 1, 3, and 10. Although the scatter of the simulated results is large, the median curve roughly follows the theoretical curve with fboost = 3. The 16–84 confidence level lies between the fboost = 1 and 10 curves at all the stages.

In the text

	Fig. 10. Joint probability distribution of ϵint and Σtot color-coded by the tidal strength. Here λtid is the value when the cloud reaches its maximum baryon mass.
In the text

Fig. 11. Cloud mass functions (top panel) and young cluster (< 5 Myr) initial mass functions (CIMFs, bottom panel) of the merger system at different stages. The dashed gray line denotes the power-law relation with a slope of −1.5, −2, and −2.5, respectively. A starburst induced by the merger leads to the formation of massive clusters, especially during the second SFR peak. The maximum cloud masses are (1.8,  23,  11,  90)×103 M⊙ for the Approach, 1st peak, Plateau, and 2nd peak stages, respectively. The maximum cluster masses are (0.93,  3.4,  1.8,  26)×103 M⊙ for these four stages.

In the text

	Fig. 12. Most massive young cluster mass vs. galaxy-wide SFR (Mcluster, max–Ṁ★) relation. The simulation data are color-coded by stage of the merger. The dashed red line represents the theoretical relation derived following Yan et al. (2017). The gray crosses represent observational data from Weidner & Kroupa (2004) and Schulz et al. (2015), which have a typical uncertainty of 0.3 dex.
In the text

Fig. 13. 2D (top panel) and 3D (bottom panel) “tuning fork” diagrams from our simulation at different phases. These diagrams reflect the deviation of the H2 depletion time centered on H2 (τH2, gas, upper branch) or recent SFR peaks (τH2, SFR, lower branch) from the galaxy-averaged value as a function of the spatial scale L. The curves of the 1st peak, Plateau, and 2nd peak are multiplied by constants of 10, 100, and 1000, respectively, to avoid overlapping. The vertical lines indicate the scale when τH2, gas exceeds τH2, SFR by 50%.

In the text

Fig. 14. CIMFs for the simulation with weakened radiative feedback. The red curve shows the result of 10 M⊙ resolution stellar radiation luminosities reduced by a factor of 0.1 across all bands. Solid black and dashed gray curves represent the results from the simulation with normal radiative feedback at 2 M⊙ and 10 M⊙, respectively. The CIMF is not sensitive to the numerical resolution but is very susceptible to the intensity of radiative feedback.

In the text