7.3. Gravitational potentials from disk density distributions¶
So far we have discussed models for disk density distribution that derive from simple models of the potential. While we have been successful in obtaining some simple models that can be used as a rough approximation to the mass distribution of a disk, all of the models discussed above have some issues. The density of the Miyamoto-Nagai potential falls of much slower than that of the exponential disks that we observe. The flattened logarithmic potential that can describe the overall potential of a disk plus halo system has significant regions with negative densities as soon as the flattening becomes more than mild and can therefore not be a physically realistic model. When we flatten the density instead of the potential (e.g., using Equation 7.17) we avoid the issue of negative densities, but this comes at the cost of having to solve the Poisson equation for flattened densities. We will tackle that problem in this section for highly flattened systems and come back to this problem for mildly flattened systems in Chapter 12.
7.3.1. Some general considerations¶
We start with some general considerations about solving for the potential for a highly flattened density, which we will represent as \(\rho(R,\phi,z) = \Sigma(R,\phi|z)\,\zeta(z)\), with \(\Sigma(R,\phi|z)\) the surface density of a layer at height \(z\) and \(\zeta(z)\) defined such that \(\int \mathrm{d}z\,\zeta(z) = 1\). For such a density, the Poisson equation is best written down in cylindrical coordinates. Using the expression for the Laplacian in cylindrical coordinates in Equation (A.14), the Poisson equation becomes \begin{equation}\label{eq-poisson-cylindrical} \frac{1}{R}\,\frac{\partial}{\partial R}\left(R\,\frac{\partial \Phi(R,\phi,z)}{\partial R}\right) + \frac{1}{R^2}\,\frac{\partial^2 \Phi(R,\phi,z)}{\partial \phi^2}+\frac{\partial^2\Phi(R,\phi,z)}{\partial z^2} = 4\pi\,G\,\Sigma(R,\phi|z)\,\zeta(z)\,. \end{equation} To make headway in solving this equation, we introduce the Fourier transform \(f_m(\phi)\) over \(\phi\) of a function \(f(\phi)\) \begin{equation}\label{eq-fourier-phi} f_m = \frac{1}{2\pi}\,\int_0^{2\pi}\mathrm{d}\phi\,f(\phi)\,e^{-im\,\phi}\,, \end{equation} with the inverse transformation given by (see Equation B.25) \begin{equation}\label{eq-fourier-phi-inv} f(\phi) = \sum_{m=-\infty}^{\infty}\,f_m\,e^{im\,\phi}\,. \end{equation} Taking the Fourier transform over \(\phi\) of Equation (7.21) leads to \begin{equation}\label{eq-poisson-cylindrical-m} \frac{1}{R}\,\frac{\partial}{\partial R}\left(R\,\frac{\partial \Phi_m(R,z)}{\partial R}\right) - \frac{m^2}{R^2}\,\Phi_m(R,z)+\frac{\partial^2\Phi_m(R,z)}{\partial z^2} = 4\pi\,G\,\Sigma_m(R|z)\,\zeta(z)\,, \end{equation} where every function with subscript \(m\) is defined by Equation (7.22). Expanding out the first term and multiplying both sides with \(R^2\) gives \begin{equation}\label{eq-poisson-cylindrical-intermed1} R^2\,\frac{\partial^2 \Phi_m(R,z)}{\partial R^2} +R\,\frac{\partial \Phi_m(R,z)}{\partial R} - m^2\,\Phi_m(R,z)+R^2\,\frac{\partial^2\Phi_m(R,z)}{\partial z^2} = 4\pi\,G\,R^2\,\Sigma_m(R|z)\,\zeta(z)\,, \end{equation} Compare this equation to the differential equation that defines Bessel functions of the first kind \(J_m(x)\) (Equation B.45), which we can re-write as \begin{equation}\label{eq-bessel-def} x^2\frac{\mathrm{d}^2 J_m(x)}{\mathrm{d} x^2} + x\,\frac{\mathrm{d}J_m(x)}{\mathrm{d}x}-m^2\,J_m(x) = -x^2\,J_m(x)\, \end{equation} The left-hand side of this equation is equal to the first three terms in Equation (7.25). Therefore, we can expect that solutions to Equation (7.24) will be given as a sum over terms of the form \(\tilde{\Phi}_m(k,z)\,J_m(k\,R)\), because for such a form, the left-hand side of Equation (7.25) becomes that of an ordinary differential equation \begin{equation} -k^2\,\tilde{\Phi}_m(k,z)\,J_m(k\,R)+J_m(k\,R)\,\frac{\partial^2\tilde{\Phi}_m(k,z)}{\partial z^2} = \ldots\,. \end{equation}
We will continue with this solution in Section 7.3.5, where we will see that \(\tilde{\Phi}_m(k,z)\) is the Hankel transform of order \(m\) and we can solve the above equation in terms of the Hankel transform of \(\Sigma_m(R|z)\,\zeta(z)\). But we will first consider the easier case of an axisymmetric, razor-thin disk and generalize this to a finite-thickness disk before coming back to the general solution.
7.3.2. Potentials for razor-thin disks¶
As we already discussed in Section 7.2.1, a razor-thin disk is a disk with density \(\rho(R,\phi,z) = \Sigma(R,\phi)\,\delta(z)\). Such a disk is characterized by its surface density \(\Sigma(R,\phi)\). We start by considering an axisymmetric razor-thin disk and follow the approach introduced by Toomre (1963) for determining the potential of a disk with surface density \(\Sigma(R)\).
The approach we follow is this. Based on the considerations above, we guess a class of potential–density pairs for razor-thin disks. Then it will turn out that the set of densities forms a complete basis for all smooth \(\Sigma(R)\) mass distributions and that there is a simple way of decomposing any \(\Sigma(R)\) into a sum over the density basis functions. The potential is then obtained as the same sum over the set of potentials corresponding to the density basis functions.
For an axisymmetric disk only the \(m=0\) component of Equation (7.24) is non-trivial and we have that \(\Sigma_0(R) = \Sigma(R)\) and therefore \(\Phi_0(R,z) = \Phi(R,z)\). For a razor-thin disk, this equation becomes \begin{equation} \frac{1}{R}\,\frac{\partial}{\partial R}\left(R\,\frac{\partial \Phi(R,z)}{\partial R}\right) +\frac{\partial^2\Phi(R,z)}{\partial z^2} = 4\pi\,G\,\Sigma(R)\,\delta(z)\,. \end{equation} It is easy to check that the function \begin{equation}\label{eq-pot-razorthin-toomrebasissolution} \tilde{\Phi}(R;k) = -2\pi\,G\,e^{-k|z|}\,J_0(kR)\,, \end{equation} with \(k\) some constant, is a solution of this equation, with the corresponding \(\tilde{\Sigma}(R;k)\) given by (e.g., Equation 7.9) \begin{equation} \tilde{\Sigma}(R;k) = k\,J_0(kR)\,. \end{equation} Therefore, if we can write the surface density \(\Sigma(R)\) of a razor-thin disk that we are interested in as a sum over contributions \(S_0(k)\,\tilde{\Sigma}(R;k)\) with amplitudes \(S_0(k)\), then the potential is the sum over contributions \(S_0(k)\,\tilde{\Phi}(R;k)\), because of the linearity of the Poisson equation. Thus, we write \begin{equation}\label{eq-invhankel-surfdens-0} \Sigma(R) = \int_0^\infty \mathrm{d}k\,J_0(kR)\,k\,S_0(k)\,, \end{equation} where we can determine the \(S_0(k)\) as \begin{equation}\label{eq-hankel-surfdens-0} S_0(k) = \int_0^\infty \mathrm{d} R' \,J_0(kR')\,R'\,\Sigma(R')\,, \end{equation} because of the Fourier-Bessel theorem (Equation B.50). \(S_0(k)\) is the zero-th order Hankel transform of \(\Sigma(R)\). The gravitational potential is then given by \begin{align}\label{eq-pot-razorthin} \Phi(R,z) & = -2\pi\,G\,\int_0^\infty\mathrm{d}k\,e^{-k|z|}\,J_0(kR)\,S_0(k)\\ & = -2\pi\,G\,\int_0^\infty\mathrm{d}k\,e^{-k|z|}\,J_0(kR)\,\int_0^\infty \mathrm{d} R' \,J_0(kR')\,R'\,\Sigma(R')\,.\nonumber \end{align}
Equation (7.33) is the general solution for the potential of a razor-thin, axisymmetric disk. The gravitational potential in the \(z=0\) mid-plane is \begin{align}\label{eq-pot-razorthin-midplane} \Phi(R,0) & = -2\pi\,G\,\int_0^\infty\mathrm{d}k\,\,J_0(kR)\,\int_0^\infty \mathrm{d} R' \,J_0(kR')\,R'\,\Sigma(R')\,. \end{align} The circular velocity in the mid-plane can be obtained from this, noting that \(\mathrm{d}\,J_0(t)/\mathrm{d}\,t = -J_1(t)\) (Equation B.48) \begin{align}\label{eq-vc-razorthin} v_c^2(R) & = 2\pi\,G\,R\,\int_0^\infty\mathrm{d}k\,\,J_1(kR)\,k\,\int_0^\infty \mathrm{d} R' \,J_0(kR')\,R'\,\Sigma(R')\,. \end{align} When we use the expressions above to compute the potential, rotation curve, and radial and vertical forces for a given surface density profile with numerical integration, care must be taken because the Bessel functions are highly oscillatory. In the next section we consider examples where we can do all integrals analytically, and so we do not have to worry about this issue. The form of the integral for the gravitational potential in Equation (7.33) is that of the Laplace transformation of \(J_0(kR)\,S_0(k)\), so razor-thin disk surface densities for which the potential has a closed form are those with known Laplace transformations.
Before continuing with a few examples, we consider the limit \(R,z \rightarrow \infty\) for a razor-thin mass distribution with a finite mass \(M\). Physically, we expect that the potential \(\Phi(R,z)\) should approach that of a point-mass with mass \(M\), because at distances \(r \gg\) the scale of the disk, the details of how the mass \(M\) is distributed should not matter. We can write Equation (7.33) using that \(\mathrm{d} M(<R') = 2\pi\,R'\Sigma(R')\,\mathrm{d}R'\) and using integration by parts as \begin{align} \Phi(R,z) & = -\,G\,\int_0^\infty\mathrm{d}k\,e^{-k|z|}\,J_0(kR)\,\int_0^\infty \mathrm{d} R' \,J_1(kR')\,k\,M(< R')\,, \end{align} because \(M(0) = 0\) and \(M\,J_0(\infty) = 0\). Swapping the order of integration and using \(J_1(kR')\,k = -\mathrm{d} J_0(kR') / \mathrm{d}\,R'\) we find that \begin{align} \Phi(R,z) & = G\,\int_0^\infty \mathrm{d} R'\,M(< R')\frac{\mathrm{d}}{\mathrm{d}R'}\left[\int_0^\infty\mathrm{d}k\,e^{-k|z|}\,J_0(kR) \,J_0(kR')\right]\label{eq-razorthin-larger-deriv0}\,. \end{align} To solve the integral in square brackets, we need the Laplace transformation of a product of two Bessel functions of the first kind of order zero, and the result is (e.g., Kausel & Irfan Baig 2012) \begin{equation}\label{eq-J0J0z-int} \int_0^\infty\mathrm{d}k\,e^{-k|z|}\,J_0(kR) \,J_0(kR') = \frac{2}{\pi\,\sqrt{(R+R')^2+z^2}}\,K\left(\frac{4RR'}{(R+R')^2+z^2}\right)\,, \end{equation} where \(K(m)\) is the complete elliptic integral of the first kind (see Equation B.41). At large distances, we can assume that \(M(<R') = M\) at all \(R'\) that contribute significantly to the integral in Equation (7.37) and we can then use the fundamental theorem of calculus to obtain the potential at large distances \begin{align} \Phi(R,z) & = G\,M\,\int_0^\infty \mathrm{d} R'\,\frac{\mathrm{d}}{\mathrm{d}R'}\left[\frac{2}{\pi\,\sqrt{(R+R')^2+z^2}}\,K\left(\frac{4RR'}{(R+R')^2+z^2}\right)\right]\nonumber\\ & = -G\,M\,\frac{2}{\pi\,\sqrt{R^2+z^2}}\,K\left(0\right) = -\frac{G\,M}{r}\,, \end{align} because \(K(0) = \pi/2\). Thus, we find that at large distances from a finite-mass, razor-thin mass distribution the potential approaches that of a point-mass of the same mass. In particular, this means that the rotation curve of any finite-mass disk will tend to the Keplerian behavior \(v_c(R) \propto R^{-1/2}\) at large distances.
7.3.3. Examples of razor-thin disks¶
7.3.3.1. The Mestel disk¶
The simplest example of using Equation (7.33) is for the Mestel disk (Mestel 1963). This is the razor-thin disk with surface density \begin{equation}\label{eq-diskgrav-mestel-surfdens} \Sigma(R) = \frac{v_c^2}{2\pi\,G\,R}\,. \end{equation} For this surface density, \(S_0(k)\) is easily computed using Equation (7.32) \begin{align} S_0(k) & = \int_0^\infty \mathrm{d} R' \,J_0(kR')\,R'\,\frac{v_c^2}{2\pi\,G\,R'}= \frac{v_c^2}{2\pi\,G\,k}\,, \end{align} because \(\int_0^\infty\mathrm{d} t\,J_m(t) = 1\) for all \(m\). The gravitational potential is then given by Equation (7.33) \begin{align} \Phi(R,z) & = -v_c^2\,\int_0^\infty\mathrm{d}k\,\frac{e^{-k|z|}}{k}\,J_0(kR)= v_c^2\,\ln\left(\sqrt{R^2+z^2}+|z|\right)\,.\label{eq-mestel-pot} \end{align} (the integral can be solved by noting that the derivative of the integral with respect to \(|z|\) is \(v_c^2\,\int_0^\infty\mathrm{d}k\,e^{-k|z|}\,J_0(kR) = 1/\sqrt{R^2+|z|^2}\) by Equation B.51). The circular velocity in the mid-plane from Equation (7.35) is given by \begin{align} v_c^2(R) & = v_c^2\,R\,\int_0^\infty\mathrm{d}k\,J_1(kR)= v_c^2\,, \end{align} which is simply constant (again using that \(\int_0^\infty\mathrm{d} t\,J_m(t) = 1\) for all \(m\)).
The Mestel disk therefore has a flat rotation curve and can be considered the infinitely thin limit of the flattened logarithmic potential, which also has a flat rotation curve. It is interesting to write the rotation curve \(v_c(R)\) in terms of the mass contained within \(R\) for the Mestel disk. The mass is given by \(M(< R) = 2\pi\,\int_0^R\,\mathrm{d}R'\,R'\,\Sigma(R')= v_c^2 R/G\), and therefore \begin{equation} v_c^2(R) = v_c^2 = \frac{G\,M(< R)}{R}\,. \end{equation} This is the same relation between \(v_c(R)\) and \(M(<R)\) as that which holds for any spherical potential. As we discussed at the beginning of this chapter, for general disk mass distributions the mass outside of \(R\) affects the rotation curve inside of \(R\), but this is not the case for the Mestel disk!
7.3.3.2. The razor-thin exponential disk¶
As we discussed in Chapter 1.2.1, the surface density profiles of disk galaxies are exponential with \(\Sigma(R) = \Sigma(0)\,e^{-R/R_d}\), which we will also find convenient to write as \(\Sigma(R) = \Sigma(0)\,e^{-\alpha R}\) below (we write the central surface density as \(\Sigma(0)\) to avoid confusion with \(\Sigma_0\), the \(m=0\) component of the surface density). For future reference, we compute the enclosed mass \(M(<R)\) \begin{equation}\label{eq-razorexp-encmass} M(< R) = 2\pi\,\Sigma(0)\,\left[R_d^2-R_d\left(R+R_d\right)\,e^{-R/R_d}\right]\,. \end{equation} For a razor-thin exponential disk, we can compute \(S_0(k)\) using Equation (7.32) \begin{align}\label{eq-S0-razorthinexp} S_0(k) &= -\Sigma(0)\,\frac{\partial}{\partial \alpha}\,\int_0^\infty \mathrm{d} R' \,J_0(kR')\,e^{-\alpha R'}= -\Sigma(0)\,\frac{\partial}{\partial \alpha}\,\left[(\alpha^2+k^2)^{-1/2}\right]\,, \end{align} using a trick introduced by Kuijken & Gilmore (1989c) and again using Equation (B.51). We could work out the derivative with respect to \(\alpha\), but will find it convenient to keep it. Using this result, we can compute the potential in the mid-plane using Equation (7.34) \begin{align} \Phi(R,0) & = 2\pi\,G\,\Sigma(0)\,\frac{\partial}{\partial \alpha}\,\int_0^\infty\mathrm{d}k\,J_0(kR)\,\left(\alpha^2+k^2\right)^{-1/2}\nonumber\\ & = \pi\,G\,\Sigma(0)\,R\,\left[I_1\left(\frac{\alpha R}{2}\right)\,K_0\left(\frac{\alpha R}{2}\right)-I_0\left(\frac{\alpha R}{2}\right)\,K_1\left(\frac{\alpha R}{2}\right)\right]\,, \end{align} where \(I_m(\cdot)\) and \(K_m(\cdot)\) are the modified Bessel functions of the first and second kind, respectively, defined in Equation (B.53). We have also used that \(\mathrm{d}\,I_0(x) / \mathrm{d}\,x = I_1(x)\) and \(\mathrm{d}\,K_0(x) / \mathrm{d}\,x = -K_1(x)\) (Equations B.55). The rotation curve is then most easily computed as the derivative of this expression for \(\Phi(R,0)\) (using that \(\mathrm{d}\,I_1(x) / \mathrm{d}\,x = I_0(x)-I_1(x)/x\) and \(\mathrm{d}\,K_1(x) / \mathrm{d}\,x = -K_0(x)-K_1(x)/x\); Equations B.56) \begin{equation} v_c^2(R) = \pi\,G\,\Sigma(0)\,\frac{R^2}{R_d}\,\left[I_0\left(\frac{R}{2R_d}\right)\,K_0\left(\frac{R}{2R_d}\right)-I_1\left(\frac{R}{2R_d}\right)\,K_1\left(\frac{R}{2R_d}\right)\right]\,, \end{equation} a result first derived by Freeman (1970). Thus, we see that using the formalism introduced in the previous section, we have been able to analytically compute the rotation curve for a realistic disk surface density. That’s pretty good!
Before showing the rotation curve of the razor-thin exponential disk, we can compute where it reaches a maximum, by finding \(\mathrm{d}\,v_c^2(R) / \mathrm{d}\,R = 0\). The derivative of \(v_c^2(R)\) is \begin{align} & \frac{\mathrm{d} v_c^2(R)}{\mathrm{d} R} = \pi\,G\,\Sigma(0)\nonumber\\ & \ \times \left[2\,\frac{R}{R_d}\,I_0\left(\frac{R}{2R_d}\right)\,K_0\left(\frac{R}{2R_d}\right)-\frac{R^2}{R_d^2}\,\left\{I_0\left(\frac{R}{2R_d}\right)\,K_1\left(\frac{R}{2R_d}\right)-I_1\left(\frac{R}{2R_d}\right)\,K_0\left(\frac{R}{2R_d}\right)\right\}\right]\,.\label{eq-diskgrav-dvcdr-exp} \end{align}
This is a rather complicated looking function, but by graphing it we can easily determine that the derivative becomes zero somewhere in the range \(2 < R/R_d < 3\). We cannot solve for this zero crossing analytically, but let's find it numerically:
[15]:
from scipy import optimize
opt= optimize.brentq(dvc2dr,2.,3.)
print("The exponential-disk rotation curve peaks at R/R_d = %.3f" % opt)
The exponential-disk rotation curve peaks at R/R_d = 2.150
we find the classical result that the rotation curve of a razor-thin exponential disk peaks at \(R \approx 2.15\,R_d\).
We now compare the rotation curve for the razor-thin exponential disk to that of the razor-thin Kuzmin disk with the same total mass and with scale length \(a\) such that the rotation curves peak at the same radius in Figure 7.11. Remember that the rotation curve of a Kuzmin disk peaks at \(R = \sqrt{2}\,a\).
[10]:
from galpy import potential
class EquivalentMassSphericalPotential(potential.Potential):
def __init__(self,amp=1.,cumulmassfunc=None,ro=None,vo=None):
potential.Potential.__init__(self,amp=amp,ro=ro,vo=vo)
self._cumulmassfunc= cumulmassfunc
def _Rforce(self,R,z,phi=0.,t=0.):
r= numpy.sqrt(R**2.+z**2.)
return -self._cumulmassfunc(r)*R/r**3.
figure(figsize=(8,5))
Rd= 1.
kzp= potential.KuzminDiskPotential(amp=10.,a=2.15/numpy.sqrt(2.)*Rd)
rp= potential.RazorThinExponentialDiskPotential(amp=10./(2.*numpy.pi),hr=Rd)
rp_spher= EquivalentMassSphericalPotential(amp=10.,\
cumulmassfunc=lambda R: Rd**2.-Rd*(R+Rd)*numpy.exp(-R/Rd))
line_rp= rp.plotRotcurve(Rrange=[0.,15.],\
label=r'$\mathrm{Exponential}$',yrange=[0.,2.5],
xlabel=r'$R/R_d$',ylabel=r'$v_c$',gcf=True)
line_rp_spher= rp_spher.plotRotcurve(Rrange=[0.,15.],\
label=r'$\mathrm{Spherical\ with\ same}\ M(<R)$',
overplot=True)
line_kuzmin= kzp.plotRotcurve(Rrange=[0.,15.],\
label=r'$\mathrm{Kuzmin}$',
overplot=True)
legend(handles=[line_rp[0],line_rp_spher[0],line_kuzmin[0]],
fontsize=18.,loc='upper right',frameon=False);
Figure 7.11: The rotation curve of the razor-thin exponential disk.
For the same mass and peak radius, the exponential disk has a higher maximum circular velocity than the Kuzmin disk. Indeed, the circular velocity for the exponential disk is larger than that for the equivalent Kuzmin disk at all radii. This is essentially because the exponential disk is more centrally concentrated than the Kuzmin disk.
Also shown in Figure 7.11 is a comparison of the rotation curve of the razor-thin exponential disk to that of the spherical mass distribution with the same enclosed mass, in this case given by Equation (7.45). We see that the razor-thin exponential disk’s rotation curve rises to a higher peak rotation velocity and then approaches the rotation curve of the equivalent-mass spherical distribution slowly from above. As for the Kuzmin disk in Figure 7.5, we see that the rotation curve of this disk mass distribution is only slightly different from that of the equivalent-mass spherical distribution, with the main difference being that the disky density gives a somewhat elevated circular velocity for the same radial profile.
7.3.3.3. A thin ring … and any razor-thin disk¶
Now that we are well warmed up in our use of Bessel functions and elliptic integrals, we can tackle the problem of computing the gravitational potential of a thin ring that we showed earlier in Figure 7.4 to demonstrate the essential differences between spherical and flattened potentials. We consider an infinitely thin, axisymmetric ring of radius \(a\) with density \begin{equation} \rho(R,z) = {M \over 2\pi\,a}\,\delta(R-a)\,\delta(z)\,. \end{equation} Therefore, the surface density is \(\Sigma(R) = M\,\delta(R-a)/(2\pi a)\). The parameter \(M\) is the total mass of the ring, as can be readily shown by integrating the surface density over all \((R,\phi)\). Using Equation (7.32), we can compute the function \(S_0(k)\) \begin{align} S_0(k) & = \int_0^\infty \mathrm{d} R' \,J_0(kR')\,R'\,\Sigma(R') = {M \over 2\pi}\,J_0(ka)\,. \end{align} So far so good. The next step is where things get more complicated. Using Equation (7.33), we can compute the potential \begin{align} \Phi(R,z) & = -GM\,\int_0^\infty\mathrm{d}k\,e^{-k|z|}\,J_0(kR)\,J_0(ka)\label{eq-thinring-almostpot} \end{align} To solve this integral, we again need the Laplace transformation of a product of two Bessel functions of the first kind of order zero and the result is (e.g., Kausel & Irfan Baig 2012) \begin{equation}\label{eq-thinring-pot} \Phi(R,z) = -{2 \over \pi}\,{GM \over \sqrt{(R+a)^2+z^2}}\,K\left({4\,R\,a \over [R+a]^2+z^2}\right)\,, \end{equation} where \(K(m)\) is again the complete elliptic integral of the first kind (see Equation B.41). Using the expression for the derivative of \(K(m)\) (Equation B.44), we can determine the forces for this potential. Taking the limit where \(R, z \gg a\), we again see that \(\Phi(R,z) \rightarrow -GM/r\), the potential of a point-mass (because \(K(0) = \pi/2\)).
We can decompose any razor-thin axisymmetric disk as an infinite collection of thin rings with infinitesimal masses \(\delta M(a)\) at radii \(a\) as \begin{align} \Sigma(R) & = \int \mathrm{d} a {\delta M(a) \over 2\pi\,a}\,\delta(R-a)\,, \end{align} and therefore \(\delta M(a) = 2\pi\,a\,\Sigma(a)\). Substituting \(\delta M(a)\) for \(M\) in Equation (7.53) and integrating over \(a\) then gives the potential of any razor-thin surface density \begin{equation}\label{eq-diskpot-pot-anyrazor} \Phi(R,z) = -4\,G\,\,\int \mathrm{d} a\,{a\,\Sigma(a)\over \sqrt{(R+a)^2+z^2}}\,K\left({4\,R\,a \over [R+a]^2+z^2}\right)\,. \end{equation} Another way of obtaining this result is by swapping the order of integration in Equation (7.33) and using the integral that we used in going from Equation (7.52) to Equation (7.53).
7.3.4. Potentials for disks with finite thickness¶
So far we have discussed solving the Poisson equation to obtain the gravitational potential and rotation curves for razor-thin disks. But real disks are not razor-thin: their aspect ratios are \(\approx10\) and the thickness of the disk has a significant effect on the orbits of stars in the disk. This will become clearer when we discuss orbits in disk mass distributions in Chapter 9.
To compute the potential for a disk with finite thickness, we can treat each layer as a razor-thin disk, obtain the contribution to the gravitational potential from each layer using Equation (7.33), and then add all of these contributions together to compute the total potential. In general this is cumbersome, because the \(S_0(k)\) factor in the integral for the potential will depend on the height of each layer. Fortunately, observations of disks demonstrate that the thickness is roughly independent of radius and therefore to a good approximation we can approximate the density of galactic disks as \begin{equation} \rho(R,\phi,z) = \Sigma(R,\phi)\,\zeta(z)\,, \end{equation} where we define \(\zeta(z)\) such that \(\int \mathrm{d}z\,\zeta(z) = 1\), such that \(\Sigma(R,\phi)\) is the surface density. If the disk is axisymmetric, this becomes \(\rho(R,z) = \Sigma(R)\,\zeta(z)\). For such a density distribution, the surface density of a layer at \(z=z'\) is given by \(\Sigma(R)\,\zeta(z')\,\mathrm{d}z\). When we only consider this layer, \(\zeta(z')\) is a constant and the potential \(\mathrm{d}\Phi(R,z)\) of this layer at \((R,z)\) is then given by Equation (7.33) \begin{equation} \mathrm{d}\Phi(R,z) = -2\pi\,G\,\mathrm{d}z'\,\zeta(z')\,\int_0^\infty\mathrm{d}k\,e^{-k|z-z'|}\,J_0(kR)\,S_0(k)\,. \end{equation} The total potential is therefore given by \begin{equation}\label{eq-pot-thickdisk} \Phi(R,z) = -2\pi\,G\,\int \mathrm{d}z'\,\zeta(z')\,\int_0^\infty\mathrm{d}k\,e^{-k|z-z'|}\,J_0(kR)\,S_0(k)\,. \end{equation} and the radial and vertical component of the force are \begin{align} F_R(R,z) & = 2\pi\,G\,\int \mathrm{d}z'\,\zeta(z')\,\int_0^\infty\mathrm{d}k\,e^{-k|z-z'|}\,J_1(kR)\,k\,S_0(k)\,,\\ F_z(R,z) & = 2\pi\,G\,\int \mathrm{d}z'\,\zeta(z')\,\mathrm{sign}(z-z')\,\int_0^\infty\mathrm{d}k\,e^{-k|z-z'|}\,J_0(kR)\,k\,S_0(k)\,. \end{align}
7.3.4.1. The double-exponential disk¶
As an example of a disk with a finite thickness, let’s look at a finite-thickness version of the razor-thin exponential disk that we discussed above. We thicken the disk by replacing the \(\delta(z)\) vertical profile by an exponential \begin{equation} \rho(R,z) = \Sigma(0)\,e^{-R/R_d}\,\frac{1}{2\,z_d}\,e^{-|z|/z_d} =\Sigma(0)\,e^{-\alpha R}\,\frac{\beta}{2}\,e^{-\beta |z|}\, \end{equation} where we \(z_d = 1/\beta\) is the scale height. Because of the two exponentials, this disk is known as the double-exponential disk. It is a very popular model for disk galaxies.
Computing the gravitational potential and forces of the double-exponential disk can be done through the equations given in the previous subsection. These can be somewhat simplified, as described in Appendix A of Kuijken & Gilmore (1989c). We won’t go into that level of detail here. Instead, we will gain some analytical insight on the dependence of the rotation curve on the thickness of the disk.
To compute the rotation curve of the double-exponential disk, we first compute \(\Phi(R,0)\). Equation (7.58) applied with \(\zeta(z') = \beta\,e^{-\beta|z|}/2\) gives \begin{align}\label{eq-phi0-dblexp-deriv1} \Phi(R,0) & = -\pi\,G\,\beta\,\int \mathrm{d}z'\,\int_0^\infty\mathrm{d}k\,e^{-(k+\beta)|z'|}\,J_0(kR)\,S_0(k)\nonumber\\ & = -2\pi\,G\,\beta\,\int_0^\infty\mathrm{d}k\,\frac{1}{\beta+k}\,J_0(kR)\,S_0(k)\,. \end{align} where \(S_0(k)\) is given by Equation (7.46). This is a tough integral to solve! In the limit of a thin disk, \(\beta \gg 1\) and we can approximate \(1/(\beta+k) \approx (1-k/\beta)/\beta\). Note that the approximation we are making is not entirely warranted: the integral extends over all \(k\), but we require \(k/\beta \ll 1\) for this to be a good approximation. However, for large \(k\), \(S_0(k) \propto 1/k\) and the amplitude of the oscillatory \(J_0(kR) \propto 1/\sqrt{k}\), so contributions to the integral at high \(k\) are small. In this approximation we have that \begin{align} \Phi(R,0) & \approx -2\pi\,G\,\int_0^\infty\mathrm{d}k\,\left(1-\frac{k}{\beta}\right)\,J_0(kR)\,S_0(k)\nonumber\\ & \approx \Phi_{\mathrm{razor-thin}}(R,0)+\frac{2\pi\,G}{\beta}\,\int_0^\infty\mathrm{d}k\,k\,J_0(kR)\,S_0(k)\nonumber\\ & \approx \Phi_{\mathrm{razor-thin}}(R,0)+\frac{2\pi\,G}{\beta}\,\Sigma(R)\,, \end{align} where to go to the second line we have substituted the razor-thin potential for the first term and to go to the third line we have used Equation (7.31). The circular velocity curve is then approximately \begin{equation}\label{eq-vc-dblexp-approx} v_c^2(R) \approx v_{c,\mathrm{razor-thin}}^2(R) -2\pi\,G\,\Sigma(0)\,\frac{z_d}{R_d}\,R\,e^{-R/R_d}\,, \end{equation} where we have used that \(\mathrm{d}\,\Sigma(R) / \mathrm{d}\,R = -\Sigma(R)/R_d\) for the exponential disk. Thus, we see that a thickened exponential disk has a lower circular velocity than the razor-thin disk with the same surface-density profile; the difference in \(v_c^2(R)\) is proportional to the aspect ratio \(z_d/R_d\) of the disk and proportional to \(R\,e^{-R/R_d}\). To get a better sense of the magnitude of the difference, we use that \(v_{c,\mathrm{razor-thin}}^2(R) \approx 2\,\pi\,G\,\Sigma(0)\,R/5\) up to a numerical factor close to 1 at \(R\approx 2R_d\) (near the peak of the rotation curve) and therefore the relative difference in the rotation curve is \begin{equation} \frac{\Delta v_c}{v_{c,\mathrm{razor-thin}}} \approx -\frac{5z_d}{2R_d}\,e^{-R/R_d} \approx -7\,\%\,\left(\frac{z_d/R_d}{0.2}\right)\,. \end{equation} Thus, we see that the rotation curve is only affected by \(\lesssim 10\,\%\) for thick, exponential disks compared to the razor-thin case. The disk thickness causes the rotation curve to be lower than that for a razor-thin disk: very flattened mass is able to sustain a higher rotation velocity than less flattened mass.
Figure 7.12 compares the rotation curve for the double-exponential disk with \(z_d/R_d = 0.2\), numerically computed using galpy, to that of the razor-thin disk and to the approximation in Equation (7.64). We see that the approximation does an excellent job of correcting the razor-thin rotation curve for the effect of disk thickness.
[12]:
from scipy import special
from galpy import potential
from matplotlib.ticker import FuncFormatter
figure(figsize=(8.5,5))
Rd= 0.2
zd= 0.2*0.2
rp= potential.RazorThinExponentialDiskPotential(amp=10.,hr=Rd)
dp= potential.DoubleExponentialDiskPotential(amp=10./2./zd,hr=Rd,hz=zd)
line_rp= rp.plotRotcurve(Rrange=[0.,3.],\
label=r'$\mathrm{Razor-thin\ exponential}$',
yrange=[0.,2.5],gcf=True,
xlabel=r'$R/R_d$',ylabel=r'$v_c$')
line_dp= dp.plotRotcurve(Rrange=[0.,3.],\
label=r'$\mathrm{Double\ exponential}\,(z_d/R_d = 0.2)$',
overplot=True)
# Also compute the approximation above
Rs= numpy.linspace(0.,3.,1001)
line_approx= plot(Rs,numpy.sqrt(numpy.pi*rp._amp*Rs**2./Rd\
*(special.i0(Rs/2./Rd)*special.k0(Rs/2/Rd)
-special.i1(Rs/2/Rd)*special.k1(Rs/2./Rd))\
-2*numpy.pi*rp._amp*Rs*zd/Rd*numpy.exp(-Rs/Rd)),
label=r'$\mathrm{Approximation}$',ls='--')
# To get the ticklabels correct, need to multiply the labels
def perRd(x,pos):
return r'$%.0f$' % (x*5.)
gca().xaxis.set_major_formatter(FuncFormatter(perRd))
legend(handles=[line_rp[0],line_dp[0],line_approx[0]],
fontsize=18.,loc='upper right',frameon=False);
Figure 7.12: The rotation curve of the double-exponential disk.
7.3.5. Gravitational potentials for non-axisymmetric density distributions¶
We are now ready to discuss the general solution of the Poisson equation for a disky, non-axisymmetric mass distribution. For simplicity, we will assume that the vertical and planar parts of the density separate as \(\rho(R,\phi,z) = \Sigma(R,\phi)\,\zeta(z)\), but everything we derive is easy to generalize to the general case: the Fourier and Hankel transforms of \(\Sigma(R,\phi)\) that we perform to obtain a solution would depend on \(z\) in the general case, but otherwise all expressions would be the same.
As discussed in Section 7.3.1, we start by taking the Fourier transform over \(\phi\), to get the Poisson equation for \(\Phi_m(R,z)\) (Equation 7.24, which we repeat here for the case where the density separates) \begin{equation} \frac{1}{R}\,\frac{\partial}{\partial R}\left(R\,\frac{\partial \Phi_m(R,z)}{\partial R}\right) - \frac{m^2}{R^2}\,\Phi_m(R,z)+\frac{\partial^2\Phi_m(R,z)}{\partial z^2} = 4\pi\,G\,\Sigma_m(R)\,\zeta(z)\,. \end{equation} Similar to what we did in Section 7.3.2, we can then easily check using the definition of the Bessel functions that the function \begin{equation} \tilde{\Phi}_m(R;k) = -2\pi\,G\,e^{-k|z|}\,J_m(kR) \end{equation} is a solution of the Poisson equation, with surface density \begin{equation} \tilde{\Sigma}(R;k) = k\,J_m(kR)\,. \end{equation} Just like for the axisymmetric case, if we can write the target surface density \(\Sigma_m(R)\) as a sum over contributions \(S_m(k)\,\tilde{\Sigma}_m(R;k)\) with amplitudes \(S_m(k)\), then the potential is the sum over contributions \(S_m(k)\,\tilde{\Phi}_m(R;k)\), because of the linearity of the Poisson equation. Thus, we write \begin{equation}\label{eq-invhankel-surfdens-m} \Sigma_m(R) = \int_0^\infty \mathrm{d}k\,J_m(kR)\,k\,S_m(k)\,, \end{equation} where we can determine the \(S_m(k)\) as \begin{equation}\label{eq-hankel-surfdens-m} S_m(k) = \int_0^\infty \mathrm{d} R' \,J_m(kR')\,R'\,\Sigma_m(R')\,, \end{equation} because of the Fourier-Bessel theorem (Equation B.50). \(S_m(k)\) is the m-th order Hankel transform of \(\Sigma_m(R)\). Like for the thick, axisymmetric disks in Equation (7.58), we obtain the total gravitational potential \(\Phi_m(R,z)\) by combining the contributions from all of the different layers and then use the inverse Fourier transform (Equation 7.23) to get \(\Phi(R,\phi,z)\). After all is said and done we arrive at \begin{align}\label{eq-pot-disk-general} \Phi(R,\phi,z) & = -G\,\sum_{m=-\infty}^{\infty}\int \mathrm{d}z'\,\zeta(z')\,\int_0^\infty\mathrm{d}k\,e^{-k|z-z'|}\,J_m(kR)\\ & \qquad \qquad \qquad \int_0^\infty \mathrm{d} R' \,J_m(kR')\,R'\,\int_0^{2\pi}\mathrm{d}\phi'\,e^{-im\,\left(\phi'-\phi\right)}\,\Sigma(R',\phi')\,.\nonumber \end{align}
Equation (7.71) is the general solution for the potential of a disk. This is clearly a complicated expression: the solution involves a quadruple integral and an infinite sum over \(m\). It is therefore best to think of this as a formal full solution that can be usefully applied in cases of higher symmetry, such as those that we discussed above. For example, it can be used to get the potential for a razor-thin, non-axisymmetric disk with \(\zeta(z) = \delta(z)\) and a simple Fourier transform in \(\phi\) (e.g., an elliptical disk which only has non-zero \(\Sigma_m\) for \(m=0\) and \(m=2\)).
.
.
.