12.2. Gravitational potentials for mildly-flattened and triaxial mass distributions¶
Calculating the gravitational potential for a general, triaxial mass distribution is a difficult problem. In this section, we first tackle the simpler problem of how to compute the potential for an oblate spheroidal (axisymmetric) shell and generalize it to a mass distribution where the density is constant on oblate spheroids—the shape obtained when rotating an ellipse around its minor axis. We then discuss the more general case where the density is constant on (triaxial) ellipsoids. The three qualitatively different geometries that we consider are shown in Figure 12.7.
[6]:
from mpl_toolkits.mplot3d import Axes3D
# we fiddle with the aspect ratio here
# to make the axes line up more tightly
# (adjusted for at the bottom)
stretch= 1.5
fig= figure(figsize=(11,11))
coeffs_sets = [(1.,1.,0.5),
(1.,1.,2.),
(1.,2.,.75)]# Coefficients in (x/a)**2 + (y/b)**2 + (z/c)**2 = 1
labels= ['oblate','prolate','triaxial']
# Spherical angle grid for parametric ellipsoids
phi= numpy.linspace(0.,2.*numpy.pi,101)
theta= numpy.linspace(0.,numpy.pi,101)
for ii,coeffs in enumerate(coeffs_sets):
ax= fig.add_subplot(1,3,ii+1,projection='3d')
# Radii in x,y,z
rx, ry, rz= coeffs
# Convert to cartesian coordinates
x= rx*numpy.outer(numpy.cos(phi),numpy.sin(theta))
y= ry*numpy.outer(numpy.sin(phi),numpy.sin(theta))
z= rz*numpy.outer(numpy.ones_like(phi),numpy.cos(theta))
# Plot as wireframe
ax.plot_wireframe(x, y, z,rstride=5,cstride=6,lw=1.,
color=rcParams['axes.prop_cycle'].by_key()['color'][ii])
maxr= max([rx,ry,rz])
for axis in 'xy':
getattr(ax,f'set_{axis}lim')((-maxr/stretch,maxr/stretch))
getattr(ax,'set_zlim')((-maxr/stretch,maxr/stretch))
ax.view_init(15+ii*10,20.)
ax._axis3don = False
annotate(rf'$\mathrm{{{labels[ii]}}}$',
(0.5,0.05),xycoords='axes fraction',
horizontalalignment='center',size=28.)
tight_layout();
Figure 12.7: The basic ellipsoidal shapes.
They are (a) oblate or flattened spheroids, (b) prolate or elongated spheroids, and (c) triaxial ellipsoids.
12.2.1. The potential of a spheroidal shell¶
We start by computing the gravitational potential of mass distributions that are axisymmetric and constant on ellipses in the \((R,z)\) plane that all have the same center and the same axis ratio. The density is then a function only of the quantity \(m\) that is related to \((R,z)\) as \begin{equation}\label{eq-oblate-prolate-mq} m^2 = R^2 + \frac{z^2}{q^2}\,, \end{equation} with \(q > 0\). When \(q < 1\), this is the equation of an ellipse with its major axis along \(z=0\) and with eccentricity \(e = \sqrt{1-q^2}\); the three-dimensional surface corresponding to this is an oblate spheroid. When \(q > 1\) we get a prolate spheroid instead, with Equation (12.8) the equation of an ellipse with its major axis along \(R=0\) and eccentricity \(e = \sqrt{1-1/q^2}\).
Oblate and prolate spheroids require a similar, but different treatment, and we start by computing the gravitational potential of a thin oblate shell of matter with total mass \(M_\mathrm{shell}\) and \(\rho \propto \delta(m-m_0)\) with \(m_0\) the location of the shell and where \(m\) and \(m_0\) are computed with the same, constant \(q < 1\). Because the Laplace equation \(\nabla^2 \Phi =0\) can be solved using the separation-of-variables technique in oblate spheroidal coordinates (see below) in which, for the correct definition of the coordinate system, we can also express the density as being a function of only a single coordinate, we switch to this coordinate system. Oblate spheroidal coordinates are defined by a focal length \(\Delta\) and form a transformation of the \((R,z)\) coordinates of the cylindrical coordinate system (the azimuthal angle \(\phi\) is unchanged) \begin{align}\label{eq-oblate-coords} R & = \Delta\, \cosh u\,\sin v\,;\quad z = \Delta\, \sinh u\,\cos v\,, \end{align} where \(u \geq 0\) and \(0 \leq v \leq \pi\) and \(\Delta\) is the constant focal length. For completeness, the full transformation between oblate spheroidal and Cartesian coordinates is given by \begin{align} x & = \Delta\, \cosh u\,\sin v\,\cos \phi\,;\quad y = \Delta\, \cosh u\,\sin v\,\sin \phi\,;\quad z = \Delta\, \sinh u\,\cos v\,. \end{align} In the \((R,z)\) plane, curves of constant \(u\) are (half-) ellipses centered on the origin with focus at \((R,z) = (\Delta,0)\) and eccentricity \(1/\cosh u\), while curves of constant \(v\) are hyperbolae which are the normals to the ellipses. This is illustrated in the left panel of Figure 12.8, which displays curves of constant \(u\) as solid lines and curves of constant \(v\) as dashed lines, with the focus \((\Delta,0)\) shown as a large dot.
[7]:
figure(figsize=(10,5))
from galpy.util import coords
delta= 1.5 # Focal length
# Oblate spheroidal
subplot(1,2,1)
# First plot curves of constant u
us= numpy.linspace(0.3,1.75,5)
vs= numpy.linspace(0.,numpy.pi,101)
for tu in us:
Rs,zs= coords.uv_to_Rz(tu*numpy.ones_like(vs),vs,
delta=delta,oblate=True)
plot(Rs,zs,'k-')
# Plot the focus
plot([delta],[0.],'o',ms=10.)
# Also plot a shell with m_0 = 2
m0= 2.
u0,v0= coords.Rz_to_uv(m0,0.,delta=delta,oblate=True)
Rs,zs= coords.uv_to_Rz(u0*numpy.ones_like(vs),vs,
delta=delta,oblate=True)
plot(Rs,zs,'-',lw=3.)
# Then plot curves of constant v
us= numpy.linspace(0.,10.,101)
vs= numpy.linspace(0.,numpy.pi,11)
for tv in vs:
Rs,zs= coords.uv_to_Rz(us,tv*numpy.ones_like(us),
delta=delta,oblate=True)
plot(Rs,zs,'k--')
xlim(0.,5.)
ylim(-2.5,2.5)
gca().set_aspect('equal')
xlabel(r'$R$')
ylabel(r'$z$')
# Prolate spheroidal
subplot(1,2,2)
# First plot curves of constant u
us= numpy.linspace(0.3,1.75,5)*1.3
vs= numpy.linspace(0.,numpy.pi,101)
for tu in us:
Rs,zs= coords.uv_to_Rz(tu*numpy.ones_like(vs),vs,
delta=delta)
plot(Rs,zs,'k-')
# Plot the foci
plot([0.,0.],[delta,-delta],'o',ms=10.)
# Also plot a shell with m_0 = 2
m0= 2.
u0,v0= coords.Rz_to_uv(m0,0.,delta=delta)
Rs,zs= coords.uv_to_Rz(u0*numpy.ones_like(vs),vs,
delta=delta)
plot(Rs,zs,'-',lw=3.)
# Then plot curves of constant v
us= numpy.linspace(0.,10.,101)
vs= numpy.linspace(0.,numpy.pi,11)
for tv in vs:
Rs,zs= coords.uv_to_Rz(us,tv*numpy.ones_like(us),
delta=delta)
plot(Rs,zs,'k--')
gca().set_aspect('equal')
xlim(0.,10.)
ylim(-5.,5.)
xlabel(r'$R$')
ylabel(r'$z$')
tight_layout(w_pad=-0.1);
Figure 12.8: Oblate (left) and prolate (right) spheroidal coordinates.
A shell with \(m = m_0\) corresponds to a shell with constant \(u = u_0\) if we set \begin{equation}\label{eq-oblate-shell-focus} \Delta^2 = m_0^2\,(1-q^2) = m_0^2\,e^2\,. \end{equation} The coordinate \(u_0\) of the shell can be found from \(\tanh u_0 = q\) or from \(\cosh u_0 = m_0/\Delta = 1/e\). As an example, the thick shell highlighted in the left panel of Figure 12.8 has \(m_0 = 2\) and \(q \approx 0.66\) with \(\Delta = 1.5\).
Thus, in the oblate spheroidal coordinate system with focus set by Equation (12.11), the shell with constant density \(\rho(m) \propto \delta(m-m_0)\) and axis ratio \(q\) corresponds to a density \(\rho(u,\phi,v) \propto \delta(u-u_0)\). Inside and outside of the shell with \(m=m_0\), the Poisson equation becomes the Laplace equation \(\nabla^2 \Phi = 0\) and we can solve for the potential of the shell by solving the Laplace equation inside and outside of the shell and connecting the two solutions at the location of the shell. As already mentioned above, the oblate spheroidal coordinate system is one of a small number of coordinate systems for which the Laplace equation can be solved using the method of separation of variables. In this method, solutions to the Laplace equation of the form \(\Phi_{mn}(u,\phi,v) = U_{mn}(u)\,\Pi_n(\phi)\,V_{mn}(v)\) can be found, the functions \(U_{mn}\), \(\Pi_n\), and \(V_{mn}\) are orthogonal (when integrating over the full volume) and form a basis for all smooth functions, and general solutions can be expressed as the sum over the basis functions \(\Phi_{mn}(u,\phi,v)\). A direct consequence of this is that if the density of the thin shell that forms the boundary condition for solving the Laplace equation is only a function of a single coordinate, in this case \(u\), then the solution of the Laplace equation can only depend on the same coordinate. Thus, the potential \(\Phi\) corresponding to the shell with \(\rho(u,\phi,v) \propto \delta(u-u_0)\) can only be a function of \(u\): \(\Phi \equiv \Phi(u)\).
Before working out what the potential \(\Phi(u)\) of the shell is, we can determine the behavior of the potential at \(r \gg m_0\) simply by considering the properties of the oblate spheroidal coordinate system. As we have discussed above, contours of constant \(u\) in the \((R,z)\) plane are ellipses with the focus at \((R,z) = (\Delta,0)\) and axis ratio \(\tanh u\). As is clear from the above figure and from the behavior of the hyperbolic tangent function at large values of its argument (\(\tanh u \rightarrow 1\) as \(u \rightarrow \infty\); e.g., \(\tanh(2.65) = 0.99\)), these ellipses becomes more and more circular as \(u\) increases. This implies that the potential far from the shell becomes spherical. Because the shell has a finite mass \(M_\mathrm{shell}\), the radial behavior of the potential approaches that of a point-mass potential, \(\Phi(r) = -GM/r\) at large \(r\). Thus, without doing any calculation, we have determined the asymptotic behavior of the potential of an oblate shell.
To work out the Laplace equation and the gravitational field in the oblate spheroidal coordinate system, we need the expressions for the gradient \(\nabla\) and Laplacian \(\nabla^2\) in this coordinate system. The gradient is given by \begin{equation} \nabla = \frac{1}{\Delta\,\sqrt{\sinh^2 u + \cos^2 v}}\,\left(\vec{\hat{e}}_u\,\frac{\partial}{\partial u} +\vec{\hat{e}}_v\,\frac{\partial}{\partial v} \right) + \frac{1}{\Delta\,\cosh u\,\sin v}\,\vec{\hat{e}}_\phi\frac{\partial}{\partial \phi}\,, \end{equation} and the Laplacian is given by \begin{align} \nabla^2 = \frac{1}{\Delta^2\,\left(\sinh^2 u + \cos^2 v\right)} & \,\left(\frac{1}{\cosh u}\,\frac{\partial}{\partial u}\,\left[\cosh u\,\frac{\partial}{\partial u}\right] + \frac{1}{\sin v}\,\frac{\partial}{\partial v}\,\left[\sin v\,\frac{\partial}{\partial v}\right]\right) \nonumber\\ & + \frac{1}{\Delta^2\,\cosh^2u\,\sin^2 v}\,\frac{\partial^2}{\partial \phi^2}\,. \end{align}
Because the potential only depends on \(u\) for a shell with constant density in \((\phi,v)\), the Laplace equation for the potential outside of the shell luckily simplifies significantly, because we only need to keep derivatives with respect to \(u\). After dropping an unnecessary pre-factor, we are left with \begin{equation} \frac{1}{\cosh u}\,\frac{\mathrm{d}}{\mathrm{d} u}\,\left[\cosh u\,\frac{\mathrm{d} \Phi}{\mathrm{d} u}\right] = 0\,. \end{equation} The general solution to this equation is \begin{equation} \Phi(u) = A\,\sin^{-1}\left(\frac{1}{\cosh u}\right) + B\,. \end{equation} Therefore, the potential of the shell is given by \begin{equation} \Phi(u) = \begin{cases}A_\mathrm{in}\,\sin^{-1}\left(\frac{1}{\cosh u}\right) + B_\mathrm{in} & u \leq u_0\\ A_\mathrm{out}\,\sin^{-1}\left(\frac{1}{\cosh u}\right) + B_\mathrm{out} & u > u_0\\\end{cases}\,. \end{equation} This general solution has four constants \((A_\mathrm{in},B_\mathrm{in},A_\mathrm{out},B_\mathrm{out})\) that we need to determine from boundary conditions at \(u=\infty\), \(u=0\), and \(u=u_0\). We set \(B_\mathrm{out}\) such that \(\lim_{u\rightarrow\infty} \Phi(u) = 0\); because \(\sin^{-1}(1/\cosh u)\) goes to zero as \(u\rightarrow\infty\), this means that we set \(B_\mathrm{out} = 0\). The \(u\) component of the gravitational field is given by \begin{equation}\label{eq-oblate-shell-force} -\nabla_u\Phi = -\frac{A_\mathrm{in}}{\Delta\,\cosh u\,\sqrt{\sinh^2 u + \cos^2v}}\,. \end{equation} At the center of the shell, \((u,v) = (0,0)\), this field is equal to \(-A_\mathrm{in}/\Delta\). However, from symmetry we know that the field at the center of the shell must vanish and therefore \(A_\mathrm{in} = 0\). Continuity of the potential at \(u=u_0\) then requires \(B_\mathrm{in} = A_\mathrm{out}\sin^{-1}\left(1/\cosh u\right)\). Similar to how we determined a relation between the surface density and the gradient of the potential for a razor-thin disk in Chapter 7.2.1 by integrating the Poisson equation over a small volume around the disk and using the divergence theorem from Equation (B.5), we can integrate the Poisson equation over a volume that encompasses the entire shell in \((\phi,v)\) in a small range of \(u\) around \(u_0\). To do this integration, we need the surface element in \((\phi,v)\) which is \(\Delta^2\,\cosh u\,\sin v\,\sqrt{\sinh^2u+\cos^2v}\). This leads to a relation between the surface integral over the gradient given in Equation (12.17) and the total mass \begin{align} 4\pi\,G\,M_\mathrm{shell} = -\int_0^{2\pi} \mathrm{d}\phi\, \int_0^\pi\mathrm{d}v\,& \Delta^2\,\cosh u_0\,\sin v\,\sqrt{\sinh^2u_0+\cos^2v}\nonumber\\ & \times \frac{A}{\Delta\,\sqrt{\sinh^2 u_0 + \cos^2v}\,\cosh u_0} \,, \end{align} which gives the following relation between \(A\) and \(M_\mathrm{shell}\): \(A = -GM_\mathrm{shell}/\Delta\). The gravitational potential of an oblate shell of mass \(M_\mathrm{shell}\) is therefore \begin{equation}\label{eq-triaxialgrav-potoblateshell} \Phi(u) = -\frac{G\,M_\mathrm{shell}}{\Delta}\,\begin{cases}\sin^{-1}\left(\frac{1}{\cosh u_0}\right) & u \leq u_0\\ \sin^{-1}\left(\frac{1}{\cosh u}\right) & u > u_0\\\end{cases}\,. \end{equation} As \(u \rightarrow \infty\), the transformations in Equations (12.9) approach \(R = \Delta \cosh u\sin v\) and \(z = \Delta \cosh u\cos v\), because \(\tanh u \rightarrow 1\). The oblate coordinate system approaches the spherical coordinate system, with \(u\) a different representation of the radial coordinate \(r\): \(\Delta \cosh u = r\). Because \(\sin^{-1}(1/\cosh u) \rightarrow 1/\cosh u\) when \(u \rightarrow \infty\), the potential then approaches \(\Phi(u) \rightarrow -GM_\mathrm{shell}/\Delta/\cosh u =-GM_\mathrm{shell}/r\), the point-mass potential.
Figure 12.9 shows the potential \(\Phi(u)\) for the shell highlighted in the left panel of Figure 12.8 as a function of \(u\). This potential is compared to the potential along \(z=0\) for a spherical shell with radius \(r_0 = m_0\) and one with radius \(r_0 = q\,m_0\) (which are the spherical shells that touch the spheroidal shell on the inside and outside).
[8]:
figure(figsize=(7,5))
xs= numpy.linspace(0.,10.,1001)
delta= 1.5
m0= 2.
u0,v0= coords.Rz_to_uv(m0,0.,delta=delta,oblate=True)
q= numpy.tanh(u0)
us= numpy.linspace(0.,u0,101)
line_oblate= plot(us,
-1/delta*numpy.arcsin(1./numpy.cosh(u0*numpy.ones_like(us))),
label=r'$\mathrm{Oblate\ shell}$',zorder=2)
us= numpy.linspace(u0,5.,101)
plot(us,-1/delta*numpy.arcsin(1./numpy.cosh(us)),
color=line_oblate[0].get_color())
# If spherical with r_0 = m_0
us= numpy.linspace(0.,u0,101)
line_sphere= plot(us,-1./delta/numpy.cosh(u0)*numpy.ones_like(us),zorder=1,
label=r'$\mathrm{Spherical\ shell\ with}\ r_0 = m_0$')
us= numpy.linspace(u0,5.,101)
plot(us,-1./delta/numpy.cosh(us),color=line_sphere[0].get_color())
# If spherical with r_0 = m_0*q
us= numpy.linspace(0.,u0,101)
line_sphere2= plot(us,-1./q/delta/numpy.cosh(u0)*numpy.ones_like(us),
zorder=0,
label=r'$\mathrm{Spherical\ shell\ with}\ r_0 = q\,m_0$')
us= numpy.linspace(u0,5.,101)
plot(us,-1./delta/q/numpy.cosh(us),color=line_sphere2[0].get_color())
ylim(-1.,0.05)
xlabel(r'$u$')
ylabel(r'$\Phi$')
legend(handles=[line_oblate[0],line_sphere[0],line_sphere2[0]],
fontsize=18.,loc='lower right',frameon=False);
Figure 12.9: The potential of an oblate shell compared to that of a spherical shell inside or outside the oblate shell.
The oblate shell creates a potential well that is somewhere between that of the enclosing spherical shells with the same mass: shallower than the inscribed spherical shell, but deeper than that of the spherical shell that is outside of the spheroidal shell. At \(u \gg u_0\), all of the potentials converge to the same \(1/r\) behavior.
The potential of a prolate shell of matter (Equation 12.8 with \(q>1\)) can be obtained using a similar procedure as that we followed for an oblate shell, but working instead in prolate spheroidal coordinates, which are defined in a similar manner to the oblate spheroidal coordinates using Equations (9.42). In the \((R,z)\) plane, curves of constant \(u\) are again (half-) ellipses centered on the origin, but they now have two foci \((R,z) = (0,\pm\Delta)\) and their major axis lies along \(R=0\); curves of constant \(v\) are again hyperbolae which are the normals to the ellipses. This is illustrated in the right panel of Figure 12.8, which displays contours of constant \(u\) and \(v\). Similar to the oblate shell, a prolate shell’s gravitational potential is constant inside of the shell. We do not derive the full expression for the potential of a prolate shell here.
12.2.2. Potentials of ellipsoidal systems¶
Just like we derived Equation (2.21) for the potential for any spherical density distribution by splitting the density into shells and adding the contribution from all of these shells, we can derive the potential for an spheroidal density distribution that is constant on similar ellipsoids (ellipses with the same axis ratio \(q\)) such that the density is a function of \(m\): \(\rho \equiv \rho(m)\). We will work this out for an oblate spheroid before generalizing to the general triaxial case. This is not quite as simple as it was for spherical densities, because the foci of each of these similar ellipsoids differ and this needs to be taken into account. The first order of business is to split the density \(\rho(m)\) into a set of infinitesimal spheroids with mass \(\mathrm{d}M(m)\), such that we can obtain the potential by summing over the potentials of all these thin shells. The volume of an oblate spheroid with major axis \(m\) and flattening \(q\) is \(V = 4\pi\,m^3\,q/3\)—which reduces to the volume of a sphere for \(q=1\)—so the amount of mass \(\mathrm{d}M(m)\) that lies between the shells with coordinates \(m\) and \(m+\mathrm{d}m\) is \begin{align}\label{eq-triaxialgrav-dM} \mathrm{d} M(m) & = 4\pi\,\rho(m)\,m^2\,q\,\mathrm{d}m\,. \end{align} We then have to add the potential from a sequence of thin shells with masses \(\mathrm{d}M(m)\) and we first re-write the potential of a thin shell from Equation (12.19) to be a function of its semi-major axis \(m\) as much as possible. For a shell with mass \(\mathrm{d} M(m)\) we have that \begin{equation}\label{eq-triaxialgrav-potoblateshell2} \Phi(u) = -\frac{G\,\mathrm{d} M(m)}{m e}\,\begin{cases}\sin^{-1}e & u \leq u_0\\ \sin^{-1}\left(\frac{1}{\cosh u}\right) & u > u_0\\\end{cases}\,. \end{equation}
To obtain the potential at a point \((R,z)\) that corresponds to the shell indexed by \(\tilde{m} = \sqrt{R^2+z^2/q^2}\), we sum over the contributions to the potential interior and exterior to this shell. The contribution from shells exterior to the shell with \(m = \tilde{m}\) is from Equations (12.20) and (12.21) \begin{equation} \Phi_{\mathrm{ext}} = -4\pi G\,{q\over e}\,\,\sin^{-1}e\,\int_{\tilde{m}}^\infty\mathrm{d}m\,m\,\rho(m)\,. \end{equation} As will become clear below, it is convenient to write this in terms of the function \(\psi(m)\) defined as \begin{equation}\label{eq-psi-triaxial} \psi(m) = 2\int_0^{m}\mathrm{d}\,m\,m\,\rho(m)\,, \end{equation} for which \(\Phi_{\mathrm{ext}}\) becomes \begin{equation}\label{eq-triaxialgrav-phiext} \Phi_{\mathrm{ext}} = -2\pi G\,{q\over e}\,\,\sin^{-1}e\,\left[\psi(\infty)-\psi(\tilde{m})\right]\,. \end{equation} Similarly, the contribution from shells interior to the shell with \(m=\tilde m\) is \begin{equation}\label{eq-triaxialgrav-phiint-1} \Phi_{\mathrm{int}} = -4\pi G\,{q\over e}\,\int_0^{\tilde{m}}\mathrm{d}m\,m\,\rho(m)\,\sin^{-1}\left(\frac{1}{\cosh u}\right)\,. \end{equation} In this equation, \(u\) is a function of \((m,R,z)\) through the requirement that \(u\) is the coordinate of \((R,z)\) in the oblate coordinate frame defined using \(\Delta = m e\) (and which is therefore different for each \(m\) in the integral). From Equations (12.9), it is clear that \begin{equation}\label{eq-triaxialgrav-umRz} {R^2 \over m^2 e^2\,\cosh^2 u} + {z^2 \over m^2 e^2\,\sinh^2 u} = \sin^2 v + \cos^2 v = 1\, \end{equation} To get rid of the inverse sine function in Equation (12.25), we can integrate by parts and use the properties of hyperbolic functions to obtain \begin{equation}\label{eq-triaxialgrav-phiint-2} \Phi_{\mathrm{int}} = -2\pi G\,{q\over e}\,\left[\left\{\psi(m)\,\sin^{-1}\left({1 \over \cosh u}\right)\right\}\Bigg|_{m=0}^{\tilde m}+\int_{m=0}^{\tilde{m}}\mathrm{d}\left(\sinh u\right)\,{\psi(m) \over 1+\sinh^2 u}\right]\,. \end{equation} From Equation (12.26), it is clear that \(m=0\) requires \(1/\cosh u = 0\) and \(\sinh u = \infty\), while \(m=\tilde{m}\) has \(1/\cosh u = e\) and \(\sinh u = q/e = \sqrt{1-e^2}/e\), such that we can simplify the first term to \(-2\pi G q \psi(\tilde{m})\,\sin^{-1}e / e\) and the integration limits become \(\sqrt{1-e^2}/e\) to \(\infty\). The integration variable \(m\) is implicitly defined by Equation (12.26) written solely in terms of \(\sinh u\) as \(R^2/\left(1+\sinh^2 u\right)+z^2/\sinh^2 u = m^2 e^2\). Adding up the contribution from shells exterior to \((R,z)\) (Equation 12.24) and that from shells interior, we therefore get the potential \begin{equation}\label{eq-triaxialgrav-phi-spheroid} \Phi(R,z) = -2\pi G\,{q\over e}\,\left[\psi(\infty)\,\sin^{-1}e-\int_{\sqrt{1-e^2}/e}^{\infty}\mathrm{d}\left(\sinh u\right)\,{\psi(m) \over 1+\sinh^2 u}\right]\,, \end{equation} Standard practice is to re-write the integral further by defining a new integration variable \(\tau\) defined as \begin{equation} \tau = a_0^2 e^2\left(\sinh^2 u - {1-e^2\over e^2}\right)\,, \end{equation} where \(a_0\) is an arbitrary non-zero constant. Then \(\tau\) is a function of \((m,R,z)\) through \begin{equation}\label{eq-triaxialgrav-mtau-Rz-oblate} \frac{R^2}{\tau+a_0^2}+\frac{z^2}{\tau+q^2\,a_0^2} = \frac{m^2}{a_0^2}\,, \end{equation} and the potential becomes \begin{equation}\label{eq-triaxialgrav-pot-oblate} \Phi(R,z) = -2\pi\,G\,\frac{q}{e}\,\left(\psi(\infty)\,\sin^{-1} e - \frac{a_0\,e}{2}\,\int_0^\infty\mathrm{d}\tau\,\frac{\psi(m)}{(\tau+a_0^2)\,\sqrt{\tau+q^2\,a_0^2}}\right)\,, \end{equation} As \((R,z) \rightarrow \infty\), this potential approaches zero (because \(\int_0^\infty\mathrm{d}\tau \big/\left[(\tau+a_0^2)\,\sqrt{\tau+q^2\,a_0^2}\right] = 2\,\sin^{-1} e / [a_0\,e]\)).
The right-hand side of Equation (12.31) is a function of \((R,z)\) through Equation (12.30). The \(\psi(m)\) function in Equation (12.23) can be computed analytically for many density profiles of interest and the potential can then be obtained efficiently using a one-dimensional quadrature. Moreover, expressions for the gravitational field, which we will give below after generalizing this result to triaxial ellipsoids, only involve the derivative of \(\psi(m)\), which is simply proportional to the density and the gravitational field can thus be computed using a one-dimensional quadrature for any spheroidal mass distribution.
Ellipsoidal systems are the triaxial generalization of the oblate and prolate spheroids that we have considered above. They have densities that are constant on similar ellipsoids defined by \begin{equation}\label{eq-triaxial-m} m^2 = \frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{c^2}\,. \end{equation} Note that this definition of the coordinate \(m\) is not the direct generalization of the definition of \(m\) in Equation (12.8) that we used for spheroidal systems, but is instead a generalization that treats each coordinate axis equivalently. The parameters \((a,b,c)\) can either have dimensions of length, making \(m\) dimensionless—the convention we will generally use—or \((a,b,c)\) can be dimensionless, in which case \(a=b=1\) and \(c=q\) reduces Equation (12.32) to Equation (12.8). For a triaxial density profile characterized by a scale parameter (e.g., the NFW profile), we can set \(a\) equal to the scale parameter.
The systems of oblate and prolate spheroidal coordinates can be generalized to the ellipsoidal coordinate system \((\lambda,\mu,\nu)\) and a density \(\rho \equiv \rho(m)\) can be expressed as a density \(\rho \equiv \rho(\lambda)\). The Laplace equation can be solved using the method of separation of variables in the general ellipsoidal coordinate system. Thus, it again immediately follows that the potential of a thin ellipsoidal shell \(\rho \propto \delta(\lambda - \lambda_0)\) is a function of \(\lambda\) only and we can derive the potential of this shell using similar methods as those used for the oblate spheroidal shell above. We do not go through this level of detail here, but Chapter 13.4.1 discusses ellipsoidal coordinates in more detail.
We can calculate the gravitational potential for a density that is stratified on similar ellipsoids, \(\rho \equiv \rho(m)\), in a similar way as how this is done for the more restricted oblate case above. The math to do this, however, is far more difficult and requires a specialized framework that is not all that useful in other contexts. Therefore, we will simply state the result (e.g., Chandrasekhar 1969) without proof. We again use the function \(\psi(m)\) defined in Equation (12.23). The potential at a point \((x,y,z)\) is then given by \begin{equation}\label{eq-triaxialgrav-pot-triaxial} \Phi(x,y,z) = -\pi\,G\,a\,b\,c\,\int_0^\infty\mathrm{d}\tau\,\frac{\psi(\infty)-\psi(m)}{\sqrt{(\tau+a^2)(\tau+b^2)(\tau+c^2)}}\,, \end{equation} where the relation between \(\tau\) and \(m\) is now given by \begin{equation} m^2 = \frac{x^2}{\tau+a^2}+\frac{y^2}{\tau+b^2}+\frac{z^2}{\tau+c^2}\,. \end{equation} In the oblate case (\(a=b > c\)), Equation (12.33) reduces to Equation (12.31) that we derived above. As an example of these formulae, we can compare the gravitational potential for oblate, spherical, and prolate versions of the Milky-Way-like NFW halo that we used in Chapter 5.4.2, along \(z=0\), in Figure 12.10.
[9]:
from galpy.potential import TriaxialNFWPotential, NFWPotential
from matplotlib.ticker import FuncFormatter
# Base spherical NFW Potential
np= NFWPotential(mvir=1,conc=8.37,
H=67.1,overdens=200.,wrtcrit=True,
ro=8.,vo=220.)
# We will need to define the oblate versions using the amp parameter
# of the NFW profile in galpy, so we need to figure that out in
# internal coordinates
base_amp= 16.*numpy.pi*np.a**3.*np.dens(np.a,0)*np._ro**3.*u.kpc**3.
# Now we set up oblate and prolate versions with the same enclosed mass profile
oblate_q= 0.5
oblate_np= TriaxialNFWPotential(amp=base_amp/oblate_q,
a=np.a,c=oblate_q,ro=8.,vo=220.)
prolate_q= 2.
prolate_np= TriaxialNFWPotential(amp=base_amp/prolate_q,
a=np.a,c=prolate_q,ro=8.,vo=220.)
rs= numpy.geomspace(0.1*u.kpc,1000.*u.kpc,101)
figure(figsize=(13,5))
subplot(1,2,1)
line_oblate= semilogx(rs,u.Quantity([oblate_np(r,0.) for r in rs]),
label=rf'$\mathrm{{Oblate\ NFW\ w/}}\ q={oblate_q}$',zorder=2)
line_sphere= plot(rs,np(rs,0.),
label=r'$\mathrm{Spherical\ NFW}$',zorder=2)
line_prolate= plot(rs,u.Quantity([prolate_np(r,0.) for r in rs]),
label=rf'$\mathrm{{Prolate\ NFW\ w/}}\ q={prolate_q}$',zorder=2)
xlabel(r'$r\,(\mathrm{kpc})$')
ylabel(r'$\Phi(r)\,(\mathrm{km}^2\,\mathrm{s}^{-2})$')
gca().xaxis.set_major_formatter(FuncFormatter(
lambda y,pos: (r'${{:.{:1d}f}}$'.format(int(\
numpy.maximum(-numpy.log10(y),0)))).format(y)))
legend(handles=[line_oblate[0],line_sphere[0],line_prolate[0]],
fontsize=18.,loc='upper left',frameon=False);
subplot(1,2,2)
line_oblate= semilogx(rs,-u.Quantity([oblate_np(r,0.) for r in rs])/oblate_np(0,0.),
label=rf'$\mathrm{{Oblate\ NFW\ w/}}\ q={oblate_q}$',zorder=2)
line_sphere= plot(rs,-np(rs,0.)/np(0,0),
label=r'$\mathrm{Spherical\ NFW}$',zorder=2)
line_prolate= plot(rs,-u.Quantity([prolate_np(r,0.) for r in rs])/prolate_np(0,0.),
label=rf'$\mathrm{{Prolate\ NFW\ w/}}\ q={prolate_q}$',zorder=2)
xlabel(r'$r\,(\mathrm{kpc})$')
ylabel(r'$\Phi(r)/|\Phi(0)|$')
gca().xaxis.set_major_formatter(FuncFormatter(
lambda y,pos: (r'${{:.{:1d}f}}$'.format(int(\
numpy.maximum(-numpy.log10(y),0)))).format(y)))
legend(handles=[line_oblate[0],line_sphere[0],line_prolate[0]],
fontsize=18.,loc='upper left',frameon=False)
tight_layout();
Figure 12.10: The potential of oblate, spherical, and prolate NFW profiles.
We have normalized the different profiles such that they contain the same mass within a sphere or spheroid with semi-major axis \(r\). Wee see that for the same enclosed mass, a flattened (oblate) profile gives a deeper potential well, while a prolate profile gives a shallower potential well. The right panel shows that once we normalize by the central potential, the shape of all three potentials is very similar.
As always, the gravitational field \(\vec{g} = -\nabla \Phi\) and is therefore given by \begin{equation}\label{eq-triaxialgrav-field-triaxial} \vec{g} = -\pi\,G\,a\,b\,c\,\int_0^\infty\mathrm{d}\tau\,\frac{\rho(m)\,\nabla m^2}{\sqrt{(\tau+a^2)(\tau+b^2)(\tau+c^2)}}\,, \end{equation} where \begin{equation} \nabla m^2 = \frac{2x}{\tau+a^2}\,\vec{\hat{e}}_{x}+\frac{2y}{\tau+b^2}\,\vec{\hat{e}}_{y}+\frac{2z}{\tau+c^2}\vec{\hat{e}}_{z}\,. \end{equation} For axisymmetric spheroids and defining \(m\) through Equation (12.8), we obtain somewhat simpler similar expressions. E.g., for the oblate spheroids of Equation (12.31), we have that \begin{equation}\label{eq-triaxialgrav-field-oblate} \vec{g} = -\pi\,G\,q\,\int_0^\infty\mathrm{d}\tau\,\frac{\rho(m)\,\nabla m^2}{(\tau+1)\,\sqrt{\tau+q^2}}\,, \end{equation} where we have set the arbitrary constant \(a_0\) to one and we now have that \begin{equation} \nabla m^2 = \frac{2R}{\tau+1}\,\vec{\hat{e}}_{R}+\frac{2z}{\tau+q^2}\vec{\hat{e}}_{z}\,. \end{equation} These expression, in fact, applies to prolate density distributions as well. Thus, it is clear that for any density distribution that is constant on similar triaxial ellipsoids, we can compute the gravitational field using a simple one-dimensional quadrature.
From the radial component of the gravitational field at \(z=0\), we can then also determine the circular velocity curve \begin{equation}\label{eq-triaxialgrav-vcirc-oblate} v^2_c(R) = 2\pi\,G\,q\,R^2\,\int_0^\infty\mathrm{d}\tau\,\frac{\rho(R/\sqrt{\tau+1})}{(\tau+1)^2\,\sqrt{\tau+q^2}}\,, \end{equation} because \(m = R/\sqrt{\tau+1}\) in this case. We can also write this as an integral over \(m\) instead of \(\tau\) as \begin{equation}\label{eq-triaxialgrav-vcirc-oblate-2} v^2_c(R) = 4\pi\,G\,q\,\int_0^R\mathrm{d}m\,m^2\,\frac{\rho(m)}{\sqrt{R^2-m^2 e^2}}\,, \end{equation} It is easy to see that this reduces to \(v_c^2 = GM(<r)/r\) for a spherical distribution (\(e = 0\) and \(q=1\)). For mildly-flattened systems with \(e \ll 1\) we have that \(me \ll R\) and we can Taylor expand the integrand in Equation (12.40) as \begin{equation}\label{eq-triaxialgrav-vcirc-oblate-approx} v^2_c(R) \approx {GM(< R) \over R} + {2\pi\,G\,q\,e^2 \over R^3}\,\int_0^R\mathrm{d}m\,m^4\,\rho(m)\,, \end{equation} where \(M(<R)= 4\pi\,q\,\int_0^R\mathrm{d}m\,m^2\,\rho(m)\) is the mass within the spheroid with semi-major axis \(R\). We can also write this as \begin{equation}\label{eq-triaxialgrav-vcirc-oblate-approx-2} v^2_c(R) \approx {GM(< R) \over R}\left[1+{e^2 \over 2q}\,X(R)\right]\,, \end{equation} where the dimensionless function \(X(R)\) is given by \begin{equation}\label{eq-triaxialgrav-X-rotcurveflattening} X(R) = {1 \over R^2}\,{\int_0^R\mathrm{d}m\,m^4\,\rho(m) \over \int_0^R \mathrm{d}m\,m^2\,\rho(m)} \,. \end{equation} This is manifestly positive for any density distribution, so we arrive at the general result that the circular velocity at a given radius increases when we flatten any spherical mass distribution by squashing it down (i.e., conserving mass from a sphere with radius \(R\) to a spheroid of semi-major axis \(R\)). To get a better sense of \(X(R)\), we can compute it for a few example density distributions. The simplest case is for a constant density, for which \begin{equation} X(R) = {3 \over 5}\,\quad (\mathrm{constant\ density})\,, \end{equation} while for any power-law \(\rho(m) \propto m^{-\alpha}\) \begin{equation} X(R) = {3-\alpha \over 5-\alpha}\,\quad (\rho[m] \propto m^{-\alpha})\,. \end{equation} For more complex distributions such as the Hernquist profile from Chapter 2.4.6 \begin{equation} X(R) = 2\,\tilde{R}^{-4}\,(1+\tilde{R})^2\,\left[ {5 \over 2}+\tilde{R}-{5+6\tilde{R}\over 2(1+\tilde{R})^2}-3 \ln(1+\tilde{R})\right]\,\quad (\mathrm{Hernquist})\,, \end{equation} or the NFW profile \begin{equation} X(R) = {(\tilde{R}-3)\,\tilde{R}^2-6\tilde{R}+6(1+\tilde{R})\ln(1+\tilde{R})\over 2\tilde{R}^2\,\left[(1+\tilde{R})\,\ln(1+\tilde{R})-\tilde{R}\right]}\,\quad (\mathrm{NFW})\,, \end{equation} where in both expressions \(\tilde{R} \equiv R/a\), with \(a\) the profile’s scale radius. A reasonable value for many galactic density profiles is \(X(R) \approx 1/3\), so we have that approximately \begin{equation}\label{eq-triaxialgrav-vcirc-oblate-approx-3} v_c(R) \approx v_{c,\mathrm{sphere}}(R)\left(1+{e^2 \over 12q}\right)\,, \end{equation} where \(v_{c,\mathrm{sphere}}(R)\) is the circular velocity if the mass were distributed spherically and \(R\) is the coordinate along the major axis. As an example of this approximation, we can compare it to the exact calculation for a flattened NFW profile for a Milky-Way-like NFW halo at three different radii (the same example halo that we used above); this comparison is shown in Figure 12.11.
[10]:
from galpy.potential import TriaxialNFWPotential, NFWPotential
# Base spherical NFW Potential
np= NFWPotential(mvir=1,conc=8.37,
H=67.1,overdens=200.,wrtcrit=True,
ro=8.,vo=220.)
# We will need to define the oblate versions using the amp parameter
# of the NFW profile in galpy, so we need to figure that out in
# internal coordinates
base_amp= 16.*numpy.pi*np.a**3.*np.dens(np.a,0)*np._ro**3.*u.kpc**3.
es= numpy.linspace(0.,1,41)
qs= numpy.sqrt(1-es**2.)
Rs= numpy.array([0.25,1.,10.])
vcs= numpy.zeros((len(es),len(Rs)))*u.km/u.s
for ii,e in enumerate(es):
# to conserve mass, need to divide base_amp by q
tp= TriaxialNFWPotential(amp=base_amp/qs[ii],a=np.a,c=qs[ii],
ro=8.,vo=220.)
vcs[ii]= [tp.vcirc(R*np.a) for R in Rs]
figure(figsize=(11,4.75))
lines= []
for ii, color in enumerate(['#1f77b4','#ff7f0e','#2ca02c']):
subplot(1,2,1)
lines.append(plot(es,vcs[:,ii],color=color,
label=rf'$R = {Rs[ii]:g}\,a$' if ii != 1 else r'$R = a$')[0])
plot(es,np.vcirc(Rs[ii]*np.a)*(1.+es**2./12./qs),'--',color=color)
subplot(1,2,2)
plot(qs,vcs[:,ii],color=color)
plot(qs,np.vcirc(Rs[ii]*np.a)*(1.+es**2./12./qs),'--',color=color)
subplot(1,2,1)
xlabel(r'$e$')
ylabel(r'$v_c$')
legend(handles=lines,
fontsize=18.,loc='upper left',frameon=False)
subplot(1,2,2)
xlabel(r'$q$')
ylabel(r'$v_c$')
from matplotlib.lines import Line2D
legend_lines= [Line2D([0],[0],color='k',),
Line2D([0],[0],color='k',ls='--')]
legend(legend_lines,
[r'$\mathrm{exact}$',r'$v_{c,\mathrm{sphere}}(R)\left(1+{e^2 \over 12q}\right)$'],
fontsize=18.,loc='upper right',frameon=False)
tight_layout();
Figure 12.11: Effect of flattening on the rotation curve of a Milky-Way-like NFW halo.
We plot the circular velocity as a function of \(e\) in the left-hand panel and of \(q\) in the right-hand panel.
For many spheroidal density profiles, we can also calculate the circular velocity using Equation (12.40) analytically. For example, for the NFW profile (see Chapter 2.4.6) we find that \begin{equation}\label{eq-triaxialgrav-vcirc-nfw} v_c^2(R) = {4\pi G\,\rho_0\,a^3\,R \over (R^2-a^2 e^2)}\,\left[{a\,(q-1)-R\over (R+a)}+{\mathrm{tanh}^{-1}\left({R\over \sqrt{R^2-a^2e^2}}\right)-\mathrm{tanh}^{-1}\left({R+a e^2 \over q\,\sqrt{R^2-a^2 e^2}}\right) \over \sqrt{1-a^2 e^2/R^2}}\right]\,, \end{equation} where \(e^2 = 1-q^2\) for both oblate and prolate spheroids (the usual definition of the eccentricity of a prolate distribution is \(e^2 = 1-1/q^2\), such that \(0 \leq e \leq 1\)). At \(R^2 < a^2e^2\), the argument of the inverse hyperbolic tangent function is imaginary, for which we can use \(\mathrm{tanh}^{-1}(z) = \tan^{-1}(iz)/i\); the square root in the denominator of this term is also imaginary in this case and the result is a real number. Codes that implement operations on complex numbers can directly use Equation (12.49). We can use this result to compare the rotation curves of oblate and prolate NFW halos if they contain the same amount of mass as a function of \(m\), using the same Milky-Way-like NFW halo as above. This is shown in Figure 12.12.
[10]:
from galpy.potential import NFWPotential
def vcirc_NFW(R,q,amp=1.,a=1.):
"""Circular velocity for a spheroidal NFW profile"""
e2= 1.-q**2.
R2a2e2= (R**2.-a**2.*e2)*(1.+0J) # Make use of complex arithmetic!
return numpy.sqrt(4.*numpy.pi*amp*a**3*q*R/R2a2e2\
*((numpy.arctanh(R/numpy.sqrt(R2a2e2))
-numpy.arctanh((R+a*e2)/q/numpy.sqrt(R2a2e2)))*R/R2a2e2**0.5
+(a*(q-1.)-R)/(R+a)))
# Base spherical NFW Potential
np= NFWPotential(mvir=1,conc=8.37,
H=67.1,overdens=200.,wrtcrit=True,
ro=8.,vo=220.)
# We will need to define the oblate/prolate versions using the amp
# parameter of the NFW profile in galpy, so we need to figure that
# out in internal coordinates
base_amp= np.dens(np.a,0,use_physical=False)*4.
qs= [1e-16,0.1,0.5,0.9999,2.,4.]
Rs= numpy.linspace(0.01,250.,101)*u.kpc
figure(figsize=(10,5))
for q in qs:
#amp = amp/q to keep mass within spheroids constant
plot(Rs,vcirc_NFW((Rs/np._ro/u.kpc).to_value(u.dimensionless_unscaled),
q,amp=base_amp/q,a=np.a)*np._vo*u.km/u.s,
label=rf'$q = {q}$' if numpy.fabs(1.-q) > 0.1 and q > 0.01
else r'$q = 0$' if numpy.fabs(1.-q) > 0.1 else r'$\mathrm{spherical}$')
xlabel(r'$R\,(\mathrm{kpc}$)')
ylabel(r'$v_c\,(\mathrm{km\,s}^{-1})$')
xlim(0.,350.)
legend(fontsize=18.,frameon=False,loc=(0.73,0.3));
Figure 12.12: Rotation curves of oblate and prolate NFW profiles.
Thus, we see the expected result that flattening the mass distribution increases the rotation curve, while stretching the mass distribution along the \(z\) axis reduces it. We can in fact flatten the NFW profile all the way into a disk as shown by the \(q=0\) curve!
Further worked out expressions for the gravitational potential and field for the triaxial generalizations of the NFW and Hernquist profiles from Chapter 2.4.6 are given by Merritt & Fridman (1996).
12.2.3. The perfect ellipsoid¶
We have shown above that it is straightforward to compute the gravitational potential, field, and rotation curve for any ellipsoidal mass distribution using a simple one-dimensional quadrature that is easy to perform numerically. Therefore, there is little reason to use less realistic, but more tractable density distributions when studying the properties of spheroidal and ellipsoidal density distributions. However, an exception has to be made for the perfect ellipsoid (de Zeeuw 1985). Not only can the gravitational potential from Equation (12.33) be worked out in terms of incomplete elliptical integrals for any combination of axis ratios \(b/a\) and \(c/a\) and the gravitational field can therefore also be obtained analytically. But the potential is also separable in ellipsoidal coordinates in such a way that the Hamilton-Jacobi equation (Chapter 3.4.3) separates and we can compute orbital actions and angles using simple one-dimensional quadratures similar to those described in Chapter 9.4.3 (the potential is of so-called Staeckel form). This implies that all orbits in the perfect ellipsoid are regular (that is, non-chaotic; see Chapter 13) and that we can easily explore all of the orbits in the perfect ellipsoid because we can use the action labels to go through them. We will make use of this in Chapter 13.4 to explore orbits in triaxial potentials.
The perfect ellipsoid is defined by the following density distribution on similar ellipsoids: \begin{equation} \rho(m) = {M \over abc\,\pi^2}\,{1 \over (1+m^2)^2}\,, \end{equation} for the dimensionless coordinate \(m\) given by Equation (12.32), where we set the profile’s scale parameter using \(a\). The parameter \(M\) is the total mass of the system. While the gravitational potential can be worked out in terms of incomplete elliptic integrals, the expression is quite cumbersome and we won’t give it here (see de Zeeuw 1985 for full details). At \(m \ll 1\), the density is cored, while at large \(m\), the density drops as \(m \propto m^{-4}\), leading to the finite mass. Interesting edge cases can be obtained by setting one or more of the shape parameters \((a,b,c)\) to zero. For example, for \(a=b\) and \(c=0\), the perfect ellipsoid becomes the Kuzmin disk that we discussed as a simple razor-thin disk model in Chapter 7.2.1.
We can derive the rotation curve for the perfect spheroid (i.e., \(a=b\)) using Equation (12.40). First we switch to the dimensionful definition of \(m\) from Equation (12.8) and write the density as \(\rho(m) = M a/(\pi^2q[a^2+m^2]^2)\), where the scale parameter \(a\) is now explicit. The circular velocity is then given by \begin{equation} v_c^2(R) = {2M a\over \pi}\,{R^2 \over (R^2+a^2e^2)^2}\,\left[\sqrt{R^2/a^2+e^2}\,\sin^{-1}\left(\sqrt{{R^2+a^2e^2 \over R^2+a^2}}\right)-q\,{R^2+a^2e^2\over R^2+a^2}\right]\,. \end{equation} When \(q\rightarrow 0\) (and \(e \rightarrow 1\)), this becomes the circular velocity of the Kuzmin disk (Equation 7.13). Computing the rotation curve for different values of \(q\), we again see the expected trend in Figure 12.13, that more flattened distribution provide more rotational support.
[11]:
def vcirc_perfect(R,q,mass=1.,a=1.):
"""Circular velocity for a spheroidal perfect ellipsoid"""
e2= (1.-q**2.)*(1.+0J) # Make use of complex arithmetic!
return numpy.sqrt(2*mass*a/numpy.pi*R**2./(R**2.+a**2.*e2)**2.
*(numpy.sqrt(R**2./a**2.+e2)
*numpy.arcsin(numpy.sqrt((R**2.+a**2.*e2)/(R**2.+a**2.)))
-q*(R**2.+a**2.*e2)/(R**2.+a**2.)))
mass= 5
a= 1.
qs= [1e-16,0.1,0.5,0.9999,2.,4.]
Rs= numpy.linspace(0.01,15*a,101)
figure(figsize=(10,6))
for q in qs:
plot(Rs/a,vcirc_perfect(Rs,q,mass=mass,a=a),
label=rf'$q = {q}$' if numpy.fabs(1.-q) > 0.1 and q > 0.01
else r'$q = 0$' if numpy.fabs(1.-q) > 0.1 else r'$\mathrm{spherical}$')
xlabel(r'$R/a$')
ylabel(r'$v_c$')
ylim(0.,1.55)
legend(fontsize=18.,frameon=False,loc=(0.73,0.45));
Figure 12.13: Rotation curves of oblate and prolate perfect ellipsoids.
.
.
.