© The Authors 2026

1 Introduction

Stars are generally thought to form via the fragmentation and collapse of giant molecular clouds. Most are born in dense and crowded stellar environments (Lada & Lada 2003) known as stellar complexes. Several decades of accumulated data indicate that star clusters and associations also form in groups within the largest stellar complexes (Efremov 1995). A subset of star clusters are formed nearly simultaneously and in close proximity (de la Fuente Marcos & de la Fuente Marcos 2008), where transient but powerful dynamical interactions can cause the star clusters to be disrupted or to merge (Fellhauer & Kroupa 2005; Mora et al. 2019). We refer to such clusters as binary (or higher-multiplicity) star clusters.

When a pair of clusters forms from the same parent molecular cloud, the clusters generally have comparable ages, distances, and metallicities (Qin et al. 2023) and remain gravitationally bound. Such systems are termed primordial binary clusters. In recent decades, extensive studies have been carried out on the dynamical evolution of primordial binary clusters in both the Magellanic Clouds and the Milky Way (Sugimoto & Makino 1989; Portegies Zwart & Rusli 2007; Priyatikanto et al. 2016; Darma et al. 2021; Li & Spurzem 2026). Sugimoto & Makino (1989) investigated the evolution of primordial binary clusters in the Magellanic Clouds through N-body simulations. Their results show that tidal effects lead to orbital circularisation, and that the spin of each cluster becomes synchronised with the orbital rotation, a process that depletes orbital angular momentum. At the same time, escaping stars carry angular momentum away from the system, with both processes contributing to the shrinkage of the orbit. Portegies Zwart & Rusli (2007) incorporated the stellar mass spectrum and stellar evolution, which had previously been neglected in simulations of primordial binary clusters, into their models in a self-consistent fashion (Makino et al. 1991). In addition to stellar escape, the early orbital evolution of binary clusters is also affected by mass loss resulting from stellar evolution. Mass loss from stellar winds and supernova explosions during the first ∼30 Myr causes the binary to expand and the orbit to become more eccentric.

The evolution of primordial binary clusters under the influence of the Galactic tidal field was investigated by de la Fuente Marcos & de la Fuente Marcos (2010) using N-body simulations. Their simulations show that the initial orbital elements (semi-major axis and eccentricity) and the mass ratio of the cluster pair can lead to different evolutionary paths, including merging, extreme tidal distortion, and ionisation. In general, close pairs merge on short timescales, secondaries in low mass-ratio pairs are rapidly disrupted, and wide pairs are quickly ionised by the background tidal field owing to mass loss. In the merging cases, the trajectories of the binary clusters become more complex and deviate from the typical spiral pattern observed in the absence of an external tidal field (Makino et al. 1991; Priyatikanto et al. 2016). In addition, Priyatikanto et al. (2016) show that even binary clusters with identical initial orbital elements can experience completely different dynamical fates, depending on their initial orientations relative to the Galaxy and their initial orbital phase. Under the influence of the external potential, the orbital evolution of binary clusters can involve orbital reversal, spiral-in motion, and vertical oscillations about the Galactic plane before the clusters reach their final configurations (merger or separation). Recently, Li & Spurzem (2026) developed a direct N-body simulation code, NbodyCP, that can generate pairs of star clusters with different ages and metallicities, providing a useful tool for investigating the formation and evolution of star clusters with multiple stellar populations.

In addition to simultaneous formation, tidal capture (van den Bergh 1996) has been proposed as a significant mechanism for the formation of binary star clusters. Clusters of non-common origin can become gravitationally bound through a close encounter followed by the loss of angular momentum (Mora et al. 2019). Several numerical simulations have investigated the dynamical processes involved in close encounters of globular clusters (GCs) in the Magellanic Clouds (Devadas Rao et al. 1987; de Oliveira et al. 1998, 2000; Gavagnin et al. 2016). Interest in the dynamics of close encounters between star clusters is largely motivated by the fact that most observed GCs comprise significant multiple stellar populations (Carretta et al. 2007; Milone et al. 2008; Carretta et al. 2009; Pancino et al. 2010). The progenitor GCs involved in a merger can possess different compositions, ages, and/or evolutionary histories (van den Bergh 1996). Therefore, the merging of clusters with different metallicities can naturally account for the metallicity spreads observed in GCs (Gavagnin et al. 2016). Similarly, close encounters between open clusters (OCs) are not uncommon in the Milky Way (Palma et al. 2025; Li & Zhu 2025; Liu et al. 2025; Zhu et al. 2025). Although several studies have investigated the dynamics of close encounters involving massive GCs, there remains a lack of dynamical studies focusing on close encounters between low-mass, more loosely structured OCs, particularly regarding the detailed mechanisms governing their merging processes. Because OCs lack the strong gravitational effects of GCs, the dynamical processes between OCs may differ significantly from those between GCs.

Therefore, in this work we investigated the dynamical interactions between low-mass OCs across a range of mass and density ratios, under different orbital configurations. This paper is organised as follows. In Sect. 2 we describe the methods and initial conditions adopted in our simulations. The numerical results are presented in Sect. 3. Our summary and a discussion are presented in Sect. 4.

2 Methods and initial conditions

The aim of this simulation is to investigate the interaction mechanisms that arise during close encounters between OCs, to identify the key structural parameters of the clusters (e.g. mass and density) that determine whether they merge or become disrupted, and to examine the physical properties of the resulting merger remnants. We employed the high-performance N-body simulation code PETAR (Wang et al. 2020) to carry out the simulations. The code can accurately handle an arbitrary fraction of multiple systems (e.g. binaries and triples) while maintaining high performance through hybrid parallelisation using Message Passing Interface (MPI), Open Multi-Processing (OpenMP), Single Instruction Multiple Data (SIMD) instructions, and Graphics Processing Unit (GPU) acceleration. McLuster is a powerful tool that can be used both to set initial conditions for N-body calculations and to generate simulated star clusters for direct study (Küpper et al. 2011). We used McLuster to generate the simulated OCs. The positions and velocities of the stars in each simulated OC follow a Plummer density profile (Plummer 1911). The Plummer density profile is defined as $$\rho (r)=\frac{3M}{4\pi R_p^3}\frac{1}{{\left(1+{\left(r/R_p\right)}^2\right)}^{5/2}},$$(1)

where M is the total mass of the cluster, and R_p is the Plummer radius, a scale parameter that sets the characteristic size of the cluster’s central regions (de la Fuente Marcos & de la Fuente Marcos 2010). Although the Plummer profile is in principle infinitely extended, the adoption of a realistic Galactic potential results in a natural truncation at the theoretical tidal radius of the cluster. The star clusters follow the canonical Kroupa (2001) initial mass function, with a slope of α₁ = −1.3 for 0.08–0.5 M_⊙ and the Salpeter (1955) slope of α₂ = −2.3 for m > 0.5 M_⊙. The lower and upper limits of the initial mass function are 0.08 M_⊙ and 100 M_⊙, respectively.

The mass and virial radius of the second OC (OC2) were kept fixed, with a total mass of 1024 M_⊙ and a virial radius (R_vir) of 1 pc. To investigate the effects of relative mass and density during close encounters between OCs, we adopted a range of mass ratios and density ratios consistent with those used by Gavagnin et al. (2016). The mass ratios considered were 1, 2, and 4, while the density ratios were 0.25, 0.5, 1, 2, and 4. Different mass and density ratios were achieved by varying the mass and R_vir of OC1. For example, in model M1ρ0.5, OC1 has the same mass as OC2 but half its density, which is accomplished by setting the R_vir of OC1 to R_vir,OC1 = 1.26 R_vir,OC2, so that the density ratio between the two clusters is 0.5. For a Plummer sphere in virial equilibrium, the relation between the R_vir and the half-mass radius (R_h) is (Aarseth & Fall 1980; de la Fuente Marcos & de la Fuente Marcos 2010) $$R_{\text{vir}}=\frac{16}{3\pi \times 1.305}{R}_{\text{h}}.$$(2)

Using this relation, the R_h of each cluster is computed and subsequently used to generate the simulated clusters.

Except for the initial orbital parameters, the N-body simulations of primordial binary clusters and close-encounter cluster pairs are essentially identical. The former are typically modelled as gravitationally bound systems on elliptical orbits, whereas the latter involve gravitationally unbound systems on hyperbolic or parabolic orbits (van den Bergh 1996). For a comprehensive analysis, both parabolic and hyperbolic orbits are considered in this work. To define these two orbits, we fixed the minimum peri-centre distance (r_peri = 5, 10, 15, 20, 25, and 30 pc) and specified the orbital eccentricity (e_p = 1 and e_h = 2). Here, the subscripts p and h denote parabolic and hyperbolic orbits, respectively. The orbital plane was chosen as the xy-plane with the x-axis in the direction of the closest approach (Devadas Rao et al. 1987). The initial conditions (at t = 0) are as follows: (i) OC2 is placed at the coordinates (x, y) = (0, 0); and (ii) OC1 is positioned at a distance of 4 × max(R_t1, R_t2), where R_t1 and R_t2 are the tidal radii of OC1 and OC2, respectively. Beyond this distance the tidal effects are negligible (Innanen et al. 1972). The initial relative velocity at this distance is obtained from the two-body formulae (see Table 1). The initial binary fraction is zero, and no binaries form during the subsequent evolution. Each set of initial conditions serves as the input for the PETAR code. At regular intervals of 1 Myr, the simulations generate output snapshots in which the essential data, including the masses and 3D positions and velocities of all particles, are recorded for later analysis. All snapshots are recorded in astrophysical units (M_⊙, pc, and Myr), with velocities expressed in pc Myr⁻¹, as adopted by PETAR. The computation was stopped when the system showed no further change in its properties (Devadas Rao et al. 1987). The present study concentrates on the tidal interaction between the two clusters. Stellar and binary evolution, as well as the influence of the Galactic tidal field, are omitted and will be addressed in forthcoming work.

3 Results

Although at the initial stage of each simulation it is known to which OC each particle belongs, this becomes increasingly ambiguous as the system evolves, particularly when the two OCs pass through pericentre. Therefore, accurately determining to which OC each particle is bound in every snapshot is essential for reliably identifying the OC’s centre of mass, mass loss, and structural evolution. To achieve this, we adopted the method introduced by Dutta Chowdhury et al. (2020). A brief summary of the implementation steps is provided below:

1. For each OC in every snapshot, the membership is initially assigned to the particles that belonged to it at t = 0, and the position and velocity of the centre of mass of this collection are then determined;

2. For each OC, this collection of particles, together with those particles that do not belong to it but whose masses are set to zero, is then used to compute the gravitational potential using the softened Newtonian potential formula;

3. The potential calculated using the softened Newtonian form for each particle is employed to obtain its binding energy in the centre-of-mass frame of each OC, as determined in the previous iteration (or in step 1 for the initial iteration). The OC membership is then updated by assigning each particle to the OC with respect to which it has the most negative binding energy. Particles that possess positive binding energy relative to all OCs are not assigned to any of them and constitute the collection of unbound particles;

4. For each OC, the position and velocity of its centre of mass are updated using the 50% most strongly bound particles from the bound-particle collection determined in step 3.

Steps 2, 3, and 4 are repeated, using each OC’s collection of bound particles as determined in the previous iteration until the relative differences between the positions and velocities of the centres of mass obtained in two successive iterations converge to better than 1 × 10⁻³ for every OC. In the above procedure, we computed the gravitational potential between particles using a softened Newtonian potential, given by $${\varphi }_i=-G{\displaystyle \sum }_{j\ne i}\frac{m_j}{\sqrt{r_{ij}^2+{\epsilon }^2}},$$(3)

where G is the gravitational constant, m_j is the mass of particle j, r_ij is the distance between particles i and j, and ε = 0.01 pc is the softening parameter, introduced to prevent divergences arising from close encounters between stars (de Oliveira et al. 1998).

3.1 General features

Across all simulations, we find that significant dynamical interactions such as merging and tidal disruption occur exclusively in the parabolic models. In the hyperbolic models, although close interactions strip some stars from the outskirts of the lower-mass cluster, the main body remains stable, and neither complete disruption nor merging occurs. Although hyperbolic encounters are the most common, as the phase space for encounters increases with the cube of the encounter velocity and with the cube of the impact parameter, the probability of a merging encounter drops sharply once the orbit becomes even mildly hyperbolic (Gavagnin et al. 2016).

Because the parameter space explored in this study is relatively broad, the detailed dynamical evolution in the merging cases is not identical. Nevertheless, the mergers share a common characteristic: during a pericentric passage, tidal forces couple most strongly to particles in the outer regions whose angular momentum is most closely aligned with the orbital angular momentum. These particles experience a tidal field that is nearly in phase with their orbital motion and therefore display an enhanced response (White 1978). Such tidal interactions heat the outer particles of the clusters, raising their energies above the binding energy of the system and causing them to escape. As a consequence, the orbits of the interacting clusters evolve and gradually decay, ultimately leading to a merger (Barnes 2016).

Varying the cluster mass ratio results in distinct collision behaviours and different levels of disruption of the companion cluster (OC2). As examples, Figs. 1 and 2 illustrate the merging processes for systems with different mass ratios. Figure 1 shows projections of model M1ρ0.25 (p10, where the number denotes r_peri) onto the x–y plane at various stages of its evolution. Tidal interactions between two extended, self-gravitating stellar systems transfer energy and momentum from their relative motion into the internal degrees of freedom (Miocchi et al. 2006). This transfer of energy begins well before the system reaches pericentre. A dynamical effect noted by van Albada & van Gorkom (1977) is also apparent in Fig. 1. When the two stellar systems first pass through one another, they experience a strong inward impulse. As a result, immediately after the encounter they contract significantly relative to their equilibrium configuration. This contraction is followed by a rebound and a rapid expansion of the outer regions of the systems, which acquire a substantial fraction of the orbital energy (White 1978; Miller 1978). The expansion occurs primarily within two conical regions surrounding the symmetry axis of the collision (see Fig. 1d).

In contrast to the smooth merging observed in equal-mass models, the dominant merger mechanism in the unequal-mass models is the accretion of the lower-mass cluster by the more massive one following its disruption. Figure 2 shows the projection of model M4ρ4 (p5) onto the x–y plane at various times during its evolution. When the more massive cluster passes through pericentre, strong tidal forces impart the largest accelerations to stars in the outskirts of the lower-mass cluster, leading to significant tidal heating and driving a rapid expansion. In addition, the characteristic timescale for equipartition of kinetic energy is much shorter for the lower-mass cluster, contributing to its more rapid internal dynamical evolution. The coupling of these two processes effectively drives the lower-mass cluster away from equilibrium, as the total potential energy per unit mass decreases rapidly (de la Fuente Marcos & de la Fuente Marcos 2010). The lower-mass cluster in proximity to the massive cluster becomes stretched by tidal forces, forming a pronounced tidal bridge. The formation of the bridges consumes part of the orbital energy and angular momentum, and the asymmetric distribution of matter within these structures exerts a torque that further diminishes the orbital energy and angular momentum of the two cores (Toomre & Toomre 1972). To illustrate this structure more clearly, all particles are plotted in the same colour in Fig. 2h. A similar structure and dynamical process is also evident in Fig. 1 of de Oliveira et al. (2000).

In the equal-mass models, the dissipation of orbital energy during the first close encounter converts the initially parabolic orbit into a nearly elliptical one, allowing the clusters to orbit each other before ultimately merging. In the unequal-mass models, however, the clusters merge rapidly after only one close interaction. Consequently, the merger timescale in the unequal-mass cases is much shorter than in the equal-mass cases (de la Fuente Marcos & de la Fuente Marcos 2010). The detailed dynamical processes are discussed in the following sections.

Fig. 1 Projection of model M1ρ0.25 (p10) onto the x–y plane at nine different evolutionary times. Black points represent OC2, located at the origin. Red points denote OC1, and grey squares indicate unbound particles.

Fig. 2 Projection of model M4ρ4 (p5) onto the x–y plane at nine different evolutionary times. To clearly illustrate the bridge structure, all particles in panel h are displayed as black points.

Fig. 3 Time evolution of the relative separation between the two clusters, r12 (top panel), the relative orbital energy, E12 (middle panel), and the orbital eccentricity, e (bottom panel), for model M1ρ0.25 (p10).

3.2 Orbital energy dissipation

Part of the relative orbital energy (E₁₂) in both parabolic and hyperbolic encounters can be lost through subsequent interactions between the two systems. If the interaction is sufficiently strong, the resulting energy loss can drive the system into a bound state with E₁₂<0 (González-García & van Albada 2005). In Fig. 3, the upper, middle, and lower panels show the time evolution of the relative separation r₁₂ between the two clusters, E₁₂, and the orbital eccentricity e, respectively, for model M1ρ0.25 (p10). E₁₂ is computed following the formulation of Dutta Chowdhury et al. (2020) and is defined as $$E_{12}=\frac{1}{2}\mu v_{12}^2+W_{12},$$(4)

where v₁₂ is the relative velocity of the two clusters, µ is the reduced bound mass, and W₁₂ is their mutual gravitational potential energy. For two mass distributions with densities ρ₁ and ρ₂, the mutual potential energy is $$W_{12}=\frac{1}{2}\int {\rho }_1(r){\text{Φ}}_2(r){\text{d}}^{3}r+\frac{1}{2}\int {\rho }_2(r){\text{Φ}}_1(r){\text{d}}^{3}r,$$(5)

where Φ₁ and Φ₂ are the gravitational potentials generated by OC1 and OC2, respectively. In contrast to Dutta Chowdhury et al. (2020), who computed W₁₂ by modelling each cluster as a Plummer sphere with its instantaneous bound mass and projected R_h, here W₁₂ was evaluated directly from the softened Newtonian potentials of the particles belonging to the two clusters. The instantaneous e was derived from the two-body orbital equations.

At the beginning of the simulation, e is close to 1 and E₁₂ is approximately zero, consistent with the characteristics of a parabolic orbit. After the pericentric passage, interactions between the clusters dissipate part of the orbital energy, converting it into internal energy. As a result, E₁₂ becomes negative and the system transitions to a gravitationally bound elliptical orbit with e<1. We note that the actual r_peri (∼5 pc) is significantly smaller than the theoretical value of 10 pc; this is due to the loss of orbital energy during the approach caused by strong tidal interactions, which lead to orbital contractions. After the first encounter, the separation between the two clusters increases to about 40 pc, after which they approach each other again and eventually merge. For comparison, the evolution of r₁₂, E₁₂, and e for model M1ρ1 (h10) is shown in Fig. 4. Throughout the evolution, E₁₂ remains positive and e stays close to 2, consistent with the initial orbital configuration. Moreover, r_peri (∼10 pc) agrees with the expected value. These results indicate that no significant interaction takes place between the two clusters.

As noted in Sect. 3.1, none of the hyperbolic models in our simulations results in a merger of the clusters, whereas all merging or disruption events occur in the parabolic models. In Table 2, we present the merger times for all parabolic models. We find that mergers become more likely as r_peri decreases. In the simulations of González-García & van Albada (2005), all runs with zero orbital energy result in a merger, except for the case in which the impact parameter exceeds the sum of the radii of the two systems. However, in our models, mergers do not necessarily occur even when the clusters partially overlap, which may be related to the relatively low masses of the systems considered.

Fig. 4 Same as Fig. 3 but for model M1ρ1 (h10).

3.3 Tidal synchronization

As discussed in previous works (Sugimoto & Makino 1989), tidal interactions cause the spin of each stellar system to tend towards synchronisation with the orbital rotation, at the expense of the orbital angular momentum, L_orb. This mechanism shows a qualitative analogy to stellar evolution theory in close binary systems (Meibom & Mathieu 2005), although important differences exist between the constituent stars in clusters and the gas particles in stellar interiors. Moreover, the stellar density distribution within a cluster is fundamentally different from the gas density distribution inside a star (Sugimoto & Makino 1989). The loss of L_orb through escaping particles, together with its transfer into the spin angular momentum of the clusters, leads to a reduction in the orbital separation. When the separation reaches a critical value, the pair rapidly merges into a single cluster.

Tidal synchronisation is also observed in our merger cases. To clearly characterise this synchronisation, we defined the net angular momentum of a star cluster as the sum of the angular momenta of its member stars: $$L={\displaystyle \sum }_i{m}_i{r}_i\times v_i.$$(6)

Here, the positions r_i and velocities u_i are measured relative to the centre of mass of the cluster. If a cluster exhibits significant rotation (velocity anisotropy), the net angular momentum L = |L| is non-zero, and the direction of the angular moment direction $\widehat{L}\equiv L/\left|L\right|v$ defines the rotation axis of the cluster (Priyatikanto et al. 2016). It should be noted that, although the initial clusters generated with McLuster are non-rotating by construction, a small but non-zero spin angular momentum is present at t = 0 (see Fig. 5). This residual spin originates from finite-N sampling noise introduced during the realisation of the isotropic velocity distribution, rather than from any physically imposed rotation. Its magnitude is typically less than a few per cent of L_orb and therefore has a negligible impact on the dynamical evolution and does not affect our conclusions. The orbital angular momentum of the cluster pair, L_orb, is given by $$L_{\text{orb}}=\mu R\times V,$$(7)

where R = r₂ − r₁ and V = u₂ − u₁ are the relative position and velocity vectors between the centres of mass of OC1 and OC2, respectively, and µ is the reduced mass of the system. The magnitude of the orbital angular momentum is then given by L_orb ≡ |L_orb|, which quantifies the angular momentum associated with the relative orbital motion of the two clusters about their common centre of mass.

Figure 5 shows the time evolution of L_orb and the spin angular momenta of the two clusters for models M1ρ0.25 (p10) and M4ρ4 (p5). In the equal-mass models (Fig. 5a), during the early stages of the simulations, internal dynamical processes within the clusters cause some particles to gain energy and escape from the system, leading to a gradual decline in L_orb. Escaping stars remove part of the angular momentum, and as the system evolves and the orbital separation decreases, tidal synchronisation and tidal friction become significant (Priyatikanto et al. 2016). As the two clusters gradually begin to interpenetrate (∼70 Myr), part of L_orb is transferred to the spin angular momentum of the clusters, leading to a gradual increase in their spin angular momentum. As discussed in Sect. 3.1, particles in the outer regions of the clusters whose angular momentum is most closely aligned with L_orb experience the strongest tidal coupling. Such tidal interactions efficiently heat these outer particles, causing some of them to escape from the system. This also explains the pronounced rise and subsequent decline of the spin angular momentum seen in Fig. 5a: while tidal interactions initially transfer L_orb to cluster rotation, the later escape of tidally heated outer stars removes angular momentum from the system.

After the two clusters pass through pericentre (∼90 Myr), dissipation of orbital energy renders the system gravitationally bound on an elliptical orbit. During the subsequent evolution, the spin angular momenta of both clusters gradually increase, and their spins become progressively synchronised with the orbital rotation. For clarify, the distributions of spin angular velocity ω_spin within a cluster are computed in cylindrical shells, together with the corresponding orbital angular velocity ω_orb. In axially symmetric systems, the natural decomposition of the velocity field is obtained by adopting cylindrical coordinates centred on the cluster centre of mass (Kim et al. 2008). For each star, we defined the cylindrical radius (R_i) as the perpendicular distance from the rotation axis, $$R_i=\left|r_i-\left(r_i\cdot \widehat{L}\right)\widehat{L}\right|.$$(8)

The global spin angular velocity of a cluster can be characterised as $${\omega }_{\text{spin}}=\frac{\left|L\right|}{I_{\Vert }},$$(9)

where $I_{\Vert }={\displaystyle {\sum }_i{m}_i{R}_i^2}$ is the moment of inertia about the rotation axis. The instantaneous orbital angular velocity is then given by $${\omega }_{\text{orb}}=\frac{|R\times V|}{|R{|}^2},$$(10)

where R and V are defined as in Eq. (7). The top panels of Fig. 6 show the ω_spin and ω_orb of model M1ρ0.25 (p10) at several evolutionary stages. The ω_spin tends to approach the ω_orb, indicating tidal synchronisation, although appreciable delays are clearly present. When L_orb becomes insufficient to supply the angular momentum required for spin synchronisation, the clusters undergo merger at a critical separation (Sugimoto & Makino 1989). We note that, after the merger, the remaining L_orb is largely transferred into the spin angular momentum of the remnant, producing a merger remnant with significant rotation (Alladin et al. 1985). The rotational properties of the merger remnant will be discussed later.

In the unequal-mass model (Fig. 5b), L_orb exhibits two pronounced drops. The first occurs when the two clusters begin to interpenetrate, as tidally heated outer stars escape and remove angular momentum from the system. The second drop coincides with the rapid disruption of the lower-mass cluster, during which its bound mass decreases sharply, leading to a further reduction in L_orb. We also find that the spin angular momentum of the merger remnant does not inherit all of the remaining L_orb. This is because the massive cluster accretes only ∼71% of the stars from the disrupted companion, while the escaping fraction (∼29%) carries away a substantial amount of angular momentum. We can see that tidal synchronisation also occurs in the unequal-mass model, although appreciable delays are present (bottom panels of Fig. 6).

Fig. 5 Time evolution of the orbital angular momentum, Lorb, and the spin angular momenta of the two clusters, Lspin,1 and Lspin,2, for models M1ρ0.25 (p10; panel a) and M4ρ4 (p5; panel b).

Fig. 6 Distributions of the spin angular velocity, ωspin (solid lines), and the orbital angular velocity, ωorb (dashed lines), as functions of the cylindrical radius (R⊥) for models M1ρ0.25 (p10; panels a–c) and M4ρ4 (p5; panels d–f) at different evolutionary times, as indicated in each panel.

Fig. 7 Rotation curves of the merger remnants for models M1ρ0.25 (p10; panel a) and M4ρ4 (p5; panel b) at different evolutionary times, as indicated. The inset shows, in logarithmic scales, a representative rotation curve of an isolated rotating cluster taken from Kim et al. (2004), illustrating the rapid decline of the rotational velocity beyond Rh.

3.4 Rotation of merger remnants

As noted above, the merger process imparts rotation to the merger remnant in the same sense as the orbital motion (de la Fuente Marcos & de la Fuente Marcos 2010). To better understand the rotational properties of the merger remnant, we computed the rotational velocity v_θ of each star. The v_θ of each star was defined as the velocity component perpendicular to the rotation axis, projected onto the equatorial plane of the merger remnant (Priyatikanto et al. 2016), and was calculated as follows: $$v_{\theta ,i}=\frac{\left(r_i\times v_i\right)\cdot \widehat{L}}{R_i}.$$(11)

Rotation can also be characterised by a rotation curve, defined as the average rotational velocity as a function of the cylindrical radius. In Fig. 7 we present the rotation curves of the merger remnants for models M1ρ0.25 (p10) and M4ρ4 (p5) at several evolutionary stages following the merger. In the equal-mass models, the merger remnants with significant rotation exhibit an asymmetric distribution of v_θ about zero (Priyatikanto et al. 2016). Within R_h, the merger remnant exhibits an approximately solid-body rotation. However, outside R_h, the mean v_θ does not show the same rapid decline as typically observed in an isolated rotating cluster (Meylan & Mayor 1986; Kim et al. 2004, 2008; Fabricius et al. 2014). The complex dynamical evolution of the binary system prior to merging is likely responsible for this behaviour. We note that beyond R_h, the rotation curve shows a clear secondary rise after an initial decline. This feature may originate from incomplete phase mixing and the retention of L_orb in the outer halo of the merger remnant.

In the unequal-mass models, the inner regions of the merger remnant also exhibit an approximately solid-body rotation. However, the extent of the rigidly rotating region is larger than that in the equal-mass case, in agreement with the results of Gavagnin et al. (2016). The rotation curve of the merger remnant shows a clear rise in the halo region. This behaviour can be attributed to particles originating from the disrupted lower-mass cluster, which are distributed in the halo of the remnant. The angular momentum associated with the lower-mass cluster is preferentially deposited in the outer parts of the remnant (Gavagnin et al. 2016). To verify this interpretation, Fig. 8 shows equatorial-plane projections of the merger remnant at different evolutionary times, with stars colour-coded according to their original cluster membership. While the massive cluster remains centrally concentrated, the stars accreted from the lower-mass cluster undergo a continuous and anisotropic expansion after the merger. These stars gradually populate the outer halo of the merger remnant, while retaining a large fraction of their original orbital angular momentum. This redistribution naturally accounts for the rise in rotational velocity observed in the outer regions of the remnant. In addition, we note a rapid expansion of the merger remnant, which may be related to its intrinsic rotation (Fabricius et al. 2014).

3.5 Structure of remnants

In this section we compare the structural properties of the initial clusters with those of the merger remnants, using models M1ρ0.25 (p10) and M4ρ4 (p5) as representative examples. To do this, we divided the remnant into a series of concentric spherical shells centred on the remnant centre. The stellar mass density in the i-th shell is then calculated as $${\rho }_i=\frac{{{\displaystyle \sum }}_j{m}_j}{V_i},$$(12)

where Σ_jm_j is the total stellar mass contained within the i-th shell and V_i is the volume of that shell. The curves shown in Fig. 9 present the resulting radial stellar mass density profiles of the merger remnants. For comparison, the initial Plummer models are shown as dotted curves. The merger remnant exhibits an extended halo with a radial density dependence of ρ ∝ r⁻⁴, similar to that observed in flattened GCs in the Large Magellanic Cloud (Elson et al. 1987), which is a direct consequence of strong gravitational interactions (Aguilar & White 1986). The absence of data points in the innermost regions of the merger remnants is caused by limited particle statistics within very small spherical shells after the merger. In these regions the number of particles per bin falls below the adopted threshold required for reliable density estimates, and the corresponding bins are therefore excluded from the profiles. Nevertheless, within the resolved radial range the merger remnants exhibit systematically lower inner densities than the initial Plummer models, indicating dynamical heating and expansion during the merger process. The R_h of the merger remnant increases by nearly an order of magnitude compared to the progenitor clusters, and continues to expand in its inner regions as the system evolves. This behaviour differs markedly from mergers that originate from initially bound orbits, in which the density distribution in the inner regions remains broadly similar to that of the initial model (Okumura et al. 1991).

The core of the merger remnant becomes more diffuse, while the outer regions exhibit pronounced expansion, a behaviour that is also reflected in the distribution of particle binding energies. Figure 10 presents the distributions of particle binding energies for the initial Plummer model and the merger remnant of models M1ρ0.25 (p10) and M4ρ4 (p5). The differences in the shapes of the histograms in Fig. 10 demonstrate that the final structure differs markedly from that of the initial model. Compared to the broad distribution of particle binding energies in the initial model, the peak of the binding-energy distribution of the merger remnant shifts towards zero. This shift is partly due to the formation of an extended halo, and partly driven by the overall expansion and dynamical heating of the remnant induced by violent relaxation during the merger process. In this sense, the merger is characterised by strong phase-space mixing and dynamical heating, which drive the remnant towards a new virial equilibrium while producing an overall expansion and an extended halo (Sugimoto & Makino 1989).

Fig. 8 Equatorial-plane projections of the merger remnant for model M4ρ4 (p5) at different evolutionary times. Stars are colour-coded according to their original cluster membership.

Fig. 9 Radial stellar mass density profiles of the merger remnants for models M1ρ0.25 (p10; panel a) and M4ρ4 (p5; panel b) at different evolutionary times. The dashed and dotted curves show the initial Plummer models of OC1 and OC2, respectively.

4 Summary and conclusions

The coeval formation scenario, in which clusters form simultaneously within the same molecular cloud and are therefore expected to share similar space velocities, ages, and metallicities, is considered to be the primary formation mechanism of binary star clusters. Such systems are referred to as primordial binary clusters. The dynamical evolution of primordial binary clusters has been extensively studied in different configurations, including stellar evolution and the influence of Galactic tidal fields (Portegies Zwart & Rusli 2007; de la Fuente Marcos & de la Fuente Marcos 2010). In these N-body simulations, primordial binary clusters are typically modelled on Keplerian orbits with different semi-major axes and eccentricities (de la Fuente Marcos & de la Fuente Marcos 2010; Darma et al. 2021). In addition to simultaneous formation, tidal capture has been proposed as an important mechanism for the formation of binary star clusters (van den Bergh 1996), whereby clusters of different origins become gravitationally bound following a close encounter and subsequent loss of angular momentum (Mora et al. 2019). The merging of clusters with different ages and metallicities can naturally explain the metallicity spread observed in GCs that host multiple stellar populations (Gavagnin et al. 2016). The dynamical processes involved in close encounters between star clusters of non-common origin have been investigated in several simulation studies, but they have primarily focused on massive GCs in the Large Magellanic Cloud (Devadas Rao et al. 1987; Gavagnin et al. 2016). Compared with the Large Magellanic Cloud, the Milky Way hosts a significantly stronger tidal field; consequently, all currently observed candidate binary clusters in the Milky Way are OCs. The most recent catalogues list several hundred candidate tidal capture binary clusters (Palma et al. 2025; Liu et al. 2025). Therefore, investigating the dynamical interactions of low-mass OCs during close encounters is important for understanding the formation and evolution of binary star clusters in the Milky Way. In fact, the dynamical processes involved in OCs may be fundamentally different from those of GCs. In Mastrobuono-Battisti et al. (2019), when clusters with masses of 10⁷ M_⊙ were replaced by lower-mass systems (10⁶ M_⊙) and the orbital parameters were left unchanged, the merging and mass exchange observed in the original models no longer occurred.

In this work, we investigated the long-term evolution of cluster encounters across a range of mass ratios and density ratios, under different orbital parameters. Across all models, significant dynamical interactions, such as merging and tidal disruption, occur exclusively in the parabolic cases. When two clusters follow parabolic orbits and experience substantial overlap at pericentre, merging or tidal disruption can occur. Cluster pairs with different mass ratios can follow entirely different evolutionary scenarios. We present the evolutionary processes of two models with different mass ratios, tracing their evolution from initial separation to merger, and analysed the physical mechanisms governing the merging process. In the equal-mass merger cases, when the clusters pass through pericentre, strong tidal interactions heat the outer particles of the clusters, raising their energies above the binding energy of the system and causing them to escape. As stars escape, the orbits of the interacting clusters evolve and gradually decay, leading to a progressive reduction in the separation between the clusters. During the interaction, a fraction of E₁₂ is transformed into the internal energy of the clusters, driving the system into a gravitationally bound state (E₁₂ < 0). At the same time, the orbit gradually circularises (with a decreasing eccentricity), and the spin of each cluster tends to synchronise with the orbital rotation. In unequal-mass merger cases, the merger process more closely resembles the tidal disruption of the companion cluster followed by its accretion, differing from the merging pathway observed in equal-mass cases. Nevertheless, the merger remnants all exhibit significant rotation. The inner regions rotate approximately as a rigid body, while the outer parts also show a rising rotational profile, which is markedly different from that of an isolated rotating cluster.

Clusters on parabolic orbits are more prone to experiencing strong interactions. This is partly because the relative velocity between the two clusters is lower, allowing more time for energy transfer, and partly because parabolic orbits have an initial orbital energy close to zero. As a result, interactions and energy exchange between the clusters can more readily drive the system into a gravitationally bound (elliptical) orbit. In the hyperbolic models, although close interactions strip some stars from the outskirts of the clusters, the main body remains stable, and neither complete disruption nor merging occurs. Although this result is consistent with previous simulation works (Devadas Rao et al. 1987; de Oliveira et al. 1998), hyperbolic orbits are more representative of the initial conditions of cluster close encounters in realistic environments, with E₁₂ > 0. For a pair of initially unbound clusters, the relative orbital energy must be positive, rather than zero as in the case of a parabolic orbit. Observational results likewise indicate that clusters undergoing close encounters are unlikely to become gravitationally bound pairs, consistent with the behaviour seen in the hyperbolic models (Piatti & Malhan 2022; Casado 2022; Zhu et al. 2025). Casado (2022) revisited all the remaining binary cluster candidates in the Galaxy that include at least one cluster older than 100 Myr using Gaia data and a careful revision of the literature. They find that the majority of the 120 pairs or groups of old systems that were revised were optical pairs or flyby encounters, with no convincing cases of true binarity.

Dutta Chowdhury et al. (2020) summarised three mechanisms that can cause an initially unbound cluster pair (E₁₂ > 0) to become gravitationally bound (E₁₂ < 0). The primary mechanism involves a close, impulsive encounter between two clusters, during which relative orbital energy is transferred into the internal energy of the clusters (Spitzer 1958). If this energy transfer is sufficiently significant, E₁₂ becomes negative, a process referred to as ‘tidal shock capture’. The second mechanism is known as ‘dissipative capture’ in which clusters experience different amounts of dynamical friction, allowing the two clusters to become gravitationally bound by chance. The final mechanism involves compressive tidal forces exerted on a cluster pair by one or more additional clusters, which can cause the pair to become gravitationally bound; this process is referred to as ‘three-body capture’. This mechanism is analogous to three-body scattering in stellar systems, where two initially unbound stars require a third star to carry away excess angular momentum in order for the pair to become gravitationally bound (Atallah et al. 2024). Therefore, studying the evolution of cluster groups may reveal interesting dynamical phenomena (Capuzzo-Dolcetta & Miocchi 2008).

Finally, we discuss the limitations of modelling cluster close encounters using two-body formulations, and outline aspects that require improvement in future work. First, when the relative velocity between clusters is derived using two-body formulae, it depends explicitly on the cluster masses. Because OCs are relatively low-mass systems, the resulting relative velocities are lower than the velocities expected for cluster encounters in the actual Galactic environment. Grinenko & Kovaleva (2025) report that the relative velocities of most OC encounters do not exceed 10–20 km s⁻¹, with a median encounter velocity of about 12.9 km s⁻¹. In addition, stellar evolution, binary evolution, and the presence of an external potential can also influence the dynamical evolution of cluster encounters. Incorporating these effects in future studies will therefore be essential for a comprehensive understanding of the formation and evolution of star clusters.

Fig. 10 Distributions of the binding energy (B) for the merger remnants in models M1ρ0.25 (p10; panel a) and M4ρ4 (p5; panel b), compared with those of the initial models OC1 and OC2. The solid black histograms show the merger remnants at the indicated times, while the dashed and dash-dotted curves represent the initial models.

Acknowledgements

This work has been supported by the National Natural Science Foundation of China (No. 12473029), Dali Expert Workstation of Rainer Spurzem, Yunnan Academician Workstation of Wang Jingxiu (202005AF150025). S.H.Z. acknowledges the support from the National Natural Science Foundation of China (grant Nos. 12173026, 12141302), the Innovation Program of Shanghai Municipal Education Commission (grant No. 2025GDZKZD04), the scientific research grants from the China Manned Space Project (grant Nos. CMS-CSST-2025-A06, CMS-CSST-2025-A07), the Program for Professor of Special Appointment (Eastern Scholar) at Shanghai Institutions of Higher Learning, and the Shuguang Program (23SG39) of Shanghai Education Development Foundation and Shanghai Municipal Education Commission.

References

All Tables

All Figures

	Fig. 1 Projection of model M1ρ0.25 (p10) onto the x–y plane at nine different evolutionary times. Black points represent OC2, located at the origin. Red points denote OC1, and grey squares indicate unbound particles.
In the text

	Fig. 2 Projection of model M4ρ4 (p5) onto the x–y plane at nine different evolutionary times. To clearly illustrate the bridge structure, all particles in panel h are displayed as black points.
In the text

	Fig. 3 Time evolution of the relative separation between the two clusters, r12 (top panel), the relative orbital energy, E12 (middle panel), and the orbital eccentricity, e (bottom panel), for model M1ρ0.25 (p10).
In the text

	Fig. 4 Same as Fig. 3 but for model M1ρ1 (h10).
In the text

	Fig. 5 Time evolution of the orbital angular momentum, Lorb, and the spin angular momenta of the two clusters, Lspin,1 and Lspin,2, for models M1ρ0.25 (p10; panel a) and M4ρ4 (p5; panel b).
In the text

	Fig. 6 Distributions of the spin angular velocity, ωspin (solid lines), and the orbital angular velocity, ωorb (dashed lines), as functions of the cylindrical radius (R⊥) for models M1ρ0.25 (p10; panels a–c) and M4ρ4 (p5; panels d–f) at different evolutionary times, as indicated in each panel.
In the text

	Fig. 7 Rotation curves of the merger remnants for models M1ρ0.25 (p10; panel a) and M4ρ4 (p5; panel b) at different evolutionary times, as indicated. The inset shows, in logarithmic scales, a representative rotation curve of an isolated rotating cluster taken from Kim et al. (2004), illustrating the rapid decline of the rotational velocity beyond Rh.
In the text

	Fig. 8 Equatorial-plane projections of the merger remnant for model M4ρ4 (p5) at different evolutionary times. Stars are colour-coded according to their original cluster membership.
In the text

	Fig. 9 Radial stellar mass density profiles of the merger remnants for models M1ρ0.25 (p10; panel a) and M4ρ4 (p5; panel b) at different evolutionary times. The dashed and dotted curves show the initial Plummer models of OC1 and OC2, respectively.
In the text

	Fig. 10 Distributions of the binding energy (B) for the merger remnants in models M1ρ0.25 (p10; panel a) and M4ρ4 (p5; panel b), compared with those of the initial models OC1 and OC2. The solid black histograms show the merger remnants at the indicated times, while the dashed and dash-dotted curves represent the initial models.
In the text