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.MWPotential2014 potential using the adiabatic approximation for the actions, we import and define the following

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

and then setup the quasiisothermaldf instance

>>> qdf= quasiisothermaldf(1./3.,0.2,0.1,1.,1.,pot=MWPotential2014,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 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.32908034635647182

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

>>> qdf.estimate_hsr(1.)
# 1.1913935820372923
>>> qdf.estimate_hsz(1.)
# 1.0506918075359255

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.021389597757156088

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.)
# 1.006386030587223

or even

>>> qdf.estimate_hz(1.,0.)
# 221852.19839263527

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.91988346380781227

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)

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.20807112565801389
>>> numpy.sqrt(qdf.sigmaz2(1.,0.))
# 0.090453510526130904

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=MWPotential2014,delta=0.45,c=True)
>>> qdfS= quasiisothermaldf(1./3.,0.2,0.1,1.,1.,pot=MWPotential2014,aA=aAS,cutcounter=True)

we can similarly calculate the tilt

>>> qdfS.tilt(1.,0.125)
# 0.10314272868452541

or about 5 degrees (the returned value has units of rad). As a function of height, we find

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

which gives

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

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

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([ 16.86790643])

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

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
# 58.708924787596544

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

>>> plot(vts,pvt/numpy.sum(pvt)/(vts[1]-vts[0]))

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],density=True,histtype='step',bins=101,range=[0.,1.5])

gives

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