Dynamical models of the Galaxy

I discuss the importance of dynamical models for exploiting survey data, focusing on the advantages of"torus"models. I summarize a number of applications of these models to the study of the Milky Way, including the determination of the peculiar Solar velocity and investigation of the Hyades moving group.


Introduction
Current studies of the structure of the Milky Way are dominated by a series of major observational programs running from ESA's Hipparcos mission [1], through major photometric surveys such as the SDSS [2], spectroscopic surveys such as RAVE [3], to ESA's scheduled Gaia mission [4], which aims to return photometric and astrometric data for 10 9 stars and low-dispersion spectra for ∼ 10 8 stars.
Turning these data-sets into a consistent picture of the current structure and the assembly of the whole Galaxy, including the dark-matter content, is an ambitious and important goal. It is likely to be impossible without sufficiently sophisticated models that can be used to interpret the data (for example compensating for the observational biases of the various surveys).
Models of the gross structure of the Galaxy have been produced with varying levels of complexity. Mass models [5,6] simply give the density distribution of the various components of the Galaxy, and thus the Galactic potential. Kinematic models, such as those produced by galaxia [7], specify the density and velocity distributions of the luminous components of the Galaxy, but do not consider the question of whether these are consistent with a steady state in any Galactic potential. The Besançon Galaxy model [8] is primarily a kinematic model (in that it is not constrained by Newton's laws of motion on large scales), with a dynamical element used to determine the vertical structure of the disc. Fully dynamical models [9] have a distribution of stars in phase-space which is made from phase-mixed orbits in a given Galactic potential, and therefore a distribution function (df) which is a function of the integrals of motion, which by Jeans theorem [10] means it is in a steady-state.

Benefits of dynamical models
A previous paper [9] provides a detailed discussion of the benefits of dynamical models. Here I just sum up the main points. a e-mail: p.mcmillan1@physics.ox.ac.uk The most obvious advantage of dynamical models is that, unlike kinematic models, they allow us to deduce the gravitational potential of the Galaxy. Existing mass models were fit to observations under the assumption of near circular orbits for various tracers, which is only suitable for a small subset of astrophysical objects, found close to the Galactic plane. Exploiting richer data-sets requires far more careful modelling. This will allow us to infer the distribution of dark matter in the Galaxy, under the assumption that the Galaxy is approximately in a steady state.
A second, complementary, advantage is that dynamical models allow us to connect objects that we can observe to objects those we can not. This means we can use observations in the solar neighbourhood to learn about the structure of the Galaxy at large. For example [11] showed that if the stellar halo were in virial equilibrium, more than half the stars of the stellar halo would be on orbits that bring them through the solar neighbourhood.
An additional advantage of dynamical models is that the associated df depends only on the three integrals of motion (either explicitly or implicitly), as opposed to the full six dimensions of phase-space.

Torus models
Dynamical modelling has been dominated by Schwarzschild modelling [12], especially for analysing the dynamics of early-type galaxies. This technique involves first integrating a number of orbits in the adopted gravitational potential and then seeking weights for these orbits such that the weighted sum of the orbits reproduces the observational data. More recently, the "made-to-measure" (M2M) technique introduced by [13] has been used to produce Nbody models which can be fitted in a broadly similar way [14,15], though in this case the particle weights (which are effectively the orbit weights) are determined "on-the-fly", rather than after the orbit has been integrated.
Both Schwarzschild and M2M models have dfs which are implicitly function of the integrals of motion, as they are constructed from phase-averaged orbits. An alternative and v φ (right) in the Solar neighbourhood, data from the Geneva-Copenhagen survey assuming the value of v ⊙ given in eq. 2 (histogram), and model from [17] which is a single df f (J) which must fit both the v R and v φ distributions (curves, solid curve is from the full model, dashed curve is thin disk contribution only). The dotted vertical line in the v φ plot is the assumed circular speed at the Sun. Clearly the v φ distributions are significantly offset from one another. The only way to bring the data and model to reasonable agreement is to apply a correction of ∼ 7 km s −1 to V ⊙ . Figure reproduced from [17] with permission. modelling strategy is to produce models which are explicitly functions of the integrals of motion. It is not possible to produce plausible models of the Galactic disc with dfs which are functions of the classical integrals, f (E, L z ), as these models have equal velocity dispersions in the radial and vertical directions, which is not the case in the Solar neighbourhood [16]. However, one can instead use the three orbital "actions", J i , as the integrals of motion and use analytic dfs, f (J), which are appropriate for the Galactic disc (e.g. [17]).
The three actions J i and three conjugate angles coordinates θ i provide canonical coordinates for six-dimensional phase space [18]. The conventional phase space coordinates w ≡ (x, v) are 2π-periodic in the angles. The actions are conserved quantities for any orbit, and the angles increase linearly with time: where the components of Ω are the orbital frequencies.
The major obstacle to using dfs of the form f (J) is that the relationship between phase space coordinates w and J, θ is only known analytically for a very limited set of gravitational potentials, none of which provides a reasonable approximation to the Galactic potential. There are two available approaches for determining the relationship between w and J: -The adiabatic approximation, in which motion in the radial and vertical directions are largely decoupled [17,19]. -Torus modelling, in which the relationship between J, θ and w in an isochrone potential (for which it is known analytically [18]) is numerically distorted to fit the potential of interest [20,21,22].
Comparison of these two approaches has shown that for most purposes they agree to reasonable accuracy up to as far as ∼ 2.5 kpc from the Galactic plane [19]. The adiabatic approximation has the advantage that it does not require specialised torus-fitting computer code, and can straightforwardly determine the value of J for a given w. Torus modelling finds all of the values of w associated with a given J, but can only find J given w as an iterative process. Torus modelling has the advantages that it can tell us about the coupling between different components of motion (e.g. the tilt of the velocity ellipsoid [19]), and that it allows us to find the angle variables θ.

The local standard of rest
As mentioned in Section 2, one major advantage of using dynamical models is that the dimensionality of the models is reduced by the assumption that the Galaxy is made up of phase mixed orbits -f (J) depends on only three actions as opposed to the six dimensions of f (w). The value of this simplification of the model is well illustrated by a recent revision in the peculiar Solar velocity relative to the local standard of rest, v ⊙ . The value of v ⊙ was found by [23] using observations by Hipparcos. The components of velocity towards the Galactic centre (U), in the direction of Galactic rotation (V) and towards the north Galactic pole (W) were analysed separately. The V-component of the Solar velocity is the most difficult to determine as asymmetric drift [18] means that average stellar velocity lags the circular velocity. Stromberg's equation was used by [23] to extrapolate from the observed populations (separated by colour) to a hypothetical population with zero velocity dispersion, which would have zero asymmetric drift. The value of v ⊙ found, U ⊙ , V ⊙ , W ⊙ = (10.00±0.36, 5.25±0.62, 7.17±0.38) km s −1 , (2) was the widely accepted value for over a decade. Figure 1 shows the distribution of stars in v φ and v R in the Solar neighbourhood, determined from Hipparcos observations and the Geneva-Copenhagen survey [24] assuming this value of v ⊙ (histogram), and the best fitting Assembling the Puzzle of the Milky Way Frequencies Ω r & Ω φ (left) and a convenient projection in action space (right) for stars in the solar neighbourhood from a simulation of the accretion of a satellite galaxy, taken t = 7.5 Gyr after the satellite was accreted. The stars form patches in frequency space with spacing ∆Ω φ between patches such that ∆Ω φ t/2π ∼ 1 (the spacing in Ω r is not as simple because a two ranges of θ r correspond to the solar neighbourhood). The same separation into discrete clumps is visible in action because Ω = Ω(J). These figures are taken from [22].
analytic df f (J) from [17] (solid curve). This gives new insight because the distributions in U and V are not considered as if they were independent, but instead it is recognised that a single dynamically consistent df must fit both. The only physically plausible way to bring the model and data v φ distributions into agreement is to apply a correction of ∼ 7 km s −1 to the value of V ⊙ (altering the circular velocity moves both distributions and has a negligible effect). This is ∼ 11σ from the widely accepted value! This contention, that the previously assumed value of V ⊙ must be in error, has since been supported by analysis of Galactic maser sources with measured parallaxes [25], and explained as an error in the application of Stromberg's equation, associated with the metallicity gradient in the disc [26].

Beyond steady-state models
Axisymmetric models with dfs of the form f (J) cannot fully describe the Galaxy, as it is not in fact in a steady state, but they are an important first step towards interpreting observations of our Galaxy, and it is likely to be very fruitful to study observational data for structures that cannot be explained by these models. We can then look for explanations of these structures as signatures of other features, such as the Galactic bar, spiral or warp, or matter that has been accreted. Some examples of using torus models to explore these signatures already exist.

Signatures of accreted satellites
The appearance in angle-action coordinates of an accreted satellite galaxy was explored by [22]. Long after phase mixing has rendered an accreted satellite indistinguishable from the background population in position, there is a strong relationship between the stars' positions and their orbital frequencies (because all the stars were once collected in the same small volume, when they were part of the satellite). This means that a sample of these stars taken from a finite volume is only found in certain small volumes in frequency space.
In Fig 2 I show figures from [22] of the frequencies and actions of stars in a finite volume about the solar position in a simulation of the disruption of a satellite galaxy. In both cases the figure shows the distribution 7.5 Gyr after the satellite was disrupted. The distribution in frequency is clearly divided into finite "patches", and the distribution in action is even more cleanly divided, because Ω = Ω(J) and it provides a more convenient projection of the distribution.
By considering the spacing between the "patches" in frequency space, along with the angles of the individual stars, [22] showed that it was even possible to determine with high accuracy the time since the satellite was disrupted. Similar techniques could even be used to determine the potential of the Galaxy (as the potential must allow the stars to all have come from the same initial satellite). Similar work has been carried out by other authors showing that using the orbital frequencies alone one can achieve some of these objectives, even for cosmological simulations with numerous accreting satellites and a non-static background potential [27].

Signatures of Lindblad resonances
A recent study [28] showed that in addition to having a relationship of the form lΩ r + mΩ φ ≃ const, stars that have recently been trapped at a resonance with a perturber also follow a relationship of the form lθ r + mθ φ ≃ const, where in both cases l and m are integers, with the perturbation having m-fold symmetry and l being −1 for an inner Lindblad resonance and +1 for an outer Lindblad resonance. The Hyades moving group, which is a very strong feature of the local velocity distribution [29], lies around a straight line in the J φ , J r plane (Fig 3), of the sort associated with the condition lΩ r +mΩ φ ≃ const (a resonance line in action space). This finding is indeed consistent with a resonance line for the both the cases l = ±1 and a range of values of m (with the exact details of the resonance lines depending quite sensitively on the details of the Galactic potential assumed).
It was claimed by [28] that the distribution of stars in angle coordinates clearly indicated that the Hyades were associated with an inner Lindblad resonance. This is because, in angle space, the stars of the Hyades moving group are associated with an overdensity in the quantity −θ r +mθ φ for various values of m, and this overdensity did not appear to significantly shift in position as a function of J r .
In [30] I compared the observed structure to torus models, and demonstrated the significant and non-intuitive impact of selection effects. In Fig 3 I show the distribution of Solar neighbourhood stars in the θ r , θ φ plane. Selection effects are responsible for the overall structure of the density distribution, most notably the high densities around θ φ = 0, θ r = 0 or ±π and the near absence of stars with θ r < 0, θ φ > 0 or θ r > 0, θ φ < 0. A less obvious selection effect brought to light by torus models is that stars with a given Fig. 3. Distribution of stars in action (upper) and angle (lower) in the Solar neighbourhood with positions and velocities given by the Geneva Copenhagen survey. The gross structures of the distributions are due to selection effects. The Hyades moving group can be seen as an overdensity in action that is spread out in J r at around J φ = 0.97, tending towards slightly lower J φ with increasing J r (consistent with a resonance line in Ω), and as an overdensity in angle around θ r = −π/2. The expected angle relation for stars at resonance would produce overdensities that tend to increased θ r with increased θ φ for inner Lindblad resonances, and to decreased θ r with increased θ φ for outer Lindblad resonances. However, selection effects reshape overdensities in angle space quite significantly. Figures are adapted from [30] value of J, found in the Solar neighbourhood are also very strongly restricted in their possible range of θ.
Using these models I was able to show that observed overdensities in θ r +mθ φ (which should correspond to outer Lindblad resonances) were not the result of selection effects (as claimed by [28]). Also, these selection effects mean that the stars associated with the resonance line in action space have to lie near certain lines in angle space, otherwise they will not be observed in the Solar neighbourhood. It is the interplay of the resonant conditions on J and θ, and the fact that the stars lie in the Solar neighbourhood that determines the distribution in both angle and action, and considering either distribution independently of the other can lead to false conclusions. The examination of torus models which included resonant components (with resonant conditions on both J and θ) indicated that the Hyades overdensity was consistent with stars trapped at either an inner or outer Lindblad resonance.