import clustertools as cts
import numpy as np
from galpy.orbit import Orbit
from galpy.potential import MWPotential2014,PlummerPotential
from galpy.util import conversion

Orbits and Tidal Tails

For cases where a cluster’s galactocentric position and velocity are known and a galpy potential (Bovy J., 2015, ApJS, 216, 29) is given, it is possible to calculate additional cluster properties based on their orbit.

To begin, lets load the final snapshot of an N-body simulation that was meant to reproduce the Galactic globular cluster Pal 5. Note one could have used load_cluster(ctype='limepy',gcname='Pal5') to generate a Pal 5 - like cluster, but it would not have any tidal tails (which will also be discussed here). The external tidal field is assumed to be equal to the MWPotential2014 model from Bovy (2015). For all types of orbit analysis, this can be replaced with a different galpy potential using the pot variable.

cluster = cts.load_cluster('snapshot',filename='pal5.dat', units='kpckms',origin='galaxy')
cts.starplot(cluster)

0

If you wish to extract the orbit of the cluster as a galpy orbit to analyze within galpy, simply initialize the orbit via:

cluster.initialize_orbit()

<galpy.orbit.Orbits.Orbit at 0x156e9ef70>

After the above step, cluster.orbit now contains the galpy orbit. You can also integrate the orbit forwards to tfinal Gyr, with the option to do the standard galpy orbit plot. For illustrative purposes, I have shown how the pot variable can be used to specify the background potential.

ts=np.linspace(0,1/conversion.time_in_Gyr(ro=cluster.orbit._ro,vo=cluster.orbit._vo),1000)
cluster.orbit.integrate(ts,MWPotential2014)
cluster.orbit.plot()

[<matplotlib.lines.Line2D at 0x145c592b0>]

If instead you wish to extract that orbital path +/- dt Gyr from the cluster’s current position, you can run the below command which also has a plot option. Note the default dt is 0.1 Gyr, but for this particular snapshot a dt of 0.3 Gyr is needed to approximately cover the length of the tails.

t,x,y,z,vx,vy,vz=cluster.orbital_path(plot=True,pot=MWPotential2014,dt=0.3)

A related feature implemented in clustertools is the ability to find the distance of each star from the stream path (dpath) and its distance along the stream path from the progenitor (dprog). When calling orbital_path_match, the timestep along the orbital path that each star is closest to is also returned (tpath).

t,dprog,dpath=cluster.orbital_path_match(pot=MWPotential2014,plot=True,dt=0.3)

Alternatively one can calculate what clustertools calls the tail path. To determine the tail path, clustertools first matches each star to the orbital path, bins stars along the orbital path, and then finds the mean position and velocity of stars in the same bin. The tail path better follows tidal tail stars which are generally offset from the orbital path due to them having non-zero velocities when they escape the cluster.

cluster.tail_path(pot=MWPotential2014,plot=True,dt=0.3,filename='../images/tail_path.png')
cluster.tail_path_match(pot=MWPotential2014,plot=True,dt=0.3,filename='../images/tail_path_match.png')

(array([-0.003, -0.003, -0.003, ...,  0.033, -0.135,  0.045]),
 array([  0.        ,   0.        ,   0.        , ...,   8.16113693,
        -22.17925342,  12.00729019]),
 array([-0.1377047 ,  0.13570313,  0.13583119, ...,  0.05525546,
        -0.3445197 , -0.15470062]))

Finally, clustertools also allows for a system to be interpolated to a past or future timestep. Note this does NOT handle any internal cluster evolution, as it simply moves all cluster stars to a new orbital position. The advanage to the orbit interpolation function is that you can specificy what constitutes a tail star and integrate their orbits separately. This function allows one to more directly compare the tidal tail systems of N-body simulations to an observed cluster by ensuring the N-body system is at the exact same position as the progenitor.

For example, lets start with an NGC 6101-like system initialized using the load_cluster, where a lowered isothermal model from LIMPY is used to generate NGC 6101 based on the best fir parameters from de Boer et al. 2019. Not we initialize the model with origin='galaxy' and units='kpckms'.

cluster = cts.load_cluster('limepy',gcname='NGC6101',origin='galaxy',units='kpckms')
cts.starplot(cluster,do_centre=True)

LOAD_CLUSTER MADE USE OF:
Gieles, M. & Zocchi, A. 2015, MNRAS, 454, 576
Vasiliev E., 2019, MNRAS, 484,2832
de Boer, T. J. L., Gieles, M., Balbinot, E., Hénault-Brunet, V., Sollima, A., Watkins, L. L., Claydon, I. 2019, MNRAS, 485, 4906

0

If we then interpolate the system 1 Gyr ahead in time, the centre of the cluster will move but the clustercentric positions and velocities of the stars themselves will stay the same

cluster.interpolate_orbit(pot=MWPotential2014,tfinal=1.0)
cts.starplot(cluster,do_centre=True)

0

Alternatively you can use interpolate_orbits which will interpolate the orbits of stars within the cluster. If origin=centre or origin=cluster, clustertools will assume it is the cluster’s potential that is being provided and orbits will be integrated within the cluster. If origin='galaxy', clustertools will assume its the galactic potential being provided. This type of interpolation is really only valid for a stellar stream, unless the combined potential of the galaxy and cluster are provided (see MovingPotential within galpy’s documentation).

cluster = cts.load_cluster('limepy',gcname='NGC6101',origin='cluster',units='pckms')
cts.starplot(cluster,do_centre=True)

LOAD_CLUSTER MADE USE OF:
Gieles, M. & Zocchi, A. 2015, MNRAS, 454, 576
Vasiliev E., 2019, MNRAS, 484,2832
de Boer, T. J. L., Gieles, M., Balbinot, E., Hénault-Brunet, V., Sollima, A., Watkins, L. L., Claydon, I. 2019, MNRAS, 485, 4906

0

To interpolate stars within the cluster, lets assume the cluster’s potential is a Plummer Sphere. Also note that tfinal is given in Myr since units=pckms

pot=PlummerPotential(cluster.mtot/conversion.mass_in_msol(ro=cts.solar_ro,vo=cts.solar_vo),b=cluster.rm/1.35/1000.0/cts.solar_ro,ro=cts.solar_ro,vo=cts.solar_vo)
cluster.interpolate_orbits(pot=pot,tfinal=1000.0)

(array([  4.08001261,  -1.73009228,   0.98027894, ..., -22.84690542,
         12.15213098,   7.30701486]),
 array([ -0.84052929,  -1.82759056,   3.99883035, ...,   4.06053046,
         32.32265448,  13.91795629]),
 array([  6.24731999,   3.08734876,  -2.38146302, ..., -23.63221055,
        -18.65189181,  -3.68864173]),
 array([  5.32189375e-01,  -3.12568732e-01,   3.62392791e-01, ...,
         -3.38536676e-01,  -4.72837211e-04,  -3.18064879e-01]),
 array([-0.12288152, -0.19749976, -0.31355414, ...,  0.30403016,
        -0.3263936 , -0.24657509]),
 array([ 0.70982236,  0.77686439, -0.14442882, ..., -0.17207202,
         0.20851848, -0.13694187]))

A quick check reveals the stellar positions have moved

cts.starplot(cluster,do_centre=True)

0