7.2. Simple gravitational potentials for flattened mass distributions¶
Determining the gravitational potential for general flattened mass distributions is difficult. Newton’s theorems, which we discussed in Chapter 2.2, do not hold for flattened systems and we can therefore not proceed in the same manner as we did for spherical mass distributions (where we could add the potential from concentric, spherical shells to build up the total potential for any spherical mass distribution). To illustrate the complication, we compare the gravitational potential from a thin spherical shell with that of a thin ring of the same total mass in the left panel of Figure 7.4 (we compare the potentials along \(z = 0.1\,R\) to avoid a singularity that exists when crossing through the ring at \(z=0\); we will discuss how to actually compute the potential and force due to a ring of matter below).
[6]:
from galpy.potential import SphericalShellPotential, RingPotential
figure(figsize=(12,4))
sp= SphericalShellPotential(amp=1.,a=1.5)
rp= RingPotential(amp=1.,a=1.5)
rs= numpy.linspace(0.1,7.,101)
subplot(1,2,1)
line_shell= galpy_plot.plot(rs,[sp(r,0.1*r) for r in rs],
yrange=[-1.,0.],
xlabel=r'$r$',
ylabel=r'$\Phi$',
gcf=True,
label=r'$\mathrm{Shell}$')
line_ring= galpy_plot.plot(rs,[rp(r,0.1*r) for r in rs],overplot=True,
label=r'$\mathrm{Ring}$')
legend(handles=[line_shell[0],line_ring[0]],
fontsize=18.,loc='lower right',frameon=False)
subplot(1,2,2)
line_shell= galpy_plot.plot(rs,[sp.Rforce(r,0.) for r in rs],
yrange=[-1.,1.],
xlabel=r'$r$',
ylabel=r'$g_r$',
gcf=True,
label=r'$\mathrm{Shell}$')
line_ring= galpy_plot.plot(rs,[rp.Rforce(r,0.) for r in rs],overplot=True,
label=r'$\mathrm{Ring}$')
legend(handles=[line_shell[0],line_ring[0]],
fontsize=18.,loc='lower right',frameon=False)
tight_layout(w_pad=0.2);
Figure 7.4: The gravitational potential (left) and force (right) of a ring and of a shell.
Unlike for the spherical shell, the potential within the ring is not constant and thus Newton’s first shell theorem does not hold (why the theorem doesn’t hold is clear from the geometric proof that we discussed in Chapter 2.2). The potential actually decreases with radius inside of the ring, implying that the gravitational force is pointing towards the closest side of the ring (why this is should again be clear from the geometric proof of Newton’s shell theorem). The right panel of Figure 7.4 shows the actual gravitational force, now at \(z=0\) to highlight the (integrable) singularity when passing through the ring. The amplitude of the force due to the ring within its plane is much higher than that of a shell and a ring of matter will therefore accelerate bodies within its plane more than a shell of the same mass. At large distances from the ring, the potential and force approach that of the spherical shell, which is the same as that of a point mass. This is a general rule for mass distributions of finite mass: at large distances their force becomes indistinguishable from that of a point of the same mass.
Thus, the potential of a ring is significantly different from that of a spherical shell. Because thin disks are more accurately approximated as consisting of rings rather than shells, their gravitational forces differ from those of spherical mass distributions in the same way that those of a ring differ from those of a shell: For a given mass, disks cause larger accelerations (and thus, e.g., velocities of a circular orbit) and the mass outside a given radius affects the gravitational acceleration felt by a body. We will discuss a general solution to determining the potential for any mass distribution below, but we will see that it is quite complicated, even in the limit that the distribution is axisymmetric and that the mass distribution is flattened into a zero thickness disk.
An easier method for obtaining gravitational potential models of flattened mass distributions is to posit a potential \(\Phi(R,z)\) (for axisymmetric potentials) and then compute the corresponding density using the Poisson equation as \(\rho = \nabla^2 \Phi / (4\pi\,G)\). By not solving for the potential for a given density, but rather solving for the density corresponding to a given potential, we get around having to solve a partial differential equation. Naturally, this simplification comes at a cost: (i) we cannot choose the density distribution, but rather have to work with whatever we compute with \(\nabla^2 \Phi / (4\pi\,G)\) (this may not be a good enough model for the system that we are studying) and (ii) there is no guarantee that the density resulting from \(\rho = \nabla^2 \Phi / (4\pi\,G)\) is everywhere positive, so we may end up with a model that is unphysical in certain regions. Despite these issues, there are several useful, simple flattened gravitational potential models that are in wide-spread use and that we discuss in this section.
7.2.1. Razor-thin disk: The Kuzmin model¶
A simple flattened axisymmetric potential is that given by the following expression \begin{equation}\label{eq-kuzmin-pot} \Phi(R,z) = -\frac{G\,M}{\sqrt{R^2+(|z|+a)^2}}\,, \end{equation} with \(a \geq 0\) a parameter of the model. This potential resembles a point-mass potential \(\Phi(R,z) = -G\,M/\sqrt{R^2+z^2}\): for any \(z > 0\), Equation (7.5) is equivalent to the potential of a point-mass located at \((R,z) = (0,-a)\), while for any \(z < 0\), Equation (7.5) is equivalent to the potential of a point-mass located at \((R,z) = (0,a)\). Because the potential at any point \((R,z)\) with \(z \neq 0\) is equivalent to that of a point mass centered on a different point (which therefore has zero density at \((R,z)\)), the density corresponding to the potential in Equation (7.5) at any point with \(z \neq 0\) must be zero. This reasoning does not hold for points with \(z=0\): For those points, the potential is equivalent to a Plummer sphere, which has non-zero density everywhere (but the density at \(z=0\) is not that of a Plummer sphere, because the equivalence only holds for \(z=0\) exactly, and not for \(z=\epsilon\) and therefore not for the derivatives of the potential). The density corresponding to the potential in Equation (7.5) is therefore of the form \begin{equation}\label{eq-razorthin-dens} \rho(R,\phi,z) = \Sigma(R,\phi)\,\delta(z)\,, \end{equation} where \(\delta(z)\) is the Dirac delta function (Appendix B.3.2). This represents a zero-thickness disk or a razor-thin disk. The function \(\Sigma(R,\phi)\) is the surface density, which has units of mass over length squared. For the potential in Equation (7.5) the density will clearly be axisymmetric and, thus, \(\Sigma(R,\phi) = \Sigma(R)\).
To determine the relation between the surface density \(\Sigma(R)\) and the potential for any razor-thin disk, we integrate the Poisson equation over a volume \(\mathcal{V}\) \begin{equation} \int_{\mathcal{V}} \mathrm{d} \vec{x} \,\nabla^2 \Phi(\vec{x}) = 4\pi\,G\,\int_{\mathcal{V}}\mathrm{d}\vec{x}\,\rho(\vec{x})\,. \end{equation} Plugging in Equation (7.6), switching to cylindrical coordinates, and using the divergence theorem (Equation B.5) to write the volume integral on the left-hand side as an integral over the surface \(\mathcal{S}\) of the volume \(\mathcal{V}\), this becomes \begin{equation} \int_{\mathcal{S}} \mathrm{d}\mathcal{S} \cdot \nabla \Phi(R,z) = 4\pi\,G\,\int R\,\mathrm{d} R\,\mathrm{d} \phi \,\Sigma(R)\,, \end{equation} were \(\mathrm{d}\mathcal{S}\cdot\) indicates integration over the surface of the component of \(\nabla \Phi(R,z)\) perpendicular to the surface. We now choose the volume \(\mathcal{V}\) such that it is a small wedge centered on \(z=0\) that includes a part of the \(z=0\) plane that is so small that we can assume that \(\Sigma(R)\) does not vary over it and that extends to the same height above and below the \(z=0\) plane and lies between lines of constant \(\phi\). Then the contribution to the surface integrand at constant \(R\) cancel each other, the contributions at constant \(\phi\) are zero because the potential is axisymmetric, and the contribution from the portions above and below the plane, at constant \(z\), are equal to each other and multiplied by the surface area of the wedge, \(R\,\mathrm{d} R\,\mathrm{d} \phi\). In the limit that the thickness of the volume goes to zero, this gives \begin{equation}\label{eq-razorthin-potderiv-surfdens} \pm\frac{\partial \Phi(R,z)}{\partial z}\Bigg|_{z=0\pm} = 2\pi\,G\,\Sigma(R)\,, \end{equation} where \(z=0\pm\) indicates that the limit is taken from the positive or negative \(z\) side. This equation holds in general for any razor-thin, axisymmetric disk.
For the potential in Equation (7.5), applying Equation (7.9) gives \begin{equation}\label{eq-kuzmin-surfdens} \Sigma(R) = \frac{M\,a}{2\pi}\,\frac{1}{\left(R^2+a^2\right)^{3/2}}\,. \end{equation} The potential–density pair in Equations (7.5) and (7.10) was introduced by Kuzmin (1956) and later independently re-discovered by Toomre (1963) as part of a more general family of potential–density pairs. It is known as the Kuzmin model or Toomre’s model 1.
The Kuzmin model demonstrates the success of the first approach to finding gravitational-potential models for flattened mass distributions: By starting from the simple form for the potential in Equation (7.5), we have ended up with a model corresponding to a razor-thin disk. This model could therefore be a good first approximation to the gravitational potential of a disk. It also shows that gravitational potentials for disks do not have to be particularly complicated. The form in Equation (7.5) is of a similar complexity to most of the spherical potentials that we discussed in Chapter 2.4.
For disk potentials, we can compute the rotation curve in the mid-plane \(z=0\) by balancing the centripetal acceleration \(a_R = v_c^2/R\) with the radial force as in \begin{equation}\label{eq-vc-midplane-general-firststep} a_R(R) = \frac{v_c^2(R)}{R} = -F_R(R,z=0) = \frac{\partial \Phi(R,z)}{\partial R}\Bigg|_{R=R;\, z=0}\,. \end{equation} Therefore, the rotation curve in the mid-plane is given by \begin{equation}\label{eq-vc-midplane-general} v_c(R) = \sqrt{R\,\frac{\partial \Phi(R,z)}{\partial R}\Bigg|_{R=R;\, z=0}}\,. \end{equation} The rotation curve for the Kuzmin disk is then given by \begin{equation}\label{eq-vc-kuzmin} v_c(R) = \sqrt{\frac{G\,M\,R^2}{\left(R^2+a^2\right)^{3/2}}}\,. \end{equation}
We can compare this rotation curve to the one we would get if the relation \(v_c^2(R) = G\,M(<R)/R\) that holds for spherical potentials would also hold for disks. For this we need to cumulative mass profile, which for the Kuzmin disk is \(M(< R) = M\,\left(1-a/\sqrt{R^2+a^2}\right)\). To plot the rotation curve of a spherical potential with this cumulative mass distribution, we define a new class of potentials in galpy. This can be easily done in the galpy framework. Because we are only interested (for now) in using this potential to get \(v_c(R)\), all we need to implement is \(F_R(R,z=0)\) (but for generality’s sake, we will properly implement \(F_R(R,z)\)):
[7]:
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.
The rotation curve as a function of \(R/a\) for the Kuzmin model together with the equivalent-mass spherical model is shown in Figure 7.5.
[8]:
from galpy import potential
figure(figsize=(8,5))
a= 1.
kzp= potential.KuzminDiskPotential(amp=5.,a=a)
kzp_spher= EquivalentMassSphericalPotential(amp=5.,\
cumulmassfunc=lambda R: 1.-a/numpy.sqrt(R**2+a**2))
line_kuzmin= kzp.plotRotcurve(Rrange=[0.,15.],\
label=r'$\mathrm{Kuzmin}$',yrange=[0.,1.75],
xlabel=r'$R/a$',ylabel=r'$v_c$',gcf=True)
line_kuzmin_spher= kzp_spher.plotRotcurve(Rrange=[0.,15.],\
label=r'$\mathrm{Spherical\ with\ same}\ M(<R)$',
overplot=True)
legend(handles=[line_kuzmin[0],line_kuzmin_spher[0]],
fontsize=18.,loc='upper right',frameon=False);
Figure 7.5: The rotation curve of the Kuzmin disk.
At \(R/a \ll 1\) the rotation curve of the Kuzmin disk in Equation (7.13) goes as \(v_c(R) \propto R\); this is similar to the behavior in a homogeneous-density sphere, because the surface density is constant at \(R/a \ll 1\) (see Equation 7.10). At \(R/a \gg 1\), the rotation curve drops as \(v_c(R) \propto R^{-1/2}\), similar to the rotation curve of a point-mass. Because the rotation curve goes from rising linearly to falling, the rotation curve reaches a peak. From Equation (7.13) it is straightforward to compute the location of the peak, which is \begin{equation} R_{\mathrm{peak}} = \sqrt{2}\,a\,. \end{equation} These features are all obvious in Figure 7.5. The shape of the Kuzmin rotation curve is a typical disk rotation curve for \(a \approx R_d\), the scale length of the disk. Comparing the Kuzmin rotation curve to that of the spherical mass distribution with the same enclosed-mass profile, we see that the rotation curves are quite similar even though the shape of the density distributions is completely different. This is an excellent demonstration of a general rule for disk rotation curves: they are largely set by the overall radial behavior of the mass distribution and only mildly dependent on its shape. This is why many useful insights in galactic dynamics can be obtained by approximating mass distributions as easy-to-work-with spherical distributions.
7.2.2. Thickened disk: the Miyamoto-Nagai model¶
The Kuzmin model provides a simple, convenient model of a razor-thin disk. But we would also like to have simple models for disks with a finite thickness. The Milky Way disk is razor-thin to a good approximation with a ratio of the scale length \(R_d\) to the scale height \(z_d\) of the mass profile of \(R_d / z_d \approx 10\), but at the next level of approximation taking into account the non-zero \(z_d\) is important, especially when looking at the properties of orbits in the disk mass distribution.
To get a thickened-disk model starting from the Kuzmin model, we can use a similar approach as that which we used in Chapter 2.4: There we saw several examples of extended spherical mass distributions that were obtained by different types of smoothing of the \(\Phi(r) \propto 1/r\) behavior of the point-mass potential. The smoothing was done by replacing \(1/r\) by expressions such as \(1/\sqrt{r^2+b^2}\) (Plummer), \(1/(b+\sqrt{r^2+b^2})\) (isochrone), and \(1/(r+a)\) (Hernquist). Similarly, we can obtain a thickened Kuzmin model by replacing \(|z|\) in Equation (7.5) by \(\sqrt{z^2+b^2}\), with \(b\) a new parameter of the model. This gives the potential \begin{equation}\label{eq-pot-mn} \Phi(R,z) = -\frac{G\,M}{\sqrt{R^2+(\sqrt{z^2+b^2}+a)^2}}\,. \end{equation} Besides reducing to the razor-thin Kuzmin disk when \(b \rightarrow 0\), this potential reduces to that of a Plummer sphere when \(a\rightarrow0\), a spherical density distribution. Therefore it is clear that the ratio \(b/a\) controls how flattened the mass distribution is: \(b/a \gg 1\) gives an approximately spherical distribution with a density profile similar to that of a Plummer sphere, while \(b/a \ll 1\) corresponds to a significantly flattened distribution, with a surface-density profile similar to that of the Kuzmin model. The potential in Equation (7.15) is a Miyamoto-Nagai potential and it remains a very popular simplified model for the gravitational potential of disks (Miyamoto & Nagai 1975).
The density corresponding to the potential in Equation (7.15) can be computed using the Poisson equation, which gives the following complicated expression \begin{equation}\label{eq-dens-mn} \rho(R,z) = \left(\frac{b^2\,M}{4\pi}\right)\,\frac{a\,R^2+\left(3\,\sqrt{z^2+b^2}+a\right)\,\left(\sqrt{z^2+b^2}+a\right)^2}{\left(R^2+(\sqrt{z^2+b^2}+a)^2\right)^{5/2}\,\left(z^2+b^2\right)^{3/2}}\,. \end{equation}
As an example of the density and potential flattening, let’s look at a Miyamoto-Nagai potential with \(b/a = 0.2\). Contours of the density in in the \((R,z)\) plane are displayed in Figure 7.6.
[9]:
from galpy import potential
figure(figsize=(10,4))
mnp= potential.MiyamotoNagaiPotential(amp=10.,a=1.,b=0.2)
mnp.plotDensity(justcontours=True,rmin=-1.5,rmax=1.5,nrs=100,nzs=100,gcf=True)
gca().set_aspect('equal')
xlabel(r'$R/a$')
ylabel(r'$z/a$');
Figure 7.6: The density of a Miyamoto-Nagai disk.
We see that the density is indeed highly flattened (the aspect ratio of this figure is such that a spherical density would appear circular). Contours of the potential are shown in Figure 7.7.
[10]:
from galpy import potential
figure(figsize=(10,4))
mnp= potential.MiyamotoNagaiPotential(amp=10.,a=1.,b=0.2)
mnp.plot(justcontours=True,rmin=-1.5,rmax=1.5,nrs=100,nzs=100,gcf=True)
gca().set_aspect('equal')
xlabel(r'$R/a$')
ylabel(r'$z/a$');
Figure 7.7: The potential of a Miyamoto-Nagai disk.
The potential is also flattened, but much less so than the density. This is a general rule for flattened potentials: the flattening (as a deviation from sphericity) in the potential is always much smaller than that in the density. The Kuzmin model above provides another extreme example of this. There the density is contained in an infinitely thin disk, while the potential is extended in \(z\).
The radial behavior of the Miyamoto-Nagai density in Equation (7.16) is quite different from the observed exponential profile of galactic disks. At \(z=0\), the density in Equation (7.16) behaves as \(R^{-3}\) for \(R \gg a\) and the large-radius density fall-off is therefore simply a power-law of \(R\) rather than exponential. This means in particular that the density of the Miyamoto-Nagai model remains too high at large radii (e.g., at \(R \gg 15\,\mathrm{kpc}\) for the Milky Way) compared to realistic models of galactic disks. However, because the model also has a large flattening (\(b/a \ll 1\)) the total mass at large radii in the disk component at large radii is still much smaller than that in the dark halo at similar radii. Therefore, unless one is looking at orbits exactly in the mid-plane of the disk at radii much beyond the disk, the Miyamoto-Nagai model provides a reasonable model for a disk.
Because of its simplicity, the Miyamoto-Nagai model is perhaps the most often used model for a galactic disk in orbit modeling. This is because computing the gravitational potential for more realistic density distribution, such as that of an exponential disk, is significantly more complex and comes at large computational cost (as we will see below). However, it is important to keep the limitations of the Miyamoto-Nagai disk in mind when using it as an approximation for a galactic disk, in particular its problematic large \(R\) behavior.
We can compute the rotation curve for the Miyamoto-Nagai potential using Equation (7.12). In Figure 7.8, we compare the rotation curves of a Kuzmin potential and a Miyamoto-Nagai potential with the same mass \(M\) and the same scale length \(a\), but with a flattening \(b/a = 0.2\).
[9]:
from galpy import potential
figure(figsize=(8,5))
kzp= potential.KuzminDiskPotential(amp=10.,a=1.)
mnp= potential.MiyamotoNagaiPotential(amp=10.,a=1.,b=0.2)
line_kuzmin= kzp.plotRotcurve(Rrange=[0.,15.],\
label=r'$\mathrm{Kuzmin}$',yrange=[0.,2.5],
xlabel=r'$R/a$',ylabel=r'$v_c$',gcf=True)
line_mn= mnp.plotRotcurve(Rrange=[0.,15.],\
label=r'$\mathrm{Miyamoto-Nagai}\,(b/a = 0.2)$',
overplot=True)
legend(handles=[line_kuzmin[0],line_mn[0]],
fontsize=18.,loc='upper right',frameon=False);
Figure 7.8: The rotation curve of the Miyamoto-Nagai disk.
We see that the rotation curve of the Miyamoto-Nagai model rises somewhat less quickly than that of the Kuzmin potential and reaches a lower peak rotation velocity at a slightly larger radius. At \(R > 2a\), the Miyamoto-Nagai rotation curve quickly approaches that of the Kuzmin model from below.
7.2.3. Models with flat rotation curves: flattened logarithmic potential¶
The Miyamoto-Nagai model above is a good and simple model for the potential of a galactic disk, but to represent the total mass distribution of a galaxy it needs to be combined with a bulge component and a dark-halo component. Because the Miyamoto-Nagai model itself has three parameters (mass \(M\), scale length \(a\), and scale height \(b\)) and each of the additional components requires at least two parameters (a mass and scale parameter; see the examples of spherical potentials in Section 2.4), the total number of parameters to represent a galaxy’s mass distribution becomes quite large. In general this is a good thing: the large number of parameters gives us much flexibility in the types of mass distributions that we can build from these simple building blocks. But sometimes we want to represent a galaxy with the fewest parameters as possible, trading flexibility in the mass model for an even simpler parameterization.
In Chapter 2.4.5 we discussed the logarithmic potential: a spherical potential \(\Phi(r) = v_c^2\, \ln r\) for which \(v_c(r) = \mathrm{constant}\). Because the logarithmic potential has a flat rotation curve, which is the rotation curve that many galaxies are observed to have, it is a simple, one-parameter model for the total mass distribution of a galaxy. However, it is a spherical model and therefore does not provide a good model for a disk galaxy with a flat rotation curve: Even when the dark-matter halo is approximately spherical, the total mass distribution that combines the halo and disk (and possible other components) is significantly flattened for large disk galaxies. Therefore, we would like to construct a flattened version of the spherical logarithmic potential to represent the total mass distribution of disk galaxies.
A general approach to flattening any spherical potential is to replace every occurrence of \(r\) with \begin{equation}\label{eq-flatten-m} m^2 = R^2 + \frac{z^2}{q^2}\,. \end{equation} The parameter \(q\) is known as the flattening. This can be done either in the expression for the potential or in that of the density. However, it cannot be done in both at the same time without breaking the Poisson equation: you have to choose whether to flatten the potential or the density. The contours in the \((R,z)\) plane of a spherical potential \(\Phi(r)\) flattened to \(\Phi_{\mathrm{flat}}(R,z) = \Phi(m)\) through Equation (7.17) will be ellipses with axis ratio \(q\); when the density \(\rho(r)\) is flattened to \(\rho_{\mathrm{flat}}(R,z) = \rho(m)\) the contours of the density are ellipses with axis ratio \(q\). It is easier to flatten the potential through Equation (7.17) and to calculate the resulting density using the Poisson equation as \(\rho = \nabla^2 \Phi / (4\pi\,G)\) than to flatten the density and solving the Poisson equation for the potential. Therefore, flattening the potential is the more popular approach. But as always when positing a form for the potential rather than the density, this approach comes at the cost of not controlling the resulting density and having to work with whatever density comes out of \(\rho = \nabla^2 \Phi / (4\pi\,G)\).
When we apply Equation (7.17) to the spherical logarithmic potential, we get the flattened logarithmic potential \begin{equation}\label{eq-pot-flatlog} \Phi(R,z) = \frac{v_c^2}{2}\,\ln \left(R^2 + \frac{z^2}{q^2}\right)\,. \end{equation} Like the spherical logarithmic potential, this flattened version has a flat rotation curve in the mid-plane: \(v_c(R) = v_c\).
The potential contours in \((R,z)\) of a flattened logarithmic potential with \(q=0.7\) are shown in Figure 7.9.
[12]:
from galpy import potential
figure(figsize=(10,4))
lp= potential.LogarithmicHaloPotential(amp=1.,q=0.7)
lp.plot(justcontours=True,rmin=-1.5,rmax=1.5,nrs=100,nzs=100,gcf=True)
gca().set_aspect('equal')
xlabel(r'$R$')
ylabel(r'$z$');
Figure 7.9: The density of a mildly-flattened logarithmic potential.
These contours appear only mildly flattened (compare, e.g., to the potential contours of the Miyamoto-Nagai model with \(b/a = 0.2\) in Figure 7.7 above), as would be appropriate for a disk plus a spherical dark-matter halo. However, when we look at the density corresponding to this potential in Figure 7.10, we notice something funny.
[13]:
from galpy import potential
# Code to create a colormap that plots NaN as gray, all else white
from matplotlib import cm
from matplotlib.colors import LinearSegmentedColormap
cdict1 = {'red': ((0.0, 1.0, 1.0),
(1.0, 1.0, 1.0)),
'green': ((0.0, 1.0, 1.0),
(1.0, 1.0, 1.0)),
'blue': ((0.0, 1.0, 1.0),
(1.0, 1.0, 1.0))}
all_white= LinearSegmentedColormap('BlueRed1', cdict1)
all_white.set_bad('0.5')
figure(figsize=(10,4))
p= lp.plotDensity(justcontours=False,rmin=-1.5,rmax=1.5,
nrs=100,nzs=100,log=True,aspect=1.,gcf=True)
p.set_cmap(all_white)
Rs= numpy.linspace(-0.3,0.3,101)
line, = plot(Rs,Rs/numpy.sqrt(1./0.7**2.-2.),lw=4.)
plot(Rs,-Rs/numpy.sqrt(1./0.7**2.-2.),color=line.get_color(),lw=4.)
gca().set_aspect('equal')
xlabel(r'$R$')
ylabel(r'$z$');
Figure 7.10: The density of a strongly-flattened logarithmic potential.
Figure 7.10 shows contours of the logarithm of the density in the \((R,z)\) plane and we have set the colormap such that all values are white, except for where the density is negative. We see that there is a large vertical region around \(R \approx 0\) where the density is negative; this region is bounded by the blue curves (see below). In the positive-density region, the density is indeed flattened, although it has a somewhat odd dumbbell shape (compare again to the Miyamoto-Nagai density contours in Figure 7.6 above).
We see that even with a mild flattening of \(q=0.7\), the flattened logarithmic potential defined in Equation (7.18) has regions with negative density. In those regions, the potential can therefore not be considered a physical model. This is not to say that the model with \(q=0.7\) is completely useless: as long as we stay far away from the negative density region, the potential can still provide a reasonable galactic model.
Let’s take a closer look at the density for the flattened logarithmic potential. The density is given by \begin{equation} \rho(R,z) = \frac{v_c^2}{4\pi\,G\,q^2}\,\frac{R^2+2\,z^2-z^2/q^2}{\left(R^2+z^2/q^2\right)^2}\,. \end{equation} This density is negative when \(R^2 < z^2\,(1/q^2-2)\). When \(q \geq 1/\sqrt{2}\) we have that \((1/q^2-2) \leq 0\) and the density is nowhere negative. When \(q < 1/\sqrt{2}\), the density is negative in regions with \(z^2 > R^2/(1/q^2-2)\); this includes the entire \(z\) range at \(R=0\) and is bounded by the curves \(z = \pm R/\sqrt{1/q^2-2}\). These curves are displayed as the blue lines in the plot of the density above.
If we approximate the density contours in \((R,z)\) with ellipses (clearly not the greatest approximation if one looks at the density contours in Figure 7.10), the axis ratio of these ellipses is \(q_\rho\) with \begin{equation}\label{eq-diskpot-densityflattening-vs-potentialflattening} q_\rho^2 = q^4\,\left(2-\frac{1}{q^2}\right)\,. \end{equation} For mildly flattened distributions, we have \(1-q \approx (1-q_\rho)/3\). That is, the deviation from sphericity in the density is about three times larger than that in the potential. This is a useful rule of thumb to relate the flattening in the density and that in the potential for any potential, but don’t expect exact results from this, deviations from this particular rule can be tens of percent!
.
.
.