Dynamical modeling of tidal streams

galpy contains tools to model the dynamics of tidal streams, making extensive use of action-angle variables and also containing a simple particle-spray method. Both of these are introduced in the following sections.

Modeling streams in action-angle coordinates with streamdf

Action-angle coordinates are ideal for modeling tidal streams in potentials that allow for the calculation of action-angle coordinates (any static potential in the absence of chaos, which describes at least regions such as the inner Milky-Way halo well). As an example of modeling streams in action-angle coordinates, we can investigate the dynamics of the following tidal stream (that of Bovy 2014; 2014ApJ…795…95B). This movie shows the disruption of a cluster on a GD-1-like orbit around the Milky Way:

The blue line is the orbit of the progenitor cluster and the black points are cluster members. The disruption is shown in an approximate orbital plane and the movie is comoving with the progenitor cluster.

Streams can be represented by simple dynamical models in action-angle coordinates. In action-angle coordinates, stream members are stripped from the progenitor cluster onto orbits specified by a set of actions \((J_R,J_\phi,J_Z)\), which remain constant after the stars have been stripped. This is shown in the following movie, which shows the generation of the stream in action space

The color-coding gives the angular momentum \(J_\phi\) and the black dot shows the progenitor orbit. These actions were calculated using galpy.actionAngle.actionAngleIsochroneApprox. The points move slightly because of small errors in the action calculation (these are correlated, so the cloud of points moves coherently because of calculation errors). The same movie that also shows the actions of stars in the cluster can be found here. This shows that the actions of stars in the cluster are not conserved (because the self-gravity of the cluster is important), but that the actions of stream members freeze once they are stripped. The angle difference between stars in a stream and the progenitor increases linearly with time, which is shown in the following movie:

where the radial and vertical angle difference with respect to the progenitor (co-moving at \((\theta_R,\theta_\phi,\theta_Z) = (\pi,\pi,\pi)\)) is shown for each snapshot (the color-coding gives \(\theta_\phi\)).

One last movie provides further insight in how a stream evolves over time. The following movie shows the evolution of the stream in the two dimensional plane of frequency and angle along the stream (that is, both are projections of the three dimensional frequencies or angles onto the angle direction along the stream). The points are color-coded by the time at which they were removed from the progenitor cluster.

It is clear that disruption happens in bursts (at pericenter passages) and that the initial frequency distribution at the time of removal does not change (much) with time. However, stars removed at larger frequency difference move away from the cluster faster, such that the end of the stream is primarily made up of stars with large frequency differences with respect to the progenitor. This leads to a gradient in the typical orbit in the stream, and the stream is on average not on a single orbit.

Action-angle modeling with streamdf in galpy

In galpy we can model streams using action-angle coordinateu with the tools in galpy.df.streamdf. We setup a streamdf instance by specifying the host gravitational potential pot=, an actionAngle method (typically galpy.actionAngle.actionAngleIsochroneApprox), a galpy.orbit.Orbit instance with the position of the progenitor, a parameter related to the velocity dispersion of the progenitor, and the time since disruption began. We first import all of the necessary modules

>>> from galpy.df import streamdf
>>> from galpy.orbit import Orbit
>>> from galpy.potential import LogarithmicHaloPotential
>>> from galpy.actionAngle import actionAngleIsochroneApprox
>>> from galpy.util import conversion #for unit conversions

setup the potential and actionAngle instances

>>> lp= LogarithmicHaloPotential(normalize=1.,q=0.9)
>>> aAI= actionAngleIsochroneApprox(pot=lp,b=0.8)

define a progenitor Orbit instance

>>> obs= Orbit([1.56148083,0.35081535,-1.15481504,0.88719443,-0.47713334,0.12019596])

and instantiate the streamdf model

>>> sigv= 0.365 #km/s
>>> sdf= streamdf(sigv/220.,progenitor=obs,pot=lp,aA=aAI,leading=True,nTrackChunks=11,tdisrupt=4.5/conversion.time_in_Gyr(220.,8.))

for a leading stream. This runs in about half a minute on a 2011 Macbook Air.

Bovy (2014) discusses how the calculation of the track needs to be iterated for potentials where there is a large offset between the track and a single orbit. One can increase the default number of iterations by specifying nTrackIterations= in the streamdf initialization (the default is set based on the angle between the track’s frequency vector and the progenitor orbit’s frequency vector; you can access the number of iterations used as sdf.nTrackIterations). To check whether the track is calculated accurately, one can use the following

>>> sdf.plotCompareTrackAAModel()

which in this case gives

This displays the stream model’s track in frequency offset (y axis) versus angle offset (x axis) as the solid line; this is the track that the model should have if it is calculated correctly. The points are the frequency and angle offset calculated from the calculated track’s \((\mathbf{x},\mathbf{v})\). For a properly computed track these should line up, as they do in this figure. If they do not line up, increasing nTrackIterations is necessary.

We can calculate some simple properties of the stream, such as the ratio of the largest and second-to-largest eigenvalue of the Hessian \(\partial \mathbf{\Omega} / \partial \mathbf{J}\)

>>> sdf.freqEigvalRatio(isotropic=True)
# 34.450028399901434

or the model’s ratio of the largest and second-to-largest eigenvalue of the model frequency variance matrix

>>> sdf.freqEigvalRatio()
# 29.625538344985291

The fact that this ratio is so large means that an approximately one dimensional stream will form.

Similarly, we can calculate the angle between the frequency vector of the progenitor and of the model mean frequency vector

>>> sdf.misalignment()
# 0.0086441947505973005

which returns this angle in radians. We can also calculate the angle between the frequency vector of the progenitor and the principal eigenvector of \(\partial \mathbf{\Omega} / \partial \mathbf{J}\)

>>> sdf.misalignment(isotropic=True)
# 0.02238411611147997

(the reason these are obtained by specifying isotropic=True is that these would be the ratio of the eigenvalues or the angle if we assumed that the disrupted materials action distribution were isotropic).

Calculating the average stream location (track)

We can display the stream track in various coordinate systems as follows

>>> sdf.plotTrack(d1='r',d2='z',interp=True,color='k',spread=2,overplot=False,lw=2.,scaleToPhysical=True)

which gives

which shows the track in Galactocentric R and Z coordinates as well as an estimate of the spread around the track as the dash-dotted line. We can overplot the points along the track along which the \((\mathbf{x},\mathbf{v}) \rightarrow (\mathbf{\Omega},\boldsymbol{\theta})\) transformation and the track position is explicitly calculated, by turning off the interpolation

>>> sdf.plotTrack(d1='r',d2='z',interp=False,color='k',spread=0,overplot=True,ls='none',marker='o',scaleToPhysical=True)

which gives

We can also overplot the orbit of the progenitor

>>> sdf.plotProgenitor(d1='r',d2='z',color='r',overplot=True,ls='--',scaleToPhysical=True)

to give

We can do the same in other coordinate systems, for example X and Z (as in Figure 1 of Bovy 2014)

>>> sdf.plotTrack(d1='x',d2='z',interp=True,color='k',spread=2,overplot=False,lw=2.,scaleToPhysical=True)
>>> sdf.plotTrack(d1='x',d2='z',interp=False,color='k',spread=0,overplot=True,ls='none',marker='o',scaleToPhysical=True)
>>> sdf.plotProgenitor(d1='x',d2='z',color='r',overplot=True,ls='--',scaleToPhysical=True)
>>> xlim(12.,14.5); ylim(-3.5,7.6)

which gives

or we can calculate the track in observable coordinates, e.g.,

>>> sdf.plotTrack(d1='ll',d2='dist',interp=True,color='k',spread=2,overplot=False,lw=2.)
>>> sdf.plotTrack(d1='ll',d2='dist',interp=False,color='k',spread=0,overplot=True,ls='none',marker='o')
>>> sdf.plotProgenitor(d1='ll',d2='dist',color='r',overplot=True,ls='--')
>>> xlim(155.,255.); ylim(7.5,14.8)

which displays

Coordinate transformations to physical coordinates are done using parameters set when initializing the sdf instance. See the help for ?streamdf for a complete list of initialization parameters.

Mock stream data generation

We can also easily generate mock data from the stream model. This uses streamdf.sample. For example,

>>> RvR= sdf.sample(n=1000)

which returns the sampled points as a set \((R,v_R,v_T,Z,v_Z,\phi)\) in natural galpy coordinates. We can plot these and compare them to the track location

>>> sdf.plotTrack(d1='r',d2='z',interp=True,color='b',spread=2,overplot=False,lw=2.,scaleToPhysical=True)
>>> plot(RvR[0]*8.,RvR[3]*8.,'k.',ms=2.) #multiply by the physical distance scale
>>> xlim(12.,16.5); ylim(2.,7.6)

which gives

Similarly, we can generate mock data in observable coordinates

>>> lb= sdf.sample(n=1000,lb=True)

and plot it

>>> sdf.plotTrack(d1='ll',d2='dist',interp=True,color='b',spread=2,overplot=False,lw=2.)
>>> plot(lb[0],lb[2],'k.',ms=2.)
>>> xlim(155.,235.); ylim(7.5,10.8)

which displays

We can also just generate mock stream data in frequency-angle coordinates

>>> mockaA= sdf.sample(n=1000,returnaAdt=True)

which returns a tuple with three components: an array with shape [3,N] of frequency vectors \((\Omega_R,\Omega_\phi,\Omega_Z)\), an array with shape [3,N] of angle vectors \((\theta_R,\theta_\phi,\theta_Z)\) and \(t_s\), the stripping time. We can plot the vertical versus the radial frequency

>>> plot(mockaA[0][0],mockaA[0][2],'k.',ms=2.)

or we can plot the magnitude of the angle offset as a function of stripping time

>>> plot(mockaA[2],numpy.sqrt(numpy.sum((mockaA[1]-numpy.tile(sdf._progenitor_angle,(1000,1)).T)**2.,axis=0)),'k.',ms=2.)

Evaluating and marginalizing the full PDF

We can also evaluate the stream PDF, the probability of a \((\mathbf{x},\mathbf{v})\) phase-space position in the stream. We can evaluate the PDF, for example, at the location of the progenitor

>>> sdf(obs.R(),obs.vR(),obs.vT(),obs.z(),obs.vz(),obs.phi())
# array([-33.16985861])

which returns the natural log of the PDF. If we go slightly higher in Z and slightly smaller in R, the PDF becomes zero

>>> sdf(obs.R()-0.1,obs.vR(),obs.vT(),obs.z()+0.1,obs.vz(),obs.phi())
# array([-inf])

because this phase-space position cannot be reached by a leading stream star. We can also marginalize the PDF over unobserved directions. For example, similar to Figure 10 in Bovy (2014), we can evaluate the PDF \(p(X|Z)\) near a point on the track, say near Z =2 kpc (=0.25 in natural units. We first find the approximate Gaussian PDF near this point, calculated from the stream track and dispersion (see above)

>>> meanp, varp= sdf.gaussApprox([None,None,2./8.,None,None,None])

where the input is a array with entries [X,Y,Z,vX,vY,vZ] and we substitute None for directions that we want to establish the approximate PDF for. So the above expression returns an approximation to \(p(X,Y,v_X,v_Y,v_Z|Z)\). This approximation allows us to get a sense of where the PDF peaks and what its width is

>>> meanp[0]*8.
# 14.267559400127833
>>> numpy.sqrt(varp[0,0])*8.
# 0.04152968631186698

We can now evaluate the PDF \(p(X|Z)\) as a function of X near the peak

>>> xs= numpy.linspace(-3.*numpy.sqrt(varp[0,0]),3.*numpy.sqrt(varp[0,0]),21)+meanp[0]
>>> logps= numpy.array([sdf.callMarg([x,None,2./8.,None,None,None]) for x in xs])
>>> ps= numpy.exp(logps)

and we normalize the PDF

>>> ps/= numpy.sum(ps)*(xs[1]-xs[0])*8.

and plot it together with the Gaussian approximation

>>> plot(xs*8.,ps)
>>> plot(xs*8.,1./numpy.sqrt(2.*numpy.pi)/numpy.sqrt(varp[0,0])/8.*numpy.exp(-0.5*(xs-meanp[0])**2./varp[0,0]))

which gives

Sometimes it is hard to automatically determine the closest point on the calculated track if only one phase-space coordinate is given. For example, this happens when evaluating \(p(Z|X)\) for X > 13 kpc here, where there are two branches of the track in Z (see the figure of the track above). In that case, we can determine the closest track point on one of the branches by hand and then provide this closest point as the basis of PDF calculations. The following example shows how this is done for the upper Z branch at X = 13.5 kpc, which is near Z =5 kpc (Figure 10 in Bovy 2014).

>>> cindx= sdf.find_closest_trackpoint(13.5/8.,None,5.32/8.,None,None,None,xy=True)

gives the index of the closest point on the calculated track. This index can then be given as an argument for the PDF functions:

>>> meanp, varp= meanp, varp= sdf.gaussApprox([13.5/8.,None,None,None,None,None],cindx=cindx)

computes the approximate \(p(Y,Z,v_X,v_Y,v_Z|X)\) near the upper Z branch. In Z, this PDF has mean and dispersion

>>> meanp[1]*8.
# 5.4005530328542077
>>> numpy.sqrt(varp[1,1])*8.
# 0.05796023309510244

We can then evaluate \(p(Z|X)\) for the upper branch as

>>> zs= numpy.linspace(-3.*numpy.sqrt(varp[1,1]),3.*numpy.sqrt(varp[1,1]),21)+meanp[1]
>>> logps= numpy.array([sdf.callMarg([13.5/8.,None,z,None,None,None],cindx=cindx) for z in zs])
>>> ps= numpy.exp(logps)
>>> ps/= numpy.sum(ps)*(zs[1]-zs[0])*8.

and we can again plot this and the approximation

>>> plot(zs*8.,ps)
>>> plot(zs*8.,1./numpy.sqrt(2.*numpy.pi)/numpy.sqrt(varp[1,1])/8.*numpy.exp(-0.5*(zs-meanp[1])**2./varp[1,1]))

which gives

The approximate PDF in this case is very close to the correct PDF. When supplying the closest track point, care needs to be taken that this really is the closest track point. Otherwise the approximate PDF will not be quite correct.

Modeling gaps in streams using action-angle coordinates

galpy also contains tools to model the effect of impacts due to dark-matter subhalos on streams (see Sanders, Bovy, & Erkal 2015). This is implemented as a subclass streamgapdf of streamdf, because they share many of the same methods. Setting up a streamgapdf object requires the same arguments and keywords as setting up a streamdf instance (to specify the smooth underlying stream model and the Galactic potential) as well as parameters that specify the impact (impact parameter and velocity, location and time of closest approach, mass and structure of the subhalo, and helper keywords that specify how the impact should be calculated). An example used in the paper (but not that with the modifications in Sec. 6.1) is as follows. Imports:

>>> from galpy.df import streamdf, streamgapdf
>>> from galpy.orbit import Orbit
>>> from galpy.potential import LogarithmicHaloPotential
>>> from galpy.actionAngle import actionAngleIsochroneApprox
>>> from galpy.util import conversion

Parameters for the smooth stream and the potential:

>>> lp= LogarithmicHaloPotential(normalize=1.,q=0.9)
>>> aAI= actionAngleIsochroneApprox(pot=lp,b=0.8)
>>> prog_unp_peri= Orbit([2.6556151742081835,
                          0.2183747276300308,
                          0.67876510797240575,
                          -2.0143395648974671,
                          -0.3273737682604374,
                          0.24218273922966019])
>>> V0, R0= 220., 8.
>>> sigv= 0.365*(10./2.)**(1./3.) # km/s
>>> tdisrupt= 10.88/conversion.time_in_Gyr(V0,R0)

and the parameters of the impact

>>> GM= 10.**-2./conversion.mass_in_1010msol(V0,R0)
>>> rs= 0.625/R0
>>> impactb= 0.
>>> subhalovel= numpy.array([6.82200571,132.7700529,149.4174464])/V0
>>> timpact= 0.88/conversion.time_in_Gyr(V0,R0)
>>> impact_angle= -2.34

The setup is then

>>> sdf_sanders15= streamgapdf(sigv/V0,progenitor=prog_unp_peri,pot=lp,aA=aAI,
                               leading=False,nTrackChunks=26,
                               nTrackIterations=1,
                               sigMeanOffset=4.5,
                               tdisrupt=tdisrupt,
                               Vnorm=V0,Rnorm=R0,
                               impactb=impactb,
                               subhalovel=subhalovel,
                               timpact=timpact,
                               impact_angle=impact_angle,
                               GM=GM,rs=rs)

The streamgapdf implementation is currently not entirely complete (for example, one cannot yet evaluate the full phase-space PDF), but the model can be sampled as in the paper above. To compare the perturbed model to the unperturbed model, we also set up an unperturbed model of the same stream

>>> sdf_sanders15_unp= streamdf(sigv/V0,progenitor=prog_unp_peri,pot=lp,aA=aAI,
                               leading=False,nTrackChunks=26,
                               nTrackIterations=1,
                               sigMeanOffset=4.5,
                               tdisrupt=tdisrupt,
                               Vnorm=V0,Rnorm=R0)

We can then sample positions and velocities for the perturbed and unperturbed preduction for the same particle by using the same random seed:

>>> numpy.random.seed(1)
>>> xv_mock_per= sdf_sanders15.sample(n=100000,xy=True).T
>>> numpy.random.seed(1) # should give same points
>>> xv_mock_unp= sdf_sanders15_unp.sample(n=100000,xy=True).T

and we can plot the offset due to the perturbation, for example,

>>> plot(xv_mock_unp[:,0]*R0,(xv_mock_per[:,0]-xv_mock_unp[:,0])*R0,'k,')

for the difference in \(X\) as a function of unperturbed \(X\):

or

>>> plot(xv_mock_unp[:,0]*R0,(xv_mock_per[:,4]-xv_mock_unp[:,4])*V0,'k,')

for the difference in \(v_Y\) as a function of unperturbed \(X\):

Particle-spray modeling of streams with streamspraydf

galpy also contains an implementation of the particle-spray method for generating tidal streams, roughly following the parametrization of Fardal et al. (2015). Full details on the galpy implementation are given in Qian et al. (2022). Here, we give a simple example of the method.

Like in the streamdf example above, we use the same orbit, potential, and cluster mass as in Bovy (2014). We setup the orbit of the progenitor and the gravitational potential (modeled as a simple LogarithmicHaloPotential):

>>> from galpy.potential import LogarithmicHaloPotential
>>> from galpy.orbit import Orbit
>>> o= Orbit([1.56148083,0.35081535,-1.15481504,0.88719443,-0.47713334,0.12019596])
>>> lp= LogarithmicHaloPotential(normalize=1.,q=0.9)

Then, we setup streamspraydf models for the leading and trailing arm of the stream:

>>> from astropy import units
>>> from galpy.df import streamspraydf
>>> spdf= streamspraydf(2*10.**4.*units.Msun,progenitor=o,pot=lp,tdisrupt=4.5*units.Gyr)
>>> spdft= streamspraydf(2*10.**4.*units.Msun,progenitor=o,pot=lp,leading=False,tdisrupt=4.5*units.Gyr)

To sample a set of 300 stars in both arms, we do

>>> orbs,dt= spdf.sample(n=300,returndt=True,integrate=True)
>>> orbts,dt= spdft.sample(n=300,returndt=True,integrate=True)

which returns a galpy.orbit.Orbit instance with all 300 stars. We can plot these in \(Z\) versus \(X\) and compare to Fig. 1 in Bovy (2014). First, we also integrate the orbit of the progenitor forward and backward in time for a brief period to show its location in the area of the stream:

>>> ts= numpy.linspace(0.,3.,301)
>>> o.integrate(ts,lp)
>>> of= o()
>>> of.integrate(-ts,lp)

Then we plot

>>> o.plot(d1='x',d2='z',color='k',xrange=[0.,2.],yrange=[-0.1,1.45])
>>> of.plot(d1='x',d2='z',overplot=True,color='k')
>>> plot(orbs.x(),orbs.z(),'r.')
>>> plot(orbts.x(),orbts.z(),'b.')

which gives

We can also compare to the track for this stream as predicted by streamdf. For this, we first setup a similar streamdf model (they are not exactly the same, as streamdf uses a velocity dispersion to set the progenitor’s mass, while streamspraydf uses the mass directly); see the streamdf documentation for a full explanation of this code:

>>> from galpy.actionAngle import actionAngleIsochroneApprox
>>> from galpy.df import streamdf
>>> aAIA= actionAngleIsochroneApprox(b=0.8,pot=lp)
>>> sigv= 0.365 #km/s
>>> sdf= streamdf(sigv/220.,progenitor=o(),pot=lp,aA=aAIA,leading=True,
                  nTrackChunks=11,tdisrupt=4.5*units.Gyr)
>>> sdft= streamdf(sigv/220.,progenitor=o(),pot=lp,aA=aAIA,leading=False,
                   nTrackChunks=11,tdisrupt=4.5*units.Gyr)

Then, we can overplot the track predicted by streamdf:

>>> o.plot(d1='x',d2='z',color='k',xrange=[0.,2.],yrange=[-0.1,1.45])
>>> of.plot(d1='x',d2='z',overplot=True,color='k')
>>> plot(orbs.x(),orbs.z(),'r.',alpha=0.1)
>>> plot(orbts.x(),orbts.z(),'b.',alpha=0.1)
>>> sdf.plotTrack(d1='x',d2='z',interp=True,color='r',overplot=True,lw=1.)
>>> sdft.plotTrack(d1='x',d2='z',interp=True,color='b',overplot=True,lw=1.)

This gives then

We see that the track from streamdf agrees very well with the location of the points sampled from streamspraydf.

The streamspraydf sample function can also return the points at the time of stripping, that is, not integrated to the present time (when using integrate=False); this can be useful for visualizing where stars get stripped from the progenitor. When initializing streamspraydf, you can also specify a different potential for computing the tidal radius and velocity distribution of the tidal debris, which can be useful when the overall potential contains pieces that are irrelevant for computing the tidal radius and that don’t allow the tidal radius to be computed (using the rtpot= option). If you want to generate a stream around a moving object, for example, a stream created within a satellite galaxy of the Milky Way, you can specify an orbit for the center of the satellite (center=) and the stream will be generated around this center rather than around the center of the total potential (this was used in Qian et al. 2022); the center orbit is integrated in centerpot, which can also differ from the potential that stream stars are integrated in (e.g., the stream stars may feel the potential from the satellite itself and/or the satellite could be experiencing dynamical friction which the stream stars do not feel).