Two-dimensional disk distribution functions
=============================================
galpy contains various disk distribution functions, both in two and
three dimensions. This section introduces the two-dimensional
distribution functions, useful for studying the dynamics of stars that
stay relatively close to the mid-plane of a galaxy. The vertical
motions of these stars may be approximated as being entirely decoupled
from the motion in the plane.
.. _usagediskdfs:
Types of disk distribution functions
---------------------------------------
galpy contains the following distribution functions for razor-thin
disks: ``galpy.df.dehnendf``, ``galpy.df.shudf``, and
``galpy.df.schwarzschilddf``. These are the distribution functions of
Dehnen (`1999AJ....118.1201D
`_), Shu
(`1969ApJ...158..505S
`_), and
Schwarzschild (the Shu DF in the epicycle approximation, see Binney &
Tremaine 2008). Everything shown below for ``dehnendf`` can also be
done for ``shudf`` and ``schwarzschilddf``. The Schwarzschild DF is
primarily included as an educational tool; it is *not* a true
steady-state DF, because it uses the approximate energy from the
epicycle approximation rather than the true energy, and is fully
superseded by the Shu DF, which *is* a good steady-state DF.
These disk distribution functions are functions of the energy and the
angular momentum alone. They can be evaluated for orbits, or for a
given energy and angular momentum. At this point, only power-law
rotation curves are supported. A ``dehnendf`` instance is initialized
as follows
>>> from galpy.df import dehnendf
>>> dfc= dehnendf(beta=0.)
This initializes a ``dehnendf`` instance based on an exponential
surface-mass profile with scale-length 1/3 and an exponential
radial-velocity-dispersion profile with scale-length 1 and a value of
0.2 at R=1. Different parameters for these profiles can be provided as
an initialization keyword. For example,
>>> dfc= dehnendf(beta=0.,profileParams=(1./4.,1.,0.2))
initializes the distribution function with a radial scale length of
1/4 instead.
We can show that these distribution functions have an asymmetric drift
built-in by evaluating the DF at R=1. We first create a set of
orbit-instances and then evaluate the DF at them
>>> from galpy.orbit import Orbit
>>> os= [Orbit([1.,0.,1.+-0.9+1.8/1000*ii]) for ii in range(1001)]
>>> dfro= [dfc(o) for o in os]
>>> plot([1.+-0.9+1.8/1000*ii for ii in range(1001)],dfro)
.. image:: images/basic-df-asydrift.png
Or we can plot the two-dimensional density at R=1.
>>> dfro= [[dfc(Orbit([1.,-0.7+1.4/200*jj,1.-0.6+1.2/200*ii])) for jj in range(201)]for ii in range(201)]
>>> dfro= numpy.array(dfro)
>>> from galpy.util.plot import dens2d
>>> dens2d(dfro,origin='lower',cmap='gist_yarg',contours=True,xrange=[-0.7,0.7],yrange=[0.4,1.6],xlabel=r'$v_R$',ylabel=r'$v_T$')
.. image:: images/basic-df-2d.png
.. _dftwod-moments:
Evaluating moments of the DF
-----------------------------
galpy can evaluate various moments of the disk distribution
functions. For example, we can calculate the mean velocities (for the
DF with a scale length of 1/3 above)
>>> dfc.meanvT(1.)
# 0.91715276979447324
>>> dfc.meanvR(1.)
# 0.0
and the velocity dispersions
>>> numpy.sqrt(dfc.sigmaR2(1.))
# 0.19321086259083936
>>> numpy.sqrt(dfc.sigmaT2(1.))
# 0.15084122011271159
and their ratio
>>> dfc.sigmaR2(1.)/dfc.sigmaT2(1.)
# 1.6406766813028968
In the limit of zero velocity dispersion (the epicycle approximation)
this ratio should be equal to 2, which we can check as follows
>>> dfccold= dehnendf(beta=0.,profileParams=(1./3.,1.,0.02))
>>> dfccold.sigmaR2(1.)/dfccold.sigmaT2(1.)
# 1.9947493895454664
We can also calculate higher order moments
>>> dfc.skewvT(1.)
# -0.48617143862047763
>>> dfc.kurtosisvT(1.)
# 0.13338978593181494
>>> dfc.kurtosisvR(1.)
# -0.12159407676394096
We already saw above that the velocity dispersion at R=1 is not
exactly equal to the input velocity dispersion (0.19321086259083936
vs. 0.2). Similarly, the whole surface-density and velocity-dispersion
profiles are not quite equal to the exponential input profiles. We can
calculate the resulting surface-mass density profile using
``surfacemass``, ``sigmaR2``, and ``sigma2surfacemass``. The latter
calculates the product of the velocity dispersion squared and the
surface-mass density. E.g.,
>>> dfc.surfacemass(1.)
# 0.050820867101511534
We can plot the surface-mass density as follows
>>> Rs= numpy.linspace(0.01,5.,151)
>>> out= [dfc.surfacemass(r) for r in Rs]
>>> plot(Rs, out)
.. image:: images/diskdf-surfacemass.png
or
>>> plot(Rs,numpy.log(out))
.. image:: images/diskdf-logsurfacemass.png
which shows the exponential behavior expected for an exponential
disk. We can compare this to the input surface-mass density
>>> input_out= [dfc.targetSurfacemass(r) for r in Rs]
>>> plot(Rs,numpy.log(input_out)-numpy.log(out))
.. image:: images/diskdf-surfacemassdiff.png
which shows that there are significant differences between the desired
surface-mass density and the actual surface-mass density. We can do
the same for the velocity-dispersion profile
>>> out= [dfc.sigmaR2(r) for r in Rs]
>>> input_out= [dfc.targetSigma2(r) for r in Rs]
>>> plot(Rs,numpy.log(input_out)-numpy.log(out))
.. image:: images/diskdf-sigma2diff.png
That the input surface-density and velocity-dispersion profiles are
not the same as the output profiles, means that estimates of DF
properties based on these profiles will not be quite
correct. Obviously this is the case for the surface-density and
velocity-dispersion profiles themselves, which have to be explicitly
calculated by integration over the DF rather than by evaluating the
input profiles. This also means that estimates of the asymmetric drift
based on the input profiles will be wrong. We can calculate the
asymmetric drift at R=1 using the asymmetric drift equation derived
from the Jeans equation (eq. 4.228 in Binney & Tremaine 2008), using
the input surface-density and velocity dispersion profiles
>>> dfc.asymmetricdrift(1.)
# 0.090000000000000024
which should be equal to the circular velocity minus the mean rotational
velocity
>>> 1.-dfc.meanvT(1.)
# 0.082847230205526756
These are not the same in part because of the difference between the
input and output surface-density and velocity-dispersion profiles (and
because the ``asymmetricdrift`` method assumes that the ratio of the
velocity dispersions squared is two using the epicycle approximation;
see above).
Using corrected disk distribution functions
-----------------------------------------------
As shown above, for a given surface-mass density and velocity
dispersion profile, the two-dimensional disk distribution functions
only do a poor job of reproducing the desired profiles. We can correct
this by calculating a set of *corrections* to the input profiles such
that the output profiles more closely resemble the desired profiles
(see `1999AJ....118.1201D
`_). galpy supports
the calculation of these corrections, and comes with some
pre-calculated corrections (these can be found `here
`_). For
example, the following initializes a ``dehnendf`` with corrections up
to 20th order (the default)
>>> dfc= dehnendf(beta=0.,correct=True)
The following figure shows the difference between the actual
surface-mass density profile and the desired profile for 1, 2, 3, 4,
5, 10, 15, and 20 iterations
.. image:: images/testSurfacemassCorrections_sigma0_0.5.png
and the same for the velocity-dispersion profile
.. image:: images/testSigmaCorrections_sigma0_0.5.png
galpy will automatically save any new corrections that you calculate.
All of the methods for an uncorrected disk DF can be used for the
corrected DFs as well. For example, the velocity dispersion is now
>>> numpy.sqrt(dfc.sigmaR2(1.))
# 0.19999985069451526
and the mean rotation velocity is
>>> dfc.meanvT(1.)
# 0.90355161181498711
and (correct) asymmetric drift
>>> 1.-dfc.meanvT(1.)
# 0.09644838818501289
That this still does not agree with the simple ``dfc.asymmetricdrift``
estimate is because of the latter's using the epicycle approximation
for the ratio of the velocity dispersions.
Oort constants and functions
------------------------------
galpy also contains methods to calculate the Oort functions for
two-dimensional disk distribution functions. These are known as the
*Oort constants* when measured in the solar neighborhood. They are
combinations of the mean velocities and derivatives thereof. galpy
calculates these by direct integration over the DF and derivatives of
the DF. Thus, we can calculate
>>> dfc= dehnendf(beta=0.)
>>> dfc.oortA(1.)
# 0.43190780889218749
>>> dfc.oortB(1.)
# -0.48524496090228575
The *K* and *C* Oort constants are zero for axisymmetric DFs
>>> dfc.oortC(1.)
# 0.0
>>> dfc.oortK(1.)
# 0.0
In the epicycle approximation, for a flat rotation curve *A* =- *B* =
0.5. The explicit calculates of *A* and *B* for warm DFs quantify how
good (or bad) this approximation is
>>> dfc.oortA(1.)+dfc.oortB(1.)
# -0.053337152010098254
For the cold DF from above the approximation is much better
>>> dfccold= dehnendf(beta=0.,profileParams=(1./3.,1.,0.02))
>>> dfccold.oortA(1.), dfccold.oortB(1.)
# (0.49917556666144003, -0.49992824742490816)
Sampling data from the DF
--------------------------
We can sample from the disk distribution functions using
``sample``. ``sample`` can return either an energy--angular-momentum
pair, or a full orbit initialization. We can sample 4000 orbits for
example as (could take two minutes)
>>> o= dfc.sample(n=4000,returnOrbit=True,nphi=1)
We can then plot the histogram of the sampled radii and compare it to the input surface-mass density profile
>>> Rs= [e.R() for e in o]
>>> hists, bins, edges= hist(Rs,range=[0,2],normed=True,bins=30)
>>> xs= numpy.array([(bins[ii+1]+bins[ii])/2. for ii in range(len(bins)-1)])
>>> plot(xs, xs*exp(-xs*3.)*9.,'r-')
E.g.,
.. image:: images/basic-df-sampleR.png
We can also plot the spatial distribution of the sampled disk
>>> xs= [e.x() for e in o]
>>> ys= [e.y() for e in o]
>>> figure()
>>> plot(xs,ys,',')
E.g.,
.. image:: images/basic-df-samplexy.png
We can also sample points in a specific radial range (might take a few
minutes)
>>> o= dfc.sample(n=1000,returnOrbit=True,nphi=1,rrange=[0.8,1.2])
and we can plot the distribution of tangential velocities
>>> vTs= [e.vxvv[2] for e in o]
>>> hists, bins, edges= hist(vTs,range=[.5,1.5],normed=True,bins=30)
>>> xs= numpy.array([(bins[ii+1]+bins[ii])/2. for ii in range(len(bins)-1)])
>>> dfro= [dfc(Orbit([1.,0.,x]))/9./numpy.exp(-3.) for x in xs]
>>> plot(xs,dfro,'r-')
.. image:: images/basic-df-samplevT.png
The agreement between the sampled distribution and the theoretical
curve is not as good because the sampled distribution has a finite
radial range. If we sample 10,000 points in ``rrange=[0.95,1.05]`` the
agreement is better (this takes a long time):
.. image:: images/basic-df-samplevTmore.png
We can also directly sample velocities at a given radius rather than
in a range of radii. Doing this for a correct DF gives
>>> dfc= dehnendf(beta=0.,correct=True)
>>> vrvt= dfc.sampleVRVT(1.,n=10000)
>>> hists, bins, edges= hist(vrvt[:,1],range=[.5,1.5],normed=True,bins=101)
>>> xs= numpy.array([(bins[ii+1]+bins[ii])/2. for ii in range(len(bins)-1)])
>>> dfro= [dfc(Orbit([1.,0.,x])) for x in xs]
>>> plot(xs,dfro/numpy.sum(dfro)/(xs[1]-xs[0]),'r-')
.. image:: images/basic-df-samplevTatR.png
galpy further has support for sampling along a given line of sight in
the disk, which is useful for interpreting surveys consisting of a
finite number of pointings. For example, we can sampled distances
along a given line of sight
>>> ds= dfc.sampledSurfacemassLOS(30./180.*numpy.pi,n=10000)
which is very fast. We can histogram these
>>> hists, bins, edges= hist(ds,range=[0.,3.5],normed=True,bins=101)
and compare it to the predicted distribution, which we can calculate as
>>> xs= numpy.array([(bins[ii+1]+bins[ii])/2. for ii in range(len(bins)-1)])
>>> fd= numpy.array([dfc.surfacemassLOS(d,30.) for d in xs])
>>> plot(xs,fd/numpy.sum(fd)/(xs[1]-xs[0]),'r-')
which shows very good agreement with the sampled distances
.. image:: images/basic-df-sampled.png
galpy can further sample full 4D phase--space coordinates along a
given line of sight through ``dfc.sampleLOS``.
Non-axisymmetric, time-dependent disk distribution functions
-------------------------------------------------------------
``galpy`` also supports the evaluation of non-axisymmetric,
time-dependent two-dimensional DFs. These specific DFs are constructed
by assuming an initial axisymmetric steady state, described by a DF of
the family discussed above, that is then acted upon by a
non-axisymmetric, time-dependent perturbation. The DF at a given time
and phase-space position is evaluated by integrating the orbit
backwards in time in the non-axisymmetric potential until the time of
the initial DF is reached. From Liouville's theorem, which states that
phase-space volume is conserved along the orbit, we then know that we
can evaluate the non-axisymmetric DF today as the initial DF at the
initial point on the orbit. This procedure was first used by `Dehnen
(2000) `_.
This is implemented in ``galpy`` as ``galpy.df.evolveddiskdf``. Such a
DF is setup by specifying the initial DF, the non-axisymmetric
potential, and the time of the initial state. For example, we can look
at the effect of an elliptical perturbation to the potential like that
described by `Kuijken & Tremaine
`_. To do this, we
set up an elliptical perturbation to a logarithmic potential that is
grown slowly to minimize non-adiabatic effects
>>> from galpy.potential import LogarithmicHaloPotential, EllipticalDiskPotential
>>> lp= LogarithmicHaloPotential(normalize=1.)
>>> ep= EllipticalDiskPotential(twophio=0.05,phib=0.,p=0.,tform=-150.,tsteady=125.)
This perturbation starts to be grown at ``tform=-150`` over a time
period of ``tsteady=125`` time units. We will consider the effect of
this perturbation on a very cold disk (velocity dispersion
:math:`\sigma_R = 0.0125\,v_c`) and a warm disk (:math:`\sigma_R =
0.15\,v_c`). We set up these two initial DFs
>>> idfcold= dehnendf(beta=0.,profileParams=(1./3.,1.,0.0125))
>>> idfwarm= dehnendf(beta=0.,profileParams=(1./3.,1.,0.15))
and then set up the ``evolveddiskdf``
>>> from galpy.df import evolveddiskdf
>>> edfcold= evolveddiskdf(idfcold,[lp,ep],to=-150.)
>>> edfwarm= evolveddiskdf(idfwarm,[lp,ep],to=-150.)
where we specify that the initial state is at ``to=-150``.
We can now use these ``evolveddiskdf`` instances in much the same way
as ``diskdf`` instances. One difference is that there is much more
support for evaluating the DF on a grid (to help speed up the rather
slow computations involved). Thus, we can evaluate the mean radial
velocity at ``R=0.9``, ``phi=22.5`` degree, and ``t=0`` by using a grid
>>> mvrcold, gridcold= edfcold.meanvR(0.9,phi=22.5,deg=True,t=0.,grid=True,returnGrid=True,gridpoints=51,nsigma=6.)
>>> mvrwarm, gridwarm= edfwarm.meanvR(0.9,phi=22.5,deg=True,t=0.,grid=True,returnGrid=True,gridpoints=51)
>>> print(mvrcold, mvrwarm)
# -0.0358753028951 -0.0294763627935
The cold response agrees well with the analytical calculation, which
predicts that this is :math:`-0.05/\sqrt{2}`:
>>> print(mvrcold+0.05/sqrt(2.))
# -0.000519963835811
The warm response is slightly smaller in amplitude
>>> print(mvrwarm/mvrcold)
# 0.821633837619
although the numerical uncertainty in ``mvrwarm`` is large, because
the grid is not sufficiently fine.
We can then re-use this grid in calculations of other moments of
the DF, e.g.,
>>> print(edfcold.meanvT(0.9,phi=22.5,deg=True,t=0.,grid=gridcold))
# 0.965058551359
>>> print(edfwarm.meanvT(0.9,phi=22.5,deg=True,t=0.,grid=gridwarm))
# 0.915397094614
which returns the mean rotational velocity, and
>>> print(edfcold.vertexdev(0.9,phi=22.5,deg=True,t=0.,grid=gridcold))
# 0.0560531474616
>>> print(edfwarm.vertexdev(0.9,phi=22.5,deg=True,t=0.,grid=gridwarm))
# 0.0739164830253
which gives the vertex deviation in rad. The reason we have to
calculate the grid out to ``6nsigma`` for the cold response is that
the response is much bigger than the velocity dispersion of the
population. This velocity dispersion is used to automatically to set
the grid edges, but sometimes has to be adjusted to contain the full
DF.
``evolveddiskdf`` can also calculate the Oort functions, by directly
calculating the spatial derivatives of the DF. These can also be calculated on a grid, such that we can do
>>> oortacold, gridcold, gridrcold, gridphicold= edfcold.oortA(0.9,phi=22.5,deg=True,t=0.,returnGrids=True,gridpoints=51,derivGridpoints=51,grid=True,derivphiGrid=True,derivRGrid=True,nsigma=6.)
>>> oortawarm, gridwarm, gridrwarm, gridphiwarm= edfwarm.oortA(0.9,phi=22.5,deg=True,t=0.,returnGrids=True,gridpoints=51,derivGridpoints=51,grid=True,derivphiGrid=True,derivRGrid=True)
>>> print(oortacold, oortawarm)
# 0.575494559999 0.526389833249
It is clear that these are quite different. The cold calculation is
again close to the analytical prediction, which says that :math:`A =
A_{\mathrm{axi}}+0.05/(2\sqrt{2})` where :math:`A_{\mathrm{axi}} =
1/(2\times0.9)` in this case:
>>> print(oortacold-(0.5/0.9+0.05/2./sqrt(2.)))
# 0.0022613349141670236
These grids can then be re-used for the other Oort functions, for
example,
>>> print(edfcold.oortB(0.9,phi=22.5,deg=True,t=0.,grid=gridcold,derivphiGrid=gridphicold,derivRGrid=gridrcold))
# -0.574674310521
>>> print(edfwarm.oortB(0.9,phi=22.5,deg=True,t=0.,grid=gridwarm,derivphiGrid=gridphiwarm,derivRGrid=gridrwarm))
# -0.555546911144
and similar for ``oortC`` and ``oortK``. These warm results should
again be considered for illustration only, as the grid is not
sufficiently fine to have a small numerical error.
The grids that have been calculated can also be plotted to show the
full velocity DF. For example,
>>> gridcold.plot()
gives
.. image:: images/basic-df-nonaxi-cold.png
which demonstrates that the DF is basically the initial DF that has been displaced (by a significant amount compared to the velocity dispersion). The warm velocityd distribution is given by
>>> gridwarm.plot()
which returns
.. image:: images/basic-df-nonaxi-warm.png
The shift of the smooth DF here is much smaller than the velocity
dispersion.
Example: The Hercules stream in the Solar neighborhood as a result of the Galactic bar
---------------------------------------------------------------------------------------
.. _hercules:
We can combine the orbit integration capabilities of galpy with the
provided distribution functions and see the effect of the Galactic bar
on stellar velocities. By backward integrating orbits starting at the
Solar position in a potential that includes the Galactic bar we can
evaluate what the velocity distribution is that we should see today if
the Galactic bar stirred up a steady-state disk. For this we
initialize a flat rotation curve potential and Dehnen's bar potential
>>> from galpy.potential import LogarithmicHaloPotential, DehnenBarPotential
>>> lp= LogarithmicHaloPotential(normalize=1.)
>>> dp= DehnenBarPotential()
The Dehnen bar potential is initialized to start bar formation four bar
periods before the present day and to have completely formed the bar two
bar periods ago. We can integrate back to the time before
bar-formation:
>>> ts= numpy.linspace(0,dp.tform(),1000)
where ``dp.tform()`` is the time of bar-formation (in the usual
time-coordinates).
We initialize orbits on a grid in velocity space and integrate them
>>> ins= Orbit(numpy.array([[[1.,-0.7+1.4/100*jj,1.-0.6+1.2/100*ii,0.] for jj in range(101)] for ii in range(101)]))
>>> ins.integrate(ts,[lp,dp])
We can then evaluate the weight of these orbits by assuming that the
disk was in a steady-state before bar-formation with a Dehnen
distribution function. We evaluate the Dehnen distribution function at
``dp.tform()`` for each of the orbits (evaluating the distribution
function only works for an Orbit with a single object, so we need to
unpack the Orbit instance that contains all orbits)
>>> dfc= dehnendf(beta=0.,correct=True)
>>> out= [[dfc(o(dp.tform())) for o in j] for j in ins]
>>> out= numpy.array(out)
This gives
>>> from galpy.util.plot import dens2d
>>> dens2d(out,origin='lower',cmap='gist_yarg',contours=True,xrange=[-0.7,0.7],yrange=[0.4,1.6],xlabel=r'$v_R$',ylabel=r'$v_T$')
.. image:: images/diskdf-dehnenhercules.png
Now that ``galpy`` contains the ``evolveddiskdf`` described above,
this whole calculation is encapsulated in this module and can be done
much more easily as
>>> edf= evolveddiskdf(dfc,[lp,dp],to=dp.tform())
>>> mvr, grid= edf.meanvR(1.,grid=True,gridpoints=101,returnGrid=True)
The gridded DF can be accessed as ``grid.df``, which we can plot as before
>>> dens2d(grid.df.T,origin='lower',cmap='gist_yarg',contours=True,xrange=[grid.vRgrid[0],grid.vRgrid[-1]],yrange=[grid.vTgrid[0],grid.vTgrid[-1]],xlabel=r'$v_R$',ylabel=r'$v_T$')
.. image:: images/diskdf-dehnenhercules-edf.png
For more information see `2000AJ....119..800D
`_ and
`2010ApJ...725.1676B
`_. Note that the
x-axis in the Figure above is defined as minus the x-axis in these
papers.