6.6. The Wolf mass estimator¶
Before discussing how we can distinguish between a cored and a cuspy inner dark-matter profile using Jeans analyses of the dSph data, let’s understand one key result mentioned above better: that the enclosed mass at the half-mass radius is robust against variations in the assumed anisotropy \(\beta\). It turns out that we can construct a very robust mass estimator for spherical systems using this result. We present here a slightly simplified derivation of the resulting mass estimator first derived by Wolf et al. (2010).
We start by re-writing Equation (6.18) in the form \(\Sigma(R)\,\sigma_\mathrm{los}^2(R) = \int_{R^2}^\infty\mathrm{d}r^2\,\frac{1}{\sqrt{r^2-R^2}}\,f(r^2)\) to make it amenable to inversion through an Abel transform. The result of a somewhat tedious calculation is \begin{align} \Sigma(R) \,\sigma_\mathrm{los}^2(R) & = \int_{R^2}^\infty\mathrm{d}r^2\,\frac{1}{\sqrt{r^2-R^2}}\, \left\{[\nu\,\sigma_r^2\,(1-\beta)](r^2) + \int_{r^2}^\infty\mathrm{d} r'^2\,\frac{\left[\beta\,\nu\,\sigma_r^2\right](r'^2)}{2r'^2}\right\}\,. \end{align} This can be inverted through an Abel transform to give (see Appendix B.4.1) \begin{align}\label{eq-wolf-nusig-abeltransformed} \nu(r)\,\sigma_r^2(r)\,[1-\beta(r)] + \int_{r}^\infty\,\mathrm{d}r'\,\frac{\beta(r')\,\nu(r')\,\sigma_r^2(r')}{r'} & = -\frac{1}{\pi}\,\int_{r^2}^\infty\mathrm{d}R\,\frac{1}{\sqrt{R^2-r^2}}\,\frac{\mathrm{d}[\Sigma(R)\,\sigma_\mathrm{los}^2(R)]}{\mathrm{d}R}\,. \end{align} The right-hand side of this equation depends solely on observable quantities, the two-dimensional projected stellar surface density \(\Sigma(R)\) and the line-of-sight velocity dispersion \(\sigma_\mathrm{los}\). The left-hand side of the equation therefore also needs to be independent of the underlying model. That is, the unknown functions need to combine to give a model-free value to match the right-hand side. The three-dimensional density \(\nu(r)\) is directly related to \(\Sigma(R)\) through Equation (6.17), and the model functions that can vary in Equation (6.20) are therefore only \(\sigma_r(r)\) and \(\beta(r)\). Assuming that an isotropic solution exists, we can equate the left-hand side for non-zero \(\beta\) with that for zero \(\beta\) \begin{equation}\label{eq-wolf-nusig-diffbeta} \nu(r)\,\sigma_r^2(r;\beta=0)= \nu(r)\,\sigma_r^2(r; \beta)\,[1-\beta(r)]+ \int_{r}^\infty\,\mathrm{d}r'\,\frac{\beta(r')\,\nu(r')\,\sigma_r^2(r'; \beta)}{r'} \,, \end{equation} where we have written \(\sigma_r(r)\) as \(\sigma_r(r;\beta=0)\) or \(\sigma_r(r;\beta)\) to indicate that the necessary radial dispersion to match the observable right-hand side of Equation (6.20) depends on \(\beta\). Then we take the derivative of Equation (6.21) with respect to \(\ln r\) and use the spherical Jeans Equation (5.48) to write this equation in terms of the mass derived for different \(\beta\) \begin{equation} M(r;\beta)-M(r;\beta=0) = -\frac{\beta(r)\,r\,\sigma_r^2(r)}{G}\,\left(\frac{\mathrm{d} \ln [\nu\,\sigma_r^2]}{\mathrm{d} \ln r}+ \frac{\mathrm{d} \ln \beta}{\mathrm{d} \ln r} + 3\right)\,. \end{equation} We can then solve for a radius \(r_\mathrm{eq}\) where the difference in the mass estimated using different \(\beta\) is zero. Assuming that \(\beta \neq 0\) this means that we need \begin{equation}\label{eq-wolf-solve-req} \frac{\mathrm{d} \ln [\nu\,\sigma_r^2]}{\mathrm{d} \ln r}\Bigg|_{r_\mathrm{eq}}+ \frac{\mathrm{d} \ln \beta}{\mathrm{d} \ln r}\Bigg|_{r_\mathrm{eq}} + 3 = 0\,. \end{equation}
For realistic systems, the stellar density varies by orders of magnitude, while changes in \(\beta\) are \(\mathcal{O}(1)\), and we can therefore set \(\mathrm{d} \ln \beta/\mathrm{d} \ln r \approx 0\). As we saw above for the dSphs, for many observed systems we have that \(\sigma_\mathrm{los}(R) \approx \mathrm{constant}\), which is difficult to achieve with a strongly varying \(\sigma_r(r)\), so we can also set \(\mathrm{d} \ln \sigma_r^2/\mathrm{d} \ln r \approx 0\) (see Wolf et al. 2010 for more detailed arguments for this assumption). Equation (6.23) then becomes \begin{equation}\label{eq-wolf-req-3} \frac{\mathrm{d} \ln \nu}{\mathrm{d} \ln r}\Bigg|_{r_\mathrm{eq}} = -3\,. \end{equation} This means that the radius where the enclosed mass does not depend strongly on \(\beta\) is that radius \(r_3\) where the logarithmic slope of the density profile is \(-3\).
To obtain the mass enclosed within \(r_3\), we go back to Equation (6.20), set \(\beta=0\) (because the mass does not depend on \(\beta\)) and take the derivative of both sides of the equation with respect to \(\ln r\) \begin{equation} \frac{\mathrm{d}\{\nu(r)\,\sigma_r^2(r)\}}{\mathrm{d} \ln r} = -\frac{\mathrm{d}}{\mathrm{d} \ln r} \left(\frac{1}{\pi}\,\int_{r^2}^\infty\mathrm{d}R\,\frac{1}{\sqrt{R^2-r^2}}\,\frac{\mathrm{d}[\Sigma(R)\,\sigma_\mathrm{los}^2(R)]}{\mathrm{d}R} \right)\,, \end{equation} where both sides are evaluated at \(r_\mathrm{eq}\). For the case where \(\sigma_\mathrm{los} \approx \mathrm{constant}\), we can bring \(\sigma_\mathrm{los}^2\) outside of the integral and use Equation (6.17) to replace the integral with \(\nu(r)\). Using the spherical Jeans Equation (5.48), this then becomes \begin{equation} M(< r_\mathrm{eq}) = -\frac{\sigma^2_\mathrm{los}\,r_{\mathrm{eq}}}{G}\,\frac{\mathrm{d} \ln \nu(r)}{\mathrm{d} \ln r}\Bigg|_{r_\mathrm{eq}} \,. \end{equation} Using Equation (6.24), we find \begin{equation} M(< r_3) = 3\,G^{-1}\,\sigma^2_\mathrm{los}\,r_3\,. \end{equation}
Because \(r_3\) is defined by a property of the deprojected stellar density profile, it is not a directly observable quantity. For most realistic stellar density profiles (e.g., the Plummer and King profiles, and common Sérsic profiles), it happens to be the case that \(r_3 \approx r_{1/2}\), with \(r_{1/2}\) the three-dimensional radius that contains half of the mass (the half-mass radius). For the same profiles, it is also the case that \(r_{1/2} \approx 4 R_e/3\), with \(R_e\) the two-dimensional radius that contains half of the mass. In terms of these we have the final Wolf mass estimator \begin{align}\label{eq-wolf-rhalf} M(< r_{1/2}) & = 3\,G^{-1}\,\sigma^2_\mathrm{los}\,r_{1/2}\\ & = 4\,G^{-1}\,\sigma^2_\mathrm{los}\,R_e \,.\label{eq-wolf-re} \end{align} When using this estimator, it is important to be clear on the meaning of the half-mass radius. This is the radius that contains half of the mass in the tracer population. If most of the mass is in a different component, e.g., dark matter, then the half-mass radius of all of the mass can be much larger than the half-mass radius of the tracer population.
To see in practice how the enclosed mass near the half-mass radius is insensitive to the assumed anisotropy, we can model the line-of-sight velocity dispersion for the Carina dSphs from Walker et al. (2009b) that was shown above. For simplicity, we take these data to be that \(\sigma_\mathrm{los} = 6\,\mathrm{km\,s}^{-1}\) at five radii from \(r=50\,\mathrm{pc}\) to \(950\,\mathrm{pc}\). We use the Plummer profile for the stars given by Walker et al. (2009b), which has \(M = 3\times2.4\times 10^5\,M_\odot\)—assuming a mass-to-light ratio of 3—and \(b = 137\,\mathrm{pc}\) (for a Plummer profile \(b = R_e\), the two-dimensional half-mass radius). We then apply the spherical Jeans equation for different values of \(\beta\), fitting the gravitational potential as a general two-power density model from Chapter 2.4.6 with four free parameters: the amplitude, scale radius, and inner- and outer slopes of \(\mathrm{d} \ln \rho / \mathrm{d} r\). This fit is done for four values of of the constant anisotropy in the following code:
[7]:
from galpy.potential import PlummerPotential, TwoPowerSphericalPotential
from galpy.df import jeans
from scipy import optimize
Re= 137.*u.pc
stellar_pot= PlummerPotential(amp=3.*2.4*1e5*u.Msun,b=Re)
def setup_mass_profile(params):
mass_pot= TwoPowerSphericalPotential(amp=10.**params[0]*u.Msun,
a=10**params[1]*u.pc,
alpha=params[2],
beta=params[3])
mass_pot.turn_physical_off()
return mass_pot
def predict_sigmalos(rs,mass_pot,beta):
return numpy.array([jeans.sigmalos(mass_pot,r,
dens=lambda r: stellar_pot.dens(r,0,
use_physical=False),
surfdens=lambda R: stellar_pot.surfdens(R,\
1000,use_physical=False),
beta=beta)*220.
for r in rs])
def fit_mass_profile(beta):
def sigmalos_chi2(params):
# Check that p[2] = inner slope of the mass density is >0 and <3
if params[2] < 0 or params[2] >= 3: return numpy.inf
if params[3] <= 2: return numpy.inf
if 10**params[1] < 50. or 10**params[1] > 3000.: return numpy.inf
# Initialize mass model
mass_pot= setup_mass_profile(params)
# Predict sigma_los(r) at 5 different radii from 100 to 1000 pc
sigmalos_pred= predict_sigmalos(numpy.linspace(50.,950.,5)*u.pc,
mass_pot,beta)
return numpy.sum((sigmalos_pred-6.)**2.)
# We only do a rough fit with 50 evaluations to get a relatively quick result
bestfit= optimize.minimize(sigmalos_chi2,[7.,2.3,1,3],method='powell',
options={'xtol':0.05,'ftol':0.1,
'maxiter':10,'maxfev':50}).x
return setup_mass_profile(bestfit)
betas= [0.5,0,-0.5,-3.0]
bf_mass_pots= []
for beta in betas:
bf_mass_pots.append(fit_mass_profile(beta))
We then show the resulting \(\sigma_\mathrm{los}(r)\) in Figure 6.4 to demonstrate that all four values of \(\beta\) are able to produce a flat \(\sigma_\mathrm{los}(r)\) over Carina’s inner kpc. We also plot the resulting best-fit enclosed-mass profile for the four fits.
[8]:
figure(figsize=(12,5))
subplot(1,3,1)
for ii,beta in enumerate(betas):
prs= numpy.linspace(50.,1000.,11)*u.pc
plot(prs.to_value(u.pc),
predict_sigmalos(prs,bf_mass_pots[ii],beta),
label=rf'$\beta = {beta:+.1f}$')
xlabel(r'$r\,(\mathrm{pc})$')
ylabel(r'$\sigma_\mathrm{los}\,(\mathrm{km\,s}^{-1})$')
ylim(0.,15.)
legend(frameon=False,fontsize=18.)
subplot(1,3,2)
for ii,beta in enumerate(betas):
prs= numpy.linspace(0.01,1.,61)*u.kpc
semilogy(prs.to_value(u.pc),
[bf_mass_pots[ii].mass(r,ro=8.,vo=220.).to_value(u.Msun)
for r in prs],
label=rf'$\beta = {beta:+.1f}$')
semilogy(prs.to_value(u.pc),
[stellar_pot.mass(r,ro=8.,vo=220.).to_value(u.Msun) for r in prs],
ls=':',color='k')
axvline(4.*Re.to_value(u.pc)/3.,ls='--',color='0.7',lw=2.,zorder=0)
import astropy.constants as apyconst
axhline((4.*(6.*u.km/u.s)**2*Re/apyconst.G).to_value(u.Msun),
ls='--',color='0.7',lw=2.,zorder=0)
text(225.,5e7,r'$r_{1/2}$',fontsize=18.,horizontalalignment='left')
text(600.,4e5,r'$\mathrm{stars}$',fontsize=18.,horizontalalignment='left')
xlabel(r'$r\,(\mathrm{pc})$')
ylabel(r'$M(<r)/M_\odot$')
ylim(1e5,1e8)
legend(frameon=True,fontsize=18.,facecolor='w',edgecolor='None')
subplot(1,3,3)
for ii,beta in enumerate(betas):
prs= numpy.linspace(0.01,1.,61)*u.kpc
plot(prs.to_value(u.pc),
[bf_mass_pots[ii].mass(r,ro=8.,vo=220.).to_value(u.Msun)
/bf_mass_pots[1].mass(r,ro=8.,vo=220.).to_value(u.Msun)
for r in prs],
label=rf'$\beta = {beta:+.1f}$')
axvline(4.*Re.to_value(u.pc)/3.,ls='--',color='0.7',lw=2.,zorder=0)
xlabel(r'$r\,(\mathrm{pc})$')
ylabel(r'$M(<r;\beta)/M(<r;\beta=0)$')
tight_layout();
Figure 6.4: The Wolf mass estimator.
Because the enclosed mass spans a few orders of magnitude over the inner kpc, we also display the enclosed-mass profile relative to that obtained for \(\beta = 0\) in the right-most panel. From that panel, it becomes clear that all four values of \(\beta\) do indeed produce a similar value for the enclosed mass near the half-mass radius (although in this case, the value where all anisotropies agree is a little higher still). We also show the contribution to the enclosed-mass profile from the stars, which demonstrates that Carina, like all dSphs, is indeed a very dark-matter dominated system, with \(M/L \approx 150\) within 1 kpc.
.
.
.