# Final Stages of N-body Star Cluster Encounters

## 1 Introduction

The study of star cluster pairs in the Magellanic Clouds has been explored extensively in the past years (Bhatia & Hatzidimitriou 1988; Bhatia & McGillivray 1988; Bica, Clariá & Dottori 1992; Rodrigues et al. 1994, hereafter RRSDB; Bica & Schmitt 1995; de Oliveira et al. 1998, hereafter Paper I; Bica et al. 1999). The evolution of physical pairs can provide fundamental insight into the past history of cluster formation in the Magellanic Clouds.

In the last decade, N-body simulations of stellar system encounters have been the main tool to investigate the dynamical processes that can occur in such interactions, like mergers and tidal disruption. Many of the first simulations worked with equal mass encounters (White 1978; Lin & Tremaine 1983; Barnes 1988). Different mass encounters have been carried out by Rao, Ramamani & Alladin (1987), Barnes & Hut (1986, 1989), RRSDB and Paper I. These studies indicate that tidal disruption and merger are two important processes in the dynamical evolution of a binary stellar system. However, there remain important issues not yet well tested, as the long term stages of such interactions, and the fraction and how stars are stripped from the cluster by the parent galaxy tidal field.

Some isolated clusters in the Magellanic Clouds present interesting structures with higher ellipticities than Galactic globular clusters (van den Bergh & Morbey 1984). Several explanations have been proposed like the presence of sub clumps merging with the main cluster which could produce the impression of high ellipticities in these clusters (Elson 1991). Another possibility is that old mergers in advanced evolutionary stages could show such ellipticities (Sugimoto & Makino 1989).

This work is a continuation of Paper I where we analyzed the morphology of selected cluster pairs in the Magellanic Clouds and have compared them to those obtained from numerical simulations of star cluster encounters. A preliminary discussion was given in de Oliveira et al.  (1999). In the present paper we study the long term behaviour of some of our numerical models (up to 1 Gyr), and compare the final stages of these simulations (isodensity maps, ellipticities, isophotal twisting) with the structure of isolated clusters.

In Section 2 we describe the method and the conditions employed in the present simulations. In Section 3 we describe the Magellanic Cloud cluster images and the procedure to derive the isodensity maps and ellipticity measures. The numerical results and discussions are presented in Section 4. In Section 5 we give the main conclusions.

## 2 The Method and the Initial Conditions

We performed the simulations using TREECODE (Herquist 1987) in the CRAY YMP-2E computer of Centro Nacional de Supercomputação of the Universidade Federal do Rio Grande do Sul (CESUP-UFRGS). For a complete description of the method see Paper I.

### 2.1 The Initial Conditions

In Paper I, cluster N-body models were generated for 16384, 4096 and 512 equal mass stars, corresponding to total masses of $`10^5\mathrm{M}_{},10^4\mathrm{M}_{}`$ and $`10^3\mathrm{M}_{}`$, respectively. The cutoff radii of the Plummer clusters are 20, 15 and 8 parsecs, respectively. These values are comparable to diameters found in Magellanic Clouds clusters (Bhatia et al. 1991). We generated two different concentration degrees for the 512 particle cluster, one with half-mass radius 1/2 of the cutoff radius and another with 1/4 (Table 1). The interaction models are characterized by a pericentre distance $`p`$ (distance of closest approach) and the initial relative velocity $`V_i`$. These initial velocities are comparable to the observed random velocities in the LMC disk (Freeman, Illingworth & Oemler 1983). The initial conditions (t=0) used in Paper I were: (i) the large cluster (in the case of equal mass clusters, the more concentrated one) is located at coordinates (X,Y)=(0,0); (ii) the small cluster lies at a distance $`r_0=5\times r_t`$ (where $`r_t`$ is the tidal radius). Beyond this distance the tidal effects are negligible. The initial relative velocity at this distance was obtained from the two-body formulae. The position and velocities of the particles during the encounters were computed in the centre-of-mass frame of the entire system. At various times, the essential data containing the positions and velocities of all the particles were stored for later analysis. The computation was stopped when disruption of the cluster occurred or, in the case of open orbits, when the relative separation of the clusters reached distance $`d>2r_0`$. We point out that the softening parameter adopted for the interaction is the smallest value of the two cluster models (column 7 of Table 1). Indeed, with reference to the AB cluster interaction (which has the largest $`ϵ`$ difference) we carried out a simulation where cluster B is allowed to evolve isolated adopting the cluster A softening value for $`\mathrm{\Delta }t`$200 Myr ($``$ 2 relaxion times). Analysis of the time evolution of cluster B structure shows that it is essentially the same for $`ϵ`$ = 0.30 as that obtained for the cluster with $`ϵ`$ = 0.47. Consequently in the AB cluster interaction significant structural changes are not expected by the adoption of the smaller $`ϵ`$ value.

In the present Paper, we proceeded the simulation from this point on for some models with closed orbit where disruption occurred (see Table 2). We used the same procedures as in Paper I. These simulations were ran up to 990 Myr for the models presented here.

We illustrate in Fig. 1 the time evolution of an elliptic orbit encounter involving two clusters with 16384 and 4096 particles respectively (model E9AB10, see Table 2). In this Fig. we can see the small cluster being deformed by the massive one with the formation of a bridge ($``$50 Myr), subsequently complete disruption occurs ($``$135 Myr); at the end of this simulation the clusters coalesce into a single one, with the stars of the disrupted cluster forming a halo around the final cluster, and some of them being ejected.

## 3 Isodensities for LMC Elliptical Clusters

In the final stages of our simulations the pairs coalesced into a single cluster with a distinct structure as compared to the original ones. We noticed that in viewing angles close to the original orbital plane of encounter, the final single cluster presented an ellipticity larger than those of the initial clusters (our models have a spherical symmetry). In this paper we have selected some LMC clusters with morphologies resembling those of the present models (Sect. 4.4). We checked many clusters with important ellipticity, as indicated in previous studies (Geisler & Hodge 1980, Zepka & Dottori 1987, Kontizas et al. 1989, 1990). In particular the selected clusters present increasing ellipticity outwards. In Table 3 we show these clusters with their V magnitudes, SWB types derived from interpreted UBV photometry (Bica et al. 1996, Girardi & Bertelli 1998) and corresponding ages. Also there are available ages from colour-magnitude diagrams: (i) NGC1783 has 700- 1100 Myr depending on adopted distance modulus (Mould et al. 1989). Based on the same data, using the $`\mathrm{\Delta }`$V turnoff/clump method, Geisler et al. (1997) obtained 1300 Myr; (ii) NGC1831 has an age of 400 Myr according to Hodge (1984). CCD data provided 500-700 Myr (Mateo 1988) and 350-550 Myr (Vallenari et al. 1992), depending on adopted models; (iii) NGC2156 has age $``$60 Myr according to Hodge (1983); (iv) NGC1978 has an age $``$2000 Myr (Olszwewski 1984, Bomans et al. 1995). Geisler et al. (1997) derived 2000 Myr by means of the $`\mathrm{\Delta }V`$ turnoff/clump method. There is good agreement of ages derived from the CMD studies and the integrated colours (Table 3). NGC1783 and NGC1831 have ages comparable to the long term stages of our models, while NGC2156 falls somewhat short, and NGC1978 has an age $``$ a factor 2 larger.

The images of this selection were obtained from the Digitized Sky Survey (DSS). The plates are from the SERC Southern Sky Survey and include IIIa-J long (3600s), V band medium (1200s) and V band short (300s) exposures. The PDS pixel values correspond to photographic density measures from the original plates, and are not calibrated.

The digitized images, similary to those generated by the models, were treated with the IRAF<sup>1</sup><sup>1</sup>1IRAF is distributed by the National Optical Astronomy Observatories, which is operated by the Association of Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation, U.S.A. package at the Instituto de Física — UFRGS, applying a 2-d Gaussian filter to smoothen out individual stars, and to create isodensity maps. The ellipticity measurement was made with the IRAF package task Ellipse which fits elliptical isophotes to data images with the use of a Fourier expansion:

$$y=\underset{n=0}{\overset{4}{}}A{}_{n}{}^{}sin(nE)+B_ncos(nE)$$
(1)   
where E is the position angle. The amplitudes ($`A_3,B_3,A_4,B_4`$), divided by the semi-major axis length and local intensity gradient, measure the isophotal deviations from perfect ellipse. Note that the IRAF routine requires that one indicates a first guess for the object major axis position angle. Subsequently the routine iteractively fits the best solution for each isophote. We stored the coefficient $`B_4`$ which tells whether the ellipse is disk (a positive $`B_4`$ parameter) or box-shaped (a negative $`B_4`$ parameter), Bender et al. 1989. We also stored the semi-major axis position angle variation. It should be noted that in general ellipticity does not present simple behaviour. Ellipticity varies with radius (Zepka & Dottori 1987, Kontizas et el. 1989,1990) and so it is difficult to assign one to a cluster. We measured elliptical isophotes starting at $`R_h`$ and stopping at the last possible non diverging ellipse. In the models, this last ellipse occurred around the radius containing 90% of the cluster mass.

## 4 Discussion

We considered models in which the less massive cluster (perturbed) is allowed to move in elliptic orbits around the more massive one (perturber). The perturbed cluster is assumed to move in an orbit of eccentricity $`e=0.6,0.7,0.9`$. The values of the pericentre and the eccentricity are meaningful only if the clusters are assumed to be point masses moving in Keplerian orbits. However, soft potential orbits are not conical and our use of these definitions is meaningful only to a good approximation, but not strictly (White 1978).

The collision parameters of the simulations are given in Table 2. The model designation (column 1 of Table 2) contains information on the encounter conditions. E refers to an elliptic orbit encounter; the number following the letter E is related to the orbital eccentricity (column 5 of Table 2); A, B, C and D refer to the model type (Table 1); the last number is the orbit pericentre (column 6 of Table 2).

From the analysis of the time evolution of each model as illustrated in Fig. 1 for E9AB10, it is possible to observe the trend to form bridges in the beginning of the encounters (see also Paper I). As time goes by the smaller cluster is completely disrupted by the massive one, resulting in a single cluster at the end of the simulation.

### 4.1 The Structure of the Final Stages

In Fig. 2, we illustrate the same encounter as in Fig. 1, but in a different plane projection (parallel to the orbital plane). We can observe the disruption of the smaller cluster, with its stars occupying orbits around the massive cluster. As these stars remain preferentially orbiting in the same plane of the original encounter, they introduce an elliptical shape to the final single cluster, when seen in a favourable plane projection.

In Fig. 3 we show isopleths, ellipticity, position angle and $`4^{th}`$ cosine coefficient (hereafter $`B_4`$ ) for the model E9AB10 in a XZ projection. These measures were made for four different evolution times, the first panel being the initial massive cluster A before interacting (t=0 Myr). We can observe a radial variation of ellipticity during the interaction, which has also been observed for all other models. It also can be seen a general trend of decreasing ellipticity toward the inner parts of the models. Zepka & Dottori (1987) and Kontizas et al. (1989) have observed internal variations of ellipticity in a number of LMC clusters, finding a general tendency of increasing ellipticity towards the inner parts of these clusters. However, they reported a few exceptions ($`5\%`$ of the data) which in turn have counterpart in our models, and could be the result of an interaction. Indeed there are about 300 cluster pairs among a total population of 7847 extended objects in the SMC, inter-Cloud region and LMC (Bica & Schmitt 1995 and Bica et al. 1998). At least 50% of the pairs are expected to be interacting (Bhatia & Hatzidimitriou 1988) and we need a favourable projection (nearly edge-on orbits) to observe elliptical structure like those of the models.

In Fig. 3 there is little time variation of ellipticity during the lifetime of our model. This variation occurs mainly in the outer parts of the cluster, probably by the continuous escaping stars of the disrupted cluster.

There occur small position angle variations (Fig. 3), which indicate that the elliptical shape of the cluster stands out in the plane of the interacting encounter. The $`B_4`$ coefficient clearly has a positive peak, which is a sign of the presence of a disc component structure in the resulting model. This suggests that the elliptical shape in our final clusters is due to a rotation disk formed mainly by stars of the disrupted cluster.

In order to check whether mass or size could affect the behaviour of our models, we did the same analyses for models with different mass and concentration degree, but with same orbital encounter. In Fig. 4 we present the same measures of Fig. 3 for other two models at the final stage and we see similar ellipticity variations in radius. In order to compare absolute ellipticity variations between models, and also between different time stages in the same model, it is required to define ellipticity at a common radial distance. As we have the main shape changes due to variations in the outermost parts of these clusters we decided to measure the ellipticity as the mean between 3$`R_h`$ and the outermost ellipse measured.

In Table 4 we give the mean ellipticity between 3$`R_h`$ and the outermost ellipse $`e_{out}`$ for all the models in different time stages. We conclude that in model E9AB10 $`e_{out}`$ shows a small tendency of increasing with time; differences among the models are found mainly in the shape of e versus radial distance curve.

When comparing ellipticity between different models at the same appropriate time (Table 4), we observe a small reduction of $`e_{out}`$ from a more massive model (E9AB10) to a less massive model (E9BC10). This relation (ellipticity) versus (mass) has been suggested for the LMC clusters (van den Bergh & Morbey 1984) where high mass clusters have higher angular momentum or have more difficulty in shedding it than do low mass clusters. However, when comparing a more massive model (E9AB10) with a less massive one (E9BD10) — where we have in the last a different concentration degree for the disrupted cluster — we observe no significant variation in $`e_{out}`$. Thus, the initial concentration degree affects the final shape of the resulting cluster. This suggests that a more concentraded cluster has a deeper potential well, making more difficult stripping its stars by the more massive cluster. On the other hand, stripping stars in the less concentrated cluster is more efficient, resulting in a more pronounced halo expansion. This is supported by the results in Table 5. This halo expansion when seen in an edge-on view seems to be an important factor to stablish the shape of the final merger.

We also compared the possible ellipticity variation between models with equal mass and concentration degree, but different initial orbital parameters. We observed no significant differences in $`e_{out}`$ between the models.

### 4.2 Rotational Velocity Field

In Fig 5 we show the line-of-sight rotational velocity field ($`V_y`$) for the model E9AB10 along the major axis at 990 Myr seen in a XZ projection. The merger consists approximately of a rigid rotation core and an outer halo ( $`r20`$ pc) with a Keplerian fall. The rotational velocity along the major axis has a peak value of about 1 $`Kms^1`$ at a radius of $`2030pc`$. It is interesting to see that the rigid rotation extends exactly to the radius of the initial massive cluster before the encounter. This indicates that the main contribution to the velocity field for $`r>20`$ pc is due to stars of the small cluster. This becomes evident when we plot the two star cluster member sets for a merged final model (Fig. 6) in a XZ projection. In this Fig. we plot the model EAB10 at T=990 Myr and a reference box (side of 40 pc) centred at the more massive cluster. It is clear that the main contribution to the final composite system for $`r>20`$ pc is due to stars of the small cluster. Velocity distributions obtained with high resolution spectra with large telescopes might reveal such rotation curves, while in turn could be a signature of evolved stages of cluster merging.

### 4.3  Mass Loss

In Table 5 we give the radii containing 10%, 50% and 90% of the mass for the massive cluster in each interaction (model A or B) for $`t=0`$ and for the final cluster stage. We can observe that $`R_{50}`$ varies little for the initial and final cluster stage when comparing models with same initial mass, suggesting that the massive one has its original stellar content little affected by the collision. However model E9AB10 shows a contraction of the core, showing that the initial mass may affects the cluster final structure. The radius containing 90% of the mass is larger for all our models, which shows that the final cluster outer halo is formed mainly by stars from the disrupted cluster. The latter result is in agreement with those presented in Paper I for the radial mass distribution. There we showed that the maximum expansion for the disrupted cluster in E models always occurs in the plane of the encounter and is a minimum in the z direction. When we compare the projected density distribution between the models before interacting and at the end of the simulation we see that they are similar within the region $`r<10`$ pc (Fig 7). However, the final cluster extends with a dependence of $`\rho r^3`$.

Together with changes in the internal structure of our final single cluster, the present simulations suggests that a fraction of cluster stars may be ejected to the field due to the encounter. In order to quantify this fraction, we calculated the total energy per particle, identifying those that had positive total energy which are expected to escape the system.

In Table 6 we have for each model the final fraction of stars that escape the cluster (in units of solar mass). The mass fraction that is ejected to the field is very small, considering the ensemble of the pair ($``$ 3% ), so the contribution to the field stars by a cluster pair encounter is not very significant. It should be noted that when we study the merger of the cluster pair we are not taking into account the tidal field of the parent galaxy. This should truncate the outer halo of the final cluster, which would contribute with an increasing fraction of stars ejected to the field. Adopting a tidal radius $`r_{tidal}65`$pc (distance moduli of 18.5 for LMC, Westerlund 1990) for observed cluster’s halo (Elson et al., 1987), we can estimate how much mass is beyond this limit. For the model E9AB10, we estimate $`M50\%`$ of the mass of the disrupted cluster (4.5% of the total mass) is beyond this limit. So, it can be concluded that although the total mass loss for the pair is not very significant, the assumption of a maximum radius for the final cluster, implies that half of the disrupted cluster stars feed the field.

### 4.4 Comparison of Simulations with Magellanic Cloud Clusters

Isopleth maps of projected planes at a given time $`t`$ of a suitable model can be compared to the observed isodensity maps of selected Magellanic Cloud clusters to infer on their dynamics. In Paper I we have used this method to search for evidences of interacting pairs. In the present study, we use model isopleths together with other measured parameters in order to compare our simulation final stages with possible observational evolved products of cluster pair interaction. Some examples are given in Figs. 8910 and 11, where we show some isolated LMC clusters isophotes, and radial dependence of ellipticity, position angle and $`B_4`$ coefficient.

The examples shown here have an ellipticity radial variation with a trend to increase towards the outer parts. NGC 1783 has an ellipticity curve compatible with our model E9BD10 when seen in a XZ projection at 990 Myr (Fig. 4, right). The cluster age (Sect. 3) is fully compatible with the model evolutionary time. This scenario can be explained by the result of a cluster pair encounter with the subsequent merger of the pair in the early history of NGC1783.

NGC1831 has many resemblances with the models. The ellipticity (Fig. 9) increases outwards, but in the external isophote the value decreases again. This is observed in the model E9BC10 at 990 Myr (Fig. 4, left). The cluster presents a considerable isophotal twisting like in the model.

As pointed out in Sect. 3, the age of NGC2156 falls short of the long term stages of the models. But notice that asymmetrical early stages (Fig. 2) at $``$ 25 Myr , if seen favourably in a given direction might present ellipticity variations. The negative B4 coefficient behaviour suggests that if the cluster is indeed a merger, disruption has not yet occurred to create a disc shape.

NGC1978 is an interesting intermediate age cluster due to the pronounced ellipticity (Fig. 11), which is already presented in the innermost available isophote. No isophotal twisting occurs in this case. Recently Kravtsov (1999) searched for special variations in the colour-magnitude diagram of NGC 1978. They found evidence of some variations which they attributed to a possible metallicity spread. Bomans et al. (1995) did not find any evidence of age variation within the cluster. A merger scenario does not require differences in age and/or metallicity. Indeed two clusters born in the same association are expected to have a close orbit encounter, which is the more favourable case for interactions (Paper I).

The $`B_4`$ coefficient has different behaviours in the 4 clusters, and in some cases has varying values in a given cluster. We might be witnessing clusters with either discoidal (positive values) or boxy (negative values) shapes. The average value of $`B_4`$ for NGC1783 suggests an overall boxy shape. NGC2156 is also boxy, but with a systematic trend to become disc-shaped in the outer parts. NGC1831 shows a predominant positive $`B_4`$ suggesting a disc-shape. Finally, NGC1978 on the average is positive thus indicating a predominantly disc-shape. However the cluster has two distinct $`B_4`$ value regions, 20$`r`$50 is definitively disc while r $`>`$ 50 should be classified as boxy. Behaviours like this have counterparts in elliptical galaxies (see e.g. the sample of Goudfrooij et al. 1994). The $`B_4`$ coefficient is a promising tool to explore the possibility of merger in galaxies or star clusters.

## 5 Conclusions

We used N-body simulations to study the final morphology and structure of star cluster pair encounters. We also compared these morphologies with those of elliptical isolated LMC clusters. The main conclusions may be summarized as follows:

1. Close orbit encounters of cluster pairs can lead to a single final cluster with a distinct structure of the original ones. When seen in a favourable plane projection, models show an elliptical shape comparable to those of some observed isolated LMC clusters. This suggests that tidal encounters could be a mechanism to explain the ellipticity of several clusters in the Magellanic Clouds.
2. Evolved stages appear to be stable after $``$ 200 Myr, suggesting that the resulting ellipticity is not transient. The simulations indicate that the initial concentration degree and mass affect the final shape of the resulting cluster. In the final stages all models present halo expansion, and some present core contraction.
3. The models show mass loss for the composite system.The fraction of stars ejected to the field by the encounter, is not so significant ($`3\%`$ with respect to the sum of both clusters). However it can represent up to $`50\%`$ of the stars of the disrupted cluster, if we assume a tidal radius of $`r_{tidal}65`$ pc for the final merged cluster.
4. The velocity distribution of some final stage models presents characteristic velocity patterns in favourable plane projections, which could be used as an observational constraint to test the present scenarios.

Finally, we call attention that the merging of spheroidal systems can produce disc-shaped products. Thus not only boxy-shaped systems might be related to mergers.

ACKNOWLEDGMENTS

We thank Dr. Hernquist for allowing us to use TREECODE and the CESUP-UFRGS for alloted time in the CRAY YMP-2E computer. We are also grateful to an anonymous referee for interesting remarks. We acknowledge support from the Brazilian Institutions CNPq, CAPES and FINEP. The images in this study are based on photographic data obtained using the UK Schmidt Telescope, which was operated by the Royal Observatory Edinburgh, with funding from the UK Science and Engineering Research Council, until 1988 June, and thereafter by the Anglo-Australian Observatory. Original plate material is copyright by the Royal Observatory Edinburgh and the Anglo-Australian Observatory. The plates were processed into the present compressed digital form with their permission. The Digitized Sky Survey was produced at the Space Telescope Science Institute under US Government grant NAG W-2166