THE TILT OF THE HALO VELOCITY ELLIPSOID AND THE SHAPE OF THE MILKY WAY HALO

We construct our sample of subdwarfs using data from the sixth SDSS data release (Adelman-McCarthy et al. 2008), in particular utilizing the light-motion catalog of Bramich et al. (2008). This is built from the multiepoch, multiband (u, g, r, i, z) photometry available for one of the SDSS equatorial stripes (Stripe 82) and covers ∼250 deg² in the right ascension range 207 < α < 33 and in the declination range |δ| < 126. A full description of our subdwarf sample is given elsewhere (Smith et al. 2009). Here, we give a brief outline of the selection procedure.

The reduced proper motion is defined as

$$\begin{equation} H_{\rm r} = r + 5\, \log _{10} \left(\frac{\mu }{\mbox{\,mas yr$^{-1}$}} \right) - 10, \end{equation} \tag{ 2 }$$

where μ is the proper motion and r is the apparent magnitude in the r band, corrected for extinction using the maps of Schlegel et al. (1998). The reduced proper motion is useful because it is independent of distance, namely

$$\begin{equation} H_{\rm r} = M_{\rm r} + 5\, \log _{10} \left(\frac{v_{\tan }}{4.74\; \mbox{km s$^{-1}$}} \right), \end{equation} \tag{ 3 }$$

where M_r is the absolute magnitude in the r band and v_tan is the tangential component of the velocity with respect to the line of sight between the Sun and the star. The reduced proper motion diagram is a plot of H_r versus color g−i, in which the populations of disk dwarf stars, white dwarfs, and the halo subdwarfs are all separated (see, e.g., Vidrih et al. 2007, who have already constructed a reduced proper motion diagram for Stripe 82 to isolate ultracool and halo white dwarfs).

Here, we are interested in selecting a clean sample of halo subdwarfs, and so we apply two pre-selection cuts in order to reduce the contamination.⁴ First, we only use stars that pass a quality cut such that uncertainties in the proper motion are smaller than 4 masyr^-1. Secondly, the r magnitude of the star should be brighter than r = 19.5. This gives us a sample size of 370,000. The latter cut allows us to remove interlopers which may be at large distances and hence have small H_r despite having μ ∼ 0. Neither cut introduces any kinematic bias. Note that we avoid cutting on μ directly in order to simplify our calculation of the detection efficiency.

In Figure 1, we show the reduced proper motion diagram. Owing to larger tangential velocities, the halo subdwarfs are clearly differentiated from the slower moving disk dwarfs. We estimate the location of the subdwarf boundary as

$$\begin{eqnarray} H_{\rm r} < 2.85 (g-i) + 11.8\qquad & \mbox{for }& (g-i) \le 2 \nonumber \\ H_{\rm r} < 5.63 (g-i) + 6.24\qquad & \mbox{for }& (g-i) > 2 \nonumber \\ H_{\rm r} > 2.85 (g-i) + 15.0\qquad & \mbox{for }& (g-i) \le 1.3 \nonumber \\ H_{\rm r} > 5.63 (g-i) + 11.386\qquad & \mbox{for }& (g-i) > 1.3, \end{eqnarray} \tag{ 4 }$$

and further clean the sample by rejecting all objects, which are within ΔH_r = 0.5 of the boundary. These cuts result in a sample of 27,000 stars. We then cross-match with the output from the SDSS SEGUE spectral parameter pipeline (Lee et al. 2008). This provides us with the radial velocities for members of the subdwarf sample (with a median error of 4.6 kms^-1), although it significantly reduces the sample size to 2210 stars. The SEGUE spectroscopic target selection is a complicated function. Suffice it to say that, despite its complexity, it is believed to be free of any significant kinematic biases.

Zoom In Zoom Out Reset image size

In order to recover the full six-dimensional phase-space information, we must determine the distance to each subdwarf. Ivezić et al. (2008) have already constructed a photometric parallax relation from SDSS observations of clusters, obtaining an intrinsic scatter of ∼0.2mag (Ivezić et al. 2008; Jurić et al. 2008; Sesar et al. 2008). Their parallax relation is also a function of [Fe/H], which we obtain from the spectral parameter pipeline (with median error of 0.14 dex). Overall, the errors obtained from the Ivezić et al. (2008) relation are impressively small,⁵ aided by the high precision photometry from Bramich et al. (2008)—for our sample, the median relative distance error is 10.1%.

In conjunction with the good accuracy on the Bramich et al. (2008) proper motions (less than 4 masyr^-1), we are able to construct an unprecedentedly large sample of halo subdwarfs with accurate positions and kinematics. The median error on each of our velocity components is ∼30–40 kms^-1. Restricting ourselves to subdwarfs with heliocentric distances less than 5 kpc gives a final sample of 1754. These stars lie at Galactocentric cylindrical polar radii between 7 and 10 kpc, and at depths of 4.5 kpc or less below the Galactic plane.

A drawback to our sample is that it is not kinematically unbiased, as the cut in the reduced proper motion diagram is implicitly a function of the kinematics. Therefore, we have to model and understand the effect of this if we are to investigate the distributions of velocities in our sample. Note that the kinematic bias comes solely from our cut on the reduced proper motion H_r, which actually selects stars via their tangential velocities rather than their proper motions. This makes the task of quantifying the bias significantly easier since we do not need to make any assumptions about the underlying distance distribution (i.e., luminosity function).

We calculate our detection efficiency as follows. For each subdwarf in our final sample, we take the sky coordinates and create a mock sample of 50,000 fake stars. Then, for each mock star, we select M_r and (g − i) at random from our observed distributions. Note that for each realization, both the magnitude and color are assigned from one star, that is, we do not assign an absolute magnitude from one star and a color from another. We then select kinematics for each mock star using the Galactocentric halo velocity distributions from Kepley et al. (2007) but with no net rotation (Allende Prieto et al. 2006). To calculate v_tan from the Galactocentric velocities requires us to assign a distance to each mock star, which we do at random from the observed distribution. This means that the efficiency does have a dependence on the distance distribution. However, this dependence is very mild since Equation (3) is a function of the tangential velocity rather than the proper motion. The efficiency is then given by the fraction of mock stars that pass our reduced proper motion cut.

In order to check whether our results are dependent on the assumed halo velocity distribution, we repeat the calculations using the values from Kepley et al. (2007), but now assuming a rotational velocity of the halo of ∼20 kms^-1 (in the direction of disk rotation). We find that uncertainties in the efficiency make little difference to the final determination of the tilts.

With the efficiencies in hand, we can now calculate the misalignment of the velocity ellipsoid of the SDSS subdwarfs. To do this, we transform our subdwarf velocities into Galactocentric spherical polar coordinates: v_r is the radial velocity with respect to the center of the Galaxy, v_θ is the zenithal component measured from the North Galactic Pole, and v_ϕ is the azimuthal component measured such that the Galactic rotation has negative v_ϕ.

The misalignment from the spherical polar coordinate surfaces can then be described by the correlation coefficients and the tilt angles using the following formula:

$$\begin{equation} {\rm Corr}[v_i,v_j] = \frac{ \sigma _{ij}^2 }{ \big(\sigma _{ii}^2\sigma _{jj}^2\big)^{1/2} } \end{equation} \tag{ 5 }$$

and

$$\begin{equation} \tan (2\alpha _{ij}) = \frac{ 2\sigma _{ij}^2 }{ \sigma _{ii}^2 - \sigma _{jj}^2 }. \end{equation} \tag{ 6 }$$

Here, the tilt angle corresponds to the angle between the i-axis and the major axis of the ellipse formed by projecting the three-dimensional velocity ellipsoid onto the ij-plane (see, e.g., Binney & Merrifield 1998, or the Appendix of this paper). We use the tilt angles to specify the orientation of the velocity ellipsoid, instead of alternatives such as the Euler angles because the former are a natural extension of the familiar two-dimensional case and much easier to visualize than other options. See also recent examples of using the tilt angles in a similar context by Dehnen & Binney (1998) and Siebert et al. (2008).

The measured sample covariances are due to both the true underlying covariances and the correlated measurement uncertainties, that is,

$$\begin{equation} {\rm Cov}_{\rm m}[v_i, v_j]= \sigma _{ij}^2 + {\rm Cov}[\delta v_i, \delta v_j], \end{equation} \tag{ 7 }$$

where Cov_m[v_i, v_j] is the covariance as measured from the sample and Cov[δv_i, δv_j] is the covariance of the error distributions. To account for the detection efficiency, we also calculate the sample covariances weighted by the inverse efficiency, that is,

$$\begin{equation} {\rm Cov}_{\rm m}[v_i,v_j] = \frac{1}{W} \sum _{k=1}^N w_k\, (v_{i,k} - \langle v_{i,k}\rangle)\, (v_{j,k} - \langle v_{j,k} \rangle), \end{equation} \tag{ 8 }$$

where

$$\begin{eqnarray} &&W =1-\sum _{k=1}^N w_k^2,\nonumber \end{eqnarray}$$

and the summation is over the number of subdwarfs in our final sample. The weights w_k are proportional to reciprocal of the efficiency ∝⁻¹_k (normalized to unity such that ∑^N_k=1w_k = 1).

As indicated by Equation (8), when the covariances are calculated we subtract the mean velocities, i.e., our results are not dependent on the assumed reference frame. However, as can be seen from Smith et al. (2009), our sample displays no significant net motion.

Let us first look at the tilt angles in the sample ignoring any stars with height |z| < 1 kpc, so as to minimize (although not eliminate) any effect due to the Galactic disk potential. This gives us a sample of 1532 subdwarfs. We find that the correlations and tilts are

$$\begin{eqnarray} {\rm Corr} [v_{\rm r}, v_\theta ] &= 0.078\;\pm \;0.029,\quad & \alpha _{{\rm r}\theta }= 3\mbox{$.\!\!^\circ $}4\;\pm \;1\mbox{$.\!\!^\circ $}3 \nonumber \\ {\rm Corr} [v_{\rm r}, v_\phi ] &= -0.028\;\pm \; 0.039,\quad & \alpha _{{\rm r}\phi }= -2\mbox{$.\!\!^\circ $}2\;\pm \;3\mbox{$.\!\!^\circ $}3 \\ {\rm Corr} [v_\phi, v_\theta ] &= -0.087\;\pm \; 0.047,\quad & \alpha _{\phi \theta }= -37\mbox{$.\!\!^\circ $}4\;\pm \;20\mbox{$.\!\!^\circ $}4, \nonumber \end{eqnarray} \tag{ 9 }$$

where the errors are obtained using the bootstrap technique.⁶ We find no evidence of any clear tilt in the α_rϕ and α_ϕθ terms. However, the tilt angle α_rθ is measured to be nonzero at about the 3σ level, though it is still very small. The good alignment of the velocity ellipsoid in spherical polars is also apparent from the velocity distributions in the (v_r, v_θ) and (v_r, v_ϕ) planes illustrated in Figure 2, in which the dashed lines show the orientation of the tilts.

Zoom In Zoom Out Reset image size

Note that the tilt angle α_ϕθ is not well constrained since the dispersions in the two components are similar and the covariance is small. Hence, the calculation of α_ϕθ involves the division of one small number by another small number, giving a large error. However, this term is not so important for the purpose of constraining the shape of the Galactic potential (see the following section).

The area of sky covered by the Stripe 82 catalog is 249 deg². It is reasonable to ask whether the velocity dispersion ellipsoid locally could have a different behavior than that averaged over a larger region, especially as this effect has been observed recently in numerical simulations (Zemp et al. 2009). To study this, we investigate the variation in tilt as a function of height from the plane for |z| ⩽ 4.5 kpc in Figure 3. Neither of the tilt angles α_rθ and α_rϕ shows any significant trends. Although there are one or two isolated points at which the tilt angles lie at ∼1σ from the mean, they are localized in longitude, indicative probably of the kinematic substructure and streams which are expected in the stellar halo in hierarchical assembly models. The nature of the substructure in the SDSS subdwarfs is addressed in detail elsewhere (Smith et al. 2009). We have also carried out this test by splitting our sample into three according to right ascension, and find that the values of the tilt angles are consistent with those derived for the full sample.

Zoom In Zoom Out Reset image size

The alignment of the velocity ellipsoid of halo stars in the spherical polar coordinate has substantial implications for the overall potential of the Galaxy. This constraint does not come from the Jeans equations, which merely require that the momentum flux balances the gravitational forces (see Evans et al. 2009 or An & Evans 2009 for recent applications). Rather, the constraint comes from the deeper requirement that a phase-space distribution function must exist. The theorem has been known for some time, although its widespread applicability has been obscured by the fact that Eddington (1915) and Chandrasekhar (1939) introduced unnecessary assumptions in its proof, as realized first by Lynden-Bell (1962).

Let us first note that there are a number of trivial ways that permit the velocity dispersion tensor to be aligned in spherical polar coordinates. The simplest is to ask for the distribution function f to depend on energy E alone, in which case the velocity dispersion tensor is everywhere isotropic. Or, we could insist that the distribution function is given as $f= f(E,|\mbox{\boldmath $L$}|)$ in a spherical potential orf = f(E, L_z) in an axisymmetric potential. Here, $\mbox{\boldmath $L$}$ is the angular momentum, whilst L_z is the component of $\mbox{\boldmath $L$}$ that is parallel to the symmetry axis. All these options are not available to us because of the well-known and long-established triaxiality of the velocity dispersion of halo stars with σ²_rr > σ²_ϕϕ > σ²_θθ—see, for example, recent determinations by Kepley et al. (2007) or Smith et al. (2009).

To generate the observed triaxial anisotropy of the velocity dispersion tensor, the phase-space distribution function must depend on at least three independent integrals of motion I_i, one of which may be chosen to be the same as the Hamiltonian, I₁ = E(v²_r, v²_θ, v²_ϕ; r, θ, ϕ). Here, it is assumed that the reference frame is chosen such that there is no net bulk streaming motion, that is, 〈v_r〉 = 〈v_θ〉 = 〈v_ϕ〉 = 0. Therefore, σ²_rr = 〈v²_r〉, σ²_rθ = 〈v_rv_θ〉 and so on. Accordingly, if the cross-terms 〈v_rv_θ〉 and 〈v_rv_ϕ〉 vanish everywhere, then any additional integrals of motion that isolate the distribution function must be even functions of v_r. This is equivalent to the statement that the v_r-dependence of the integrals is only through the square of the radial velocity component v²_r, and independent of the sign of v_r. Then, the second isolating integral I₂ can be recast as a globally defined function independent of v²_r using the energy integral, that is, I₂ = I₂(E; v_θ, v_ϕ; r, θ, ϕ). However, since I₂ is an integral of motion, its Poisson bracket with the Hamiltonian must vanish, from which it follows that I₂ must also be independent of r as well (Lynden-Bell 1962). Hence, the integral can always be cast in the form I₂ = I₂(E; v_θ, v_ϕ; θ, ϕ) and so is completely independent of both r and (its conjugate momentum) v_r.

This implies that the radial coordinate in the Hamilton–Jacobi equation must separate, and so the potential has to be of the form

$$\begin{equation} \psi (r,\theta,\phi) = \psi _0(r) + {\Psi (\theta,\phi) \over r^2}, \end{equation} \tag{ 10 }$$

where ψ₀ and Ψ are arbitrary functions of the indicated arguments. In fact, by exactly the same line of reasoning applied to the third isolating integral I₃, if 〈v_θv_ϕ〉 also vanishes, then the potential has to have the form

$$\begin{equation} \psi (r,\theta,\phi) = \psi _0(r) + \frac{\xi (\theta)}{r^2} + \frac{\zeta (\phi)}{r^2 \sin ^2\!\theta }, \end{equation} \tag{ 11 }$$

although this is a stronger result than we will need here.

If the potential has the form (10), then Poisson's equation implies that the total density of stars and dark matter is

$$\begin{equation} \rho (r,\theta,\phi) = \rho _0(r) + \frac{\Omega (\theta,\phi)}{4\pi Gr^4}, \end{equation} \tag{ 12 }$$

where

$$\begin{equation*} \rho _0=\frac{1}{4\pi Gr^2}\frac{d}{dr}\left(r^2\frac{d\psi _0}{dr}\right) \end{equation*}$$

$$\begin{equation*} \Omega =\frac{1}{\sin \theta }\frac{\partial }{\partial \theta } \left(\sin \theta \,\frac{\partial \Psi }{\partial \theta }\right) +\frac{1}{\sin ^2\!\theta }\frac{\partial ^2\Psi }{\partial \phi ^2}. \end{equation*}$$

That is to say, the dipole potential is necessarily associated with an astrophysically unrealistic density cusp diverging as r⁻⁴ unless Ω = 0. Consequently, we have Ψ = 0.⁷ This leaves us with the theorem that

• If a steady state stellar population has a nondegenerate (i.e., triaxial) velocity dispersion tensor whose eigenvectors are everywhere aligned in spherical polar coordinates, then the underlying gravitational potential must be spherically symmetric.

This theorem is implicit in Lynden-Bell (1962), whereas restricted versions were known to Eddington (1915) and Chandrasekhar (1939). In fact, from the preceding proof, the crux of the result lies in the existence of the second integral that forces the separation of the radial part of the Hamilton–Jacobi equation, and the presence of the third integral is of a secondary importance. Hence, the theorem can be relaxed somewhat:

• If the potential is nonsingular, it is a sufficient condition for a spherical symmetry that one of the nondegenerate eigenvectors of the velocity dispersion tensor is aligned radially everywhere.

This follows, as spherical symmetry is guaranteed for nonsingular potentials once the radial part of the Hamilton–Jacobi equation can be separated off.

Of course, spherical alignment of the velocity dispersion tensor does not imply that the stellar density itself is spherically symmetric. In fact, many investigators have already found that σ²_θθ ≠ σ²_ϕϕ for samples of halo stars (e.g., Gould 2003; Kepley et al. 2007; Smith et al. 2009), which implies that the stellar halo has a density distribution that is flattened or triaxial even though the gravity field is nearly spherical. This is the case, for example, in the models presented by White (1985), Arnold (1990), and Evans et al. (1997), in which the phase-space distribution function depends on energy and the all three components of angular momentum, $f = f(E, \mbox{\boldmath $L$})$.

The tilt angles α_rϕ and α_ϕθ are consistent with zero, but our results give a very small, but nonzero measurement of the tilt angle α_rθ (at the 3σ level). In an exactly spherical potential, the 〈v_rv_θ〉 cross-term and hence α_rθ must vanish. However, even if the dark halo is spherical, the Galactic potential is not spherical due to the influence of the bulge and the disk. Of course, at the distances probed by our sample of SDSS subdwarfs, it is the halo that dominates the gravitational potential, whilst the disk gives the main perturbation.

The main consequence of the disk is to convert the planar rosette orbits of a spherical potential into the short-axis tube orbits of a mildly oblate potential. However, in typical axisymmetric Galactic potentials, the key 1:1 resonance (in cylindrical coordinates R and z) also occurs at Galactocentric radii between 2 and 10 kpc, causing the axial orbits to become unstable to out-of-plane perturbations and siring the family of banana or saucer orbits (see, e.g., Pfenniger 1984; Miralda-Escudé & Schwarzschild 1989; Schwarzschild 1993; Evans 1994; Binney & Tremaine 2008, Section 3.7.3). It is to these two orbital families—the short-axis tubes and the bananas—that our SDSS subdwarfs will belong.

We can confirm this by computing orbits in the Milky Way potential of Fellhauer et al. (2006), which comprises a spherical isothermal halo and a Miyamoto-Nagai disk. Our initial conditions were chosen as follows: since the potential is axisymmetric, we start our orbits off from a location x = 8 kpc, y = 0 kpc, and with z distributed uniformly from −1 to −4 kpc (which is a reasonable approximation of our observed distribution). Our initial velocities were chosen according to the trivariate ellipsoidal Gaussian σ_r = 142 kms^-1, σ_ϕ = 81 kms^-1, and σ_θ = 77 kms^-1, that is, according to the values derived for our observed sample of halo subdwarfs by Smith et al. (2009). This is of course not a bona fide distribution function, but just a convenient sampling function to scan phase space. The orbits were calculated using a fourth-order Runge–Kutta method with a time step of 0.1 Myr, which results in a fractional energy change of better than 10⁻⁶ over a typical orbital period of ∼200 Myr. The orbits so produced are indeed predominantly short-axis tubes (∼80%), but with a reasonable mixture of banana orbits (∼20%). Some examples are illustrated in the panels of Figure 4.

Zoom In Zoom Out Reset image size

It is immediately apparent on geometrical grounds that the banana orbits can only give positive contributions to 〈v_rv_θ〉. By contrast, the short-axis tubes can yield both positive and negative contributions. Hence, the effect of the Galactic disk is to create a family of orbits for which the tilt term σ²_rθ cannot exactly vanish.

In fact, it is easy to bolster this geometric argument with a quantitative calculation. Using our sampling function, we generated 96 initial conditions and integrated the orbits for 1 Tyr, recording the velocities each time a test particle returned to its initial position (defined in practice to within a spherical volume of 300 pc radius). The enormous timescale is to maximize the number of times a particle returns to its initial position, so that σ²_rθ for the orbit can be calculated as accurately as possible. Then, we average the resulting values of σ²_rθ for every crossing and for all of our orbits. This gives 〈v_rv_θ〉 = 690 km²s^-2, which is very comparable to our observed covariance of 893 ± 335 km²s^-2.

However, if we divide our orbits into banana and short-axis tubes (banana orbits are defined as those which are asymmetric about z = 0 or the Galactic plane) and recalculate the average σ²_rθ, then we find that the 19 banana orbits have 〈v_rv_θ〉 = 1740 km²s^-2, while the remaining 77 short-axis tube orbits have 〈v_rv_θ〉 = −140 km²s^-2. The banana orbits indeed contribute almost all of the effect. As a check, if we repeat this procedure in a spherical potential (i.e., without the Galactic disk), we indeed find that 〈v_rv_θ〉 is consistent with zero, as it should be.

A very mild misalignment in the velocity ellipsoid of the SDSS subdwarfs with respect to spherical polar coordinates is therefore naturally explained by the effects of the Galactic disk.