Galactic bar

Integrate the globular cluster system as point masses where the galaxy now includes a galactic bar. Thus, we no are no longer in a time-static and axis-symmetric potential

import tstrippy
import astropy.units as u
import astropy.constants as const
import astropy.coordinates as coord
import numpy as np
import matplotlib.pyplot as plt

Load the units

unitbasis = tstrippy.Parsers.potential_parameters.unitbasis
unitT=u.Unit(unitbasis['time'])
unitV=u.Unit(unitbasis['velocity'])
unitD=u.Unit(unitbasis['distance'])
unitM=u.Unit(unitbasis['mass'])
unitG=u.Unit(unitbasis['G'])
G = const.G.to(unitG).value

Load Galactic Parameters

potentialname='pouliasis2017pii'
MWparams = tstrippy.Parsers.potential_parameters.pouliasis2017pii()
MWrefframe= tstrippy.Parsers.potential_parameters.MWreferenceframe()

Add the bar parameters

This is an artificial bar used in Ferrone et al 2023. It is not intended to perfectly model the bar, but rather provide an idea of what a bar can do to stellar streams. It uses a Long-Murali bar. This bar is a needle that has been convolved with a soften kernal to give it a triaxial-shape.

# mass and size
abar = 4 * unitD
bbar = 1 * unitD
cbar = 0.5 * unitD
Mbar = 990.0*2.32*1e7 * unitM
barname = "longmuralibar"
barparams = [MWparams[0],Mbar.value,abar.value,bbar.value,cbar.value]
# oreitnation and bar pattern speed
theta0= 25 * (np.pi/180)
omega =  28  * unitV / unitD
omega = -omega.value
barpolycoeff=[theta0,omega]

Reduce the mass of the disks

Since we are adding a bar, we need to reduce the mass of the disks

MWparams[5] = 1120.0 * 2.32*10**7
MWparams[8] = 1190.0 * 2.32*10**7

Set integration time

Do this in years so that you can assure the number of timesteps is an integer, and then convert the total time T and the timestep dt into integration units.

integrationtime,dt = 5e9,1e5
Ntimestep=int(integrationtime/dt)
integrationtime,dt=integrationtime*u.yr,dt*u.yr
dt=dt.to(unitT)
currenttime = 0.0*unitT
integrationparams = [currenttime,dt,Ntimestep]

Load the Globular Clusters

Store them in an astropy object to store them to galactocentric coordinates

GCdata  =   tstrippy.Parsers.baumgardtMWGCs().data
skycoordinates=coord.SkyCoord(
    ra=GCdata['RA'],
    dec=GCdata['DEC'],
    distance=GCdata['Rsun'],
    pm_ra_cosdec=GCdata['mualpha'],
    pm_dec=GCdata['mu_delta'],
    radial_velocity=GCdata['RV'],)
galacticcoordinates = skycoordinates.transform_to(MWrefframe)
x,y,z=galacticcoordinates.cartesian.xyz.to(unitD).value
vx,vy,vz=galacticcoordinates.velocity.d_xyz.to(unitV).value
# grab the number of objects
Nobj=len(x)

Initialize the integrator

organize the input arguments

staticgalaxy = [potentialname,MWparams]
initialkinematics = [x,y,z,vx,vy,vz]
integrationparams = [currenttime.value,dt.value,Ntimestep]
galacticbar = [barname,barparams,barpolycoeff]
leapfropintime = [Ntimestep,Nobj]

input the arguments and compute the integral

tstrippy.integrator.setstaticgalaxy(*staticgalaxy)
tstrippy.integrator.setinitialkinematics(*initialkinematics)
tstrippy.integrator.setintegrationparameters(*integrationparams)
# tstrippy.integrator.setbackwardorbit()
tstrippy.integrator.initgalacticbar(*galacticbar)
xt,yt,zt,vxt,vyt,vzt=tstrippy.integrator.leapfrogintime(*leapfropintime)
timestamps = tstrippy.integrator.timestamps.copy()
tstrippy.integrator.deallocate()

Analyize the change in energy

def transform_to_bar_frame(x,y,z,barpolycoeff,thetime):
    """
    Transform the coordinates to the bar frame
    """
    theta0,omega=barpolycoeff
    theta=theta0+omega*thetime
    xbar=x*np.cos(theta)+y*np.sin(theta)
    ybar=-x*np.sin(theta)+y*np.cos(theta)
    zbar=z
    return xbar,ybar,zbar

E0 = np.zeros(Nobj)
_,_,_,phi0=tstrippy.potentials.pouliasis2017pii(MWparams,xt[:,0],yt[:,0],zt[:,0])
_,_,_,phiF=tstrippy.potentials.pouliasis2017pii(MWparams,xt[:,-1],yt[:,-1],zt[:,-1])
T0 = 0.5*(vxt[:,0]**2+vyt[:,0]**2+vzt[:,0]**2)
TF = 0.5*(vxt[:,-1]**2+vyt[:,-1]**2+vzt[:,-1]**2)

# grab the bar, who changes in time
xbar0,ybar0,zbar0=transform_to_bar_frame(xt[:,0],yt[:,0],zt[:,0],barpolycoeff,timestamps[0])
xbarF,ybarF,zbarF=transform_to_bar_frame(xt[:,-1],yt[:,-1],zt[:,-1],barpolycoeff,timestamps[-1])
_,_,_,phibar0=tstrippy.potentials.longmuralibar(barparams,xbar0,ybar0,zbar0)
_,_,_,phibarF=tstrippy.potentials.longmuralibar(barparams,xbarF,ybarF,zbarF)

E0 = T0 + phi0 + phibar0
EF = TF + phiF + phibarF
Lz0 = xt[:,0]*vyt[:,0]-yt[:,0]*vxt[:,0]
LzF = xt[:,-1]*vyt[:,-1]-yt[:,-1]*vxt[:,-1]

import matplotlib
colors = matplotlib.cm.rainbow(np.linspace(0,1,Nobj))
# make the font size bigger and also us Latex font
plt.rc('text', usetex=True)
plt.rc('font', size=20)
plt.rc('axes', titlesize=20)
plt.rc('axes', labelsize=20)
plt.rc('xtick', labelsize=20)
plt.rc('ytick', labelsize=20)
plt.rc('legend', fontsize=20)
plt.rc('figure', titlesize=20)

fig,axis=plt.subplots(1,1,figsize=(10,7))
AXIS = {
    "xlabel": r"$L_z$ [kpc km/s]",
    "ylabel": r"$E$ [km$^{2}$ s$^{-2}$]",
    "title": "Integration time: {:d} Gyr".format(int(integrationtime.to(u.Gyr).value)),
}
for i in range(Nobj):
    if i ==0:
        axis.scatter(Lz0[i],E0[i],color=colors[i],label="Initial")
        axis.scatter(LzF[i],EF[i],color=colors[i],marker='x',label="Final")
    else:
        axis.scatter(Lz0[i],E0[i],color=colors[i])
        axis.scatter(LzF[i],EF[i],color=colors[i],marker='x')
    axis.plot([Lz0[i],LzF[i]],[E0[i],EF[i]],color=colors[i],alpha=0.5,zorder=0)
axis.set(**AXIS)
axis.tick_params(axis='y', rotation=45)
axis.legend()
fig.tight_layout()

View the most relatively changed clusters

dE = np.abs(E0 - EF)
dLz = np.abs(Lz0 - LzF)
meanE = np.mean(dE)
meanLz = np.mean(dLz)
metric = np.sqrt((dE/meanE)**2 + (dLz/meanLz)**2)
top4=np.argsort(metric)[-5:-1]
# targetGC = "NGC5139"
for i in top4:
    print("GC: {:s} has metric: {:.2f}".format(GCdata['Cluster'][i],metric[i]))

GC: Gran2 has metric: 5.78
GC: NGC5946 has metric: 7.36
GC: NGC6681 has metric: 9.54
GC: NGC4372 has metric: 10.93

nstamps = len(timestamps)
colors = matplotlib.cm.rainbow(np.linspace(0,1,nstamps))

fig= plt.figure(figsize=(12,10))
gridspec=fig.add_gridspec(2,3,width_ratios=[1,1,0.05],height_ratios=[1,1])
axes=[[],[]]
axes[0].append(fig.add_subplot(gridspec[0,0]))
axes[0].append(fig.add_subplot(gridspec[0,1]))
axes[1].append(fig.add_subplot(gridspec[1,0]))
axes[1].append(fig.add_subplot(gridspec[1,1]))
caxis=fig.add_subplot(gridspec[:,2])

im=axes[0][0].scatter(xt[top4[0],:],yt[top4[0],:],c=timestamps,marker='o',s=1,cmap=plt.cm.rainbow)
axes[0][0].set_aspect("equal")
axes[0][0].set_title("{:s}".format(GCdata['Cluster'][top4[0]]));

axes[0][0].set_ylabel("Y [kpc]");

im=axes[0][1].scatter(xt[top4[1],:],yt[top4[1],:],c=timestamps,marker='o',s=1,cmap=plt.cm.rainbow)
axes[0][1].set_aspect("equal")
axes[0][1].set_title("{:s}".format(GCdata['Cluster'][top4[1]]));



im=axes[1][0].scatter(xt[top4[2],:],yt[top4[2],:],c=timestamps,marker='o',s=1,cmap=plt.cm.rainbow)
axes[1][0].set_aspect("equal")
axes[1][0].set_title("{:s}".format(GCdata['Cluster'][top4[2]]));
axes[1][0].set_xlabel("X [kpc]");
axes[1][0].set_ylabel("Y [kpc]");

axes[1][1].scatter(xt[top4[3],:],yt[top4[3],:],c=timestamps,marker='o',s=1,cmap=plt.cm.rainbow)
axes[1][1].set_aspect("equal")
axes[1][1].set_title("{:s}".format(GCdata['Cluster'][top4[3]]));
axes[1][1].set_xlabel("X [kpc]");


cbar=fig.colorbar(im,cax=caxis,orientation='vertical')
cbar.set_label("Time [~Gyr]")
fig.tight_layout()

Co-rotation frame

fig= plt.figure(figsize=(12,10))
gridspec=fig.add_gridspec(2,3,width_ratios=[1,1,0.05],height_ratios=[1,1])
axes=[[],[]]
axes[0].append(fig.add_subplot(gridspec[0,0]))
axes[0].append(fig.add_subplot(gridspec[0,1]))
axes[1].append(fig.add_subplot(gridspec[1,0]))
axes[1].append(fig.add_subplot(gridspec[1,1]))
caxis=fig.add_subplot(gridspec[:,2])

for i in range(2):
    for j in range(2):
        index = i*2+j
        xbar,ybar,zbar=transform_to_bar_frame(xt[top4[index],:],yt[top4[index],:],zt[top4[index],:],barpolycoeff,timestamps)
        im=axes[i][j].scatter(xbar,ybar,c=timestamps,marker='o',s=1,cmap=plt.cm.rainbow)
        axes[i][j].set_aspect("equal")
        axes[i][j].set_title("{:s}".format(GCdata['Cluster'][top4[index]]));

cbar=fig.colorbar(im,cax=caxis,orientation='vertical')
axes[1][1].set_xlabel("X [kpc]");
axes[1][0].set_xlabel("X [kpc]");
axes[1][0].set_ylabel("Y [kpc]");
axes[0][1].set_ylabel("Y [kpc]");
cbar.set_label("Time [~Gyr]")
fig.suptitle("Bar Frame of most perturbed GCs")
fig.tight_layout()