10.2. Distribution functions for razor-thin disks¶
Now that we have discussed the Jeans equations that relate the moments of steady-state distribution functions for disks, we present some examples of full steady-state distribution functions in this section. We start by describing distribution functions for razor-thin disks. These distribution functions are all functions \(f(E,L_z)\) of the planar part of the energy, which we will simply write as \(E\) here, and the \(z\) component of the angular momentum \(L_z\). As discussed at the beginning of this chapter, there is an implied multiplication of \(f(E,L_z)\) with \(\delta(I_3)\) that we will not carry around, but that causes the distribution function to be zero for any non-zero \((z,v_z)\). Because the distribution function is zero for non-zero \((z,v_z)\), we will simply work in the mid-plane in this section and drop the vertical dependence of any quantity, e.g., the potential will be \(\Phi(R) \equiv \Phi(R,z=0)\).
When we discussed steady-state distribution functions for spherical systems, we spend much time on self-consistent models where the gravitational potential is sourced by the density obtained by integrating the distribution function over velocity—that is, the distribution function’s implied density satisfies the Poisson equation with the total gravitational potential. As we saw in Chapter 8, galactic disks are embedded in massive dark-matter halos and at any point in the disk a significant fraction of the gravitational force is contributed by the dark matter. Moreover, stellar disks like the Milky Way are made up of an almost continuous sequence of stellar populations (e.g., of different ages) that taken together generate the gravitational potential of the stellar disk. We will find it much easier to create distribution function models for simpler components of the full stellar disk (e.g., a single age population) than for the full stellar disk at once. Thus, when dealing with disks, self-consistency is not as important a consideration when discussing steady-state disk distribution functions and most of the models that we will discuss are not self-consistent. In particular, in most examples below, we will write down equilibrium disk models for exponential disks embedded in a gravitational potential with a flat rotation curve, which as we know from Chapter 7.3.3.2 is not a self-consistent setup, but it is a realistic description of the dynamics of stellar disks.
10.2.1. The distribution function of a cold, rotating disk¶
We now turn to the question of how we can determine a razor-thin distribution function for a rotating disk with a given surface density \(\Sigma(R)\) in a given axisymmetric potential \(\Phi(R)\), which is not self-consistently generated by \(\Sigma(R)\). We start by building a disk distribution function for a cold disk with zero velocity dispersion and then warm it up by introducing non-zero velocity dispersion in Section 10.2.2 below. These distribution functions are all functions of the \(E\) and \(L_z\). The derivation in this and the next section follows that of Dehnen (1999b).
From the Stromberg asymmetric drift Equation (10.8) it is clear that an axisymmetric disk with vanishing velocity dispersion \(\overline{v_R^2} = 0\) must be rotating at the circular velocity at each radius. Thus, we can build the full disk out of circular orbits and the form of the distribution function must be \begin{equation}\label{eq-cold-df-form} f(E,L_z) = F(L_z)\,\delta(E-E_c[L_z])\,, \end{equation} where \(E_c[L_z]\) is the energy of a circular orbit with angular momentum \(L_z\) and \(F(L_z)\) is a function that we need to determine by matching \(\Sigma(R)\). Because we only need to consider orbits that are infinitesimally close to a circular orbit in this distribution function, we can work out the energy in the epicycle approximation. The energy difference in the epicycle approximation is given by \begin{align} E - E_c[L_z] & = {v_R^2 \over 2} + \Phi_{\mathrm{eff}}(R;L_z) - \Phi_{\mathrm{eff}}(R_g;L_z)\\ & \approx \frac{v_R^2}{2} + \frac{\kappa^2\,X^2\,\cos^2(\kappa t + \alpha)}{2}\,,\label{eq-epicycle-EEc} \end{align} where \(X\) is the amplitude of the epicyclic motion and \(\alpha\) is the phase (see Equations 9.10 and 9.19). In Chapter 9.3.2 we demonstrated that \(v_\phi - v_c = -2B\,X\,\cos(\kappa t + \alpha)\) where \(v_c\) is the circular velocity at \(R\) (up to terms that are second order in \(X\); Equation 9.31). substituting that \(-4B = \kappa^2/\Omega\), we can therefore also write the energy difference as \begin{align} E - E_c[L_z] & = \frac{v_R^2}{2} + \frac{4\,\Omega^2\,(v_\phi - v_c)^2}{2\,\kappa^2}= \frac{1}{2}\,\left[v_R^2 + \gamma^2\,(v_\phi - v_c)^2\right]\,,\label{eq-EEc-epicycle} \end{align} where we have substituted \(\gamma = 2\Omega/\kappa\). Thus, we can write the distribution function in Equation (10.10) in terms of \((R,\phi,v_R,v_\phi)\) as \begin{equation} f(R,\phi,v_R,v_\phi) = 2\,F(R\,v_\phi)\,\delta\left(\frac{1}{2}\,\left[v_R^2 + \gamma^2\,(v_\phi - v_c)^2\right]\right)\,. \end{equation} The factor of two accounts for the \(E < E_c[L_z]\) part of the distribution function which we have dropped when going to the epicycle approximation (this is necessary to keep the delta function properly normalized). The surface density corresponding to this distribution function is \begin{align} \Sigma'(R) & = 2\,\int_{-\infty}^\infty\mathrm{d}v_\phi\,\int_{-\infty}^\infty\mathrm{d}v_R\,F(R\,v_\phi)\,\delta\left(\frac{1}{2}\,\left[v_R^2 + \gamma^2\,(v_\phi - v_c)^2\right]\right)\nonumber\\ & = \frac{2}{\gamma}\, \int_{0}^{2\pi}\mathrm{d}\psi\,\int_{0}^\infty\mathrm{d}v\,v\,F\left(R\,\left[\frac{v\,\sin\psi}{\gamma}+v_c\right]\right)\,\delta\left(\frac{v^2}{2}\right)\nonumber\\ & = \frac{2\pi}{\gamma}\,F(R\,v_c)\,. \end{align} In this calculation, we introduced polar velocity coordinates \((v,\psi)\) with \(v_R = v\,\cos\psi\) and \(\gamma(v_\phi-v_c) = v\sin \psi\) to perform the integral over the delta function and we used that \(\int_0^\infty \mathrm{d}x\,f(x)\,\delta(x) = f(0)/2\) (Equation B.30). For \(\Sigma'(R)\) to equal \(\Sigma(R)\), we therefore need to choose the function \(F(L_z)\) as \(F(L_z) = \gamma\Sigma(R_L)/(2\pi)\), where \(R_L = R_g(L_z)\) is the guiding-center radius of a circular orbit with angular momentum \(L_z\) (see Equation 9.7; while we generally denote the guiding-center radius as \(R_g\), here we use \(R_L\) to clearly distinguish the radius definition’s dependence on \(L_z\) from an alternative radius definition in terms of \(E\) in the Dehnen distribution function below, where we will denote the equivalent radius as \(R_E\)). Note that in this expression, \(\gamma\) in general also depends on \(R_L\). The distribution function for a cold disk with surface density \(\Sigma(R)\) is therefore \begin{equation}\label{eq-colddisk-df} f(E,L_z) = \frac{\gamma(R_L)\,\Sigma(R_L)}{2\pi}\,\delta(E-E_c[L_z])\,. \end{equation}
10.2.2. Distribution functions for warm, rotating disks¶
Now that we have derived a general expression for the distribution function of a cold, rotating disk, we ask how we can generalize it to have a given velocity dispersion profile \(\sigma_R(R)\). For a distribution function that depends on \((E,L_z)\) there is no unique solution to this problem. Just like in the spherical case, the Jeans equations do not close even if we derive higher and higher-order versions of it. Specifying \(\Sigma(R)\) and \(\sigma_R(R)\) still leaves freedom in the higher-order moments of the velocity distribution that is not fixed by the assumption of a steady state. To model the distribution function of a population of stars in a disk, we therefore need to use additional assumptions about their velocity distribution to specify a unique model.
Below we discuss three examples of how to generalize the cold disk distribution function in Equation (10.16). All three of these examples generalize the expression by replacing the delta function by a function \(K[\cdot]\) that has a finite width \(\sigma_R^2(R)\): \(\delta(E-E_c[L_z]) \rightarrow K[(E-E_c[L_z])/\sigma_R^2(R)]\): whereas for the cold distribution function, the only populated orbits at a given angular momentum \(L_z\) are the circular orbits with \(E = E_c[L_z]\), in these warm distribution functions, orbits at a given \(L_z\) are populated according to \(K[(E-E_c[L_z])/\sigma_R^2(R)]\). This general approach of creating warm disk distribution functions by warming up the cold distribution function is the main way in which even more complex and realistic disk distribution functions than we discuss in this chapter are obtained (e.g. Binney 2010; Binney & McMillan 2011).
From the derivation of the cold-disk distribution function, it should be clear that simply replacing the delta function with a finite-width kernel \(K[\cdot]\) will not produce \(\Sigma'(R) = \Sigma(R)\) using the prefactor in Equation (10.16): the derivation of this prefactor depends on the kernel being a delta function, and only if both \(\Sigma(R)\) and \(\sigma_R(R)\) are independent of \(R\) would \(\Sigma'(R) = \Sigma(R)\). Similarly, if we compute the actual velocity dispersion \(\sigma'_R(R)\) in the radial direction as the square root of the second moment in \(v_R\) of the distribution function, we find that \(\sigma'_R \neq \sigma_R\). We therefore have to either (i) derive a new prefactor \(F(L_z)\) such that \(\Sigma'(R) = \Sigma(R)\) or (ii) live with the fact that we only have \(\Sigma'(R) \approx \Sigma(R)\) without having exact equality between the two. We will largely take path (ii) here and not worry too much about exactly matching the “target” \(\Sigma(R)\) and \(\sigma_R(R)\) and the distribution function’s actual \(\Sigma'(R)\) and \(\sigma'_R(R)\).
10.2.2.1. The Schwarzschild distribution function¶
The first and simplest example of a warmed-up disk distribution function replaces the delta function in Equation (10.16) with an exponential function as in Equation (10.20) below. The Schwarzschild distribution function results from approximating the argument of the exponential using the epicycle approximation in Equation (10.13) and is given by \begin{equation} f(E,L_z) = \frac{\gamma(R_L)\,\Sigma(R_L)}{2\pi\,\sigma_R^2(R_L)}\,\exp\left(-\frac{v_R^2 + \gamma^2\,(v_\phi - v_c)^2}{2\sigma_R^2(R_L)}\right)\,. \end{equation} This is a distribution function that is only approximately in a steady state, because the energy difference \(E-E_c[L_z]\) is only equal to \([v_R^2 + \gamma^2\,(v_\phi - v_c)^2]/2\) for nearly circular orbits. The surface density \(\Sigma'(R)\) corresponding to the Schwarzschild distribution function is \begin{align} \Sigma'(R)& = \int_{-\infty}^{\infty}\mathrm{d}v_\phi\,\frac{\gamma(R_L)\,\Sigma(R_L)}{\sqrt{2\pi}\,\sigma_R(R_L)}\,\exp\left(-\frac{\gamma^2\,(v_\phi - v_c)^2}{2\sigma_R^2(R_L)}\right)\,.\label{eq-schwarzschilddf-Sigmap} \end{align} Because of the \(R_L\) dependence of \(\sigma_R\) in the integrand, for general \(\sigma_R(R)\) this integral needs to be calculated numerically. However, when \(\sigma_R/v_c \ll 1\) then we can treat all functions that depend on \(R_L\) as being constant over the range where the integrand is significantly non-zero and we find that \(\Sigma'(R) = \Sigma(R)\). Thus, for \(\sigma_R/v_c \ll 1\) where only close-to-circular orbits are populated significantly, we have that \(\Sigma'(R)\approx\Sigma(R)\).
Following a similar argument, we have that for close-to-circular orbits \(\overline{v_\phi} \approx v_c(R)\), \(\sigma'_R(R) \approx \sigma_R(R)\), and \(\sigma'_\phi(R) \approx \sigma'_R(R)/\gamma\) (in agreement with \(\sigma'_\phi(R)/\sigma'_R(R) = \sqrt{B/(B-A)}\) in the epicycle approximation). In fact, the entire velocity distribution implied by the Schwarzschild distribution function for \(\sigma_R/v_c \ll 1\) is nearly Gaussian. When the velocity dispersion \(\sigma_R(R)\) is larger, the \(R_L = L_z/v_c\) dependence of the \(\Sigma\), \(\gamma\), and \(\sigma_R\) functions in the distribution function cause the distribution to no longer be Gaussian.
As an example of the Schwarzschild distribution function, we consider the family defined by \begin{equation}\label{eq-schwarzschilddf-examplefamily} \Sigma(R) = \Sigma_0\,e^{-R/3}\,;\quad \sigma_R(R) = \sigma_0\,e^{-(R-1)}\,, \end{equation} for different values of the velocity dispersion parameter \(\sigma_0\) at \(R=1\). That is, both the surface density profile \(\Sigma(R)\) and the velocity-dispersion profile \(\sigma_R(R)\) are exponential, the former with a scale length of \(1/3\) and the latter with a scale length of \(1\) (this is a not-unreasonable model for stars near the Sun). As the background potential, we use a potential with a flat rotation curve: \(\Phi(R) = v_c^2\,\ln R\) in the plane. First we numerically compute \(\Sigma'(R)\) and compare it to \(\Sigma(R)\) for a model with a relatively high \(\sigma_0/v_c = 0.15\) in the left panel of Figure 10.1 and we do the same for \(\sigma'_R(R)\) in the right panel, comparing it to \(\sigma_R(R)\).
[9]:
from galpy.df import schwarzschilddf, dehnendf
df_warm= schwarzschilddf(profileParams=(1./3.,1.0,0.15),beta=0.)
Rs= numpy.linspace(0.1,2.2,101)
figure(figsize=(12,4))
subplot(1,2,1)
sigmaprime= numpy.array([df_warm.surfacemass(R) for R in Rs])
line_sprime_schw= semilogy(Rs,sigmaprime,zorder=1,
label=r"$\Sigma_{\mathrm{Schwarzschild}}'$")
sigma= numpy.array([df_warm.targetSurfacemass(R) for R in Rs])
line_sdens= semilogy(Rs,sigma,zorder=0,
label=r"$\Sigma$")
subplot(1,2,2)
sigma2prime= numpy.array([df_warm.sigma2(R) for R in Rs])
line_s2prime_schw= plot(Rs,numpy.sqrt(sigma2prime),zorder=1,
label=r"$\sigma_{R',\mathrm{Schwarzschild}}/v_c$")
sigma2= numpy.array([df_warm.targetSigma2(R) for R in Rs])
line_sigma2= plot(Rs,numpy.sqrt(sigma2),zorder=0,
label=r"$\sigma_R/v_c$")
# Also add Dehnen
df_warm_dehnen= dehnendf(profileParams=(1./3.,1.0,0.15),beta=0.)
subplot(1,2,1)
sigmaprime= numpy.array([df_warm_dehnen.surfacemass(R) for R in Rs])
line_sprime_dehnen= semilogy(Rs,sigmaprime,zorder=1,
label=r"$\Sigma_\mathrm{Dehnen}'$")
xlabel(r'$R$')
xlim(0.,2.3)
legend(handles=[line_sprime_schw[0],line_sprime_dehnen[0],line_sdens[0]],
fontsize=18.,loc='upper right',frameon=False)
subplot(1,2,2)
sigma2prime= numpy.array([df_warm_dehnen.sigma2(R) for R in Rs])
line_s2prime_dehnen= plot(Rs,numpy.sqrt(sigma2prime),zorder=1,
label=r"$\sigma_{R,\mathrm{Dehnen}}'/v_c$")
xlabel(r'$R$')
xlim(0.,2.3)
ylim(0.,0.4)
legend(handles=[line_s2prime_schw[0],line_s2prime_dehnen[0],line_sigma2[0]],
fontsize=18.,loc='upper right',frameon=False);
Figure 10.1: Actual (primed) and target (unprimed) surface density (left) and velocity dispersion (right) profiles in Schwarzschild and Dehnen distribution functions.
We see that \(\Sigma'(R)\) is indeed approximately equal to \(\Sigma(R)\). The agreement between \(\sigma'_R(R)\) and \(\sigma_R(R)\) is less good and \(\sigma'_R(R)\) is higher by about 10% than \(\sigma_R(R)\) over much of the range in \(R\). Nevertheless, the overall shape and amplitude of both \(\Sigma'(R)\) and \(\sigma'_R(R)\) are in relatively good agreement with \(\Sigma(R)\) and \(\sigma_R(R)\).
The distribution of azimuthal velocities for different values of \(\sigma_0\) at \(R=1\) is given in Figure 10.2 (with arbitrary overall normalization of each distribution).
[10]:
from galpy.df import shudf, dehnendf, schwarzschilddf
df_cold_shu= shudf(profileParams=(1./3.,1.0,0.02),beta=0.)
df_warmer_shu= shudf(profileParams=(1./3.,1.0,0.07),beta=0.)
df_warm_shu= shudf(profileParams=(1./3.,1.0,0.15),beta=0.)
df_cold_dehnen= dehnendf(profileParams=(1./3.,1.0,0.02),beta=0.)
df_warmer_dehnen= dehnendf(profileParams=(1./3.,1.0,0.07),beta=0.)
df_warm_dehnen= dehnendf(profileParams=(1./3.,1.0,0.15),beta=0.)
df_cold_schwarz= schwarzschilddf(profileParams=(1./3.,1.0,0.02),beta=0.)
df_warmer_schwarz= schwarzschilddf(profileParams=(1./3.,1.0,0.07),beta=0.)
df_warm_schwarz= schwarzschilddf(profileParams=(1./3.,1.0,0.15),beta=0.)
vts= numpy.linspace(0.4,1.6,201)
figure(figsize=(7.75,4))
# Cold
ycold= numpy.array([df_cold_dehnen(numpy.array([1.,0.,vt])) for vt in vts])
norm= (numpy.amax(ycold)+numpy.random.uniform()/10.)
line_cold= plot(vts,ycold/norm,
label=r'$\sigma_0 = 0.02\,v_c$')
ycold= numpy.array([df_cold_shu(numpy.array([1.,0.,vt])) for vt in vts])
plot(vts,ycold/norm,color=line_cold[0].get_color(),ls='--',zorder=0)
ycold= numpy.array([df_cold_schwarz(numpy.array([1.,0.,vt])) for vt in vts])
plot(vts,ycold/norm,color=line_cold[0].get_color(),ls=':',zorder=0)
# Warmer
ywarmer= numpy.array([df_warmer_dehnen(numpy.array([1.,0.,vt])) for vt in vts])
norm= (numpy.amax(ywarmer)+numpy.random.uniform()/10.)
line_warmer= plot(vts,ywarmer/norm,
label=r'$\sigma_0 = 0.07\,v_c$')
ywarmer= numpy.array([df_warmer_shu(numpy.array([1.,0.,vt])) for vt in vts])
plot(vts,ywarmer/norm,color=line_warmer[0].get_color(),ls='--',zorder=0)
ywarmer= numpy.array([df_warmer_schwarz(numpy.array([1.,0.,vt])) for vt in vts])
plot(vts,ywarmer/norm,color=line_cold[0].get_color(),ls=':',zorder=0)
# Warm
ywarm= numpy.array([df_warm_dehnen(numpy.array([1.,0.,vt])) for vt in vts])
norm= (numpy.amax(ywarm)+numpy.random.uniform()/10.)
line_warm= plot(vts,ywarm/norm,
label=r'$\sigma_0 = 0.15\,v_c$')
ywarm= numpy.array([df_warm_shu(numpy.array([1.,0.,vt])) for vt in vts])
plot(vts,ywarm/norm,color=line_warm[0].get_color(),ls='--',zorder=0)
ywarm= numpy.array([df_warm_schwarz(numpy.array([1.,0.,vt])) for vt in vts])
plot(vts,ywarm/norm,color=line_cold[0].get_color(),ls=':',zorder=0)
xlabel(r'$v_\phi/v_c$')
# Schwarzschild and Shu legend entries
import matplotlib.lines as mlines
dehnen_line= mlines.Line2D([], [], color='k',marker=None,ls='-',
label=r'$\mathrm{Dehnen}$')
shu_line= mlines.Line2D([], [], color='k',marker=None,ls='--',
label=r'$\mathrm{Shu}$')
schwarzschild_line= mlines.Line2D([], [], color='k',marker=None,ls=':',
label=r'$\mathrm{Schwarzschild}$')
gca().add_artist(legend(handles=[dehnen_line,shu_line,schwarzschild_line],
loc='upper right',fontsize=18.,frameon=False))
legend(handles=[line_cold[0],line_warmer[0],line_warm[0]],
fontsize=18.,loc='upper left',frameon=False);
Figure 10.2: Tangential velocity distributions in Dehnen, Shu, and Schwarzschild distribution functions with different velocity dispersions.
The coldest distribution function has \(\sigma_0/v_c = 0.02\) and therefore consists of stars only on nearly circular orbits. The distribution of azimuthal velocities is strongly peaked at \(v_\phi = v_c\) and the distribution is close to symmetric and Gaussian. When we increase \(\sigma_0/v_c\), we see that the distribution becomes wider, because of the higher velocity dispersion, but also becomes asymmetric. Especially in the case where \(\sigma_0/v_c = 0.15\) the velocity distribution is highly asymmetric. This is another manifestation of the asymmetric drift. The mean azimuthal velocity \(\overline{v_\phi} < v_c\) for these examples, because both \(\Sigma(R)\) and \(\sigma_R(R)\) are declining functions of \(R\) (see Equation 10.8). Physically what happens is that at a given radius \(R\), stars with \(v_\phi > v_c\) must have guiding-center radii \(R_g > R\), because of the conservation of angular momentum: \(R\,v_\phi = R_g\,v_c\) or \(v_\phi/v_c = R_g/R > 1\) (remember that we are working in a disk with a flat rotation curve, \(v_c = \mathrm{constant}\)). Similarly, stars with \(v_\phi < v_c\) must have \(R_g < R\). Now, because the surface density declines with radius, there are more stars with \(v_\phi < v_c\) at \(R\) than there are stars with \(v_\phi > v_c\) and the average \(\overline{v_\phi} < v_c\). The decline of the velocity dispersion increases this effect: that \(\sigma_R(R_g < R) > \sigma_R(R_g > R)\) means that stars from \(R_g < R\) have \(|R_g-R|\) that is typically larger than \(|R_g-R|\) from stars with \(R_g > R\), and thus, that \(|v_\phi-v_c|\) is larger for stars coming from \(R_g < R\) than it is for stars coming from \(R_g > R\), further pulling the mean \(\overline{v_\phi}\) to smaller values. Thus, there are more stars coming from inner radii and these are typically coming from further inwards than there are stars coming from outer radii, leading to \(\overline{v_\phi} < v_c\) and skewing the \(v_\phi\) distribution.
Lest your take-away from this discussion is that \(\overline{v_\phi} < v_c\) always holds, it is important to keep in mind that it does not have to be the case that \(\overline{v_\phi} < v_c\); the sign of the \((\overline{v_\phi} - v_c)\) difference is set by the gradient of \(\Sigma'\) and \(\sigma'_R\) (see Equation 10.8). While the overall stellar distribution in a disk galaxy is decreasing outwards and thus for the disk as a whole we expect \(\overline{v_\phi} < v_c\), sub-populations of stars (e.g., young stars) might have surface density profiles that increase outwards. In that case, the asymmetric drift could work the other way around and produce \(\overline{v_\phi} > v_c\) when measuring \(\overline{v_\phi}\) for that population of stars. As an example, we show the distribution of \(v_\phi\) for the same model as above, but using \(\Sigma(R) = \Sigma_0\,e^{R_0/3}\), in Figure 10.3 (and keeping the slowly declining \(\sigma_R\)).
[11]:
from galpy.df import shudf, dehnendf, schwarzschilddf
df_cold_shu= shudf(profileParams=(-1./3.,1.0,0.02),beta=0.)
df_warmer_shu= shudf(profileParams=(-1./3.,1.0,0.07),beta=0.)
df_warm_shu= shudf(profileParams=(-1./3.,1.0,0.15),beta=0.)
df_cold_dehnen= dehnendf(profileParams=(-1./3.,1.0,0.02),beta=0.)
df_warmer_dehnen= dehnendf(profileParams=(-1./3.,1.0,0.07),beta=0.)
df_warm_dehnen= dehnendf(profileParams=(-1./3.,1.0,0.15),beta=0.)
df_cold_schwarz= schwarzschilddf(profileParams=(-1./3.,1.0,0.02),beta=0.)
df_warmer_schwarz= schwarzschilddf(profileParams=(-1./3.,1.0,0.07),beta=0.)
df_warm_schwarz= schwarzschilddf(profileParams=(-1./3.,1.0,0.15),beta=0.)
vts= numpy.linspace(0.4,1.6,201)
figure(figsize=(7.75,4))
# Cold
ycold= numpy.array([df_cold_dehnen(numpy.array([1.,0.,vt])) for vt in vts])
norm= (numpy.amax(ycold)+numpy.random.uniform()/10.)
line_cold= plot(vts,ycold/norm,
label=r'$\sigma_0 = 0.02\,v_c$')
ycold= numpy.array([df_cold_shu(numpy.array([1.,0.,vt])) for vt in vts])
plot(vts,ycold/norm,color=line_cold[0].get_color(),ls='--',zorder=0)
ycold= numpy.array([df_cold_schwarz(numpy.array([1.,0.,vt])) for vt in vts])
plot(vts,ycold/norm,color=line_cold[0].get_color(),ls=':',zorder=0)
# Warmer
ywarmer= numpy.array([df_warmer_dehnen(numpy.array([1.,0.,vt])) for vt in vts])
norm= (numpy.amax(ywarmer)+numpy.random.uniform()/10.)
line_warmer= plot(vts,ywarmer/norm,
label=r'$\sigma_0 = 0.07\,v_c$')
ywarmer= numpy.array([df_warmer_shu(numpy.array([1.,0.,vt])) for vt in vts])
plot(vts,ywarmer/norm,color=line_warmer[0].get_color(),ls='--',zorder=0)
ywarmer= numpy.array([df_warmer_schwarz(numpy.array([1.,0.,vt])) for vt in vts])
plot(vts,ywarmer/norm,color=line_cold[0].get_color(),ls=':',zorder=0)
# Warm
ywarm= numpy.array([df_warm_dehnen(numpy.array([1.,0.,vt])) for vt in vts])
norm= (numpy.amax(ywarm)+numpy.random.uniform()/10.)
line_warm= plot(vts,ywarm/norm,
label=r'$\sigma_0 = 0.15\,v_c$')
ywarm= numpy.array([df_warm_shu(numpy.array([1.,0.,vt])) for vt in vts])
plot(vts,ywarm/norm,color=line_warm[0].get_color(),ls='--',zorder=0)
ywarm= numpy.array([df_warm_schwarz(numpy.array([1.,0.,vt])) for vt in vts])
plot(vts,ywarm/norm,color=line_cold[0].get_color(),ls=':',zorder=0)
xlabel(r'$v_\phi/v_c$')
# Schwarzschild and Shu legend entries
import matplotlib.lines as mlines
dehnen_line= mlines.Line2D([], [], color='k',marker=None,ls='-',
label=r'$\mathrm{Dehnen}$')
shu_line= mlines.Line2D([], [], color='k',marker=None,ls='--',
label=r'$\mathrm{Shu}$')
schwarzschild_line= mlines.Line2D([], [], color='k',marker=None,ls=':',
label=r'$\mathrm{Schwarzschild}$')
gca().add_artist(legend(handles=[dehnen_line,shu_line,schwarzschild_line],
loc='upper right',fontsize=18.,frameon=False))
legend(handles=[line_cold[0],line_warmer[0],line_warm[0]],
fontsize=18.,loc='upper left',frameon=False);
Figure 10.3: Like Figure 10.2, but for distribution functions with increasing surface density with \(R\).
Because the velocity dispersion still decreases with increasing \(R\), there is now a competition between the two effects discussed in the previous paragraph. But on the whole, the surface density wins out and the velocity distribution has a mean \(\overline{v_\phi} > v_c\) in this case. This effect does in fact happen for low-metallicity stars in the solar neighborhood, which have negative asymmetric drifts.
10.2.2.2. The Shu distribution function¶
The Schwarzschild distribution function is only approximately a steady-state distribution function, because it computes the energy in the epicycle approximation. A true steady-state razor-thin disk distribution function can be obtained simply by computing the energy without using any approximation: \begin{equation}\label{eq-df_shu} f(E,L_z) = \frac{\gamma(R_L)\,\Sigma(R_L)}{2\pi\,\sigma_R^2(R_L)}\exp\left(-\frac{E-E_c[L_z]}{\sigma_R^2(R_L)}\right)\,. \end{equation} This is the Shu distribution function (Shu 1969). As a function of \((R,\phi,v_R,v_\phi)\), we can write this as \begin{equation} f(R,\phi,v_R,v_\phi) = \frac{\gamma(R_L)\,\Sigma(R_L)}{2\pi\,\sigma_R^2(R_L)}\exp\left(-\frac{1}{2}\frac{v_R^2+v_\phi^2-v^2_c(R_L)+2\,[\Phi(R)-\Phi(R_L)]}{\sigma_R^2(R_L)}\right)\,. \end{equation} From this form, it is clear that the distribution of \(v_R\) at \((R,\phi,v_\phi)\) is again a Gaussian with dispersion \(\sigma_R(R_L)\). We can calculate \(\Sigma'(R)\) and \(\sigma'_R(R)\) using similar expressions to that in Equation (10.18); the integration again needs to be performed numerically.
The Shu distribution function reduces to the Schwarzschild distribution function in the limit that \(\sigma_R \rightarrow 0\), but for larger \(\sigma_R\) it deviates from the Schwarzschild distribution. For the same family of \(\Sigma(R)\) and \(\sigma_R(R)\) in Equations (10.19), we compare the distribution of \(v_\phi\) for the two families in Figure 10.2. It is clear for small velocity dispersion that the two distribution functions give essentially identical velocity distributions, but for the model with the largest \(\sigma_0\), the Schwarzschild distribution function gives a slightly wider velocity distribution and one that is slightly shifted toward smaller \(v_\phi\). The difference in the mean \(\overline{v_\phi}\) is nevertheless only about \(v_c/100\), so the difference is quite small. Of course in locations where the velocity dispersion is larger (e.g., \(R<1\) in the above example) the difference will be larger.
10.2.2.3. The Dehnen distribution function¶
Dehnen (1999b) introduced a new disk distribution function by starting from a different form of the cold-disk distribution function. We can write the cold-disk distribution function alternatively as \begin{equation}\label{eq-cold-df-formalt} f(E,L_z) = F(E)\,\delta\left(\Omega(E)\,[L_c(E)-L_z]\right)\,, \end{equation} where \(L_c(E)\) is the angular momentum of a circular orbit with energy \(E\) and \(\Omega(E)\) is the azimuthal frequency of that orbit. The reason that this is a better starting point than Equation (10.10) is that the argument of the delta function is a better approximation of \(\kappa\,J_R\) (where \(J_R\) is the radial action) than \(E-E_c(L_z)\) (Dehnen 1999a). The combination \(\kappa\,J_R\) is in many ways the natural quantity to use when warming up a cold disk (Dehnen 1999b). Following a similar calculation of \(\Sigma'(R)\) as in Section 10.2.1 above, we find that we can match a desired surface density \(\Sigma(R)\) by setting \(F(E) = \gamma(R_E)\Sigma(R_E)/(2\pi)\), where now \(R_E\) is the radius of a circular orbit with energy \(E\). An alternative form for the cold-disk distribution function in Equation (10.16) is therefore \begin{equation}\label{eq-colddisk-df-alt} f(E,L_z) = \frac{\gamma(R_E)\,\Sigma(R_E)}{2\pi}\,\delta\left(\Omega(E)\,[L_c(E)-L_z]\right)\,. \end{equation} Equation (10.23) is the same distribution function as the cold-disk distribution function in Equation (10.16).
We can now build a disk distribution function with non-zero velocity dispersion by replacing the delta-function kernel in Equation (10.23) with an exponential, just like how we created the Schwarzschild and Shu distribution functions starting from Equation (10.16). This gives the Dehnen disk distribution function \begin{equation}\label{eq-df_dehnen} f(E,L_z) = \frac{\gamma(R_E)\,\Sigma(R_E)}{2\pi\,\sigma_R^2(R_E)}\exp\left(-\frac{\Omega(E)\,[L_c(E)-L_z]}{\sigma_R^2(R_E)}\right)\,. \end{equation} Dehnen (1999b) argues that this distribution function more faithfully models the process through which stellar populations in a disk attain a non-zero velocity dispersion than the Shu distribution function does.
Let’s compare the Dehnen, Shu, and Schwarzschild distribution functions for the same model family \(\Sigma(R)\) and \(\sigma_R(R)\) as in Equations (10.19). The distribution of \(v_\phi\) at \(R=1\) is shown in Figure 10.2. We see that all three distribution functions give essentially the same distribution of \(v_\phi\) for small \(\sigma_0\), but for the largest \(\sigma_0\), the Dehnen distribution function gives a narrower distribution than the other models. The difference between the Shu and Dehnen distribution functions is mostly due to the fact that the Dehnen distribution function has a smaller velocity dispersion \(\sigma'_R\) than the Shu distribution function for the same input \(\sigma_R\).
The Dehnen model produces a better match between \(\Sigma'(R)\) and \(\Sigma(R)\) and between \(\sigma'_R(R)\) and \(\sigma_R(R)\). For example, for the warmest Dehnen model in the example above, Figure 10.1 compares \(\Sigma(R)\) to \(\Sigma'(R)\) and \(\sigma_R(R)\) to \(\sigma'_R(R)\). We see again that the match is very good and much better than in the case of the Schwarzschild distribution function. Thus, the Dehnen distribution function more faithfully represents the desired surface density and velocity dispersion profile and this is typically the case over a range of models for \(\Sigma(R)\) and \(\sigma_R(R)\) and for different gravitational-potential models. Dehnen (1999b) also provides an algorithm that can be applied to all warm disk distribution functions discussed here to iteratively change the input surface density and velocity dispersion profiles to attain a better match between \(\Sigma'(R)\) and \(\Sigma(R)\) and between \(\sigma'_R(R)\) and \(\sigma_R(R)\).
A procedure for sampling phase-space points from the Dehnen and Shu distribution functions can be found in Dehnen (1999b).
.
.
.