Evolution of star clusters on eccentric orbits: semi-analytical approach

TL;DR

This paper extends a semi-analytical model to simulate the evolution of star clusters on eccentric orbits, closely matching detailed N-body simulations while being computationally more efficient.

Contribution

The authors adapt and extend the EMACSS code to accurately model star cluster evolution on eccentric orbits, enabling faster analysis of large cluster samples.

Findings

Semi-analytical models closely match N-body simulation results.

Models effectively track cluster evolution from pre- to post-collapse.

Approach significantly reduces computational time compared to N-body methods.

Abstract

We study the dynamical evolution of star clusters on eccentric orbits using a semi-analytical approach. In particular we adapt and extend the equations of EMACSS code, introduced by Gieles et al. (2014), to work with eccentric orbits. We follow the evolution of star clusters in terms of mass, half-mass radius, core radius, Jacobi radius and the total energy over their dissolution time. Moreover, we compare the results of our semi-analytical models against $N$-body computations of clusters with various initial half-mass radius, number of stars and orbital eccentricity to cover both tidally filling and under-filling systems. The evolution profiles of clusters obtained by our semi-analytical approach closely follow those of $N$-body simulations in different evolutionary phases of star clusters, from pre-collapse to post-collapse. Given that the average runtime of our semi-analytical models is significantly less than that of $N$-body models, our approach makes it feasible to study the evolution of large samples of globular clusters on eccentric orbits.

Tables

Table 1. Table 1: Initial parameters and inputs of each performed semi-analytical models and N-body runs. The evolution of each model is shown in the figures which is identified in the first column. Columns 2 and 3 are the initial number of stars and initial half-mass radius, respectively. The orbital parameters are given in columns 4 to 6. The last two columns are the initial Jacobi radius and filling factor that are analytically calculated. See the text for more details.

Figures	$N_0$	$r_{h0}$	$R_P$	$R_A$	$e$	$r_{J0}$	${ℛ}_{hJ0}$
#		[pc]	[kpc]	[kpc]		[pc]
1	5 000	3	10	30	0.5	19.23	0.16
1	8 000	3	10	30	0.5	22.50	0.13
1	16 000	3	10	30	0.5	28.36	0.11
2	4 000	0.3	10	30	0.5	17.87	0.02
2	4 000	1	10	30	0.5	17.87	0.03
2	4 000	3	10	30	0.5	17.87	0.17
2	4 000	6	10	30	0.5	17.87	0.34
3	8 000	3	30	30	0	71.13	0.04
3	8 000	3	16.1	30	0.3	37.10	0.08
3	8 000	3	5.3	30	0.7	11.66	0.26
4	8 000	3	15	15	0	35.57	0.08
4	8 000	3	8.1	15	0.3	18.60	0.16
4	8 000	3	5	15	0.5	11.20	0.26
4	8 000	3	2.65	15	0.7	5.86	0.51

Equations

$$E=-\kappa \frac{G{M}^2}{r_h},$$

$${\tau }_{rh}=0.138\frac{N^{1/2}{r}_h^{3/2}}{\sqrt{Gm}\ln (\gamma N)},$$

$$ϵ\equiv -(\dot{E}/E){\tau }_{rh},$$

$$\xi \equiv -(\dot{N}/N){\tau }_{rh}=-(\dot{M}/M){\tau }_{rh},$$

$$\mu \equiv ({\dot{r}}_h/r_h){\tau }_{rh},$$

$$\lambda \equiv (\dot{\kappa }/\kappa ){\tau }_{rh},$$

$$\displaystyle\delta\equiv(\dot{r}_{c}/r_{c})\tau_{rh}.\$$

$$ϵ=-\lambda +\mu +2\xi .$$

$$R_G=a(1-e\cos \eta ),$$

$$\displaystyle t=(T/2\pi)[(\eta-\eta_{0})-e\sin(\eta-\eta_{0})],\$$

$$r_{J}=\bigg{(}\dfrac{M}{M_{G}}\bigg{)}^{1/3}\bigg{(}\dfrac{aR_{G}^{4}}{R_{A}R_{P}+2aR_{G}}\bigg{)}^{1/3}.$$

$$\mathcal{R}_{ch}=\bigg{(}\dfrac{N_{2}}{N}+\dfrac{N_{2}}{N_{3}}\bigg{)}^{2/3}$$

$${\rho }_c=$$

$$ϵ=({\mathcal{R}}_{hJ}/\kappa )\xi$$

$$ϵ=\zeta \simeq 0.1$$

$$\xi =\mathcal{F}{\xi }_1(1-\mathcal{P})+(3/5)\zeta [(f+(1-f)\mathcal{F}]\mathcal{P}.$$

$$\mathcal{F}=$$

$$\mathcal{P}(e=0)=\bigg{(}\dfrac{\mathcal{R}_{vJ}}{\mathcal{R}_{vJ1}}\bigg{)}^{z}\bigg{[}\dfrac{N\ln(\gamma N_{1})}{N_{1}\ln(\gamma N)}\bigg{]}^{1-x}.$$

$$\delta =$$

$${\tau }_{rc}=\frac{{\sigma }_c^3}{15.4G^{2}m{\rho }_c\ln (\gamma N)},$$

$$\mu =$$

$$\kappa ({\mathcal{R}}_{ch})={\kappa }_1+({\kappa }_0-{\kappa }_1)\mathrm{erf}({\mathcal{R}}_{ch}/{\mathcal{R}}_{ch0}).$$

$$\mathcal{K}=\frac{2({\kappa }_0-{\kappa }_1){\mathcal{R}}_{ch}\exp (-{\mathcal{R}}_{ch}^2/{\mathcal{R}}_{ch0}^2)}{\sqrt{\pi }\kappa {\mathcal{R}}_{ch0}}.$$

$$\lambda =$$

$$T_{diss}(a,e)=T_{diss}(a,0)(1-e^2)(1-c{e}^2),$$

$$\mathcal{P}(e)=\bigg{(}\dfrac{\mathcal{R}_{vJ}}{\mathcal{R}_{vJ1}}\bigg{)}^{z}\bigg{[}\dfrac{N\ln(\gamma N_{1})}{N_{1}\ln(\gamma N)}\bigg{]}^{1-x}\mathcal{Q}(e).$$

$$\mathcal{Q}(e)=[(1-e^2)(1-c{e}^2){]}^{-1}.$$

Full text

††pagerange: Evolution of star clusters on eccentric orbits: semi-analytical approach–Evolution of star clusters on eccentric orbits: semi-analytical approach††pubyear: 2017

1 Introduction

The evolution of star clusters is driven by internal and external processes. Two-body relaxation, stellar evolution and binary formation are the main internal drivers and the galactic tidal field is the most significant external factor to the disruption of star clusters (Spitzer, 1987). Due to the complexity of star clusters, numerical and computer codes are often used to study their evolution. Among all available methods, direct $N$-body codes such as NBODY6/7, which includes many details (e.g. the dynamical formation of higher order systems, stellar evolution of single and binary stars as well as a 3D external tidal field) generally produce the most accurate results (Aarseth, 1999; Heggie & Hut, 2003).

Although the runtime of $N$-body simulations have decreased by the advent of GRAPE computer series (Makino, 2001) and Graphics Processing Units (GPUs, Nitadori & Aarseth 2012) and the direct modeling of some globular clusters (GCs) is now possible (e.g. Zonoozi et al. 2011, 2014, 2017), clusters with $N\gtrsim 10^{5}$ stars are still computationally expensive to simulate especially in the presence of primordial binaries (Wang et al., 2015, 2016). The other notable methods for simulating star clusters that are based on the statistical modeling of GCs include different techniques like Fokker-Planck, Monte Carlo and gaseous model simulations (Spurzem et al., 2005). Monte Carlo codes (Hénon, 1971a, b; Giersz, 1998) such as MOCCA (Hypki & Giersz, 2013) have been shown to follow the outcome of $N$-body simulations closely, albeit they need to be calibrated using $N$-body simulations first (Giersz, Heggie & Hurley, 2008). The main advantage of such Monte Carlo codes is that their runtime scales with the number of stars ($N$), hence they are faster than direct $N$-body simulations whose runtime scales as $\mathcal{O}(N^{2})$.

In terms of analytical attempts to study the evolution of star clusters, one can mention the pioneering analytical work of Hénon who used a Fokker-Planck equation for gravitational encounters, for tidally limited (Hénon, 1961) and isolated clusters (Hénon, 1965). Hénon showed that a single-mass star cluster in a tidal field loses its mass in a self-similar manner and an isolated cluster keeps its mass and expands.

Alexander & Gieles (2012) introduced emacss, a fast C++ code, to simulate the evolution of single-mass star clusters on circular orbits affected by two-body relaxation and an external tidal field. In the first version of the code, emacss solved two coupled differential equations for number of stars and half-mass radius and predicted the evolution of these parameters accurately in comparison with $N$-body results. In the second version of emacss (Gieles et al., 2014, hereafter G14), the core evolution of a single-mass star cluster was added and the pre-collapse phase or the unbalanced evolution which is an important phase in the evolution of star clusters in a strong tidal field is included. In the third version (Alexander et al., 2014), the effect of stellar evolution given an initial mass function was added to emacss.

emacss is useful to study the evolution of a large sample of star cluster since it is much faster than both Monte Carlo and $N$-body simulations. However, in its current form it is only limited to circular orbits. The Milky Way GCs with determined orbital properties have eccentric orbits (Dinescu, Girard & van Altena, 1999; Casetti-Dinescu et al., 2007, 2013). The effect of eccentric orbits over the life-time of star clusters is not negligible (Küpper et al., 2015b). The escape rate of stars from a star cluster depends on its galactocentric distance. While a cluster passes through its perigalacticon, the stars gain enough energy to go to outer orbits due to tidal heating and shocking and can leave the cluster (Webb et al., 2013, 2014). Baumgardt & Makino (2003) evolved a large set of $N$-body simulations for star clusters both on circular and eccentric orbits. They showed that the dissolution time of star clusters on eccentric orbits is approximately proportional to those on circular orbit with the radius equal to the apogalactic distances as $T_{diss}(e)\sim T_{diss}(R_{A},0)(1-e)$, where, $e$ is the orbital eccentricity and $T_{diss}(R_{A},0)$ is the lifetime of a cluster moving on a circular orbit with a radius equal to the apogalactic radius of the eccentric orbit. However, recently, Xu Cai et al. (2016) updated this relation using a series of N-body simulations and found that the dissolution time of a star cluster on an eccentric orbit is proportional to a circular orbit with the radii equal to the semi-major axis as $T_{diss}(e)\sim T_{diss}(a,0)(1-e^{2})(1-ce^{2})$ where $c\simeq 0.5$ and $T_{diss}(a,0)$ is the lifetime of a cluster moving on a circular orbit with a radius equal to the semi-major axis of the eccentric orbits.

In the present paper, we aim to extend the work of Alexander et al. (2014) to study star clusters moving on eccentric orbits. In particular we achieve this by adapting the equations of emacss to work with eccentric orbits. In order to check the accuracy of our method, we compare our results with $N$-body models.

The present paper is organized as follows. In Section 2 we describe the equations that we have used to model the evolution of star clusters on eccentric orbits. We compare the results of our semi-analytical models against $N$-body simulations in Section 3. Finally, in Section 5 we conclude and summarize our work.

2 Overview of semi-analytical approach

In this section, we set out the theoretical framework that is used to describe the evolution of the mass, half-mass radius, and core radius of star clusters. We start with introducing the key parameters and equations that determine the evolution of star clusters, the time-scales and definitions that is used in this paper and overview the basic physical details of escape considered by our prescription.

2.1 Equations of Evolution

We assume a single-mass star cluster with total mass $M$ including $N=M/m$ stars, where $m$ is the mass of each star and is constant for all stars. In the present paper we have not considered the effect of an initial mass function, stellar evolution or primordial binaries. This enable us to isolate the influence of eccentric orbits compared to circular orbits on the evolution of star clusters. We assume that star clusters move on eccentric orbits around a point-mass host galaxy with mass $M_{G}$, located at galactocentric distance $R_{G}(t)$ which changes as a function of time.

The total energy of a star cluster in virial equilibrium is given by (Spitzer, 1987)

[TABLE]

where $G$ is the gravitational constant, $r_{h}$ is the half-mass radius and $\kappa$ is the form factor that depends on the density profile of the star cluster and can be written as $\kappa=r_{h}/(4r_{v})$, where $r_{v}$ is the virial radius (G14).

The theoretical starting point of this paper is based on the diffusion of energy in a star cluster. In the lack of energy sources such as stellar evolution and primordial binaries, the single-mass star cluster contracts in the early phase of its evolution. Therefore the inner region of cluster heats up and loses energy to the outer halo. This process continues to the extent that the core radius ($r_{c}$) becomes very small and the core density grows significantly (Takahashi, 1995). This process is called core collapse and it stops when the first hard binary forms. The evolution phase before the core collapse is called pre-collapse phase or unbalanced evolution. After the core collapse, the energy flows from the core into the halo with a constant rate by the two-body relaxation (Alexander & Gieles, 2012). This evolution phase is called post-collapse or balanced evolution. The flux of energy through the star cluster changes its characteristic parameters, i.e. mass, half-mass radius, core radius and density profile. The time-scale for dynamical evolution of single-mass cluster is determined by half-mass relaxation time, $\tau_{rh}$, as (Spitzer & Hart, 1971)

[TABLE]

where $\ln({\gamma N})$ is the Coulomb logarithm with $\gamma\approx 0.11$ for single-mass star clusters (Giersz & Heggie, 1994). Following G14, we can define five dimensionless parameters for the evolution of each of the characteristic parameters of star clusters per $\tau_{rh}$ as

[TABLE]

Note that the equality in equation (4) is the direct result of the constant mean mass of stars. In equation (3), $\epsilon$ represents the energy production rate. Mass loss from the star cluster can be obtained from $\xi$ in equation (4) which remains positive during the life-time of the shrinking star cluster. According to equations (5) and (7), evolution of the half-mass radius and the core radius ($r_{c}$) can be evaluated from $\mu$ and $\delta$, respectively. Taking the time derivative on both sides of equation (1) and using equations (3) to (6), one can relate four of the dimensionless parameters as

[TABLE]

Note that $\delta$ does not appear in equation (8). This is due to the fact that the total energy does not depend on the core radius directly (see equation (1)). The evolution of the core radius on eccentric orbit is one of the goals of this paper and we can evaluate it from the equations described in section 2.2.

The apogalactic ($R_{A}$) and perigalactic ($R_{P}$) radii of the star cluster orbit are related to eccentricity $e$, and semi-major axis $a$, by $R_{A}=a(1+e)$ and $R_{P}=a(1-e)$, respectively. The galactocentric distance of the star cluster and elapsed time are related via the following equations (Binney & Tremaine, 2008)

[TABLE]

where $\eta$ is the eccentric anomaly. Equation (10) is the Kepler’s equation. In it $T=2\pi\sqrt{a^{3}/(GM_{G})}$ is the orbital period and $\eta_{0}$ is the initial phase, i.e. if $\eta_{0}=0$, the star cluster begins its evolution at perigalacticon. We use the Newton-Raphson method to numerically solve the Kepler’s equation and determine $R_{G}(t)$ as a function of time.

The external tidal field puts a limiting boundary for the star cluster so that if stars gain enough energy to pass this boundary, they become unbound from the star cluster and scatter into the host galaxy. This boundary is specified by the Jacobi radius $r_{J}$. For a star cluster on an eccentric orbit around a point-mass galaxy, the instantaneous Jacobi radius is given by Ernst & Just (2013)

[TABLE]

In this equation $R_{G}$ is time-dependent and can be calculated using equation (9) and equation (10).

It should be noted that a strict theoretical definition of the Jacobi radius exists only for circular orbits and for a cluster on an eccentric orbit there is no clear definition of the Jacobi radius. As a result, equation (11) is not accurate for the case of the eccentric orbits (See Sec. 4 for more discussion).

2.2 Review of the required formulas and initial conditions

To determine the time-evolution of the characteristic parameters of star clusters, we solve the five coupled differential equations (3) to (7). To achieve this we need additional expressions to relate the five dimensionless parameters ($\epsilon,\xi,\mu,\lambda,\delta$) to other parameters of the star clusters. Here, we review the additional formulas needed to relate the five coupled differential equations (3) to (7) as given by G14. For clarity we follow the notation of G14.

First of all, It is useful to define three ratios as $\mathcal{R}_{hJ}\equiv r_{h}/r_{J}$, $\mathcal{R}_{vJ}\equiv r_{v}/r_{J}$, $\mathcal{R}_{ch}\equiv r_{c}/r_{h}$. The ratio $\mathcal{R}_{hJ}$ (filling factor) denotes the strength of the galactic tidal field, i.e. the star clusters with initial $\mathcal{R}_{hJ}\gtrsim 0.1$ are initially Roche volume filling clusters and those with initial $\mathcal{R}_{hJ}<0.1$ are initially Roche volume under-filling clusters (Alexander & Gieles, 2013). However following G14, we use $\mathcal{R}_{vJ}$ instead of $\mathcal{R}_{hJ}$ in mass-loss rate (equation (20)).

We differentiate between the pre-collapse phase (unbalanced evolution, hereafter UE) and the post-collapse phase (balanced evolution, hereafter BE). Therefore, we present two sets of equations for each of the five dimensionless parameters ($\epsilon,\xi,\mu,\lambda,\delta$). The transition from UE to BE occurs at the core collapse time. The core collapse occurs when $\mathcal{R}_{ch}$ in balanced phase reaches the value that is given by (Heggie & Hut, 2003, G14)

[TABLE]

where $N_{2}$ and $N_{3}$ are two free parameters and are chosen to be $N_{2}=12$ and $N_{3}=15\,000$. During the pre-collapse phase, the core radius decreases and the core density, $\rho_{c}$, grows. In the post-collapse phase, the core radius increases and the core density reduces as $r_{c}^{-2}$. The relation of the core density can be written as

[TABLE]

In equation (13), the relation between the core radius and density is a power-law relation whose exponent $\alpha$, has been determined in several studies, e.g. $\alpha=2.21$ (Heggie & Stevenson, 1988), $\alpha=2.23$ (Takahashi, 1995), $\alpha=2.26$ (Baumgardt et al., 2003) and $\alpha=2.20$ (G14). We adopt $\alpha=2.2$ and $\rho_{c0}=0.055M_{0}r_{v0}^{-0.8}$. In equation (14), $\rho_{h}$ is the half-mass density and is defined as $\rho_{h}\equiv 3M/(8\pi r_{h}^{3})$.

The energy production rate, $\epsilon$, in the UE is time-dependent and increases as stars escape. In the BE this rate is approximately constant or based on Alexander & Gieles (2012) notation $\zeta\simeq 0.1$. Thus we can express the energy production rate as

[TABLE]

In equation (15), $\xi$ is the escape rate of stars from the star cluster which can be obtained as (G14)

[TABLE]

The second term on the right hand side of equation (17) accounts for the escape rate in a tidal field and the first term denotes the mass-loss of isolated cluster with $\xi_{1}=0.0142$. In equation (17), $f=0.3$ is the ratio of escape rate of a Roche volume filling cluster to that in BE, and $\mathcal{F}$ causes a smooth transition from UE to BE and is given by

[TABLE]

where $\mathcal{R}_{ch}^{min}$ is the minimum value of $\mathcal{R}_{ch}$ that occurs at core collapse time and is obtained by equation (12).

Due to the fact that the mechanism of escape rate changes in the case of eccentric orbit, we need to alter the $\mathcal{P}$-parameter in equation (17). Alexander & Gieles (2012) defined this parameter for circular orbit as

[TABLE]

The parenthesis part in equation (20) is the result of integration over the Maxwellian velocity distribution to obtain the fraction of stars with enough escape energy (Gieles & Baumgardt, 2008) and the bracket part is proportional to the time delay of the stars which has gained enough energy to escape from the star cluster (Baumgardt, 2001). The values of the free parameters in equation (20) are chosen by G14 to be $z=1.61$, $x=0.75$, $\mathcal{R}_{vJ1}=0.145$. $N_{1}$ is the scaling value of $N$ that is tuned by G14 to be $N_{1}=15\,000$ to match the N-body simulations for circular orbits. The values of $N_{1}$, $\mathcal{R}_{vJ1}$, and consequently $\mathcal{P}$ for the eccentric orbits are different from those in G14 for circular orbit. In section 2.3 we introduce a new method for fixing the $\mathcal{P}-$parameter as a function of eccentricity (based on the findings by Baumgardt & Makino (2003) and Xu Cai et al. (2016) such that the semi-analytical models, can reproduce the N-body simulation for different eccentricities.

For the evolution of the core radius we have

[TABLE]

In equation (21), $\delta_{1}=-0.09$, $\delta_{2}=-0.002$ and $\tau_{rc}$ is the core relaxation time and is defined as (Spitzer & Hart, 1971)

[TABLE]

where $\sigma_{c}^{2}=(8/3)\pi G\rho_{c}r_{c}^{2}$. Note that we can derive equation (22) by combining equations (5), (7) and (12).

Letting $\mathcal{K}\equiv\mathrm{d}\ln\kappa/\mathrm{d}\ln\mathcal{R}_{ch}$, the evolution of half-mass radius can be described as (G14)

[TABLE]

The combination of equations (8) and (15) leads to equation (24). Moreover, one can derive equation (25) using equations (8), (16) and (22).

Finally, we should determine the evolution of form factor. G14 proposed a relation between form factor and $\mathcal{R}_{ch}$ by fitting $N$-body results. They showed that the form factor in this relation is not necessarily equal to $\kappa$. Labeled as $\kappa(\mathcal{R}_{ch})$, the fitting function is

[TABLE]

Here, the three free parameters have different values in each phase: $\kappa_{0}=r_{h0}/(4r_{v0})$, $\kappa_{1}=0.295$ and $\mathcal{R}_{ch0}=0.100$ for UE and $\kappa_{0}=0.200$, $\kappa_{1}=0.265$ and $\mathcal{R}_{ch0}=0.220$ for BE.

Using the logarithmic derivative of $\mathcal{K}$ and equation (26), we have

[TABLE]

Using the definition of $\mathcal{K}$ and equations (5) to (7), we can derive the evolution rate of the form factor as

[TABLE]

We have written a code which numerically solves the coupled differential equations for an eccentric orbit111The code can be obtained by contacting the authors. using fourth-order Runge-Kutta with a constant time step of $\Delta t=1$ Myr. For each model the input parameters are $N_{0}$ (initial number of stars), $r_{h0}$ (initial half-mass radius), $\eta_{0}$ (initial phase), $R_{A}$ and $R_{P}$. The initial values of mass and form factor can be given as $M_{0}=mN_{0}$ and $\kappa_{0}=r_{h0}/(4r_{v0})$, respectively. The initial values of core and virial radius for a Plummer density profile are $r_{c0}\simeq 0.4r_{h0}$ and $r_{v0}\simeq 1.3r_{h0}$, respectively (Heggie & Hut, 2003, Table 8.1).

2.3 Determination of $\mathcal{P}-$parameter for eccentric orbits

Apart from the initial conditions, the lifetime of the cluster is proportional to ${\mathcal{R}^{z}_{vJ1}}[N_{1}/\ln(\gamma N_{1})]^{1-x}$ (see e.g., equation 23 in Alexander & Gieles 2012 or a simplified lifetime equation A10 in Gieles, Heggie & Zhao 2011). According to Xu Cai et al. 2016, there is a relation between the lifetime of a star cluster ($T_{diss}(a,e)$) on eccentric orbit with half major axis $a$, and the lifetime of a star cluster on circular orbit with radius $a$, $T_{diss}(a,0)$, as

[TABLE]

where $e$ is the eccentricity and $c\simeq 0.5$.

Therefore, one can obtain a formula for $\mathcal{P}$-parameter for a cluster on an eccentric orbit as a function of eccentricity and the corresponding parameter for a cluster evolving on a circular orbit ($\mathcal{P}(0)$) by multiplying $\mathcal{P}(0)$ by $\mathcal{Q}(e)$.

[TABLE]

where, $\mathcal{Q}(e)$ is a function of eccentricity as

[TABLE]

Using the above formulae for all models with different initial masses, radii, and eccentricities, one can see the shape of the $M(t)$, $r_{h}(t)$, and $r_{c}(t)$ curves of semi-analytical models are in general agreement with the N-body simulation (see Figs. 1 to 4). Our results also indicate that the scaling of the star cluster’s lifetime as a function of eccentricity as obtained by Xu Cai et al. (2016) is very successful to estimate the lifetime of star clusters on eccentric orbits and to reproduce the N-body simulation for different eccentricities.

3 Comparison with $N$-body models

In this section the evolution of star clusters resulted from our semi-analytical models are compared against those of N-body simulations. We used the GPU-enabled version of the collisional fourth-order Hermite $N$-body code nbody6 (Nitadori & Aarseth, 2012). We chose a Plummer model (Plummer, 1911) in virial equilibrium as the initial density profile for our models. All stars have equal mass $m=0.5M_{\odot}$. The stars with distances more than $2r_{J}$ were assumed to be unbound from the cluster and the simulations were terminated when less than one percent of stars remained in the clusters. The modeled star clusters are on eccentric orbits around a point-mass galaxy with mass $M_{G}=10^{11}M_{\bigodot}$. All sets of models are summarized in Table 1.

3.1 Models with $e=0.5$ and different initial masses and half-mass radii

First, we performed three simulations with $N_{0}=$ 5000, 8000, 16000 and the same half-mass radius $r_{h0}=3$ pc. Moreover, four additional simulations are performed with different half-mass radii of $r_{h0}=0.3,1,3,6$ pc and the same number of stars $N_{0}=$ 4000. We assume that the apogalactic and perigalactic distances of the orbits are $R_{A}=30$ kpc and $R_{P}=10$ kpc and the evolution starts at the perigalacticon. Thus the eccentricity of all orbits is $e=0.5$ which corresponds to the mean value for the globular star clusters in the Milky Way (Dinescu, Girard & van Altena, 1999; Casetti-Dinescu et al., 2007, 2013).

We first model three star clusters with different $N_{0}$ and equal half-mass radii $r_{h0}=3$ pc. The evolution of mass, half-mass radius, core radius, total energy, and form factor are shown in Fig. 1. As can be seen, the mass of these models decreases with a staircase pattern as expected from Xu Cai et al. (2016). This pattern is a result of the variation of the galactocentric distance with time which consequently causes the Jacobi radius to oscillate. In other words, when passing through the perigalacticon (apogalacticon), the star cluster experiences the strongest (weakest) tidal field on its orbit and its Jacobi radius reaches its minimum (maximum) value. As a result, the mass loss rate is larger at perigalacticon than that at the apogalacticon and the staircase shape is built. The evolution of mass which is obtained by our semi-analytical model is in good agreement with $N$-body results.

According to the evolution of $r_{h}$, three evolutionary phases are recognized: (i) in the pre-collapse phase, because of the absence of energy source in core, $r_{h}$ decreases and the amplitude of the oscillation is lower. (ii) in the post collapse phase, the flux of energy from the core to the outer regions of the cluster (as a result of two-body relaxation) leads to the expansion of cluster and $r_{h}$ increases (expansion-dominated regime). (iii) if the star cluster resides in a tidal field, its expansion is eventually halted at some points and the evaporation of stars from cluster begins until complete dissolution (evaporation-dominated regime in post-collapse phase).

Although the oscillation amplitude of our models is lower than that of $N$-body simulations, in both approaches the amplitude of the oscillations of $r_{h}$ in the post-collapse phase is higher than those in the pre-collapse phase. Moreover, the amplitude of the oscillations is even higher in the latter part of evolution.

According to Fig. 1, core radius of star clusters at core collapse moment is smaller than that at $t=0$ but it increases in the post-collapse phase. The core radius oscillates in $N$-body simulations. Our semi-analytical models do not reproduce these oscillations, since the oscillations come from the escape rate and the core evolution rate ($\delta$) does not depend on it (see equation (21)). The evolution of $r_{c}$ is in good agreement with $N$-body simulations in the post-collapse phase.

The total energy increases in both pre- and post-collapse phases but with different rates and the oscillation in energy occurs only in the pre-collapse phase in our semi-analytical models (equation (15)). The form factor increases in the pre-collapse phase and then reaches a constant value in the post-collapse phase and at the end of the cluster life time decreases. The amplitude of the Jacobi radius decreases in all models which is related to the mass-loss rate and remarkably follow the $N$-body results.

Fig. 2 shows the evolution of mass, half-mass radius, core radius, total energy, and form factor of four star clusters with different $r_{h0}$ and the same number of stars ($N_{0}=4\,000$). These figures show how our semi-analytical method works in different tidal regimes. From left to right, we model four star cluster with different initial half-mass radii. For these models the value of $\mathcal{R}_{hJ}=r_{h}/r_{J}$ increases from 0.015 to 0.30 covering both weak and strong tidal fields. According to the classification of star clusters in a tidal field in Alexander & Gieles (2013), the first two models are Roche volume under-filling clusters and the other models are Roche volume filling clusters. The evolution of mass in these models show that the mass-loss rate in the early stages of evolution of Roche volume filling clusters in both $N$-body simulations and our semi-analytical models are higher. For instance, in the very weak tidal field, cluster loses around $~{}40\%$ of its mass during the first orbital period. The Roche-volume under-filling clusters reach core collapse after just a few Myrs, whereas the Roche-volume filling clusters spend a longer time in the pre-collapse phase (few orbital periods). In the strong tidal field, the star cluster spends roughly half of its lifetime in the pre-collapse phase. The dissolution time of these models increase with $r_{h0}$. The form factor ranges between 0.19 and 0.26 approximately in all models. In the Roche volume filling clusters the form factor remains nearly constant during the pre-collapse phase. This is because the core and half-mass radii decrease at the same rate and at the instant of the core collapse, $\kappa$ rises.

3.2 Models with different eccentricities

In section 3.1 all models are evolving on eccentric orbit with the specific eccentricity, $e=0.5$, and corresponding value of $N_{1}=2\,025$ was fixed for this eccentricity to match with the N-body simulation. In order to show that our semi-analytical model can reproduce the N-body simulation for different eccentricities we performed a new set of simulations with eccentricity $e=0,0.3,0.5,0.7$ (where eccentricity is defined as $e=(R_{A}-R_{P})/(R_{A}+R_{P})$ and two different value of apocentric distance of $R_{A}=$ 15 and 30 kpc to cover both tidally filling and under-filling systems.

The results of our simulations for clusters with orbital eccentricities of $e=$0, 0.3, 0.7 and the same apogalactic distances of $R_{A}=30$ kpc are illustrated in Fig. 3. We also repeat the calculations for a lower value of apocentric distance of $R_{A}=15$ kpc to see how the higher filling factor affects the results (Fig. 4).

During the perigalactic passage, the star clusters with high eccentricities experience a stronger tidal field which leads to an increase in their escape rate. Moreover, the oscillation amplitude of the Jacobi radii is larger for star clusters with higher eccentricities. Thus the slope of decreasing mass is higher in high eccentricity cases and the staircase pattern of mass evolution is more clear. We conclude that the lifetime of the star clusters increases with decreasing the eccentricity which is in agreement with the scaling formula for lifetime suggested by Xu Cai et al. (2016).

Comparing two models with the same eccentricity of $e=0.7$ but different pericentric and apocentric distances which corresponds to different filling factor, one can see easily that the more tidally filling cluster (Table 1) can not reproduce the N-body simulation. Therefore, our semi-analytical formalism is valid for both tidally filling and under-filling systems and violates only for the very high-eccentricity filling regime. A cluster with high eccentricity ($e=0.7$) spends less than half of its lifetime in pre-collapse phase but for the other two cases with the lower eccentricity of 0.5 and 0.3, the post-collapse phase dominates the pre-collapse phase and the half-mass radius grows up more in the first half part of post-collapse phase and consequently leads to a longer lifetime.

4 Caveats: Tidal and limiting radii

As stated before tidal field is no longer static for a cluster with a non-circular orbit. This makes estimating the Jacobi radius for a GC on an eccentric orbit to be challenging. The variation of tidal field on a cluster moving on an eccentric orbit can be rapid at pericentre and results in a change in the shape of the Jacobi surface from a static approximation (Renaud, Gieles & Boily, 2011) as well as an injection of additional energy into the star cluster (Weinberg, 1994a, b, c; Küpper et al., 2010; Küpper, Lane & Heggie, 2012). We therefore model eccentric orbits in an approximate manner.

One approach in this regard is assuming that for a GC on an eccentric orbit, its tidal radius is imposed at perigalacticon ($R_{p}$) where the tidal field of the host galaxy is the strongest and the Jacobi radius has a minimum. That is, derivative of $r_{J}$ is zero at this radius (von Hoerner, 1957; King, 1962). This assumption is based on the fact that the internal relaxation time of the cluster is greater than its orbital period for almost all observed GCs. 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. In other words, the satellites are truncated during pericentre passages to the size indicated by the pericentric tidal radius (King, 1962; Inannen, Harris & Webbink, 1983). However, this is indeed true for collisionless systems such as dwarf satellite galaxies that are orbiting around the host giant galaxy, N-body simulations of collisional systems show that after the pericentre passage the satellite expands again and hence, the King (1962) conjecture that the satellites are trimmed at the pericentre and then remain unchanged is not valid (Gajda & Łokas, 2016). If the cluster moves on an eccentric orbit, stars outside the tidal radius will likely become unbound (temporarily) at perigalacticon, where the tidal radius reaches its minimum. When the cluster moves away from perigalacticon and the instantaneous tidal radius of the cluster increases again, some stars are able to be re-captured by the cluster (Küpper, Lane & Heggie, 2012). Therefore, the Jacobi radius of a cluster will be greater than the perigalactic tidal radius, and there are no clear relationship between limiting radii and perigalactic distance (see e.g., Odenkirchen et al. 1997).

Another approach is using the orbit-averaged tidal radius (Brosche, Odenkirchen & Geffert, 1999). This approach is already tested by several authors. For example, Küpper et al. (2010) and Küpper, Lane & Heggie (2012) found that the time averaged mean tidal radius of the cluster and not the perigalactic tidal radius is better approximation to reproduce the structure of tidal tails. Moreover, N-body simulations by Madrid, Hurley & Sippel (2012) found that the half-mass radius of a GC is more likely imposed at $R_{A}$ than $R_{P}$.

Alternatively, one can use the instantaneous Jacobi (tidal) radius which corresponds to the distance of $L_{1}/L_{2}$ Lagrange points, as if the satellite was on a circular orbit of radius equal to its current distance from the host. Webb et al. (2013) explored the influence of the orbital eccentricity on the tidal radius using the direct $N$-body experiments and showed that the limiting radius is not imposed at $R_{P}$, instead, it traces the "instantaneous tidal radius" of the cluster at any point in the orbit. It is shown that the assumption of the instantaneous Jacobi radius is good approximation for satellite galaxies orbiting around the main host galaxy (e.g., see figure 5 of Gajda & Łokas 2016). Here in this work we follow this approach and will show that with this assumption the semi-analytical models follow the N-body simulations remarkably well.

Using the mass and galactocentric distance of the model clusters at each time step, we calculate the instantaneous Jacobi radius of each model which is not constant. The oscillation of $R_{G}(t)$ leads to the oscillation of $r_{J}$ and then the Jacobi radius is not constant and the instantaneous Jacobi radius increases and decreases periodically along the orbit between apogalactic and perigalactic radii. In order to check the validity of our assumptions on the instantaneous Jacobi radius made in the semi-analytical models, the evolution of different parameters in semi-analytical models is compared to a realistic $N$-body models.

In order to put the criteria for determining whether a star is bound or unbound we invoke a distance cutoff such that the cluster-centric distance of stars must be greater than the cluster instantaneous tidal radius for it to be unbound. This criterion has been used frequently in $N$-body modeling of star clusters (e.g. Takahashi & Baumgardt 2012; Haghi et al. 2015). It has also been suggested that a star’s velocity plays a role in determining the unbound stars (e.g. Küpper et al. 2010; Küpper, Lane & Heggie 2012). Given the different definition of bound stars, the evolution of bound stars displays different oscillating and staircase pattern, during the pericentre passage. These artificial behaviours that emerge due to the definition of bound for clusters, illustrate that it is not possible to have a universal definition of the bound stars for clusters on eccentric orbits. Webb et al. (2013) found that this additional criteria only effected a small percentage of simulated stars and did not change any of the $N$-body results. Moreover, Xu Cai et al. (2016) showed that the differences between the evolution of cluster for different definition of bound are small. Therefore, for both N-body models and semi-analytical calculations we only assumed the Jacobi radius as a distance cutoff to determine the star to be considered unbound. These assumptions are tested in Section 3 (see Figs 1 to 4) where we show that the $N$-body results are followed by our semi-analytical models and hence, our definition of bound, i.e., being within the instantaneous Jacobi radius, is a good approximation.

5 Summary

Although modeling a few low-density medium-sized GCs of the Milky Way with $N=10^{5}$ stars, over their entire lifetime is now possible by a direct $N$-body approach, modeling the majority of the Galactic GCs with a realistic density distribution and large number of stars (i.e., $N\geq 5\times 10^{5}$) over a Hubble time is still beyond our computational ability.

Therefore, the faster codes such as mocca and emacss allow us to model the large number of Galactic GCs with realistic initial density, size and number of stars. But, both methods are limited for clusters moving on circular orbits.

Following the approach introduced in emacss code, we extend the model to include the dynamical evolution of single-mass star clusters in more realistic elliptical orbits. We calculated the evolutionary equations of emacss code so that the eccentric orbits can be now handled properly. Our approach addresses both the pre-collapse and post-collapse evolutionary phases of star clusters over their entire life-time.

We compared the evolution of single-mass star clusters using the semi-analytical approach with the outcome of $N$-body simulations and showed that the evolution of all parameters are in good agreement with each other. Our models include different initial sizes, number of stars and orbital eccentricities such that they cover both Roche-volume filling and under-filling systems with focusing merely on eccentric orbits.

The main advantage of our adopted approach is the runtime of the models (seconds to minutes) which is several orders of magnitude shorter than both Monte Carlo (hours to days) and $N$-body methods (days to months).

We intend to use our method to study the evolution of GC systems and model the observed properties of star cluster populations in a time-dependent galaxy potential.

Acknowledgements

We would like to thank the anonymous referee for his/her useful comments and suggestions which improved the quality of this work. We thank Holger Baumgardt and Antonio Solima for helpful discussions. Pouria Khalaj is a team member of the StarFormMapper project which has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 687528.