2.4. Examples of spherical potentials¶
Let’s now consider some example spherical potentials that illustrate the concepts above. This is by no means an exhaustive list of spherical potentials, but they are the main ones that we will use in the rest of this book and some of the main ones that often come up in research. More examples are discussed in the exercises.
2.4.1. Point-mass potential¶
The most basic spherical potential is that of a point mass. As we demonstrated in Section 2.1 above, the potential of a point mass \(M\) is \begin{equation}\label{eq-pot-kepler} \Phi(r) = -\frac{G\,M}{r}\,. \end{equation} Because of Newton’s second shell theorem, this is also the potential outside the boundary for any spherical mass distribution with a finite extent. Potentials of any shape with a total finite mass should also approach the point-mass potential as \(r \rightarrow \infty\).
The circular velocity is given by \begin{equation}\label{eq-vc-kepler} v_c(r) = \sqrt{\frac{G\,M}{r}}\,. \end{equation} The point-mass potential of the Sun is the dominant gravitational potential that sets planetary trajectories in the solar system (perturbations from other planets are very small). Because Equation (2.31) is equivalent to Kepler’s third law, first proposed by Kepler, the point-mass potential and orbits within it are often referred to as Keplerian. For example, a disk of mass \(m \ll M\) orbiting a point mass is Keplerian and is in Keplerian rotation. Keplerian orbits have special properties that we discuss in more detail in the next chapters.
In galactic dynamics, point-mass potentials are typically only used to represent the gravitational field from (supermassive) black holes or as an approximation. Point-mass potentials are not typically used in \(N\)-body simulations of galaxies as we will discuss in Chapter 12.4: the gravitational potential of individual particles in \(N\)-body simulations is smoothed to, for example, a Plummer potential (see below) to avoid the force divergence for (unphysical) close encounters; these encounters are unphysical because galaxies are collisionless. \(N\)-body simulations of stellar clusters do use the point-mass potentials for all of their constituent stars.
2.4.2. The homogeneous sphere¶
The potential of a homogeneous density sphere, \(\rho(r) = \rho_0 =\) constant for all \(r < R\) and zero at \(r \geq R\), can be calculated using (2.21); the result is \begin{equation}\label{eq-pot-homogeneous} \Phi(r) = \begin{cases} \frac{2\pi\,G\,\rho_0}{3}\,(r^2-3\,R^2)\,,& \qquad (r < R)\,,\\ -\frac{4\pi\,G\,\rho_0\,R^3}{3\,r}\,,& \qquad (r \geq R)\,.\end{cases} \end{equation} At \(r<R\), this quadratic potential is that of a harmonic oscillator; we only consider \(r < R\) in the following, because at \(r \geq R\) the behavior is the same as that of the point-mass potential above due to Newton’s second shell theorem. The circular velocity is \begin{equation} v_c(r) = \sqrt{\frac{4\pi\,G\,\rho_0}{3}}\,r\,, \end{equation} and the dynamical time using the definition from (2.27) is therefore \begin{equation}\label{eq-tdyn-homogeneous} t_\mathrm{dyn} = \sqrt{\frac{3\pi}{G\,\rho_0}}\,. \end{equation} The dynamical time in Equation (2.34) is the same as the one we derived in Equation (2.28) when substituting \(\rho_0\) for the average density. In a homogeneous sphere, the dynamical time is constant with \(r\). If instead of considering a body on a circular orbit, we look at a body moving perfectly radially, its equation of motion according to Newton’s second law is \begin{equation} \frac{\mathrm{d}^2 r}{\mathrm{d} t^2} = -\frac{4\pi\,G\,\rho_0}{3}\,r\,, \end{equation} which is the equation of motion of a simple harmonic oscillator with frequency \(\omega = \sqrt{4\pi\,G\,\rho_0/3}\). The period of this radial oscillation is \begin{equation} T = {2\pi \over \omega} = \sqrt{\frac{3\pi}{G\,\rho_0}} = t_\mathrm{dyn}\,, \end{equation} which is again constant. This result provides further justification for the definition of the dynamical time in Equation (2.27): no matter whether a body is on a circular orbit or a purely radial orbit, the period of the orbit is \(\approx t_\mathrm{dyn}\).
Because \(v_c(r) \propto r\) for a homogeneous sphere, a disk of bodies rotates with an angular rotation rate of \(\Omega(r) = v_c(r)/r =\) constant. Because the rotation rate is constant, such a disk rotates as a solid body and this setup is referred to as solid-body rotation.
2.4.3. Plummer sphere¶
Next we will look at a few generalizations of the Kepler point-mass potential above that often crop up in the study of galaxies. One of the most important of these is the Plummer model (Plummer 1911), which smooths the Kepler potential for a mass \(M\) by replacing the radius \(r\) in the denominator of Equation (2.30) by \(\sqrt{r^2+b^2}\) where \(b\) is a constant, the Plummer scale length, \begin{equation}\label{eq-pot-plummer} \Phi(r) = -\frac{G\,M}{\sqrt{r^2+b^2}}\,. \end{equation} This potential approaches a Kepler potential as \(b \rightarrow 0\), or, put in a way that is more practically useful, \(r \gg b\). The Plummer potential was originally introduced to represent globular clusters, but is no longer a preferred model for these systems. It remains highly useful because it has many convenient analytical properties and provides a crude model for, e.g., dark matter halos, that is easy to work with; it is also used to represent the stellar density in small dwarf galaxies. Finally, it also provides a simple kernel used to smooth the gravitational field of point particles in \(N\)-body simulations.
It is straightforward to derive the circular velocity from Equation (2.37). The density profile using the Poisson equation is \begin{equation}\label{eq-pot-plummer-density} \rho(r) = \frac{3M\,b^2}{4\pi}\,\frac{1}{\left(r^2+b^2\right)^{5/2}}\,. \end{equation}
2.4.4. Isochrone potential¶
Similar to the Plummer sphere, the isochrone potential is obtained by smoothing the potential of a point mass (Henon 1959). For the isochrone potential, this is done by replacing \(r\rightarrow b+\sqrt{r^2+b^2}\) in the denominator of the gravitational potential (see Equation 2.30), where \(b\) is again a constant \begin{equation}\label{eq-pot-isochrone} \Phi(r) = -\frac{G\,M}{b+\sqrt{r^2+b^2}}\,. \end{equation} Similar to the Plummer sphere, this potential approaches a Kepler potential when \(b \rightarrow 0\), or when \(r \gg b\). We can determine how fast the isochrone model approaches the Kepler potential by expanding the potential around \(b/r\) for \(r \gg b\) \begin{align} \Phi_{r\gg b}(r) & = -\frac{G\,M}{r\left[b/r+\sqrt{1+\left(b/r\right)^2}\right]} \approx -\frac{G\,M}{r}\,\left(1-\frac{b}{r}\right)\,. \end{align} Because the first term is the potential of a point mass, which has zero density at \(r \neq 0\), the density at large radii is given by that corresponding to the second term, \(\Phi(r) \propto r^{-2}\). From the Poisson equation, we know that the density is a second derivative of the potential, and must therefore tend to \(\rho_{r\gg b}(r) \propto r^{-4}\). This \(1/r^4\) behavior is also that of a popular model for dark-matter halos, the Hernquist model. We will discuss power-law potentials and the Hernquist model in more detail below.
When \(r/b \rightarrow 0\), we can expand the potential around \(r/b\) \begin{align} \Phi_{r\ll b}(r) & = -\frac{G\,M}{b\left[1+\sqrt{\left(r/b\right)^2+1}\right]} \approx -\frac{G\,M}{2b}\,\left(1-\frac{1}{4}\,\frac{r^2}{b^2}\right) = \frac{G\,M}{8\,b^3}\,r^2+\mathrm{constant}\,. \end{align} Up to an irrelevant constant, this is the potential of a homogeneous sphere with density \(\rho_0 = 3\,M/(16\pi\,b^3)\). This immediately tells us that the central density of the isochrone mass distribution is this \(\rho_0\). Thus, the isochrone potential smoothly interpolates between the potential of a point mass and that of a homogeneous density distribution, essentially the two extremes of the types of densities seen in stellar systems.
The circular velocity of the isochrone potential is given by \begin{equation}\label{eq-vc-isochrone} v_c^2(r) = \frac{G\,M\,r^2}{(b+a)^2\,a}\,, \end{equation} where \(a = \sqrt{r^2+b^2}\).
It is clear that the isochrone model can represent a wide range of potentials. Because it interpolates between the two extremes of density distributions, over a narrow radial range it can essentially represent all spherical potentials. But the main reason for the importance of the isochrone potential is that it is the most realistic model for galaxies in which all of the orbits can be solved analytically (Saha 1991). That is, just like for the Kepler potential, we can work out the full orbit of a body in terms of elementary functions, without requiring any numerical quadrature or numerically solving the differential equation of motion. Naturally, the analytic solution is more complicated than that for a Kepler potential (see Chapter 4.4). That all of the orbits are analytic is extremely useful in many settings in galactic dynamics.
Let’s compare the rotation curves of the spherical potentials that we have discussed so far by plotting them with galpy in Figure 2.5. We normalize the point-mass, Plummer, and isochrone potentials such that they all have \(GM = 1\), and the homogeneous sphere such that \(v_c = 1\) for \(r = 1\). We set the scale \(b\) of the isochrone and Plummer potentials to \(b=1\) and we set the size of the homogeneous sphere to be larger than the radial range that we consider such that its density is constant everywhere that we consider here.
[9]:
from galpy import potential
figure(figsize=(10,6))
# Define each potential
kp= potential.KeplerPotential(amp=1.)
hp= potential.HomogeneousSpherePotential(normalize=1.,R=20.)
pp= potential.PlummerPotential(amp=1.,b=1.)
ip= potential.IsochronePotential(amp=1.,b=1.)
# Plot the rotation curve for each
line_kepler= potential.plotRotcurve(kp,Rrange=[0.,10.],\
label=r'$\mathrm{Point\ mass}$',yrange=[0.,1.5],
xlabel=r'$r$',ylabel=r'$v_c$',gcf=True)
line_homog= potential.plotRotcurve(hp,Rrange=[0.,10.],\
label=r'$\mathrm{Homogeneous\ sphere}$',overplot=True)
line_plummer= potential.plotRotcurve(pp,Rrange=[0.,10.],\
label=r'$\mathrm{Plummer}\,(b=1)$',overplot=True)
line_isochrone= potential.plotRotcurve(ip,Rrange=[0.,10.],\
label=r'$\mathrm{Isochrone}\,(b=1)$',overplot=True)
legend(handles=[line_kepler[0],line_homog[0],
line_plummer[0],line_isochrone[0]],
fontsize=18.,loc='upper right',frameon=False);
Figure 2.5: Rotation curves of basic spherical potentials.
We see the expected behavior in these curves: the point-mass \(v_c(r)\) drops as \(r^{-1/2}\), while the other potentials all have \(v_c(r) \propto r\) at \(r \ll b\) (and everywhere for the homogeneous sphere, which quickly runs out of the shown range in \(v_c\)), because they all have an approximate constant-density core. At \(r > b\) the Plummer and isochrone potentials approach the Keplerian curve, because they all have the same total mass. The Plummer potential approaches the Keplerian curve quicker than the isochrone potential for the same \(b\).
2.4.5. Power-law models¶
An important class of models are those in which the density is a power-law of radius \begin{equation} \rho(r) = \rho_0\,\left(\frac{r_0}{r}\right)^\alpha\,, \end{equation} where \(\rho_0\) is the density at \(r=r_0\). The enclosed mass of this density profile is \begin{equation} M(< r) = \frac{4\pi\,\rho_0\,r_0^3}{3-\alpha}\,\left(\frac{r_0}{r}\right)^{\alpha-3}\,, \end{equation} which is finite at \(r=\infty\) only when \(\alpha > 3\). Using Equation (2.23), we find that the potential is \begin{equation} \Phi(r) = -\frac{4\pi\,G\,\rho_0\,r_0^2}{(3-\alpha)(2-\alpha)}\,\left(\frac{r_0}{r}\right)^{\alpha-2}\,,\qquad (\alpha \neq 2)\,, \end{equation} and \begin{equation}\label{eq-pot-powlaw-log} \Phi(r) = \frac{4\pi\,G\,\rho_0\,r_0^2}{(3-\alpha)}\,\ln r\,,\qquad (\alpha = 2)\,. \end{equation}
The circular velocity is given by \begin{equation} v^2_c(r) = \frac{4\pi\,G\,\rho_0\,r_0^2}{3-\alpha}\,\left(\frac{r_0}{r}\right)^{\alpha-2}\,. \end{equation} For \(\alpha = 2\), the circular velocity is constant with \(v_c =\sqrt{4\pi\,G\,\rho_0\,r_0^2}\) and we can write the gravitational potential in Equation (2.46) in terms of \(v_c\) \begin{equation}\label{eq-pot-powlaw-log-vc} \Phi(r) = v_c^2\,\ln r\,,\qquad (\alpha = 2)\,. \end{equation} This is an important potential that is, for obvious reasons, often referred to as the logarithmic potential. Because it gives rise to a constant circular velocity, it is a very simple and often-used model for the potential of galaxies with a flat rotation curve. When we study self-consistent dynamical equilibrium states in Chapter 5, we will see that for the logarithmic potential the velocity dispersion is constant with \(r\) and it is therefore also often referred to as the singular isothermal sphere, with the “singular” modifier due to the \(1/r^2\) density divergence as \(r\rightarrow 0\).
2.4.6. Two-power density models¶
2.4.6.1. Generalities¶
A final important class of spherical potentials for galaxies are those that are described by a different power law in their inner and outer parts, with a smooth transition region between the two. A popular set of models with this property is that with the following density law \begin{equation}\label{eq-twopower-dens} \rho(r) = \frac{\rho_0\,a^\alpha}{r^\alpha\,(1+r/a)^{\beta-\alpha}}\,. \end{equation} In the small and large \(r\) limit, these models behave as power law models \begin{equation} \rho(r) = \begin{cases} \rho_0\,\left(\frac{a}{r}\right)^\alpha, & r \ll a\,,\\ \rho_0\,\left(\frac{a}{r}\right)^\beta, & r \gg a\,. \end{cases} \end{equation} The inner slope is therefore \(\alpha\) and the outer slope is \(\beta\), with the transition occurring around \(r \approx a\), where \(a\) is known as the scale radius. For \(\beta=4\), these models have relatively simple properties and they are known as \(\gamma\) models (where \(\gamma\) is what we call \(\alpha\) above); Dehnen (1993) and Tremaine et al. (1994) provide more information on this set of models. We will mainly consider two popular instances of two-power density models: the Hernquist model (Hernquist 1990), which has \(\alpha = 1, \beta =4\), and the Navarro-Frenk-White (NFW; Navarro et al. 1997) model, with \(\alpha = 1, \beta = 3\). The enclosed mass is given by \begin{equation}\label{eq-massenc-hernnfw} M(< r) = 4\pi\,\rho_0\,a^3\,\times \begin{cases} \frac{(r/a)^2}{2\,(1+r/a)^2}, & \mathrm{(Hernquist)}\,\\ \ln(1+r/a)-\frac{r/a}{1+r/a}, & \mathrm{(NFW)}\,. \end{cases} \end{equation} Using Equation (2.23), this gives the gravitational potential \begin{equation}\label{eq-potential-hernnfw} \Phi(r) = -4\pi\,G\,\rho_0\,a^2\,\times \begin{cases} \frac{1}{2\,(1+r/a)}, & \mathrm{(Hernquist)}\,\\ a\,\ln(1+r/a)/r, & \mathrm{(NFW)}\,. \end{cases} \end{equation}
The total mass of an NFW potential does not converge as we increase \(r\), but Equation (2.51) demonstrates that the total mass of a Hernquist potential is \(M = 2\pi\,\rho_0\,a^3\). The Hernquist potential written in terms of this mass \(M\) is \begin{equation}\label{eq-sphergrav-hernquist-potential} \Phi(r) = -\frac{G\,M}{r+a}\,. \end{equation} The Hernquist potential is therefore also a simple smoothing of the Kepler point-mass potential (like the Plummer model), in this case by replacing \(r \rightarrow r+a\). As for the Plummer model, this often means that calculations involving a Hernquist potential can be analytically simplified. Another member of this family (and of the \(\gamma\) models) is that with \(\alpha=2\) and \(\beta=4\), the Jaffe profile, which is often used to represent very cuspy profiles because of its convenient analytical properties (Jaffe 1983).
The NFW profile is found to describe the density profile of dark matter halos well in cosmological simulations (Navarro et al. 1996) and remains highly important for this reason. Hernquist profiles have somewhat more tractable properties and are also a popular model for dark matter halos and for elliptical galaxies and galactic bulges, because they can approximately represent the \(R^{1/4}\) de Vaucouleurs surface-brightness profile of ellipticals and bulges (see Chapter 12.1). The outskirts of dark-matter halos at the present time have not yet fully formed and it has also been argued that the long-term state of dark-matter halos more closely resembles a Hernquist profile than an NFW profile (Busha et al. 2005).
Let’s compare the rotation curves of Hernquist and NFW potentials with the same scale radius \(a\). First, we compare the rotation curves for the case where both potentials correspond to the same enclosed mass at \(r= 12\,a\) (this is approximately the virial radius for Milky-Way-mass dark-matter halos) in the left panel of Figure 2.6. We do this by computing the enclosed mass using galpy's .mass function for each potential set up with an amplitude of unity and then setting the new amplitude such that both potentials have the same mass:
[10]:
from galpy import potential
from galpy.util import plot as galpy_plot
figure(figsize=(12,4))
subplot(1,2,1)
# Define each potential, equal mass at 12 scale radii
hp_4scale= potential.HernquistPotential(amp=1.,a=1.)
hp= potential.HernquistPotential(amp=10./hp_4scale.mass(12.),a=1.)
nfp_4scale= potential.NFWPotential(amp=1.,a=1.)
nfp= potential.NFWPotential(amp=10./nfp_4scale.mass(12.),a=1.)
# Plot the rotation curve for each
line_hernquist= potential.plotRotcurve(hp,Rrange=[0.,15.],\
label=r'$\mathrm{Hernquist}$',yrange=[0.,2.5],
xlabel=r'$r/a$',ylabel=r'$v_c$',gcf=True)
line_nfw= potential.plotRotcurve(nfp,Rrange=[0.,15.],\
label=r'$\mathrm{NFW}$',overplot=True)
legend(handles=[line_hernquist[0],line_nfw[0]],
fontsize=18.,loc='upper right',frameon=False)
galpy_plot.text(r'$\mathrm{Equal\ mass\ at}\ R=12\,a$',
top_left=True,size=18.)
subplot(1,2,2)
# Define each potential, equal inner profile
hp= potential.HernquistPotential(amp=20.,a=1.)
nfp= potential.NFWPotential(amp=20.,a=1.)
# Plot the rotation curve for each
line_hernquist= potential.plotRotcurve(hp,Rrange=[0.,15.],\
label=r'$\mathrm{Hernquist}$',yrange=[0.,2.5],
xlabel=r'$r/a$',ylabel=r'$v_c$',gcf=True)
line_nfw= potential.plotRotcurve(nfp,Rrange=[0.,15.],\
label=r'$\mathrm{NFW}$',overplot=True)
legend(handles=[line_hernquist[0],line_nfw[0]],
fontsize=18.,loc='upper right',frameon=False)
galpy_plot.text(r'$\mathrm{Equal\ inner\ density\ profile}$',
top_left=True,size=18.)
tight_layout(w_pad=-0.5)
Figure 2.6: Rotation curves of Hernquist and NFW potentials.
Because both potentials are defined to have the same mass at \(r=12\,a\), the circular velocity at these radii is the same. We see that the rotation curve of the Hernquist potential reaches a much higher maximum value than the NFW potential with the same mass, because the Hernquist sphere’s mass is much more concentrated than that of the NFW sphere. At \(r > 12\,a\), the Hernquist \(v_c(r)\) declines more steeply with radius than the NFW \(v_c(r)\), because of the \(r^{-4}\) behavior of the density at \(r \gg a\) for the Hernquist profile versus the \(r^{-3}\) behavior for the NFW density.
Next, we compare the rotation curves if we normalize the Hernquist and NFW potentials such that they have the same inner density profile \(\rho(r) \propto \rho_0\,r^{-1}\) at \(r \ll a\) in the right panel of Figure 2.6. The inner rotation curves are now the same, because the enclosed mass profile is the same at \(r \ll a\), but because the density of the NFW profile is shallower than that of the Hernquist profile, the NFW rotation curve reaches a significantly higher value than that of the Hernquist profile.
2.4.6.2. NFW profiles for dark-matter halos¶
NFW profiles are in widespread use as models of dark-matter halos and there are various other ways to parameterize these beside the \((\rho_0,a)\) parameterization used in Equation (2.49) (which is, in fact, not commonly used). The most important of these is the parameterization using the virial mass and concentration. As we discussed above, NFW models do not have a finite mass, so any mass parameter has to be defined within a specified radius. The virial mass is the mass within the virial radius. We will discuss the gravitational collapse and virialization of dark-matter halos in detail in Chapter 17.3.1. One thing we will derive there is that after formation, the virialized part of dark-matter halos has a radial extent—the virial radius—within which the mean density is a constant factor \(\Delta_v\) times the Universe’s critical density (the density for which the Universe is spatially flat, \(\rho_\mathrm{crit} = 3 H^2/[8\pi G]\), with \(H\) the Hubble parameter; see Chapter 17.1.1). This factor is often set to 200, although technically this factor depends on a halo’s formation redshift and it is only \(\approx 200\) at \(z \gtrsim 2\) and \(\approx 100\) at the present.
To relate the virial radius to the \((\rho_0,a)\) parameterization used in Equation (2.49), we first compute the mean density within a given radius for the NFW profile using Equation (2.51) \begin{equation}\label{eq-meandens-nfw} \bar{\rho}(< r) = {3\,\rho_0\,a^3 \over r^3}\,\left(\,\ln[1+r/a]-\frac{r/a}{1+r/a}\right)\,. \end{equation} At the virial radius \(r_\mathrm{vir}\), this mean density should equal \(\Delta_v \,\rho_\mathrm{crit}\) \begin{equation}\label{eq-meandens-nfw-atrvir} {3\,\rho_0 \over c^3}\,\left(\,\ln[1+c]-\frac{c}{1+c}\right) = \Delta_v\,\rho_\mathrm{crit}\,, \end{equation} where \(c\) is the concentration, \(c = r_\mathrm{vir}/a\), the second parameter used in the virial-mass–concentration parameterization. The expression in parentheses is often abbreviated as \(f(c)\) \begin{equation} f(c) = \ln[1+c]-\frac{c}{1+c}\,. \end{equation} Because of the definition of the virial radius, we have the following relation between virial mass \(M_\mathrm{vir}\) and virial radius \begin{equation}\label{eq-nfw-mvir-rvir} M_\mathrm{vir} = {4\pi r_\mathrm{vir}^3 \over 3}\,\Delta_v\,\rho_\mathrm{crit}\,. \end{equation} For example, for the Milky Way with \(M_\mathrm{vir} \approx 10^{12}\,M_\odot\), we have \(r_\mathrm{vir} \approx 250\,\mathrm{kpc}\). Thus, given a pair \((M_\mathrm{vir},c)\) and a value of \(\Delta_v\), we can set the amplitude parameter \(\rho_0\) using \begin{equation}\label{eq-nfw-rhoo-fromconc} \rho_0 = {\Delta_v\,\rho_\mathrm{crit} \over 3}\,{c^3 \over f(c)}\,, \end{equation} and \(\rho_0\) therefore only depends on the concentration and is independent of the virial mass. Similarly, we obtain \(a = r_\mathrm{vir}/c\) or \begin{equation}\label{eq-nfw-a-frommvirc} a = {1 \over c}\,\left(3\,M_\mathrm{vir} \over 4\pi\,\Delta_v\,\rho_\mathrm{crit}\right)^{1/3}\,. \end{equation}
Conversely, to determine \(M_\mathrm{vir}\) and \(c\) for a given \(\rho_0\) and \(a\), we need to numerically solve Equation (2.55) for \(c\), compute \(r_\mathrm{vir} = a\,c\), and obtain \(M_\mathrm{vir}\) from Equation (2.57). In terms of \((M_\mathrm{vir},c,r_\mathrm{vir})\), the enclosed mass profile is given by \begin{equation}\label{eq-nfw-encmass-invirial} M(< r) = M_\mathrm{vir}\,{f(c r/r_\mathrm{vir}) \over f(c)}\,. \end{equation}
Another popular way of parameterizing NFW halos is by the maximum of the circular-velocity curve, \(v_\mathrm{max}\), and the radius \(r_\mathrm{max}\) at which this maximum occurs. This is an especially useful parameterization to use when analyzing \(N\)-body simulations of dark-matter halos, because these quantities are easy to measure from the simulation. To relate this to the other parameterizations, we first derive the rotation curve of an NFW profile by combining Equations (2.24) and (2.60) \begin{align} v^2_c(r) & = {G M_\mathrm{vir}\,\over a\,f(c)}\,{f(r/a) \over r/a}\,.\label{eq-nfw-vc2} \end{align} Only the second fraction depends on \(r\) and to find the radius of the maximum, we only need to consider this factor; we find the maximum by setting the derivative of this expression to zero, using \(x \equiv r/a\) as the variable \begin{equation} 0 = {\mathrm{d} \over \mathrm{d} x}\,\left( {f(x) \over x}\right) = -{\ln(1+x) \over x^2} + {2 x+1 \over x(1+x)^2}\,. \end{equation}
and we can find the exact location of \(x_\mathrm{max}\) by finding the root of this function:
[8]:
from scipy import optimize
print("""The NFW rotation curve peaks at r/a = """
f"""{optimize.brentq(lambda x: -numpy.log(1+x)/x+(2.*x+1)/(1.+x)**2.,0.1,10):.16f}""")
The NFW rotation curve peaks at r/a = 2.1625815870646115
Thus, we have that
\begin{equation} r_\mathrm{max} = 2.1625815870646115\times a = {2.1625815870646115 \over c}\,\left(3\,M_\mathrm{vir} \over 4\pi\,\Delta_v\,\rho_\mathrm{crit}\right)^{1/3}\,, \end{equation} Note that the relation between \(r_\mathrm{max}\) and \(a\) does not depend on the second parameter of the NFW profile (\(\rho_0\), \(M_\mathrm{vir}\), or \(v_\mathrm{max}\)). To get \(v_\mathrm{max}\), we then substitute \(r_\mathrm{max}\) into Equation (2.61) and we get \begin{equation} v^2_\mathrm{max} = 0.21621659550187314\times {G M_\mathrm{vir}\over r_\mathrm{vir}}\,{c\over f(c)} = 0.21621659550187314\times 4\pi G \rho_0 \, a^2\,, \end{equation} where \(0.21621659550187314= f(2.1625815870646115)/2.1625815870646115\). Because \(\sqrt{G M_\mathrm{vir}/ r_\mathrm{vir}}\) is the circular velocity \(v_{c,\mathrm{vir}}\) at the virial radius, we can also write this as \begin{equation} v_\mathrm{max} \approx 1.21 \times v_{c,\mathrm{vir}} \,{\sqrt{c/f(c)} \over 2.6}\,, \end{equation} where 2.6 is a typical value for \(\sqrt{c/f(c)}\) (which varies by only \(\approx\pm 10\,\%\) for relevant concentrations).
Large numerical simulations of the formation of dark-matter halos in Universes with similar cosmological parameters as our own have shown that the formed dark-matter halos populate only a narrow sequence in the space of \((M_\mathrm{vir},c)\) following a relation between \(M_\mathrm{vir}\) and \(c\) (Navarro et al. 1996, Navarro et al. 1997). Such relations are known as concentration–mass relations and they are updated every few years whenever a new set of best-fit cosmological parameters are determined, but a representative relation is the redshift zero relation from Dutton & Macciò (2014) \begin{equation}\label{eq-nfw-duttonmaccio} \log_{10} c = 0.905-0.101\,\log_{10}\left( {M_\mathrm{vir} \over 10^{12}\,h^{-1}\,M_\odot}\right)\,, \end{equation} where \(\Delta_v = 200\) and the authors use a Hubble constant of \(H_0 = 100\,h\,\mathrm{km\,s}^{-1}\,\mathrm{Mpc}^{-1} = 67.1\,\mathrm{km\,s}^{-1}\,\mathrm{Mpc}^{-1}\) following Planck Collaboration et al. (2014). Figure 2.7 shows the relation between the two parameters in all three parameterizations of the NFW profile for this concentration–mass relation from the scale of the smallest dwarf galaxies to that of large clusters of galaxies.
[18]:
from galpy.potential import NFWPotential
from matplotlib.ticker import FuncFormatter
from matplotlib.cm import autumn, winter
def dutton14_conc(m200):
return 10.**(0.905-0.101*numpy.log10(m200*0.671/10.**12/u.Msun))
m200s= numpy.geomspace(1e7,1e15,101)*u.Msun
cs= numpy.linspace(2.,28.,101)
const_cs= numpy.arange(5,30,5)
const_m200s= numpy.geomspace(1e8,1e14,7)
const_c_cmap= autumn
const_m200_cmap= winter
figure(figsize=(10,5))
subplot(1,3,1)
galpy_plot.plot(m200s.to_value(u.Msun),dutton14_conc(m200s),
semilogx=True,gcf=True,color='k',lw=3.,
xrange=[3e6,3e15],
yrange=[0.,30],
xlabel=r'$M_\mathrm{vir} / M_\odot$',
ylabel=r'$c$')
for const_c in const_cs:
plot(m200s.to_value(u.Msun),
const_c*numpy.ones_like(m200s.to_value(u.Msun)),
color=const_c_cmap((const_c-5.)/20.),lw=1.5,zorder=0)
for const_m200 in const_m200s:
plot(const_m200*numpy.ones_like(cs),
cs,
color=const_m200_cmap((numpy.log10(const_m200)-8.)/6),
lw=1.5,zorder=0)
subplot(1,3,2)
# Set up NFWPotentials for each (mvir,conc)
nps= [NFWPotential(mvir=m200.to_value(1e12*u.Msun),
conc=dutton14_conc(m200).to_value(u.dimensionless_unscaled),
H=67.1,overdens=200.,wrtcrit=True,
ro=8.,vo=220.)
for m200 in m200s]
# rho0 = 4 rho(a) from definition
galpy_plot.plot([4.*np.dens(np.a,0.).to_value(u.Msun/u.pc**3) for np in nps],
[np.a*np._ro for np in nps],
loglog=True,gcf=True,color='k',lw=3.,
xrange=[3e-4,9e-2],
yrange=[0.08,900],
xlabel=r'$\rho_0\,(M_\odot\,\mathrm{pc}^{-3})$',
ylabel=r'$a\,(\mathrm{kpc})$')
# Also for const_cs
all_const_cs_nps= []
for const_c in const_cs:
all_const_cs_nps.append([NFWPotential(mvir=m200.to_value(1e12*u.Msun),
conc=const_c,
H=67.1,overdens=200.,wrtcrit=True,
ro=8.,vo=220.)
for m200 in m200s])
plot([4.*np.dens(np.a,0.).to_value(u.Msun/u.pc**3) for np in all_const_cs_nps[-1]],
[np.a*np._ro for np in all_const_cs_nps[-1]],
color=const_c_cmap((const_c-5.)/20.),lw=1.5,zorder=0)
# Also for const_m200s
all_const_m200s_nps= []
for const_m200 in const_m200s:
all_const_m200s_nps.append([NFWPotential(mvir=const_m200/1e12,
conc=c,
H=67.1,overdens=200.,wrtcrit=True,
ro=8.,vo=220.)
for c in cs])
plot([4.*np.dens(np.a,0.).to_value(u.Msun/u.pc**3) for np in all_const_m200s_nps[-1]],
[np.a*np._ro for np in all_const_m200s_nps[-1]],
color=const_m200_cmap((numpy.log10(const_m200)-8.)/6),
lw=1.5,zorder=0)
gca().xaxis.set_major_formatter(FuncFormatter(
lambda y,pos: (r'${{:.{:1d}f}}$'.format(int(\
numpy.maximum(-numpy.log10(y),0)))).format(y)))
gca().yaxis.set_major_formatter(FuncFormatter(
lambda y,pos: (r'${{:.{:1d}f}}$'.format(int(\
numpy.maximum(-numpy.log10(y),0)))).format(y)))
subplot(1,3,3)
galpy_plot.plot([np.rmax().to_value(u.kpc) for np in nps],
[np.vmax().to_value(u.km/u.s) for np in nps],
loglog=True,gcf=True,color='k',lw=3.,
xrange=[0.15,3000.],
yrange=[0.5,3000.],
xlabel=r'$r_\mathrm{max}\,(\mathrm{kpc})$',
ylabel=r'$v_\mathrm{max}\,(\mathrm{km\,s}^{-1})$')
# Also for const_cs
for const_c,const_cs_nps in zip(const_cs,all_const_cs_nps):
plot([np.rmax().to_value(u.kpc) for np in const_cs_nps],
[np.vmax().to_value(u.km/u.s) for np in const_cs_nps],
color=const_c_cmap((const_c-5.)/20.),lw=1.5,zorder=0)
# Also for const_m200s
for const_m200,const_m200s_nps in zip(const_m200s,all_const_m200s_nps):
plot([np.rmax().to_value(u.kpc) for np in const_m200s_nps],
[np.vmax().to_value(u.km/u.s) for np in const_m200s_nps],
color=const_m200_cmap((numpy.log10(const_m200)-8.)/6),
lw=1.5,zorder=0)
gca().xaxis.set_major_formatter(FuncFormatter(
lambda y,pos: (r'${{:.{:1d}f}}$'.format(int(\
numpy.maximum(-numpy.log10(y),0)))).format(y)))
gca().yaxis.set_major_formatter(FuncFormatter(
lambda y,pos: (r'${{:.{:1d}f}}$'.format(int(\
numpy.maximum(-numpy.log10(y),0)))).format(y)))
tight_layout();
Figure 2.7: Relation between the parameters in different parameterizations of the NFW profile. The thick black curve is the cosmological concentration-mass relation from Dutton & Macciò (2014).
We see that lighter dark-matter halos are more concentrated than heavier halos, which results in an increase in the typical density \(\rho_0\) as the size \(a\) decreases. The behavior of the \(v_\mathrm{max}\) vs. \(r_\mathrm{max}\) relation mostly reflects the changing mass, which affects both \(v_\mathrm{max}\) and \(r_\mathrm{max}\) significantly. To allow the effects of changing concentration and virial-mass to be disentangled, Figure 2.7 also shows the location of a grid of halos with either concentration or virial mass held fixed (at \(c = 5, 10, \ldots\) or \(M_\mathrm{vir} = 10^8, 10^9, \ldots M_\odot\)). We see that, as expected from Equation (2.58), the typical density \(\rho_0\) only depends on concentration, while \(a\) depends on both virial mass and concentration. A wide range of concentrations only results in a small range of \(v_\mathrm{max}\) at a given \(r_\mathrm{max}\) compared to the range of \(v_\mathrm{max}\).
.
.
.