Table Of Contents

Previous topic

Action-angle coordinates

Next topic

Orbit

This Page

Three-dimensional disk distribution functions

galpy contains a fully three-dimensional disk distribution: galpy.df.quasiisothermaldf, which is an approximately isothermal distribution function expressed in terms of action–angle variables (see 2010MNRAS.401.2318B and 2011MNRAS.413.1889B). Recent research shows that this distribution function provides a good model for the DF of mono-abundance sub-populations (MAPs) of the Milky Way disk (see 2013MNRAS.434..652T and 2013ApJ...779..115B). This distribution function family requires action-angle coordinates to evaluate the DF, so galpy.df.quasiisothermaldf makes heavy use of the routines in galpy.actionAngle (in particular those in galpy.actionAngleAdiabatic and galpy.actionAngle.actionAngleStaeckel).

Setting up the DF and basic properties

The quasi-isothermal DF is defined by a gravitational potential and a set of parameters describing the radial surface-density profile and the radial and vertical velocity dispersion as a function of radius. In addition, we have to provide an instance of a galpy.actionAngle class to calculate the actions for a given position and velocity. For example, for a galpy.potential.MWPotential potential using the adiabatic approximation for the actions, we import and define the following

>>> from galpy.potential import MWPotential
>>> from galpy.actionAngle import actionAngleAdiabatic
>>> from galpy.df import quasiisothermaldf
>>> aA= actionAngleAdiabatic(pot=MWPotential,c=True)

and then setup the quasiisothermaldf instance

>>> qdf= quasiisothermaldf(1./3.,0.2,0.1,1.,1.,pot=MWPotential,aA=aA,cutcounter=True)

which sets up a DF instance with a radial scale length of R_0/3, a local radial and vertical velocity dispersion of 0.2\,V_c(R_0) and 0.1\,V_c(R_0), respectively, and a radial scale lengths of the velocity dispersions squared of R_0. cutcounter=True specifies that counter-rotating stars are explicitly excluded (normally these are just exponentially suppressed). As for the two-dimensional disk DFs, these parameters are merely input (or target) parameters; the true density and velocity dispersion profiles calculated by evaluating the relevant moments of the DF (see below) are not exactly exponential and have scale lengths and local normalizations that deviate slightly from these input parameters. We can estimate the DF’s actual radial scale length near R_0 as

>>> qdf.estimate_hr(1.)
0.33843243662586048

which is quite close to the input value of 1/3. Similarly, we can estimate the scale lengths of the dispersions squared

>>> qdf.estimate_hsr(1.)
1.1527209864858059
>>> qdf.estimate_hsz(1.)
1.0441867587783933

The vertical profile is fully specified by the velocity dispersions and radial density / dispersion profiles under the assumption of dynamical equilibrium. We can estimate the scale height of this DF at a given radius and height as follows

>>> qdf.estimate_hz(1.,0.125)
0.018715154050080292

Near the mid-plane this vertical scale height becomes very large because the vertical profile flattens, e.g.,

>>> qdf.estimate_hz(1.,0.125/100.)
0.85435378664432149

or even

>>> qdf.estimate_hz(1.,0.)
128674.27506772846

which is basically infinity.

Evaluating moments

We can evaluate various moments of the DF giving the density, mean velocities, and velocity dispersions. For example, the mean radial velocity is again everywhere zero because the potential and the DF are axisymmetric

>>> qdf.meanvR(1.,0.)
0.0

Likewise, the mean vertical velocity is everywhere zero

>>> qdf.meanvz(1.,0.)
0.0

The mean rotational velocity has a more interesting dependence on position. Near the plane, this is the same as that calculated for a similar two-dimensional disk DF (see Evaluating moments of the DF)

>>> qdf.meanvT(1.,0.)
0.9150884078276913

However, this value decreases as one moves further from the plane. The quasiisothermaldf allows us to calculate the average rotational velocity as a function of height above the plane. For example,

>>> zs= numpy.linspace(0.,0.25,21)
>>> mvts= numpy.array([qdf.meanvT(1.,z) for z in zs])

which gives

>>> plot(zs,mvts)
_images/qdf-meanvtz.png

We can also calculate the second moments of the DF. We can check whether the radial and velocity dispersions at R_0 are close to their input values

>>> numpy.sqrt(qdf.sigmaR2(1.,0.))
0.20918647082092351
>>> numpy.sqrt(qdf.sigmaz2(1.,0.))
0.092564222527283468

and they are pretty close. We can also calculate the mixed R and z moment, for example,

>>> qdf.sigmaRz(1.,0.125)
0.0

or expressed as an angle (the tilt of the velocity ellipsoid)

>>> qdf.tilt(1.,0.125)
0.0

This tilt is zero because we are using the adiabatic approximation. As this approximation assumes that the motions in the plane are decoupled from the vertical motions of stars, the mixed moment is zero. However, this approximation is invalid for stars that go far above the plane. By using the Staeckel approximation to calculate the actions, we can model this coupling better. Setting up a quasiisothermaldf instance with the Staeckel approximation

>>> from galpy.actionAngle import actionAngleStaeckel
>>> aAS= actionAngleStaeckel(pot=MWPotential,delta=0.55,c=True)
>>> qdfS= quasiisothermaldf(1./3.,0.2,0.1,1.,1.,pot=MWPotential,aA=aAS,cutcounter=True)

we can similarly calculate the tilt

>>> qdfS.tilt(1.,0.125)
5.4669442080366721

or about 5 degrees. As a function of height, we find

>>> tilts= numpy.array([qdfS.tilt(1.,z) for z in zs])
>>> plot(zs,tilts)

which gives

_images/qdf_tiltz.png

We can also calculate the density and surface density (the zero-th velocity moments). For example, the vertical density

>>> densz= numpy.array([qdf.density(1.,z) for z in zs])

and

>>> denszS= numpy.array([qdfS.density(1.,z) for z in zs])

We can compare the vertical profiles calculated using the adiabatic and Staeckel action-angle approximations

>>> semilogy(zs,densz/densz[0])
>>> semilogy(zs,denszS/denszS[0])

which gives

_images/qdf-densz.png

Similarly, we can calculate the radial profile of the surface density

>>> rs= numpy.linspace(0.5,1.5,21)
>>> surfr= numpy.array([qdf.surfacemass_z(r) for r in rs])
>>> surfrS= numpy.array([qdfS.surfacemass_z(r) for r in rs])

and compare them with each other and an exponential with scale length 1/3

>>> semilogy(rs,surfr/surfr[10])
>>> semilogy(rs,surfrS/surfrS[10])
>>> semilogy(rs,numpy.exp(-(rs-1.)/(1./3.)))

which gives

_images/qdf-densr.png

The two radial profiles are almost indistinguishable and are very close, if somewhat shallower, than the pure exponential profile.

General velocity moments, including all higher order moments, are implemented in quasiisothermaldf.vmomentdensity.

Evaluating and sampling the full probability distribution function

We can evaluate the distribution itself by calling the object, e.g.,

>>> qdf(1.,0.1,1.1,0.1,0.) #input: R,vR,vT,z,vz
array([ 10.16445158])

or as a function of rotational velocity, for example in the mid-plane

>>> vts= numpy.linspace(0.,1.5,101)
>>> pvt= numpy.array([qdfS(1.,0.,vt,0.,0.) for vt in vts])
>>> plot(vts,pvt/numpy.sum(pvt)/(vts[1]-vts[0]))

which gives

_images/qdf-callvt.png

This is, however, not the true distribution of rotational velocities at R =0 and z =0, because it is conditioned on zero radial and vertical velocities. We can calculate the distribution of rotational velocities marginalized over the radial and vertical velocities as

>>> qdfS.pvT(1.,1.,0.) #input vT,R,z
15.464330302557528

or as a function of rotational velocity

>>> pvt= numpy.array([qdfS.pvT(vt,1.,0.) for vt in vts])

overplotting this over the previous distribution gives

_images/qdf-pvt.png

which is slightly different from the conditioned distribution. Similarly, we can calculate marginalized velocity probabilities `pvR, pvz, pvRvT, pvRvz, and pvTvz. These are all multiplied with the density, such that marginalizing these over the remaining velocity component results in the density.

We can sample velocities at a given location using quasiisothermaldf.sampleV (there is currently no support for sampling locations from the density profile, although that is rather trivial):

>>> vs= qdfS.sampleV(1.,0.,n=10000)
>>> hist(vs[:,1],normed=True,histtype='step',bins=101,range=[0.,1.5])

gives

_images/qdf-pvtwsamples.png

which shows very good agreement with the green (marginalized over vR and vz) curve (as it should).