A. Coordinate systems

A.1. Generalities

\label{sec-coordsapp-general}

We use a variety of coordinate systems in these notes, which we briefly introduce here. Because most stellar systems are either close-to-spherical or have a disk-like geometry, the two main coordinate systems that we use are spherical coordinates and cylindrical coordinates. You should be familiar with spherical coordinates. Our convention will be that the spherical coordinates are \((r,\phi,\theta)\) for the (radial, azimuthal, and polar) coordinate, with \(\theta\) measured from the pole. These coordinates are given in terms of the cartesian coordinates \((x,y,z)\) as

\begin{align}\label{eq-sphere-coords-1} x &= r\,\sin \theta\, \cos \phi\,,\\ y &= r\,\sin \theta\, \sin \phi\,,\\ z &= r\,\cos \theta\,,\label{eq-sphere-coords-2} \end{align}

and

\begin{align} r &= \sqrt{x^2 + y^2 + z^2}\,,\\ \phi &= \mathrm{atan2}(y,x)\,,\\ \theta &= \pi/2 - \mathrm{atan}(z/\sqrt{x^2+y^2})\,, \end{align}

where atan2 is the arctangent function with two arguments. We will typically define velocities in the spherical coordinate system to be \((v_r,v_\phi,v_\theta) = (\dot{r},r\,\sin \theta\,\dot{\phi},r\,\dot{\theta})\), where the dot denotes the time derivative. The relation between the velocity \((v_x,v_y,v_z)\) in cartesian coordinates and that in spherical coordinates is then

\begin{align}\label{eq-vel-spher} v_x & = v_r\,\sin\theta\,\cos\phi-v_\phi\,\sin\phi+v_\theta\,\cos\theta\,\cos\phi\\ v_y & = v_r\,\sin\theta\,\sin\phi+v_\phi\,\cos\phi+v_\theta\,\cos\theta\,\sin\phi\\ v_z & = v_r\,\cos\theta\phantom{\,\sin\phi+v_\phi\,\cos\phi}-v_\theta\,\sin\theta\,. \end{align}

The volume element is given by

\begin{equation} \mathrm{d} \vec{x} = r^2\,\sin\theta\,\mathrm{d}r\,\mathrm{d}\phi\,\mathrm{d}\theta\,. \end{equation}

The gradient \(\nabla\) in spherical coordinates is given by

\begin{equation}\label{eq-gradient-spherical} \nabla = \hat{\mathbf{e}}_r\,\frac{\partial}{\partial r} + \frac{\hat{\mathbf{e}}_\theta}{r}\,\frac{\partial}{\partial \theta} + \frac{\hat{\mathbf{e}}_\phi}{r\,\sin\theta}\,\frac{\partial}{\partial \phi}\,, \end{equation}

where \((\hat{\mathbf{e}}_r,\hat{\mathbf{e}}_\phi,\hat{\mathbf{e}}_\theta)\) are the unit vectors for spherical coordinates. The Laplacian \(\nabla^2\) is given by

\begin{equation}\label{eq-laplacian-spherical} \nabla^2 = \frac{1}{r^2}\,\frac{\partial}{\partial r}\left(r^2\,\frac{\partial}{\partial r}\right) +\frac{1}{r^2\sin\theta}\,\frac{\partial}{\partial \theta}\left(\sin\theta\,\frac{\partial}{\partial \theta}\right) +\frac{1}{r^2\sin^2\theta}\,\frac{\partial^2}{\partial \phi^2}\,, \end{equation}

To describe the dynamics of disk galaxies, we will use cylindrical coordinates. Our convention is that cylindrical coordinates are \((R,\phi,z)\) for the (radial, azimuthal, and vertical) coordinate. These coordinates are given in terms of the cartesian coordinates as

\begin{align} x & = R\,\cos \phi\,,\\ y & = R\,\sin \phi\,,\\ z & = z\,, \end{align}

and

\begin{align} R &= \sqrt{x^2 + y^2}\,,\\ \phi &= \mathrm{atan2}(y,x)\,,\\ z &= z\,. \end{align}

We will typically define velocities in the cylindrical coordinate system to be \((v_R,v_\phi,v_z) = (\dot{R},R\,\dot{\phi},\dot{z})\). The relation between the velocity \((v_x,v_y,v_z)\) in cartesian coordinates and that in spherical coordinates is then

\begin{align}\label{eq-vel-cyl} v_x & = v_R\,\cos\phi-v_\phi\,\sin \phi\\ v_y & = v_R\,\sin\phi+v_\phi\,\cos\phi\\ v_z & = v_z\,. \end{align}

The volume element is given by

\begin{equation} \mathrm{d} \vec{x} = R\,\mathrm{d}R\,\mathrm{d}\phi\,\mathrm{d}z\,. \end{equation}

The gradient \(\nabla\) in cylindrical coordinates is given by

\begin{equation} \nabla = \hat{\mathbf{e}}_R\,\frac{\partial}{\partial R} + \frac{\hat{\mathbf{e}}_\phi}{R}\,\frac{\partial}{\partial \phi}+\hat{\mathbf{e}}_z\,\frac{\partial}{\partial z}\,, \end{equation}

where \((\hat{\mathbf{e}}_R,\hat{\mathbf{e}}_\phi,\hat{\mathbf{e}}_z)\) are the unit vectors for cylindrical coordinates. The Laplacian \(\nabla^2\) for cylindrical coordinates given by

\begin{equation}\label{eq-laplacian-cylindrical} \nabla^2 = \frac{1}{R}\,\frac{\partial}{\partial R}\left(R\,\frac{\partial}{\partial R}\right) + \frac{1}{R^2}\,\frac{\partial^2}{\partial \phi^2}+\frac{\partial^2}{\partial z^2}\,. \end{equation}

Don’t be too harsh on the author of these notes if he occasionally writes ‘Z’ instead of ‘z’…

Polar coordinates are the two-dimensional coordinate system corresponding to the \((R,\phi)\) part of cylindrical coordinates. In particular, the Laplacian in polar coordinates is simply

\begin{equation}\label{eq-laplacian-polar} \nabla^2 = \frac{1}{R}\,\frac{\partial}{\partial R}\left(R\,\frac{\partial}{\partial R}\right) + \frac{1}{R^2}\,\frac{\partial^2}{\partial \phi^2}\,. \end{equation}

A.2. Positions in the Milky Way

\label{sec-mwpos}

To describe the dynamics of the Milky Way, we require various coordinate frames. The basic coordinate frame that astronomical measurements are reported in is the equatorial system, which is a spherical coordinate system centered in the Earth with a longitudinal angle called right ascension (RA) and a latitudinal angle called declination (Dec). The radial direction is the distance from Earth, but note that for the purpose of Galactic astronomy the difference between distances measured from the Sun and the Earth is negligible.

A more convenient coordinate system for studying the Milky Way is that of Galactic coordinates. This is simply a rotation of the equatorial system, such that the resulting system is aligned with the Milky Way. Specifically, Galactic coordinates are still centered on the Earth and the radial direction remains distance. The latitudinal angle (\(b\)) is such that an angle of zero is the (approximate position of the) Galactic midplane as seen from Earth and the longitudinal angle (\(l\)) is such that an angle of zero at latitude zero points towards the (approximate position of the) Galactic center. Note that we do not follow the convention for spherical coordinates agreed to above here, because \(b=0\) at the midplane rather than the pole (welcome to astronomy!). We won’t go into detail of how to transform from equatorial coordinates to Galactic coordinates; take a look at these notes for a detailed discussion of the definition of the equatorial and Galactic coordinate systems and how to transform coordinates between them. In practice, if you are given coordinates in the equatorial coordinate system, you can transform them to Galactic coordinates using astropy. For example, to transform (RA,Dec) = \((120^\circ,30^\circ)\) to Galactic coordinates \((l,b)\) we setup a SkyCoord object for the given (RA,Dec) and ask for this position in Galactic coordinates (the SkyCoord class does the transformation internally):

import astropy.units as u
import astropy.coordinates as apycoords
ra= 120.*u.deg
dec= 30.*u.deg
c= apycoords.SkyCoord(ra=ra,dec=dec,frame='icrs')
print(c.galactic)

<SkyCoord (Galactic): (l, b) in deg
    ( 191.27460843,  27.07427141)>

There have been various definitions of equatorial coordinates over the last few decades, but unless you are told otherwise, coordinates that you are handed are probably in the ICRS frame. RA and Dec can be specified in a large variety of ways (no, not just in degrees as would be handy, but in things like “hours, minutes, seconds” and such); astropy can deal with most of these. SkyCoord also works for array, so this works for example

ra= 120.*u.deg*numpy.ones(3)
dec= 30.*u.deg*numpy.ones(3)
c= apycoords.SkyCoord(ra=ra,dec=dec,frame='icrs')
print(c.galactic)

<SkyCoord (Galactic): (l, b) in deg
    [( 191.27460843,  27.07427141), ( 191.27460843,  27.07427141),
     ( 191.27460843,  27.07427141)]>

The cartesian frame associated with Galactic coordinates has coordinates \((X,Y,Z)\) and velocities traditionally given as \((U,V,W)\). The coordinates are simply computed as

\begin{align} X & = D\,\cos l\,\cos b\,,\\ Y & = D\,\sin l\,\cos b\,,\\ Z & = D\,\sin b\,, \end{align}

where \(D\) is the distance. Because \(l = 0\) corresponds to the Galactic center, \(X\) is positive going towards the Galactic center and negative going away from it. Confusingly, \((X,Y,Z)\) is a left-handed coordinate system; \(Y\) is positive in the direction of Galactic rotation (which is perpendicular to the Sun–Galactic center line). \(Z\) is positive towards the North Galactic Pole (NGP). We discuss velocities further below. astropy can also compute these cartesian coordinates. To do this, you also need to specify the distance in the SkyCoord instantiation and then ask for the output in cartesian Galactic coordinates:

ra= 120.*u.deg
dec= 30.*u.deg
distance= 1.2*u.kpc
c= apycoords.SkyCoord(ra=ra,dec=dec,distance=distance,frame='icrs')
print("(X,Y,Z) in (kpc,kpc,kpc)")
print("\t",c.galactic.cartesian)

(X,Y,Z) in (kpc,kpc,kpc)
         (-1.04788016, -0.20890423,  0.54617413) kpc

Galactic coordinates are useful only when we study the structure (in position, velocity, or chemistry) of the solar neighborhood. When we investigate the Milky Way on a global scale, we use Galactocentric coordinates that have the Galactic center as their zero point. Specifically, we will define Galactocentric coordinates as the left-handed coordinate system that has the Galactic center at its center, with \(x\) increasing in the disk midplane towards the location of the Sun, \(y\) increasing in the disk midplane perpendicular to \(x\) in the direction of Galactic rotation at the Sun’s position, and \(z\) increasing towards the direction of the North Galactic Pole. To convert from Galactic coordinates to these Galactocentric coordinates, we need to know the Sun’s position with respect to the Galactic center and with respect to the Galactic midplane. Neither of these quantities are perfectly known at the present time and they are uncertain enough (especially the distance to the center) that the uncertainty affects our interpretation of Galactic dynamics. Roughly, the Sun is located about \(z_\odot = 25\) pc above the Galactic midplane and about \(R_0 = 8.1\) kpc from the Galactic center. Note that because the Sun is not located in the Galactic midplane, the \(b=0\) plane does not exactly line up with the Galactic midplane and a slight rotation (through an angle that is approximate equal to \(z_\odot/R_0\)) is necessary to go between the two coordinate systems (see here for some additional explanation of this). Again, we can use astropy’s coordinate-transformation utilities to convert positions to the Galactocentric coordinate frame. One subtlety is that astropy uses a right-handed coordinate system instead of the left-handed, where the Sun is located at negative \(x\) and \(x\) increases away from the Sun. This means that the \(x\) coordinate is flipped with respect to the left-handed frame. For the \((R_0,z_\odot)\) pair given in this paragraph, we can then transform to the right-handed Galactocentric frame as (using the same example as we have been using throughout this section):

ra= 120.*u.deg
dec= 30.*u.deg
distance= 1.2*u.kpc
c= apycoords.SkyCoord(ra=ra,dec=dec,distance=distance,frame='icrs')
gc = c.transform_to(apycoords.Galactocentric)
print("(x,y,z) in (kpc,kpc,kpc) in right-handed frame")
print("\t",gc.cartesian)

(x,y,z) in (kpc,kpc,kpc) in right-handed frame
         (-9.34605497, -0.20890406,  0.57657821) kpc

and the position in the left-handed Galactocentric frame is the same, except that the sign of \(x\) is flipped (and is therefore \(x=9.34605497\) kpc in this example. To get the cylindrical Galactocentric coordinates, we set the representation_type to be cylindrical:

ra= 120.*u.deg
dec= 30.*u.deg
distance= 1.2*u.kpc
c= apycoords.SkyCoord(ra=ra,dec=dec,distance=distance,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 = 'cylindrical'
print("(R,phi,z) in (kpc,deg,kpc)")
print("\t",gc.rho.to(u.kpc).value,gc.phi.degree,gc.z.to(u.kpc).value)

(R,phi,z) in (kpc,deg,kpc)
         9.148537291648713 -178.69155460938154 0.5744039278619298

Here the angle is again in the right-handed coordinate system. \(180^\circ-\) this angle gives the angle in the left-handed coordinate system.

Of course, we could also start from Galactic coordinates (starting from Galactic coordinates approximately equal to those above):

glon= 191.*u.deg
glat= 27.*u.deg
distance= 1.2*u.kpc
c= apycoords.SkyCoord(glon,glat,distance=distance,frame='galactic')
gc_frame= apycoords.Galactocentric(galcen_distance=8.1*u.kpc,
                                   z_sun=25.*u.pc)
gc = c.transform_to(gc_frame)
gc.representation_type = 'cylindrical'
print("(R,phi,z) in (kpc,deg,kpc)")
print("\t",gc.rho.to(u.kpc).value,gc.phi.degree,gc.z.to(u.kpc).value)

(R,phi,z) in (kpc,deg,kpc)
         9.150114090588534 -178.7224064420817 0.5730235993485637

A.3. Velocities in the Milky Way

\label{sec-mwvel}

Observed velocities are typically reported as proper motions on the sky, for the part of the motion that is in the plane of the sky, and the line of sight velocity, the velocity towards or away from us. Proper motions are measured by comparing the position of a star on the sky at two different epochs, typically spaced by at least a few years and up to decades (although astrometric satellites like Gaia measure the sky positions of celestial sources at a much higher cadence and thus have better time resolution than this). Line-of-sight velocities can be measured using the Doppler shifts of spectral lines in a spectrum taken of the light of the source. Proper motions are therefore angular velocities, typically reported in mas/yr, while line-of-sight velocities are measured as Doppler shifts, which are fractions of the speed of light and the velocity in km/s can therefore be determined. To determine the full physical velocity, the angular proper motions must be multiplied with the distance. Note that to avoid confusion with the radial velocity in the cylindrical Galactocentric frame, we will attempt to always refer to the Doppler-shifted velocity as the line-of-sight velocity, but in the literature this is also commonly known as the radial velocity.

Just like the positions of celestial sources are typically reported as RA and Dec in the equatorial system, proper motions are typically reported as proper motions \(\mu_{\alpha,*}\) and \(\mu_\delta\) in \(\alpha=\)RA and \(\delta=\)Dec, respectively. These are defined based on displacement \(\Delta \alpha\) and \(\Delta \delta\) in RA and Dec, respectively, over a period of time \(T\) as

\begin{align} \mu_{\alpha,*} & = \Delta \alpha \,\cos\delta\,/\,T\,,\\ \mu_\delta & = \Delta \delta\,/\,T\,. \end{align}

We have added an asterisk in the subscript for \(\mu_{\alpha,*}\) because of the \(\cos \delta\) factor that is present. This factor is necessary to transform the observed coordinate displacement in RA into a physical displacement. Virtually all proper-motion catalogs report \(\mu_{\alpha,*}\) even though they often will not include the asterisk. Make sure that you are sure whether or not the \(\cos \delta\) factor is included before making use of any proper motions! If you want to use proper motions to compute the present-day celestial coordinate position in RA and Dec from its position at, say, the year 2000, make sure to divide out the \(\cos\delta\) factor before applying the displacement, that is, you want the coordinate displacement \(\Delta \alpha = \mu_{\alpha,*}\,T/\cos\delta\).

Velocities can be transformed from the equatorial to the Galactic coordinate frame by applying a few simple rotation matrices. This is rather cumbersome, but explained in detail in these notes if you are interested in the specifics. Because the Galactic frame is simply a rotation of the equatorial frame, proper motions can be transformed between the two without involving the line-of-sight velocity. astropy can also transform velocities (at least in versions >= 2). Because this is a new feature it is not yet available through the SkyCoord class that we have used above (this is to allow user feedback to shape the eventual implementation through the SkyCoord class). Therefore, we need to directly instantiate one of the frames, e.g., the ICRS equatorial frame. For example, to transform the proper motion, do

ra= 120.*u.deg
dec= 30.*u.deg
pmra= 5.*u.mas/u.yr
pmdec= -3.*u.mas/u.yr
c= apycoords.ICRS(ra=ra,dec=dec,
                  pm_ra_cosdec=pmra,
                  pm_dec=pmdec)
print(c.transform_to(apycoords.Galactic))

<Galactic Coordinate: (l, b) in deg
    ( 191.27460843,  27.07427141)
 (pm_l_cosb, pm_b) in mas / yr
    ( 4.34640069,  3.88700412)>

As already discussed above, the velocities in the cartesian Galactic coordinate frame as commonly known as \((U,V,W)\). These can be calculated using astropy as

ra= 120.*u.deg
dec= 30.*u.deg
distance= 1.2*u.kpc
pmra= 5.*u.mas/u.yr
pmdec= -3.*u.mas/u.yr
vlos= 55.*u.km/u.s
c= apycoords.ICRS(ra=ra,dec=dec,
                  distance=distance,
                  pm_ra_cosdec=pmra,
                  pm_dec=pmdec,
                  radial_velocity=vlos)
cg= c.transform_to(apycoords.Galactic)
cg.representation_type= 'cartesian'
print(cg)

<Galactic Coordinate: (u, v, w) in kpc
    (-1.04788016, -0.20890423,  0.54617413)
 (U, V, W) in km / s
    (-33.32415106, -31.85479021,  44.72141953)>

As for positions in the Galactic coordinate frame, these velocities in the Galactic coordinate frame are only useful when studying the kinematics of the solar neighborhood. For example, studies of the local velocity distribution (e.g., Bovy et al. 2009) will report the distribution of the \((U,V,W)\) components. However, to study the dynamics of the Milky Way on a global scale, it is necessary to go to the Galactocentric frame again. For this transformation we require the Sun’s three-dimensional velocity with respect to the Galactic center. There are many subtleties in how this velocity is determined and it remains a poorly measured quantity, especially the component in the direction of Galactic rotation. We will come back to the measurement of the Sun’s velocity in later chapters. For now we will just state that the Sun’s velocity with respect to the Galactic center is approximately \((v_x,v_y,v_z) = (-11.1,245,7.25)\,\mathrm{km\,s}^{-1}\) in the left-handed Galactocentric frame (\(v_x = 11.1\,\mathrm{km\,s}^{-1}\) in the right-handed frame; to specify \(v_x\) when using astropy we have to use the right-handed frame). We can then transform velocities all the way to Galactocentric coordinates as

ra= 120.*u.deg
dec= 30.*u.deg
distance= 1.2*u.kpc
pmra= 5.*u.mas/u.yr
pmdec= -3.*u.mas/u.yr
vlos= 55.*u.km/u.s
c= apycoords.ICRS(ra=ra,dec=dec,
                  distance=distance,
                  pm_ra_cosdec=pmra,
                  pm_dec=pmdec,
                  radial_velocity=vlos)
# Define GC frame
v_sun = apycoords.CartesianDifferential([11.1,245.,7.25]*u.km/u.s)
gc_frame= apycoords.Galactocentric(galcen_distance=8.1*u.kpc,
                                   z_sun=25.*u.pc,
                                   galcen_v_sun=v_sun)
cg= c.transform_to(gc_frame)
cg.representation_type= 'cartesian'
print("Cartesian velocity (v_x,v_y,v_z) in km/s")
print("\t",
      cg.v_x.to(u.km/u.s).value,
      cg.v_y.to(u.km/u.s).value,
      cg.v_z.to(u.km/u.s).value)
cg.representation_type= 'cylindrical'
print("Cylindrical velocity (v_r,v_phi,v_z) in km/s=")
print("\t",
      cg.d_rho.to(u.km/u.s).value,
      (cg.d_phi*cg.rho).to(u.km/u.s,
        equivalencies=u.dimensionless_angles()).value,
      cg.d_z.to(u.km/u.s).value)

Cartesian velocity (v_x,v_y,v_z) in km/s
         -22.08605720033482 213.14517509768976 52.07396464035026
Cylindrical velocity (v_r,v_phi,v_z) in km/s=
         17.213193237684614 -213.5939268065537 52.07396464035026

Note again that the sign of \(v_x\) is flipped with respect to what it would be in the left-handed coordinate frame. Similarly, in the cylindrical representation the sign of \(v_\phi\) is flipped.