THE INFLUENCE OF ORBITAL ECCENTRICITY ON TIDAL RADII OF STAR CLUSTERS

Theoretical calculations of the radius of a star cluster, or tidal radius, usually assume that the gravitational field of the host galaxy regulates cluster size. In most previous treatments, for simplicity it is further assumed that the gravitational field in which the cluster orbits is constant, i.e., the cluster has a circular orbit in a spherically symmetric galactic potential (e.g., von Hoerner 1957; King 1962; Innanen et al. 1983; Jordán et al. 2005; Binney & Tremaine 2008; Bertin & Varri 2008). The tidal radius is then assumed to be equal to the Jacobi radius (r_J), the distance beyond which the acceleration a star feels toward the galaxy center is greater than the acceleration it feels toward the cluster center, and the star is able to escape. First-order tidal theory determines the tidal radius (von Hoerner 1957) via

$$\begin{equation} r_t=R_{{\rm gc}}\left(\frac{M_{{\rm cl}}}{2M_g}\right)^{1/3}, \end{equation} \tag{ 1 }$$

where R_gc is the galactocentric distance of the cluster, M_cl is cluster mass, and M_g is the mass of the galaxy (assumed in early studies to be a point mass).

For a cluster with a non-circular orbit, the fact that the tidal field is no longer static makes an analytic approach very difficult (see Renaud et al. 2011 for another approach). Historically, it has been assumed that for a globular cluster on an eccentric orbit, its tidal radius is imposed at perigalacticon (R_p) where the tidal field of the host galaxy is the strongest. This assumption was initially suggested by von Hoerner (1957) and later by King (1962), and follows from the fact that the internal relaxation time ($t_{r_h}$) of the cluster is greater than its orbital period for almost all observed globular clusters. Therefore, after stars outside the tidal radius at perigalacticon escape, the cluster would not be able to relax and expand before it returns to perigalacticon. Thus in Equation (1), R_gc is usually replaced with R_p to calculate the tidal radius of a cluster with an eccentric orbit (e.g., Innanen et al. 1983; Fall & Zhang 2001; Read et al. 2006; Webb et al. 2012).

However, recent studies are showing with increasingly strong evidence that the actual sizes of observed clusters are not imposed at perigalacticon in this simple way. The actual size of an observed cluster is known as its limiting radius r_L, which marks the point where the cluster density falls to zero (Binney & Tremaine 2008). Using the solved orbits of 15 Galactic globular clusters, Odenkirchen et al. (1997) demonstrated that cluster limiting radii are not solely dependent on perigalactic distance. Brosche et al. (1999) suggested that some sort of orbit-averaged tidal radius is more appropriate when predicting limiting radii. Even with the orbits of an additional 29 Milky Way globular clusters currently known (Dinescu et al. 1999; Casetti-Dinescu et al. 2007), there is still no clear relationship between limiting radii and perigalactic distance, and the conclusions of Odenkirchen et al. (1997) still hold.

With the Galactic potential as given by Johnston et al. (1995; the same Galactic potential Casetti-Dinescu et al. used to solve cluster orbits), we calculated the theoretical tidal radius r_t at perigalacticon (that is, the Jacobi radius at R = R_p) of each of the 44 Galactic clusters with known orbits. The r_t values were calculated with the formalism of Bertin & Varri (2008; outlined in Section 3). The main uncertainty in the theoretical tidal radius is due to the uncertainty in R_p quoted in Dinescu et al. (1999) and Casetti-Dinescu et al. (2007). To compare theory and observations for individual clusters, we take cluster limiting radii as determined from direct King (1966) model fits r_k to the observed cluster profiles as listed in Harris (1996) (2010 edition). We then calculate the ratio of the difference (r_k − r_t) between observed and theoretically predicted values to their average ((r_k + r_t)/2). An uncertainty of 10% was assigned to values of r_k. If theory and observations are in agreement, the ratio will be approximately zero. Clusters which have a ratio greater than zero will be ones which overfill their predicted tidal radius, while clusters with ratios less than zero will be tidally underfilling. The comparison between theory and observations is illustrated in Figure 1.

Zoom In Zoom Out Reset image size

While it is not expected that all clusters are tidally filling (Gieles et al. 2010), the fact that the majority of clusters appear to be tidally overfilling is a strong signal that the simple assumptions built into the model need investigation. The known presence of tidal tails around observed (e.g., Odenkirchen et al. 2001) and simulated (e.g., Montuori et al. 2007; Küpper et al. 2012; Lane et al. 2012) globular clusters is not sufficient to explain the cases of apparent overfilling. The observed r_k values are determined from King-model profile fits that are heavily dominated by the inner populations of stars, out to a few half-light radii. In almost all cases, the extremely low densities of stars in the tails, at or beyond the nominal tidal radius, exert little leverage on these fits. Furthermore, King (1966) models are known to underestimate cluster sizes in general as they require a sharp tidal cutoff which is not observed in all clusters (McLaughlin & van der Marel 2005). Only in the most extreme cases (e.g., Pal 5; Odenkirchen et al. 2003) will large and extended tidal tails influence model fits to the cluster surface brightness profile.

Theoretical calculations may also underestimate r_k because the assumption that the tidal radius is imposed at perigalacticon implies that the shape of the tidal field and the cluster orbit are not important. Hence factors such as tidal heating and disk shocking due to a varying symmetric galactic potential are not taken into consideration. However, despite these potential inaccuracies we still expect to see some sort of correlation between r_k and perigalactic distance, if tidal radii are imposed at perigalacticon. The Milky Way cluster data therefore suggest that something is wrong with basic tidal theory.

Recent N-body simulations by Küpper et al. (2010) find that their fitted King (1962) radius was better represented by the time-averaged mean tidal radius of the cluster and not the perigalactic tidal radius. In a later study on the structure of tidal tails, Küpper et al. (2012) found that while stars outside the tidal radius as calculated at perigalacticon will likely become unbound at perigalacticon, some are able to be recaptured by the cluster as it moves away from perigalacticon and the instantaneous tidal radius of the cluster increases. That is, the limiting radius of a cluster will be greater than the tidal radius calculated at perigalacticon. This discrepancy is expected to be amplified for clusters on very eccentric orbits, as they make quick perigalactic passes and spend the majority of their time near apogalacticon (R_a). In fact, N-body simulations by Madrid et al. (2012) suggest that the half-mass radius of a globular cluster is more likely imposed at R_a than R_p.

The purpose of this study is to explore more thoroughly the influence of orbital eccentricity on cluster size. Model N-body clusters with different initial half-mass radii are evolved from zero age to 10 Gyr over a range of orbital eccentricities in the disk of a Milky-Way-like potential. The models and their initial conditions are described in Section 2. In Section 3 we focus on the influence of orbital eccentricity on the mass (M), half-mass radius (r_m), limiting radius, and tidal radius of each model cluster over time. In Section 4 we explore the influence of initial cluster half-mass radius on our results. In Section 5 we discuss the results of all our N-body models and the influence of orbital eccentricity not only on cluster radii, but on individual stars within the cluster as well. Based on our findings, in Section 6 we suggest a correction factor that can be applied to the perigalactic tidal radius of a cluster to better match its observed limiting radius. This correction factor is then applied to the Galactic globular clusters shown in Figure 1. Finally, in Section 7 we summarize our conclusions.

We use the NBODY6 direct N-body code (Aarseth 2003) to evolve model star clusters from zero age to a Hubble time, over a range of both initial cluster half-mass radii and orbital eccentricity. All models in this study begin with 48,000 single stars and 2000 binaries.

As long as we use one particle to represent one real star, clusters of this size correspond physically to either a very massive open cluster or a low-mass globular cluster in the Milky Way. Ideally, we would like to follow more "average" globular clusters, which will typically have 200,000 stars or more. However, to this point in the history of the subject, direct N-body integrations with N-values that large have only been carried out for special, specific purposes (for example, see Hurley & Shara 2012 for a 200,000-star simulation in the tidal field of a point-mass Galactic potential and a circular orbit; Zonoozi et al. 2011 for N ⩽ 100, 000-star simulations directed at modeling Pal 14, again on a circular orbit; or Heggie & Giersz 2009 for a 105,000-star simulation of NGC 6397 over 1 Gyr). Encouragingly, however, these high-N models generally confirm the trends obtained from earlier, small-N simulations (see Hurley et al. for additional discussion).

The purpose of our suite of 50,000-star models is instead to survey a wide parameter space of initial cluster half-mass radii and orbital type. Our chosen ranges (see below) match real Milky Way star clusters moving in a realistic time-varying potential. Eventually, with advances in computational capabilities this type of survey work can be extended to higher N values (up to 10⁶) that will cover almost the entire known range of globular clusters. However, as will be seen below, the models summarized here already prove to be highly informative in revealing important physical effects of orbital eccentricity on the internal dynamical evolution in a direct way that does not rely on analytical approximations.

Since we are only concerned with the influence of orbital eccentricity on clusters of different initial half-mass radii, our choice of initial parameters such as cluster metallicity, the stellar initial mass function (IMF), and binary fraction are of little consequence as long as they remain consistent between models. However, we note that binary fractions of a few percent are typical for globular clusters (e.g., Davis et al. 2008).

The masses of single stars are drawn from a Kroupa et al. (1993) IMF between 0.1 and 30 M_☉. For binary stars, the masses of two randomly selected single stars are combined to equal the total mass of the binary, with the primary and secondary masses determined by a mass ratio randomly drawn from a uniform distribution. The initial total mass of each model is 3 × 10⁴ M_☉. The initial period of each binary is drawn from the distribution of Duquennoy & Mayor (1991) and their orbital eccentricities are assumed to follow a thermal distribution (Heggie 1975). All stars were given a metallicity of Z = 0.001. The initial positions and velocities of the stars are generated based on a Plummer density profile (Plummer 1911; Aarseth et al. 1974). We note that the Plummer model extends to an infinite radius so we impose a cutoff at ∼10 r_m to avoid rare cases of large cluster-centric distances. A description of the algorithms for stellar and binary evolution can be found in Hurley (2008a, 2008b).

The model clusters orbit in a three-dimensional Galactic potential, which consists of a point-mass bulge, a Miyamoto & Nagai (1975) disk (with a = 4.5 kpc and b = 0.5 kpc), and a logarithmic halo potential. The combined mass profiles of all three potentials result in a circular velocity of 220 km s⁻¹ at a galactocentric distance of 8.5 kpc. The bulge and disk have masses of 1.5 × 10¹⁰ and 5 × 10¹⁰ M_☉, respectively (Xue et al. 2008). Aarseth (2003) and Praagman et al. (2010) describe the incorporation of the Galactic potential into NBODY6. All models were set to orbit in the plane of the disk such that a cluster on a circular orbit experiences a static tidal field, and will not be subject to factors such as tidal heating or disk shocking.

For the purposes of our study, the only parameters which are important and change from model to model are initial cluster half-mass radius, initial cluster position, and initial cluster velocity which determine the shape of the orbit. Our first models were for clusters with orbital eccentricities of 0.25, 0.5, 0.75, and 0.9, where eccentricity is defined as e = (R_a − R_p)/(R_a + R_p). All of these have the same perigalactic distances of 6 kpc and initial half-mass radii r_{m, i} of 6 pc. These clusters are located at perigalacticon at time zero. For comparison purposes a model was simulated with a circular orbit at perigalacticon (6 kpc) and four more with circular orbits at the apogalactic distance of each eccentric cluster (10 kpc, 18 kpc, 43 kpc, and 104 kpc). These simulations allow us to directly compare the properties of a globular cluster on an eccentric orbit to clusters on circular orbits at both perigalacticon and apogalacticon.

Our set of models also included re-simulations of the cluster with an orbital eccentricity of 0.5, and the corresponding e = 0 perigalactic and apogalactic simulations, but with initial half-mass radii of 4 pc, 2 pc, 1 pc, and 0.5 pc. This range allows us to study the influence of initial cluster half-mass radius. Both sets of models are summarized in Table 1. Model names are based on orbital eccentricity (e.g., e05), circular radius at apogalacticon (e.g., r18), and initial half-mass radius (rm6). Hence a model cluster with a perigalactic distance of 6 kpc, apogalactic distance of 18 kpc (orbital eccentricity of 0.5), and an initial half-mass radius of 6 pc would be labeled e05r18rm6.

The first portion of this study will focus solely on models with an initial r_m equal to 6 pc, with the only difference between each model being their Galactic orbits. However, we first need to determine whether any given star is bound to the cluster. In a cluster-centric coordinate system, we define the x-axis as pointing away from the galactic center, the y-axis pointing in the direction of motion of the cluster, and the z-axis pointing perpendicular to the orbital plane. In this coordinate system the energy of an individual star can be written as

$$\begin{equation} E = \frac{1}{2} (\dot{x}^2 + \dot{y}^2 + \dot{z}^2) - \sum _{i=1}^{N-1} \frac{ G m_i}{\Vert r-r_i \Vert } - \frac{1}{2} \Omega ^2 (z^2 - \upsilon x^2), \end{equation} \tag{ 2 }$$

where the second term is the potential energy due to the remaining N − 1 stars in the simulation, each with mass m_i and located a distance r_i from the star. The third term is the tidal potential with Ω² equal to the orbital frequency of the cluster. Here, $\upsilon$ is a dimensionless positive coefficient defined below in Equation (6). The tidal potential, taken from Bertin & Varri (2008), results in a stretching of the cluster in the x-direction, no change in the y-direction, and a compression in the z-direction. If the resultant energy is less than zero the star is bound, otherwise it is considered unbound.

We considered additional criteria for determining whether a star is bound or unbound in addition to the energy calculation. Other studies have invoked a distance cutoff such that the stars' cluster-centric distance must be greater than the cluster perigalactic or instantaneous tidal radius for it to be unbound (e.g., Takahashi & Baumgardt 2012). It has also been suggested that a star's velocity plays a role in whether or not it can be considered unbound (e.g., Küpper et al. 2010, 2012). However, we found that these additional criteria did not change any of the results found in Section 3 as they only affected a small percentage of simulated stars. Therefore, we only require that a star's energy as given by Equation (2) to be greater than zero for the star to be considered unbound.

Figure 2 shows a model cluster at a representative time step. The tidal tails formed by escaping stars are clearly visible in our simulations. The densely populated spherical collection of stars marked in red are those that satisfy our boundedness criterion and are considered cluster members. The unbound stars that appear close to the center of the cluster are foreground stars with a large z-coordinate and are simply projected onto the cluster in the x–y plane. These tails have no effect on our determination of theoretical tidal radii or observed limiting radii as we only consider stars that are bound to the cluster.

Zoom In Zoom Out Reset image size

3.1. Mass

We first wish to study how orbital eccentricity influences the total mass, or more specifically the mass loss, of a star cluster over time, which then plays an important role in determining tidal radii (r_t∝M^1/3). The total mass of stars bound to the cluster in each model is illustrated in Figure 3. In general, mass loss is due to stellar evolution, evaporation due to two-body interactions, and tidal stripping. For clusters with circular orbits, in all cases the apogalactic cluster loses mass at a lower rate than the perigalactic case as a result of less tidal stripping.

Zoom In Zoom Out Reset image size

The e = 0.25 case loses mass at almost the same rate as if it had a circular orbit at its perigalactic distance, but the final mass is still notably larger than the perigalactic case. Since eccentric clusters spend the majority of their time away from perigalacticon, they too will be subject to less tidal stripping than an ideal cluster that spends all its time at perigalacticon. As eccentricity increases, the mass-loss profile shifts farther away from the perigalactic case and closer to a cluster with a circular orbit at apogalacticon.

At higher eccentricities the mass no longer smoothly decreases, in contrast to the circular orbit cases. Instead, periodic fluctuations are present. The minima of these fluctuations correspond to perigalactic passes, where the rapid increase in tidal field strength results in episodes of significant mass loss. These fluctuations suggest that a greater change in tidal field strength between apogalacticon and perigalacticon results in stars gaining more energy at or near perigalacticon.

Especially interesting in the lower right panel of Figure 3 is the fact that once a cluster undergoes significant mass loss during a perigalactic pass, the cluster starts to regain mass before resuming its mean mass-loss rate. In these intervals just after perigalacticon, the cluster is recapturing some of the stars which were previously unbound. Stated differently, many of the stars that were formally unbound at R_p drift away slowly enough that they are recaptured as the cluster moves back outward and its instantaneous tidal radius expands again. Furthermore, these fluctuations suggest that while a perigalactic pass does have a strong effect on an eccentric cluster, the cluster cannot be treated as if it had a circular orbit at R_p. These results are in agreement with the findings of Küpper et al. (2012) discussed earlier.

3.2. Half-mass Radius

We next consider how orbital eccentricity can influence the half-mass radius r_m of bound stars within a globular cluster. It should be noted that r_m is not the same as the half-light radius r_h (also known as the effective radius), which is a directly observable parameter. In our simulations, the half-mass radius is always slightly larger than the half-light radius.

The results of our simulations are illustrated in Figure 4. The initial increase in r_m during the first ∼2000 Myr in all cases is driven by two-body relaxation and stellar evolution mass loss. However, once the cluster is relaxed, tidal stripping becomes the dominant dynamical process.

Zoom In Zoom Out Reset image size

The r_m profiles of the apogalactic cases in the lower panels do not begin to decrease after 2000 Myr, but instead continue to increase up to 10 Gyr. As discussed in the next section this trend is due to the fact that these clusters are barely tidally filling, so can still expand and not be subject to tidal stripping.

Similar to the results of Section 3.1, for low eccentricities the r_m profile of the eccentric model cluster is comparable to the circular orbit case at perigalacticon. Increasing eccentricity brings the r_m profile closer to the apogalactic case on the average, but increasing eccentricity again results in sharp fluctuations in the r_m profile which correspond to perigalactic passes. The trends shown in Figure 4 reveal what is perhaps the most striking difference between clusters on static circular orbits (the classically assumed case) and ones on more realistic eccentric orbits. If cluster limiting radii are imposed at perigalacticon, we would expect the minima of the eccentric r_m profile to be equal to the r_m profile of a cluster orbiting at R_p. While a high-e cluster may briefly expand after a perigalactic pass, the next perigalactic pass would restore r_m to a size equal to the perigalactic case. But not even at perigalacticon does the r_m of the eccentric cluster equal the perigalactic case. Instead, the time-averaged r_m is reflective of a tidal field weaker than the field at perigalacticon, in agreement with Küpper et al. (2010).

After a perigalactic pass, Figure 4 illustrates again that the cluster is able to increase in size. Especially apparent in the lower right panel of Figure 4 is the fact that the cluster is able to increase to a size greater than its mean r_m. Inspection of our N-body models shows that this increase in size is due to a combination of recapturing some of the previously unbound stars (Küpper et al. 2012) and the stars in the inner region of the cluster gaining enough energy to move outward and repopulate the halo of the cluster. These statements are discussed in further detail in Section 5.

3.3. Tidal and Limiting Radii

The Jacobi radius represents a theoretical surface around a cluster past which a star cannot pass and still remain bound. The Jacobi radius allows for the calculation of the instantaneous tidal radius of each model cluster as a function of time. To calculate the instantaneous tidal radius of each model cluster, we require a derivation of cluster tidal radius which takes into consideration the tidal field of the host galaxy. Assuming only that the tidal field must be spherically symmetric, the theoretical tidal radius as derived by Bertin & Varri (2008) is

$$\begin{equation} r_t=\left(\frac{GM}{\Omega ^2\upsilon }\right)^{1/3}, \end{equation} \tag{ 3 }$$

where Ω, κ, and $\upsilon$ are defined as

$$\begin{equation} \Omega ^2=(d\Phi _G(R)/dR)_{R_{p}}/R_{p}, \end{equation} \tag{ 4 }$$

$$\begin{equation} \kappa ^2=3\Omega ^2+(d^2\Phi _G(R)/dR^2)_{R_{p}}, \end{equation} \tag{ 5 }$$

$$\begin{equation} \upsilon =4-\kappa ^2/\Omega ^2, \end{equation} \tag{ 6 }$$

Φ_G is the galactic potential, R_p is the perigalactic distance, Ω is the orbital frequency of the cluster, κ is the epicyclic frequency of the cluster at R_p, and $\upsilon$ is a positive dimensionless coefficient. Using the tidal field of the Milky Way discussed in Section 2 and the mass and galactocentric distance of the model clusters at each time step, we calculate the instantaneous tidal radius of each model. The results of these calculations are shown in Figure 5. Here the instantaneous r_t increases and decreases periodically along the orbit, but never quite reaches the perigalactic and apogalactic cases due to differences in mass-loss rates among all three cases.

Zoom In Zoom Out Reset image size

Next we compare the tidal radius to the limiting radius. For a simulated cluster, since we know which stars are bound or unbound, we could call the limiting radius of the cluster the distance to the farthest bound star, but this approach introduces some significant problems. First, since cluster tidal radii are calculated for stars with circular prograde orbits, any star with a retrograde and/or eccentric orbit within its cluster can remain bound beyond the nominal tidal radius (Read et al. 2006). Second, any star in the process of escaping the cluster can reach large cluster-centric distances before becoming energetically unbound. Third, the stars along the y-axis of the cluster (the direction of motion) are unaffected by the tidal potential in Equation (2), allowing them to also remain bound at larger cluster-centric distances. These three issues cause the true limiting radius of the cluster to change dramatically from time step to time step. To gain a more stable indication of cluster size, we instead focus on the x-axis of the cluster, the axis along which the tidal radius is calculated, and define the limiting radius as the average x-coordinate of all stars with ||x|| > r_t. This calculation typically involves less than 1% of the total cluster population. While this is not the true limiting radius of the cluster and will always be slightly larger than the true tidal radius, it acts as a tracer of the outer region that is less affected by individual extreme outliers. If a cluster is tidally overfilling, the limiting radius will still be significantly larger than the tidal radius. For a cluster that is tidally underfilling, the limiting radius is simply the distance to the outermost bound star.

In Figure 6 we show this empirically determined r_L for each model as a function of time. For circular orbits, on average the limiting radius of the cluster decreases smoothly as a result of mass loss. For eccentric orbits, the small fluctuations with perigalactic passes in Figures 3 and 4 are much more prevalent in Figure 6. Comparing Figure 5 to Figure 6, the fluctuations in both figures indicate that the limiting radius behaves the same as the instantaneous tidal radius.

Zoom In Zoom Out Reset image size

It should be noted that the apogalactic cases for e = 0.75 and 0.9 in Figure 6 are not smooth due to the fact that these clusters are barely tidally filling and their limiting radii are easily influenced by individual escaping stars.

Directly comparing limiting radii in Figure 6 and tidal radii in Figure 5, all circular orbits have a relatively constant ratio at approximately r_L/r_t = 1.1. Since we expect r_L to slightly overestimate r_t, a ratio of 1.1 suggests that these clusters are approximately tidally filling. For eccentric clusters the ratio is in general also 1.1, suggesting that the clusters come close to filling their instantaneous tidal radius at all times. Fluctuations in the ratio for the e = 0.75 and 0.9 cases indicate that after a perigalactic pass the cluster is slightly tidally underfilling and works to fill its instantaneous tidal radius on the way to apogalacticon. When traveling back in from apogalacticon to perigalacticon, the cluster will remain tidally filled and lose stars to tidal stripping as the instantaneous tidal radius shrinks.

Up until this point we have only considered clusters with initial half-mass radii of 6 pc. This initial half-mass radius was chosen simply to ensure that the model clusters with R_p = 6 kpc would be tidally filling. As seen in Figure 4 this produces clusters with sizes at 10 Gyr ranging from 3 to 14 pc, which are larger than most (but not all) real globular clusters. In an attempt to produce Milky-Way-like clusters which have a mean effective radius of 2.5 pc, we re-simulated the e0r6rm6, e05r18rm6, and e0r18rm6 models with initial half-mass radii of 4 pc, 2 pc, 1 pc, and 0.5 pc. The results are shown in plots of half-mass radius versus time in Figure 7.

Zoom In Zoom Out Reset image size

The rm4 clusters closely resemble the rm6 clusters, with similar periodic fluctuations with perigalactic passes and the final r_m values for the perigalactic, eccentric, and apogalactic cases. However, the rm4 clusters undergo smaller initial expansion due to two-body interactions than the rm6 clusters, and thus when outer region stars are being removed through tidal stripping, since the majority of the mass is concentrated in the inner region, the mass-loss profile is less affected.

This issue of initial size becomes even more significant in the rm2, rm1, and rm05 cases. The periodic fluctuations on the eccentric orbit are barely visible in the rm2 cluster, and non-existent in the rm1 and rm05 clusters. The rm1 cases are significantly smaller at 10 Gyr, approaching the typical ∼2–3 pc sizes that match the majority of real globular clusters. The rm05 models have completely dissolved by 10 Gyr.

These small clusters are only tidally filling in the sense that two-body interactions have pushed some bound stars to orbits that take them out to the instantaneous tidal radius. With the majority of the bound stars located in the inner regions of the cluster, tidal stripping is not the dominant form of mass loss and the influence of the galactic potential and cluster orbit are minimized. Instead, stellar evolution and two-body interactions are the dominant forms of mass loss. These clusters would be classified as "tidally unaffected" (Carballo-Bello et al. 2012). For the rm1 case, the r_m profiles of the perigalactic, eccentric, and apogalactic cases begin to split only after ∼5 Gyr, when stars have been finally pushed to the outer regions of the cluster and tidal stripping is beginning to play an important role. While this is true, it is not due to two-body interactions but instead a result of core collapse. For the rm05 case, not even core collapse can push stars to the outer region of the cluster in order for tidal stripping to occur. In fact, the rm05 clusters all have the same r_m up until the complete evaporation of the cluster at approximately 7 Gyr.

Producing Milky-Way-like globular cluster effective radii of 2–3 pc for clusters on circular orbits at 6 kpc or greater appears to require initial r_m size less than 1 pc. The N-body models reveal that either clusters originally are extremely compact and tidally unaffected, or present-day cluster orbits have changed significantly from the orbit along which the clusters originally formed. Some observational support for this view can be found in recent measurements of very young, massive clusters (e.g., Bastian et al. 2008, 2012; Portegies Zwart et al. 2010; Marks & Kroupa 2010). However, recent N-body simulations by Sippel et al. (2012) showed that r_h and r_m can be very different because of stellar mass segregation, and produce clusters with final r_h values near 3 pc despite large r_{m, i}. These issues will be explored in future studies.

A perigalactic pass has three effects on a globular cluster, which we illustrate in Figure 8 for the e = 0.9 model e09r104rm6. In this figure, we plot the energy per unit mass of individual stars (as per Equation (2)) as a function of radial distance from the cluster center at nine points in the orbit. Beginning in Panel (A) of Figure 8, for a given time between apogalacticon and perigalacticon there are a few stars that are within close proximity of the cluster but remain unbound (marked in red). As the cluster moves toward R_p and the instantaneous tidal radius shrinks (Panel (B)), more and more stars become temporarily unbound. As predicted in Section 3, even stars in the inner region of the cluster are provided with a significant increase in energy by the tidal shock and can become unbound. Just after the cluster reaches perigalacticon (Panel (C)), a large number of stars are no longer bound to the cluster. As the cluster moves away from perigalacticon (Panels (D)–(I)):

• 1.  
  some stars that became unbound escape the cluster (which causes the initial decrease in r_m and r_L);
• 2.  
  some of the stars that are unbound in Panel (C) return to energies below zero and are recaptured (see Figure 9);
• 3.  
  the tidal shock gives stars initially found in the inner region enough additional energy to move outward and fill the orbits vacated by stars which permanently escaped the cluster (see Figure 10).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is even possible for inner region stars to become temporarily or permanently unbound if they move outward at a rate faster than the instantaneous tidal radius increases.

Now that we have shown that limiting radii are not imposed at perigalacticon, it is useful to know how to calculate a meaningful number that predicts the limiting radius of a globular cluster on an eccentric orbit. As we saw in Figures 5 and 6, the limiting radius essentially traces the instantaneous tidal radius. However, the ratio between cluster limiting radius and instantaneous tidal radius undergoes small periodic variations as a function of the location of a cluster along its orbit.

Orbital phase is defined as

$$\begin{equation} F = \frac{R_{gc} - R_p}{R_{a}-R_p} \end{equation} \tag{ 7 }$$

such that the cluster has F = 0 at perigalacticon and F = 1 at apogalacticon. A cluster with a circular orbit will always have F = 0. The median limiting radius of the eccentric cluster as a function of orbital phase is then determined in order to calculate the ratio of the instantaneous limiting radius to the limiting radius of the e = 0 perigalactic case (r_L(e)/r_L(e = 0)), both normalized by mass. This ratio is plotted as a function of phase in Figure 11. It is important to note that we have ignored the first 2000 Myr of evolution for each model cluster when evaporation due to two-body relaxation is the dominant source of mass loss.

Zoom In Zoom Out Reset image size

For a given orbital eccentricity, r_L(e)/r_L(e = 0) changes almost linearly with phase F. It is interesting to note that we observed a second-order effect that the rate at which the cluster expands is lower than the rate at which it contracts. When the cluster is moving away from perigalacticon, it works to fill its expanding tidal radius. Conversely, a cluster moving toward perigalacticon would be larger as it is always tidally filling on the way inward.

The rate at which the limiting radius of the cluster increases and decreases as a function of orbital phase is a strong function of orbital eccentricity. These rates were determined explicitly by finding the slopes of each line, where the y-intercept is forced to equal one. The slopes are listed in Table 2, along with the associated uncertainty (1σ).

A smooth relationship between slope and eccentricity emerges, which can be fit with an exponential, and allows us to propose a purely analytical correction to the calculation of the tidal radius of a globular cluster. For a globular cluster on an orbit with eccentricity E, with a tidal radius at perigalacticon r_t(R_p) and located at a phase F in its orbit, its limiting radius is equal to

$$\begin{equation} {r_L(F) = r_t(R_p) (1 + aFe^{bE})}, \end{equation} \tag{ 8 }$$

where a = 0.17 ± 0.03 and b = 4.1 ± 0.2. Note that since this calculation involves a single cluster over the course of a single orbit, the mass normalization is no longer necessary as cluster mass will not have changed significantly over a fraction of one orbit.

The next step is to simulate a larger suite of model clusters ranging in initial cluster mass and half-mass radii to determine if these parameters play a role in the correction factor suggested above. However, regardless of the influence of cluster mass or initial half-mass radius, all tidally affected simulations follow the rule that the limiting radius traces the instantaneous tidal radius rather well. Thus, if full orbital information or phase F is unknown, the calculation of the instantaneous tidal radius is a reasonable estimate of the limiting radius of a cluster. For globular clusters in other galaxies, in which only their projected galactocentric distances are known, it may be possible to determine their theoretical tidal radius based on their present King radius r_k. Future work will explore this possibility.

6.1. Application to the Milky Way

For many Milky Way globular clusters, their current galactocentric distance, orbital eccentricity, and orbital phase are known. In Figure 12 the revised, fully corrected version of Figure 1 is illustrated, where r_L is now the phase-corrected value from Equation (8). We now see more tidally underfilling clusters and the scatter of points more nearly around zero. A stronger agreement between theory and observations emerges. Correcting for using a non-spherically symmetric potential in calculating tidal radii and improved methods for determining observational limiting radii will likely strengthen this comparison further.

Zoom In Zoom Out Reset image size

Globular clusters have been simulated with a range of both circular and eccentric orbits. After determining which stars are bound to the cluster at a given time, we show that while eccentric clusters undergo episodes of significant mass loss during a perigalactic pass, their time-averaged mass-loss rate reflects a tidal field less than the tidal field at perigalacticon. Additionally, it was found that clusters are able to recapture unbound stars after a perigalactic pass as their instantaneous tidal radius increases.

Second, we show that the half-mass radius of a globular cluster increases and decreases about a mean value over the course of an orbit. These fluctuations suggest that the perigalactic pass also has the effect of energizing inner region stars to larger orbits. Finally, we find that the limiting radius of a cluster traces its instantaneous tidal radius at all times.

These findings argue against the historical assumption that globular cluster tidal radii, and by extension limiting radii, are imposed at perigalacticon for clusters that do not have circular orbits. While it remains true that the half-mass relaxation time is greater than one orbital period, the cluster does not need to fully relax in order to expand. The eccentric orbit introduces an effect of tidal shocking that is not experienced by clusters in a static potential (circular orbit).

While the instantaneous tidal radius is a useful first approximation of the limiting radius, we have proposed an analytically determined correction factor that is a function of orbital eccentricity and phase. This correction leads to a much stronger agreement between the predicted limiting radii and observational King (1962) radii of Milky Way globular clusters. Future studies will explore how the correction factor depends on initial mass or initial half-mass radius and how corrected limiting radii are related to King radii.

Since the tidal field of the Milky Way is not spherically symmetric, correcting limiting radii based on eccentricity and orbital phase is not the final step. We still need to correct for orbital inclination to account for factors like disk shocking and tidal heating, which may reveal important effects for the Milky Way and other disk galaxies. However, the present results already have clear applicability to elliptical galaxies, which have more nearly spherical potentials. We are currently investigating N-body simulations in these directions. The ultimate goal is to be able to predict the limiting radius of any tidally affected globular cluster, given its orbit, galactocentric position, and the galactic potential of the host galaxy.

J.W. acknowledges funding through the A. Boyd McLay Ontario Graduate Scholarship. J.W. also thanks Juan P. Madrid and Anna Sippel for valuable discussions and correspondence. A.S. and W.E.H. acknowledge financial support through research grants from the Natural Sciences and Engineering Research Council of Canada.