← Previous section Next section →

5.2. The virial theorem

\label{sec-spherequil-virial}

The virial theorem relates the kinetic and potential energy of a system of bodies. We can derive different versions of the virial theorem, depending on whether we are dealing with a self-gravitating system of bodies (bodies moving under the influence of their mutual gravity) or a system of bodies orbiting within the potential due to an external mass distribution (for example, the system of globular clusters in the Milky Way, which orbit in the total Milky Way mass distribution and contribute almost no mass to this themselves). We can derive all different versions by starting from the virial quantity \(G\) \begin{equation}\label{eq-virial-def} G = \sum_{i=1}^N w_i\,\vec{x}_i\cdot\vec{v}_i\,, \end{equation} where \(w_i\) are a set of weights that is left arbitrary for now. Assuming that the system of \(N\) bodies is in equilibrium, we have that \(\mathrm{d} G / \mathrm{d} t = 0\) or \begin{equation} \sum_{i=1}^N w_i \dot{\vec{x}}_i\cdot\vec{v}_i = -\sum_{i=1}^N w_i\,\vec{x}_i\cdot\dot{\vec{v}}_i\,, \end{equation} which we can write in terms of the gravitational field as \begin{equation}\label{eq-virial-verygeneral} \sum_{i=1}^N w_i \vec{v}_i\cdot\vec{v}_i = -\sum_{i=1}^N w_i\,\vec{x}_i\cdot\vec{g}(\vec{x}_i)\,. \end{equation} Setting the weights equal to the masses of the bodies, \(w_i = m_i\), the left-hand side of this equation equals twice the kinetic energy of the system.

If we now consider that the bodies are orbiting within a central, external point-mass potential, \(\Phi(r) = -G\,M/r\), then we get the following virial theorem \begin{align} \sum_{i=1}^N w_i |\vec{v}_i|^2 = \sum_{i=1}^N w_i\,\frac{G\,M}{r_i}\,. \end{align} where \(r_i = |x_i|\) and the right-hand side is now minus the potential energy. The mass of the system is therefore \begin{align}\label{eq-virial-mass-kepler} M = \frac{1}{G}\,\frac{\sum_{i=1}^N w_i |\vec{v}_i|^2}{\sum_{i=1}^N w_i\,\frac{1}{r_i}}\,. \end{align} This result holds for any weighting \(w_i\), although how good of an estimator the right-hand side is for the mass will strongly depend on how the bodies are weighted (for example, setting all \(w_i\) but one equal to zero will make this a very noisy estimator! Of course, then the assumption that \(\mathrm{d} G / \mathrm{d} t = 0\) also breaks down). The simplest choice for the weights is \(w_i=1\).

Alternatively, we can work out Equation (5.13) for a set of self-gravitating bodies. In that case, \(\vec{g}(\vec{x}_i) = \sum_{j \neq i} G\,m_j\,(\vec{x}_i-\vec{x}_j)/|\vec{x}_i-\vec{x}_j|^3\). The right-hand side of Equation (5.13) can then be written as \begin{align} -\sum_{i=1}^N w_i\,\vec{x}_i\cdot\vec{g}(\vec{x}_i) & = \sum_{i=1}^N \sum_{j \neq i} w_i\,\vec{x}_i\cdot\,\frac{G\,m_j\,(\vec{x}_i-\vec{x}_j)}{|\vec{x}_i-\vec{x}_j|^3}\nonumber\\ & = \sum_{i=1}^N \sum_{j < i} w_i\,\vec{x}_i\cdot\,\frac{G\,m_j\,(\vec{x}_i-\vec{x}_j)}{|\vec{x}_i-\vec{x}_j|^3}+\sum_{i=1}^N \sum_{j > i} w_i\,\vec{x}_i\cdot\,\frac{G\,m_j\,(\vec{x}_i-\vec{x}_j)}{|\vec{x}_i-\vec{x}_j|^3}\label{eq-virial-selfgrav-simplify}\\ & = \sum_{i=1}^N \sum_{j < i}\left[(w_i\,m_j\,\vec{x}_i-w_j\,m_i\,\vec{x}_j)\cdot\,\frac{G\,(\vec{x}_i-\vec{x}_j)}{|\vec{x}_i-\vec{x}_j|^3}\right]\nonumber\,. \end{align} If we then set \(w_i = m_i\), this can be simplified to \begin{equation} -\sum_{i=1}^N w_i\,\vec{x}_i\cdot\vec{g}(\vec{x}_i) = \sum_{i=1}^N \sum_{j < i}\,\frac{G\,m_i\,m_j}{|\vec{x}_i-\vec{x}_j|}\,, \end{equation} which is again minus the potential energy of the system. The virial theorem for self-gravitating systems is therefore \begin{equation}\label{eq-spherequil-virialtheorem-selfgravity} \sum_{i=1}^N m_i \vec{v}_i\cdot\vec{v}_i = \sum_{i=1}^N \sum_{j < i}\,\frac{G\,m_i\,m_j}{|\vec{x}_i-\vec{x}_j|}\,, \end{equation} If all bodies have the same mass \(m\) and the total mass is \(M = Nm\), then we can estimate the mass as \begin{equation}\label{eq-virial-mass-selfgrav} M = \frac{1}{G}\,\frac{\sum_i |\vec{v}_i^2|}{\frac{1}{N}\,\sum_{i=1}^N\sum_{j < i}\,\frac{1}{|\vec{x}_i-\vec{x}_j|}}\,. \end{equation}

Equations (5.15) and (5.19) both make it clear that under the assumption of equilibrium, the mass of a stellar system can be obtained by balancing the spread in velocities (the kinetic energy) with an appropriate weighting of the spread in positions (the potential energy). This makes sense: if we guess a mass \(M\) that is too low, then the observed velocities would be too high and the system would in the future quickly spread out and thus not be in equilibrium; if on the other hand, we guess a mass that is too high, the velocities would appear low and the system will quickly collapse to occupy a smaller volume. For the correct mass \(M\), both the spread in velocities and in positions can be maintained over time.

The virial theorem in the form above is also known as the scalar virial theorem, to distinguish it from more detailed forms of the virial theorem like the tensor virial theorem that we will discuss in Chapter 14.1 in the context of elliptical galaxies. There, we also discuss how to generalize the virial theorem to continuous, collisionless systems. When we set \(w_i = m_i\), both the external-point-mass-potential and self-gravitating cases of the virial theorem can be written as \begin{equation}\label{eq-virial-scalar} 2K+W = 0\,, \end{equation} where \(K\) is the kinetic energy and \(W\) is the potential energy. Using that the total energy \(E = K + W\), we can also write this as \begin{equation}\label{eq-virial-scalar-energy} 2E-W = 0\,. \end{equation}

In the self-gravitating case, Equation (5.18) becomes \begin{equation}\label{eq-spherequil-virialtheorem-selfgravity-continuous} \int \mathrm{d}\vec{x}\,\rho(\vec{x})\,\overline{v^2} = -\frac{1}{2} \int\mathrm{d}\vec{x}\,\rho(\vec{x},t)\,\Phi(\vec{x},t)\,, \end{equation} in the continuous limit. We derive this form in more detail in Chapter 14.1, but it directly follows from Equation (5.18) when replacing sums with integrals, masses by infinitesimal mass elements, and remembering that \(|\vec{x}_i-\vec{x}_j|\) is symmetric in \((i,j)\).

Let’s use these virial estimators to estimate the mass of the Milky Way using globular clusters from the Harris (1996) catalog available here. First we write a function to download and read the catalog, which we save in a cache directory $HOME/.galaxiesbook/cache/harris (note that the catalog is somewhat difficult to read, which is why the harris_read function below is quite complicated…):

# harris.py: download, parse, and read the Harris
# Milky-Way globular cluster catalog
# Include all necessary imports to make this cell self-contained
import sys
import os
from pathlib import Path
import shutil
import subprocess
import tempfile
import numpy
import pandas
_MAX_DOWNLOAD_NTRIES= 2
_ERASESTR= "                                                                                "
_CACHE_BASEDIR= Path(os.getenv('HOME')) / '.galaxiesbook' / 'cache'
_CACHE_HARRIS_DIR= _CACHE_BASEDIR / 'harris'
_HARRIS_URL= 'http://physwww.mcmaster.ca/~harris/mwgc.dat'
def harris_read():
    """
    NAME:
       read
    PURPOSE:
       read and parse the Harris Milky-Way globular cluster catalog
    INPUT:
       (none)
    OUTPUT:
       pandas DataFrame with the catalog
    HISTORY:
       2017-07-26 - Started - Bovy (UofT)
    """
    filePath= _CACHE_HARRIS_DIR / 'mwgc.dat'
    _download_harris(filePath)
    with open(filePath,'r') as harrisfile:
        lines= harrisfile.readlines()
        # 1st part: 72-229, 2nd part: 252-409, 3rd part: 433-590
        ids= lines[72:229]
        metals= lines[252:409]
        vels= lines[433:590]
        # Parse into dict
        ids_dict= {'ID': [],
                   'Name': [],
                   'RA': numpy.empty(len(ids)),
                   'DEC': numpy.empty(len(ids)),
                   'L': numpy.empty(len(ids)),
                   'B': numpy.empty(len(ids)),
                   'R_Sun': numpy.empty(len(ids)),
                   'R_gc': numpy.empty(len(ids)),
                   'X': numpy.empty(len(ids)),
                   'Y': numpy.empty(len(ids)),
                   'Z': numpy.empty(len(ids)),
                   '[Fe/H]': numpy.empty(len(ids)),
                   'wt': numpy.empty(len(ids)),
                   'E(B-V)': numpy.empty(len(ids)),
                   'V_HB': numpy.empty(len(ids)),
                   '(m-M)V': numpy.empty(len(ids)),
                   'V_t': numpy.empty(len(ids)),
                   'M_V,t': numpy.empty(len(ids)),
                   'U-B': numpy.empty(len(ids)),
                   'B-V': numpy.empty(len(ids)),
                   'V-R': numpy.empty(len(ids)),
                   'V-I': numpy.empty(len(ids)),
                   'spt': ['---' for ii in range(len(ids))],
                   'ellip': numpy.empty(len(ids)),
                   'v_r': numpy.empty(len(ids)),
                   'v_r_unc': numpy.empty(len(ids)),
                   'v_LSR': numpy.empty(len(ids)),
                   'sig_v': numpy.empty(len(ids)),
                   'sig_v_unc': numpy.empty(len(ids)),
                   'c': numpy.empty(len(ids)),
                   'r_c': numpy.empty(len(ids)),
                   'r_h': numpy.empty(len(ids)),
                   'mu_V': numpy.empty(len(ids)),
                   'rho_0': numpy.empty(len(ids)),
                   'lg(tc)': numpy.empty(len(ids)),
                   'lg(th)': numpy.empty(len(ids)),
                   'core_collapsed': numpy.zeros(len(ids),dtype='bool')}
        # Parse first part
        rest_keys= ['L','B','R_Sun','R_gc','X','Y','Z']
        for ii in range(len(ids)):
            # Parse ID and Name
            ids_dict['ID'].append(ids[ii][:10].strip())
            name= ids[ii][10:25].strip()
            if name == '': ids_dict['Name'].append('---')
            else: ids_dict['Name'].append(name)
            # Parse RA and Dec, Harris (l,b) don't seem exactly right...
            hms= ids[ii][25:36].split() # hour, min, sec
            ids_dict['RA'][ii]=\
                int(hms[0])*15.+int(hms[1])/4.+float(hms[2])/240.
            hms= ids[ii][38:49].split() # deg, arcmin, arcsec
            sgn= -1 if ids[ii][38] == '-' else 1
            ids_dict['DEC'][ii]=\
                int(hms[0])+sgn*(int(hms[1])/60.+float(hms[2])/3600.)
            # Split the rest
            rest= ids[ii][49:].split()
            for r,k in zip(rest,rest_keys): ids_dict[k][ii]= float(r.strip())
        # Parse second part
        idx_arr= numpy.arange(len(ids))
        for ii in range(len(metals)):
            # Not sure whether I can assume that they are in order
            idname= metals[ii][:10].strip()
            idx= idx_arr[\
                numpy.strings.equal(ids_dict['ID'],idname)][0]
            ids_dict= _parse_float_entry(ids_dict,'[Fe/H]',idx,
                                         metals[ii][13:18].strip())
            ids_dict= _parse_float_entry(ids_dict,'wt',idx,
                                         metals[ii][19:21].strip())
            ids_dict= _parse_float_entry(ids_dict,'E(B-V)',idx,
                                         metals[ii][24:28].strip())
            ids_dict= _parse_float_entry(ids_dict,'V_HB',idx,
                                         metals[ii][29:34].strip())
            ids_dict= _parse_float_entry(ids_dict,'(m-M)V',idx,
                                         metals[ii][35:40].strip())
            ids_dict= _parse_float_entry(ids_dict,'V_t',idx,
                                         metals[ii][42:46].strip())
            ids_dict= _parse_float_entry(ids_dict,'M_V,t',idx,
                                         metals[ii][48:53].strip())
            for k,si,ei in zip(['U-B','B-V','V-R','V-I'],
                               [56,62,68,74],[60,66,72,78]):
                ids_dict= _parse_float_entry(ids_dict,k,idx,
                                             metals[ii][si:ei].strip())
            ids_dict['spt'][idx]= metals[ii][80:82].strip()
            ids_dict= _parse_float_entry(ids_dict,'ellip',idx,
                                         metals[ii][86:].strip())
        # Parse third part
        idx_arr= numpy.arange(len(ids))
        for ii in range(len(vels)):
            # Not sure whether I can assume that they are in order
            idname= vels[ii][:10].strip()
            idx= idx_arr[\
                numpy.strings.equal(ids_dict['ID'],idname)][0]
            for k,si,ei in zip(['v_r','v_r_unc','v_LSR','sig_v','sig_v_unc',
                                'c','r_c','r_h','mu_V','rho_0',
                                'lg(tc)','lg(th)'],
                               [13,21,27,37,43,50,59,66,73,80,88,93],
                               [19,25,33,41,47,54,64,70,78,85,91,98]):
                ids_dict= _parse_float_entry(ids_dict,k,idx,
                                             vels[ii][si:ei].strip())
            if vels[ii][56] == 'c': ids_dict['core_collapsed'][idx]= True
        return pandas.DataFrame(ids_dict)

def _parse_float_entry(ids_dict,key,idx,val):
    try:
        ids_dict[key][idx]= float(val)
    except ValueError:
        ids_dict[key][idx]= numpy.nan
    return ids_dict

def _download_harris(filePath):
    if not filePath.exists():
        _download_file(_HARRIS_URL,filePath)
    return None

def _download_file(downloadPath,filePath):
    sys.stdout.write('\r'+"Downloading file %s ...\r" \
                         % (filePath.name))
    sys.stdout.flush()
    try:
        # make all intermediate directories
        os.makedirs(filePath.parent)
    except OSError: pass
    # Safe way of downloading
    downloading= True
    interrupted= False
    file, tmp_savefilename= tempfile.mkstemp()
    os.close(file) #Easier this way
    ntries= 1
    while downloading:
        try:
            cmd= ['curl','-L',f'{downloadPath}',
                  '-o',f'{tmp_savefilename}',
                  '--fail','--silent','--show-error',
                  '--connect-timeout','10',
                  '--retry','3']
            subprocess.check_call(cmd)
            shutil.move(tmp_savefilename,filePath)
            downloading= False
            if interrupted:
                raise KeyboardInterrupt
        except subprocess.CalledProcessError as e:
            if not downloading: #Assume KeyboardInterrupt
                raise
            elif ntries > _MAX_NTRIES:
                raise IOError('File %s does not appear to exist on the server ...' % (filePath.name))
            elif not 'exit status 4' in str(e):
                interrupted= True
            os.remove(tmp_savefilename)
        finally:
            if Path(tmp_savefilename).exists():
                os.remove(tmp_savefilename)
        ntries+= 1
    sys.stdout.write('\r'+_ERASESTR+'\r')
    sys.stdout.flush()
    return None

We now load the catalog and convert the positions to the Galactocentric frame:

import astropy.units as u
import astropy.constants as const
import astropy.coordinates as apycoords
gcdata= harris_read()
c= apycoords.SkyCoord(ra=gcdata['RA'],
                      dec=gcdata['DEC'],
                      distance=gcdata['R_Sun'],
                      unit=(u.deg,u.deg,u.kpc),
                      frame='icrs')
gc_frame= apycoords.Galactocentric(galcen_distance=8.1*u.kpc,
                                   z_sun=25.*u.pc)
gc = c.transform_to(gc_frame)
gc.representation_type = 'spherical'

To be in a regime where the assumption of spherical symmetry is closer to being correct, we only select clusters at distances well beyond the disk, \(r > 20\,\mathrm{kpc}\). We can then estimate the mass of the Milky Way using the estimator from Equation (5.15), assuming that the globular clusters are test particles in a point-mass potential. Observations of the proper motions of globular clusters are difficult, so we will estimate \(\sum_i |\vec{v}_i|^2\) from just the line-of-sight velocities under the assumption that (a) the velocity distribution is Gaussian and isotropic and (b) the line-of-sight velocity (approximately corrected for the Sun’s motion with respect to the solar neighborhood to \(v_\mathrm{LSR}\)) is approximately equal to the Galactocentric radial velocity; then we can estimate \(\sum_i |\vec{v}_i|^2\) as \(3\sum_i v_\mathrm{LSR}^2\). Putting this together with the Galactocentric distances, we find a mass of \(\approx 5\times 10^{11}\,M_\odot\):

indx= (gc.distance > 20.*u.kpc)*(True^numpy.isnan(gcdata['v_LSR']))
M_est= (3.*numpy.sum(gcdata['v_LSR'][indx]**2.)\
        *(u.km/u.s)**2\
       /numpy.sum(1./gc.distance[indx])/const.G)
print("""For a point-mass central potential,
    with %i GCs with median distance %.0f kpc
    the mass is %.2f 10^12 Msun""" % \
    (numpy.sum(indx),
     numpy.median(gc.distance[indx]).to(u.kpc).value,
     M_est.to(10**12*u.Msun).value))

For a point-mass central potential,
    with 19 GCs with median distance 35 kpc
    the mass is 0.45 10^12 Msun

If instead we assume that the globular clusters form a self-gravitating system and similarly apply Equation (5.19), we get a mass of \(\approx 2\times 10^{11}\,M_\odot\):

mutual_dist= \
 (numpy.atleast_2d(gc.distance[indx]).T-gc.distance[indx])
mutual_dist[numpy.diag_indices(numpy.sum(indx))]= 0.
mutual_dist= numpy.triu(mutual_dist)
M_est= (3.*numpy.sum(gcdata['v_LSR'][indx]**2.)*(u.km/u.s)**2\
        /numpy.sum(1./numpy.fabs(mutual_dist)[mutual_dist > 0])\
        *numpy.sum(indx)/const.G)
print("""For a self-gravitating system,
    with %i GCs with median distance %.0f kpc
    the mass is %.2f 10^12 Msun""" % \
    (numpy.sum(indx),
     numpy.median(gc.distance[indx]).to(u.kpc).value,
     M_est.to(10**12*u.Msun).value))

For a self-gravitating system,
    with 19 GCs with median distance 35 kpc
    the mass is 0.17 10^12 Msun

We see that there is only a factor of two difference between the two estimators, so no matter what estimator we use, we always find a large mass (larger than the total mass of stars and gas by at least a factor of two). It is also clear from this that the globular clusters cannot form a self-gravitating system: for the self-gravitating estimator to make sense, each globular cluster would have to have a mass of about \(10^{10}\,M_\odot\), so their mass-to-light ratio would have to be \(M/L \gtrsim 10^4\)!

Comparing the estimate assuming a point-mass potential with the enclosed-mass profile for the Milky Way that we discussed in Chapter 1.2.1, we see that the virial estimator is not that far off from our current best knowledge (it is about 30% higher than the best estimate), even though the virial estimator is highly simplified. This demonstrates the power and simplicity of equilibrium modeling.

← Previous section Next section →