© ESO, 2013

1. Introduction

Tracking large-scale orbits within the Milky Way is an important component in the study of the dynamical properties of stars, globular clusters, or satellite galaxies. Given the object’s coordinates and velocity components, this is a straightforward task once the Galactic gravitational potential is available. Owing to a lack of observational constraints, a great deal of the Milky Way’s mass distribution – in particular of its dark matter halo – remains unknown. Therefore, different mass models are able to reproduce the observations equally well. For the purpose of numerical orbit calculations, a mathematically simple and analytically closed potential is preferred because it supports fast computations, which is especially favorable when using Monte Carlo methods for error estimation.

The Galactic gravitational potential of Allen & Santillán (1991) perfectly fulfills these criteria and has thus been widely used to calculate trajectories for, e.g., globular clusters (Allen et al. 2008; Lane et al. 2012; Odenkirchen et al. 1997), dwarf spheroidals (Lépine et al. 2011), planetary nebulae (Wu et al. 2011), white dwarfs (Pauli et al. 2003, 2006), horizontal-branch stars (Altmann & de Boer 2000; Kaempf et al. 2005), subluminous B stars (Altmann et al. 2004; Tillich et al. 2011), halo stars (Pereira et al. 2012; Schuster et al. 2012), runaway stars (Irrgang et al. 2010; Silva & Napiwotzki 2011), and hypervelocity stars (Edelmann et al. 2005; Hirsch et al. 2005). However, the model parameters can be recalibrated using new and improved observational constraints. In this context, a simplified expression for the halo component is introduced as well.

As for similar studies, standard constraints encompassing terminal velocities for the inner rotation curve, circular velocities of maser sources for the outer rotation curve, in-plane proper motion of Sgr A*, local mass/surface density, and the velocity dispersion in Baade’s window are fitted to determine the respective parameters. To constrain the mass at large radii, in particular that of the dark matter halo, the approach of Przybilla et al. (2010) is followed here by requiring the dynamically peculiar halo star SDSSJ153935.67+023909.8 (J1539+0239 for short) to be gravitationally bound to the Milky Way (Sect. 2).

Despite its simplicity, the resulting revised Allen & Santillán model (hereafter denoted Model I) is capable of meeting all of its imposed conditions. In addition to the fast computation of realistic trajectories, calibration of the potential using the star J1539+0239 allows the halo mass to be estimated in a fully consistent way (Sect. 3). To investigate the influence of the applied halo component on the results, the fitting is redone with the original halo mass distribution being replaced by a truncated, flat rotation curve model according to Wilkinson & Evans (1999) (hereafter denoted Model II, Sect. 4) and a Navarro et al. (1997) dark matter halo (hereafter denoted Model III, Sect. 5).

The spatial origin of the hypervelocity star HE 0437-5439 (Edelmann et al. 2005) is investigated as a first application example (Sect. 6). Finally, conclusions are drawn (Sect. 7).

2. Observations and fitting

In this section, all the details of the fitting process are outlined, i.e., what assumptions are made, what observations are taken into consideration, and how the fitting actually is done.

2.1. Introductory remarks

All gravitational potentials Φ(r,z) considered here are the sum of a central spherical bulge component Φ_b(R), an axisymmetric disk Φ_d(r,z), and a massive spherical dark matter halo Φ_h(R): $\mathrm{\Phi }\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathrm{=}{\mathrm{\Phi }}_{\mathrm{b}}\mathrm{\left(}\mathit{R}\mathrm{\right(}\mathit{r,z}\mathrm{\left)}\mathrm{\right)}\mathrm{+}{\mathrm{\Phi }}_{\mathrm{d}}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathrm{+}{\mathrm{\Phi }}_{\mathrm{h}}\mathrm{\left(}\mathit{R}\mathrm{\right(}\mathit{r,z}\mathrm{\left)}\mathrm{\right)}\mathit{.}$(1)Here, (r,φ,z) are cylindrical coordinates and $\mathit{R}\mathrm{=}\mathit{R}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathrm{=}\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}}$ is the spherical radius.

The corresponding density can be derived by virtue of Poisson’s equation: $\mathrm{4}\mathit{\pi G\rho }\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathrm{=}{\mathrm{\nabla }}^{\mathrm{2}}\mathrm{\Phi }\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathit{.}$(2)For the sake of simplicity, all potentials are independent of time. Consequently, rotating, nonaxisymmetric features, such as the Galactic bar or the spiral arms, are not implemented here, because they require that more complex models (see, e.g., Pichardo et al. 2003, 2004) are considered. To account for the effects of bar and spiral arms, additional systematic uncertainties are added to the observational data where necessary, see for instance Sect. 2.2.2.

Sticking to the convention of Allen & Santillán (1991), throughout this paper, the gravitational potentials are expressed in units of 100 km² s^-2, lengths in kpc, and masses in galactic mass units M_gal = 100 × 1000² kpc [m] /G [SI] M_⊙ ≈ 2.325 × 10⁷ M_⊙, yielding a gravitational constant G of unity.

2.2. Observational constraints

Similar to Dehnen & Binney (1998), observations are divided into eight different groups that are introduced and discussed in the following.

2.2.1. Solar kinematics

Analyzing kinematic observations made in the celestial reference frame in a Galactocentric picture requires a coordinate transformation from one system to the other. Since this transformation obviously depends on the Sun’s position and velocity, the interpretation of any data is affected by the values chosen for these parameters. To set up the coordinate transformation, we take J2000.0 coordinates for the Galactic center (GC) ${\mathit{\alpha }}_{\mathrm{GC}}\mathrm{=}{\mathrm{17}}^{\mathrm{h}}{\mathrm{45}}^{\mathrm{m}}\mathrm{37}\stackrel{\mathrm{s}}{.}\mathrm{224}$, ${\mathit{\delta }}_{\mathrm{GC}}\mathrm{=}\mathrm{-}\mathrm{28}°\mathrm{56}\prime \mathrm{10}\stackrel{\mathrm{″}}{.}\mathrm{23}$ and north Galactic pole (NGP) ${\mathit{\alpha }}_{\mathrm{NGP}}\mathrm{=}{\mathrm{12}}^{\mathrm{h}}{\mathrm{51}}^{\mathrm{m}}\mathrm{26}\stackrel{\mathrm{s}}{.}\mathrm{282}$, ${\mathit{\delta }}_{\mathrm{NGP}}\mathrm{=}\mathrm{+}\mathrm{27}°\mathrm{07}\prime \mathrm{42}\stackrel{\mathrm{″}}{.}\mathrm{01}$, as mentioned in the appendix of Reid & Brunthaler (2004). By monitoring stellar orbits around the central supermassive black hole Sgr A*, Gillessen et al. (2009) restricted the distance r_⊙ from the Sun to the GC to the range ${\mathit{r}}_{\mathrm{\odot }}\mathrm{=}\mathrm{8.33}\mathrm{\pm }\mathrm{0.35}\hspace{0.17em}\mathrm{kpc}\mathit{,}$(3)which is the first observational constraint.

The Sun’s velocity is a superposition of the in-plane circular motion v₀ of its local standard of rest (LSR) around the GC and its peculiar motion v_⊙ = (U,V,W)_⊙ relative to the LSR. Hereby, U is the component towards the GC, V in direction of Galactic rotation, and W perpendicular to the Galactic plane. Motivated by the McMillan & Binney (2010) discussion of a systematic offset in the V-component of maser motions versus a nonstandard value of V_⊙, the ${{\mathrm{(}}^{}\mathit{U,V,}{\mathit{W}}^{\mathrm{)}}}_{\mathrm{\odot }}$ values determined by Schönrich et al. (2010) are in this case preferable to others since they allow the use of maser sources to constrain the Galactic rotation curve (see Sect. 2.2.2) without the necessity of modeling an additional net peculiar motion. Thus, ${{v}}_{\mathrm{\odot }}\mathrm{=}\mathrm{\left(}\mathrm{11.}{\mathrm{1}}_{-0.75}^{\mathrm{+}\mathrm{0.69}}\mathit{,}\mathrm{12.2}{\mathrm{4}}_{-0.47}^{\mathrm{+}\mathrm{0.47}}\mathit{,}\mathrm{7.2}{\mathrm{5}}_{-0.36}^{\mathrm{+}\mathrm{0.37}}\mathrm{\right)}\mathrm{\pm }\mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathrm{0.5}\mathrm{\right)}\hspace{0.17em}\mathrm{km}\hspace{0.17em}{\mathrm{s}}^{-1}$ is used here.

Following the argument in McMillan (2011), the proper motion of Sgr A* along Galactic longitude l, μ_SgrA^∗ = −6.379 ± 0.026 mas yr^-1 (Reid & Brunthaler 2004) is utilized as an indicator for v₀ instead of Oort’s constants A and B. With a spatial motion of Sgr A* of 0 ± 1 km s^-1 (McMillan & Binney 2010), this proper motion can be assumed to be solely caused by the solar motion around the GC yielding as second observational constraint: ${\mathit{\mu }}_{\mathrm{Sgr}{\mathrm{A}}^{\mathrm{\ast }}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{{\mathit{r}}_{\mathrm{\odot }}}\left({\mathit{V}}_{\mathrm{\odot }}\mathrm{+}{\mathit{v}}_{\mathrm{0}}\right)\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{{\mathit{r}}_{\mathrm{\odot }}}\left({\mathit{V}}_{\mathrm{\odot }}\mathrm{+}{\mathit{v}}_{\mathrm{c}}\mathrm{\left(}{\mathit{r}}_{\mathrm{\odot }}\mathrm{\right)}\right)\mathit{.}$(4)The circular velocity v_c(r) is hereby linked to the potential via ${\mathit{v}}_{\mathrm{c}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\sqrt{\begin{array}{c}\mathrm{´}\\ \mathit{r}\end{array}\frac{\mathrm{d\Phi }\mathrm{\left(}\begin{array}{c}\mathrm{´}\\ \mathit{r}\end{array}\mathit{,}\mathrm{0}\mathrm{\right)}}{\mathrm{d}\begin{array}{c}\mathrm{´}\\ \mathit{r}\end{array}}}|}_{\begin{array}{c}\mathrm{´}\\ \mathit{r}\end{array}\mathrm{=}\mathit{r}}\mathrm{·}$(5)

2.2.2. Rotation curve

The Galactic rotation curve deduced from terminal velocities in the inner and intermediate regions and from maser sources in the outer region yield the third, fourth, and fifth groups of observational constraints considered here.

Terminal velocities:

the terminal velocity v_term(l) is the measured peak velocity of the interstellar gas along the line-of-sight at Galactic coordinates b = 0 and − 90° < l < 90°. Assuming a circularly rotating interstellar medium, v_term(l) can be computed from the circular velocity v_c at r = r_⊙sin(l) corrected for the observer’s projected motion: ${\mathit{v}}_{\mathrm{term}}\mathrm{\left(}\mathit{l}\mathrm{\right)}\mathrm{=}{\mathit{v}}_{\mathrm{c}}\left({\mathit{r}}_{\mathrm{\odot }}\sin \hspace{0.17em}\mathit{l}\right)\mathrm{-}\left({\mathit{v}}_{\mathrm{c}}\mathrm{\left(}{\mathit{r}}_{\mathrm{\odot }}\mathrm{\right)}\mathrm{-}{\mathit{V}}_{\mathrm{\odot }}\right)\sin \hspace{0.17em}\mathit{l}\mathrm{-}{\mathit{U}}_{\mathrm{\odot }}\cos \hspace{0.17em}\mathit{l}\mathit{.}$(6)For this study, unprocessed terminal velocities are extracted from surveys in H i by Burton & Gordon (1978) and in CO by Clemens (1985). Similar to Dehnen & Binney (1998), a constant uncertainty of 7 km s^-1 is added in quadrature to the stated uncertainties of the Clemens data to account for non circular motions and the inability of the model potential to reproduce, e.g., spiral arm features. Complementary to the CO measurements, the H i data probe the central region with |sin l| < 0.3, which is additionally affected by distortions due to the bar. Therefore, a rather generous uncertainty of 15 km s^-1 is quadratically added to their measurement errors.

Maser observations:

maser sources are associated with high-mass star-forming regions and are consequently good tracers of the kinematic properties of their surrounding interstellar gas. Using very long baseline interferometry techniques, their positions, parallaxes, and proper motions can be measured to very high accuracy even for distant sources giving the opportunity to sample a relatively wide range of Galactocentric radii. Together with heliocentric radial velocities from Doppler shifts, these data hence allow precise determination of the masers’ three-dimensional motions. In contrast to Reid et al. (2009), McMillan & Binney (2010) show that it is fairly justified to adopt their corresponding circular velocity components as probes for the Galactic rotation curve when taking the Schönrich et al. values for the peculiar motion of the Sun. Following this argument, parallaxes, proper motions, and heliocentric radial velocities for 30 maser sources have been compiled from the literature (Asaki et al. 2010; Hirota et al. 2008a,b; Honma et al. 2011; Niinuma et al. 2011; Reid et al. 2009; Rygl et al. 2010; Sandstrom et al. 2007; Sato et al. 2010). To account for the virial motion of the masers with respect to their high-mass star-forming regions, a constant uncertainty of 7 km s^-1 (see Reid et al. 2009) is added in quadrature to the measurement errors when computing circular velocities.

2.2.3. Mass and surface densities

The local dynamical mass density ρ_⊙ and the surface density Σ_1.1 are together the sixth group of observational constraints. Their derivation is straightforward from the model: $\begin{array}{ccc}& & {\mathit{\rho }}_{\mathrm{\odot }}\mathrm{=}{\mathit{\rho }}_{\mathrm{b}}\mathrm{\left(}{\mathit{r}}_{\mathrm{\odot }}\mathrm{\right)}\mathrm{+}{\mathit{\rho }}_{\mathrm{d}}\mathrm{\left(}{\mathit{r}}_{\mathrm{\odot }}\mathit{,}\mathrm{0}\mathrm{\right)}\mathrm{+}{\mathit{\rho }}_{\mathrm{h}}\mathrm{\left(}{\mathit{r}}_{\mathrm{\odot }}\mathrm{\right)}\\ & & {\mathrm{\Sigma }}_{\mathrm{1.1}}\mathrm{=}\underset{\mathrm{-}\mathrm{1.1}\hspace{0.17em}\mathrm{kpc}}{\overset{\mathrm{1.1}\hspace{0.17em}\mathrm{kpc}}{\mathrm{\int }}}\left({\mathit{\rho }}_{\mathrm{b}}\mathrm{\left(}{\mathit{r}}_{\mathrm{\odot }}\mathit{,z}\mathrm{\right)}\mathrm{+}{\mathit{\rho }}_{\mathrm{d}}\mathrm{\left(}{\mathit{r}}_{\mathrm{\odot }}\mathit{,z}\mathrm{\right)}\mathrm{+}{\mathit{\rho }}_{\mathrm{h}}\mathrm{\left(}{\mathit{r}}_{\mathrm{\odot }}\mathit{,z}\mathrm{\right)}\right)\mathrm{d}\mathit{z}\mathit{.}\end{array}$\Holmberg & Flynn obtained ρ_⊙ = 0.102 ± 0.010 M_⊙ pc^-3 from Hipparcos data on a volume-complete sample of A and F stars (Holmberg & Flynn 2000) and Σ_1.1 = 74 ± 6 M_⊙ pc^-2 from observations of K giant stars (Holmberg & Flynn 2004).

2.2.4. Velocity dispersion in Baade’s window

The velocity dispersion of the bulge in Baade’s window, σ_BW, is implemented as the seventh constraint probing the inner most region of the Milky Way. Following the reasoning of Dehnen & Binney (1998), it is given as ${\mathit{\sigma }}_{\mathrm{BW}}\mathrm{=}{\mathit{\sigma }}_{\mathrm{b}}\mathrm{(}\mathrm{0.0175}{\mathit{r}}_{\mathrm{\odot }}\mathit{,}\mathrm{-}\mathrm{0.068}{\mathit{r}}_{\mathrm{\odot }}\mathrm{)}\mathrm{=}\mathrm{117}\mathrm{\pm }\mathrm{15}\hspace{0.17em}\mathrm{km}\hspace{0.17em}{\mathrm{s}}^{-1}\mathit{,}$(9)where σ_b is estimated from the model potential as ${\mathit{\sigma }}_{\mathrm{b}}^{\mathrm{2}}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{b}}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}}\underset{\mathit{z}}{\overset{\mathrm{\infty }}{\mathrm{\int }}}{\mathit{\rho }}_{\mathrm{b}}\mathrm{\left(}\mathit{r,}\begin{array}{c}\mathrm{´}\\ \mathit{z}\end{array}\mathrm{\right)}\frac{\mathit{\partial }\mathrm{\Phi }\mathrm{\left(}\mathit{r,}\begin{array}{c}\mathrm{´}\\ \mathit{z}\end{array}\mathrm{\right)}}{\mathit{\partial }\begin{array}{c}\mathrm{´}\\ \mathit{z}\end{array}}\mathrm{d}\begin{array}{c}\mathrm{´}\\ \mathit{z}\end{array}\mathit{.}$(10)

2.2.5. The extreme halo star J1539+0239

The most stringent constraints on the halo mass result from the kinematically most extreme objects. Usually, the very distant satellite galaxy Leo I is considered to be the most important tracer of the halo mass. Its proper motion has recently been measured with the Hubble Space Telescope (HST) for the first time (Sohn et al. 2012) and allowed to determine a limit on the halo mass (Boylan-Kolchin et al. 2012). However, close-by halo stars may provide even more compelling results. The kinematically most extreme halo star is the blue horizontal-branch star J1539+0239 for which Przybilla et al. (2010) found a Galactic rest-frame velocity of about $\mathrm{69}{\mathrm{4}}_{-221}^{\mathrm{+}\mathrm{300}}\hspace{0.17em}\mathrm{km}\hspace{0.17em}{\mathrm{s}}^{-1}$, making it the fastest halo object known so far. Unlike other hypervelocity stars, however, this object is currently approaching the Galactic disk with its pericentric passage of the GC still to occur. As discussed in the next paragraph, this indicates that J1539+0239 is gravitationally bound to the Milky Way, yielding a so far unexploited and – due to the star’s immense space motion – significant constraint on the Galaxy’s total mass. We show in Sects. 3.2, 4.2, and 5.2 that the limits on the halo mass derived from the kinematics of this star alone are even more stringent than those from Leo I, so we add J1539+0239 to our set of constraints for the potentials and use Leo I as a consistency check. Our approach, which is based on the six-dimensional phase space information of the most extreme objects, is complementary to statistical studies of several thousand stars (see, e.g., Gnedin et al. 2010; Xue et al. 2008) that only make use of line-of-sight velocities.

To justify the assumption that J1539+0239 is bound to the Galaxy, consider the opposite, i.e., a trajectory that passes by the Milky Way only once, never to come back again. Such an orbit could be explained either by a non-Galactic origin, e.g., as a hypervelocity star ejected from another galaxy, or by an extreme dynamical event, e.g., a very close multibody encounter in a globular cluster or satellite galaxy, ejecting a star from its halo environment in a direction pointing not too far away from the GC. Although both scenarios cannot be ruled out by stellar age or space motion arguments, they are unlikely since they require very special circumstances, e.g., an extraordinarily high ejection velocity, a specific ejection direction, or the ejection event occurring at a certain time. Most of these restrictions are weaker or do not apply at all in the case of a bound orbit. For instance, the peculiar motion of J1539+0239 could still be the result of an ejection event. But because the accessible phase space is finite now, there is a non-zero probability of finding a trajectory starting from any initial condition (with matching conserved quantities) sooner or later in a phase-space state very close to J1539+0239, allowing for countless more ejection locations, directions and times. In particular, the GC no longer can be excluded as the star’s spatial origin, which opens up the possibility of a dynamical interaction with the central supermassive black hole (Hills 1988). This mechanism is known to be most powerful in terms of acceleration and is the best candidate for explaining extreme velocity objects. Alternatively, as shown by Abadi et al. (2009), the latter may result from the tidal disruption of a dwarf galaxy, a non ejection scenario that is also able to produce bound (but not unbound), high-velocity stars that approach the Galactic disk.

The assumption that J1539+0239 is gravitationally bound to the Milky Way is the eighth and last observational constraint. With E_∗ denoting the star’s total energy, i.e., the sum of kinetic and potential energy, per unit mass, the standard approach to incorporating this is to require E_∗ ≤ 0. But the Milky Way is not an isolated system, and beyond some distance, the influence of other objects such as Andromeda is no longer negligible. Since those effects are not included in the presented models, a more rigorous approach is chosen here by defining a bound state as ${\mathit{E}}_{\mathrm{\ast }}\mathrm{\le }\mathrm{\Phi }\mathrm{(}\mathrm{0}\mathit{,}\mathrm{200}\hspace{0.17em}\mathrm{kpc}\mathrm{)}\mathit{,}$(11)theoretically limiting orbits to R ~ 200 kpc. However, the actual implementation in the least-squares fit is $\begin{array}{c}\\ {\mathit{\chi }}_{\mathrm{\ast }}\mathrm{=}\frac{{\mathit{E}}_{\mathrm{\ast }}\mathrm{-}\mathrm{\Phi }\mathrm{\left(}\mathrm{0}\mathit{,}\mathrm{200}\hspace{0.17em}\mathrm{kpc}\mathrm{\right)}}{\mathrm{\Delta }{\mathit{E}}_{\mathrm{\ast }}}& \mathrm{if} {\mathit{E}}_{\mathrm{\ast }}\mathrm{-}\mathrm{\Phi }\mathrm{\left(}\mathrm{0}\mathit{,}\mathrm{200}\hspace{0.17em}\mathrm{kpc}\mathrm{\right)}\mathit{>}\mathrm{0}\mathit{,}\\ \mathrm{0}& \mathit{otherwise}\mathit{.}\end{array}$Hereby, ΔE_∗ is the error in E_∗ resulting from uncertainties in determining J1539+0239’s distance, radial velocity, and proper motions. Since ΔE_∗ is large (for instance E_∗/ΔE_∗ ~ 0.82 in the original Allen & Santillán potential), violations of Eq. (11) are generously admitted. For the same reason, the particular choice of 200 kpc is not decisive as for instance − Φ(0,200 kpc)/ΔE_∗ ~ 0.16 while − Φ(0,150 kpc)/ΔE_∗ ~ 0.21 and − Φ(0,250 kpc)/ΔE_∗ ~ 0.13 (again in the original Allen & Santillán potential).

2.3. Fitting process

The model parameters (see Sects. 3–5) are simultaneously determined via a χ²-minimization procedure using the Interactive Spectral Interpretation System (ISIS, Houck & Denicola 2000). Observations of the rotation curve, as well as of J1539+0239, are always taken in raw form from the literature, i.e., as measured in celestial coordinates, in order to exclude inconsistencies during the transformation to the Galactic reference frame originating in different values for the LSR and solar motion. Moreover, any given uncertainty found in the literature is considered here, typically via Gaussian statistics. For instance, to compare the maser measurements to model predictions, their positions and velocities are converted to Galactic coordinates, allowing computation of rotational velocities v_c and Galactrocentric distances r. The corresponding uncertainties in the observational data are then accounted for using Gaussian error propagation. Additionally, the (admittedly small) uncertainties stemming from v_⊙ are propagated in the same way. As the resulting data points in the (r,v_c)-plane have errors in both directions (Δr, Δv_c), it is not obvious at which radius to compare the model to observation, i.e., how to compute the χ² without neglecting Δr. The solution chosen here is to convert Δr to an error in v_c by estimating the effects of it on the model rotation curve according to (Δv_c)_model = (dv_c/dr)_modelΔr and adding it in quadrature to Δv_c.

Very much attention is paid to take all sources of error into consideration in order to assign the correct weighting to each observational constraint in the total χ². Unfortunately, systematic effects such as the missing modeling of spiral arms are difficult to quantify and are not sufficiently dealt with by just adding a constant error. Therefore, an unweighted fitting routine, e.g., assuming that a single maser point weighs the same as the proper motion of Sgr A* (see Eq. (4)), does not yield a satisfactory result because the fit would be dominated by the large number of data points used while fitting the rotation curve. To account for this, each contribution to the total χ² of the above-mentioned eight groups of observational constraints is divided by the number of data points N within the group: ${\mathit{\chi }}_{\mathrm{w}}^{\mathrm{2}}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathrm{8}}\left(\frac{\mathrm{1}}{{\mathit{N}}_{\mathit{i}}}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{i}}}{\mathit{\chi }}_{\mathit{ij}}^{\mathrm{2}}\right)\mathit{.}$(13)The model parameters given in Tables 1–3 are those that minimize this weighted χ².

3. Model I

Model I is a revision of the Allen & Santillán (1991) potential.

3.1. Components and characteristics

3.1.1. Bulge and disk

The potential forms Φ_b(R) and Φ_d(r,z) of bulge and disk are those proposed by Miyamoto & Nagai (1975): $\begin{array}{ccc}& & {\mathrm{\Phi }}_{\mathrm{b}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{{\mathit{M}}_{\mathrm{b}}}{\sqrt{{\mathit{R}}^{\mathrm{2}}\mathrm{+}{{\mathit{b}}_{\mathrm{b}}}^{\mathrm{2}}}}\\ & & {\mathrm{\Phi }}_{\mathrm{d}}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{{\mathit{M}}_{\mathrm{d}}}{\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\left({\mathit{a}}_{\mathrm{d}}\mathrm{+}\sqrt{{\mathit{z}}^{\mathrm{2}}\mathrm{+}{{\mathit{b}}_{\mathrm{d}}}^{\mathrm{2}}}\right)}^{\mathrm{2}}}}\mathrm{·}\end{array}$The parameters M_b/M_d determine the contribution of their components to the total potential, and b_b/a_d/b_d are scale lengths.

The related densities ρ_b and ρ_d are (see Eq. (2)) $\begin{array}{ccc}& & {\mathit{\rho }}_{\mathrm{b}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\frac{\mathrm{3}{{\mathit{b}}_{\mathrm{b}}}^{\mathrm{2}}{\mathit{M}}_{\mathrm{b}}}{\mathrm{4}\mathit{\pi }\mathrm{(}{\mathit{R}}^{\mathrm{2}}\mathrm{+}{{\mathit{b}}_{\mathrm{b}}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{5}\mathit{/}\mathrm{2}}}\\ & & {\mathit{\rho }}_{\mathrm{d}}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathrm{=}\frac{{{\mathit{b}}_{\mathrm{d}}}^{\mathrm{2}}{\mathit{M}}_{\mathrm{d}}}{\mathrm{4}\mathit{\pi }}\frac{{\mathit{a}}_{\mathrm{d}}{\mathit{r}}^{\mathrm{2}}\mathrm{+}\mathrm{(}{\mathit{a}}_{\mathrm{d}}\mathrm{+}\mathrm{3}{\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}}}^{\mathrm{)}}{\mathrm{(}{\mathit{a}}_{\mathrm{d}}\mathrm{+}{\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}}}^{\mathrm{)}}}^{\mathrm{2}}}{{\mathrm{(}{\mathit{z}}^{\mathrm{2}}\mathrm{+}{{\mathit{b}}_{\mathrm{d}}}^{\mathrm{2}}\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}{\left({\mathit{r}}^{\mathrm{2}}\mathrm{+}{\left({\mathit{a}}_{\mathrm{d}}\mathrm{+}\sqrt{\mathrm{(}{\mathit{z}}^{\mathrm{2}}\mathrm{+}{{\mathit{b}}_{\mathrm{d}}}^{\mathrm{2}}\mathrm{)}}\right)}^{\mathrm{2}}\right)}^{\mathrm{5}\mathit{/}\mathrm{2}}}\mathrm{·} \end{array}$Integrating these densities over the entire volume to obtain the total masses m_b and m_d of bulge and disk gives the expected identities m_b = M_b and m_d = M_d.

3.1.2. The dark matter halo

Pursuing mathematical simplicity, the halo potential Φ_h(R) used here differs slightly from the original one by Allen & Santillán. We assume the halo mass m_h inside a sphere of radius R is given by (see Allen & Martos 1986) $\begin{array}{c}\\ {\mathit{m}}_{\mathrm{h}}\mathrm{\left(}\mathit{<}\mathit{R}\mathrm{\right)}\mathrm{=}\frac{{\mathit{M}}_{\mathrm{h}}{\left(\frac{\mathit{R}}{{\mathit{a}}_{\mathrm{h}}}\right)}^{\mathit{\gamma }}}{\mathrm{1}\mathrm{+}{\left(\frac{\mathit{R}}{{\mathit{a}}_{\mathrm{h}}}\right)}^{\mathit{\gamma }\mathrm{-}\mathrm{1}}}& \mathrm{if} \mathit{R}\mathit{<}\mathrm{\Lambda }\mathit{,}\\ \frac{{\mathit{M}}_{\mathrm{h}}{\left(\frac{\mathrm{\Lambda }}{{\mathit{a}}_{\mathrm{h}}}\right)}^{\mathit{\gamma }}}{\mathrm{1}\mathrm{+}{\left(\frac{\mathrm{\Lambda }}{{\mathit{a}}_{\mathrm{h}}}\right)}^{\mathit{\gamma }\mathrm{-}\mathrm{1}}}\mathrm{=}\mathrm{constant}\mathit{,}& \mathit{otherwise}\mathit{.}\end{array}$Here, M_h is again a weighting factor and a_h a scale length. Clearly, this specific form is chosen to have the asymptotic behavior m_h(<R) ∝ R for large R motivated by the observed flat rotation curve. To avoid an unphysical, infinite halo mass, a cutoff parameter Λ is incorporated as well. In principle, the exponent γ is a free parameter, too. However, the shortage of halo constraints means the ability of the model to reproduce the observations is not attenuated by setting γ = 2, reducing the expression’s complexity (and avoiding singularities at the origin in the equations of motion for γ < 2).

The expression m_h(<R) is the volume integral of the halo density ρ_h over the sphere of radius R. Assuming a spherically symmetric halo density ρ_h = ρ_h(R) and using Poisson’s equation (Eq. (2)), a relation between m_h(<R) and its respective potential can be derived: $\begin{array}{ccc}{\mathit{m}}_{\mathrm{h}}\mathrm{\left(}\mathit{<}\mathit{R}\mathrm{\right)}& \mathrm{=}& \mathrm{4}\mathit{\pi }\underset{\mathrm{0}}{\overset{\mathit{R}}{\mathrm{\int }}}\begin{array}{c}\mathrm{´}\\ {\mathit{R}}^{\mathrm{2}}\end{array}{\mathit{\rho }}_{\mathrm{h}}\mathrm{\left(}\begin{array}{c}\mathrm{´}\\ \mathit{R}\end{array}\mathrm{\right)}\mathrm{d}\begin{array}{c}\mathrm{´}\\ \mathit{R}\end{array}\mathrm{=}\underset{\mathrm{0}}{\overset{\mathit{R}}{\mathrm{\int }}}\begin{array}{c}\mathrm{´}\\ {\mathit{R}}^{\mathrm{2}}\end{array}{\mathrm{\nabla }}^{\mathrm{2}}{\mathrm{\Phi }}_{\mathrm{h}}\mathrm{\left(}\begin{array}{c}\mathrm{´}\\ \mathit{R}\end{array}\mathrm{\right)}\mathrm{d}\begin{array}{c}\mathrm{´}\\ \mathit{R}\end{array}\\ & \mathrm{=}& \underset{\mathrm{0}}{\overset{\mathit{R}}{\mathrm{\int }}}\begin{array}{c}\mathrm{´}\\ {\mathit{R}}^{\mathrm{2}}\end{array}\frac{\mathrm{1}}{\begin{array}{c}\mathrm{´}\\ {\mathit{R}}^{\mathrm{2}}\end{array}}\frac{\mathrm{d}}{\mathrm{d}\begin{array}{c}\mathrm{´}\\ \mathit{R}\end{array}}\left(\begin{array}{c}\mathrm{´}\\ {\mathit{R}}^{\mathrm{2}}\end{array}\frac{\mathrm{d}}{\mathrm{d}\begin{array}{c}\mathrm{´}\\ \mathit{R}\end{array}}{\mathrm{\Phi }}_{\mathrm{h}}\mathrm{\left(}\begin{array}{c}\mathrm{´}\\ \mathit{R}\end{array}\mathrm{\right)}\right)\mathrm{d}\begin{array}{c}\mathrm{´}\\ \mathit{R}\end{array}\\ & \mathrm{=}& {\mathit{R}}^{\mathrm{2}}\frac{\mathrm{d}}{\mathrm{d}\mathit{R}}{\mathrm{\Phi }}_{\mathrm{h}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathit{.}\end{array}$(19)Inserting Eq. (18) and directly integrating the result accounting for the boundary condition Φ_h(∞) = 0 yields ${\mathrm{\Phi }}_{\mathrm{h}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\underset{\mathrm{\infty }}{\overset{\mathit{R}}{\mathrm{\int }}}\frac{{\mathit{m}}_{\mathrm{h}}\mathrm{\left(}\begin{array}{c}\mathrm{´}\\ \mathit{R}\end{array}\mathrm{\right)}}{\begin{array}{c}\mathrm{´}\\ {\mathit{R}}^{\mathrm{2}}\end{array}}\mathrm{d}\begin{array}{c}\mathrm{´}\\ \mathit{R}\end{array}\mathrm{=}\begin{array}{c}\frac{{\mathit{M}}_{\mathrm{h}}}{{\mathit{a}}_{\mathrm{h}}}\left(\frac{\mathrm{1}}{\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}}\ln \left(\frac{\mathrm{1}\mathrm{+}{\left(\frac{\mathit{R}}{{\mathit{a}}_{\mathrm{h}}}\right)}^{\mathit{\gamma }\mathrm{-}\mathrm{1}}}{\mathrm{1}\mathrm{+}{\left(\frac{\mathrm{\Lambda }}{{\mathit{a}}_{\mathrm{h}}}\right)}^{\mathit{\gamma }\mathrm{-}\mathrm{1}}}\right)\mathrm{-}\frac{{\left(\frac{\mathrm{\Lambda }}{{\mathit{a}}_{\mathrm{h}}}\right)}^{\mathit{\gamma }\mathrm{-}\mathrm{1}}}{\mathrm{1}\mathrm{+}{\left(\frac{\mathrm{\Lambda }}{{\mathit{a}}_{\mathrm{h}}}\right)}^{\mathit{\gamma }\mathrm{-}\mathrm{1}}}\right) \mathrm{if} \mathit{R}\mathit{<}\mathrm{\Lambda }\mathit{,}\\ \mathrm{-}\frac{{\mathit{M}}_{\mathrm{h}}}{\mathit{R}}\frac{{\left(\frac{\mathrm{\Lambda }}{{\mathit{a}}_{\mathrm{h}}}\right)}^{\mathit{\gamma }}}{\mathrm{1}\mathrm{+}{\left(\frac{\mathrm{\Lambda }}{{\mathit{a}}_{\mathrm{h}}}\right)}^{\mathit{\gamma }\mathrm{-}\mathrm{1}}}\mathit{,}\mathit{otherwise}\mathit{.}\end{array}$For Λ = 100 kpc and γ = 2.02, this is equivalent to the expression given in Allen & Santillán (1991) obtained via integration by parts. The advantage of the above representation is that the dependent variable R appears only once and linearly in the argument of the logarithm (since γ = 2).

The corresponding density ρ_h is $\begin{array}{c}\\ {\mathit{\rho }}_{\mathrm{h}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\frac{{\mathit{M}}_{\mathrm{h}}}{\mathrm{4}\mathit{\pi }{\mathit{a}}_{\mathrm{h}}}\frac{{\left(\frac{\mathit{R}}{{\mathit{a}}_{\mathrm{h}}}\right)}^{\mathit{\gamma }\mathrm{-}\mathrm{1}}\left({\left(\frac{\mathit{R}}{{\mathit{a}}_{\mathrm{h}}}\right)}^{\mathit{\gamma }\mathrm{-}\mathrm{1}}\mathrm{+}\mathit{\gamma }\right)}{{\mathit{R}}^{\mathrm{2}}{\left(\mathrm{1}\mathrm{+}{\left(\frac{\mathit{R}}{{\mathit{a}}_{\mathrm{h}}}\right)}^{\mathit{\gamma }\mathrm{-}\mathrm{1}}\right)}^{\mathrm{2}}}& \mathrm{if} \mathit{R}\mathit{<}\mathrm{\Lambda }\mathit{,}\\ \mathrm{0}& \mathit{otherwise}\mathit{.}\end{array}$

3.2. Results

The properties of the best-fitting Model I are visualized in Figs. 1 and 2 and summarized in Table 1. Given the model’s simplicity, the overall agreement with the imposed constraints is very good. In particular, r_⊙, μ_SgrA^∗, ρ_⊙, Σ_1.1, and σ_BW are almost perfectly reproduced. The rotation curve, on the other hand, reveals systematic discrepancies between model and observation originating mainly in the inadequate treatment of the Galactic bar and the spiral arms. Moreover, most of the maser sources still seem to lag behind the Galactic rotation curve. Despite these shortcomings, the unweighted χ² per degree of freedom (d.o.f.) is about 1.5, thereby confirming the good match.

To see how well the individual parameters are confined by the observational constraints, 90%-confidence limits were computed from the weighted ${\mathit{\chi }}_{\mathrm{w}}^{\mathrm{2}}$ statistics after multiplying it with a factor yielding ${\mathit{\chi }}_{\mathrm{w}}^{\mathrm{2}}\mathit{/}\mathrm{d}\mathit{.}\mathrm{o}\mathit{.}\mathrm{f}\mathit{.}\mathrm{=}\mathrm{1}$ at the best fit, i.e., at the minimum ${\mathit{\chi }}_{\mathrm{w}}^{\mathrm{2}}$. The resulting intervals are given in Table 1, too. As illustrated in Fig. 3, there is a strong correlation between M_h and a_h that explains the large uncertainties of these two parameters. This degeneracy is a consequence of the missing constraints on the halo mass distribution. The available observations are only able to constrain the halo mass inside a sphere of radius 14.1 kpc (see Fig. 3), which, however, can be reproduced by numerous different combinations of M_h-a_h pairs. A parametrization of the halo using only one independent variable would thus be sufficient to obtain a good fit to the data. From Fig. 3, it is also clear that the degeneracy in the two halo parameters does not follow curves of constant halo mass inside 50, 100, or 200 kpc, which explains the relatively large uncertainties of the total masses inside these radii as listed in Table 1.

Przybilla et al. (2010) were the first to exploit the star J1539+0239 as a probe for the mass of the Galaxy and its main contributor, the dark matter halo, and they derived ${\mathit{m}}_{\mathrm{h}}\mathrm{\ge }\mathrm{1.}{\mathrm{7}}_{-1.1}^{\mathrm{+}\mathrm{2.3}}\mathrm{\times }{\mathrm{10}}^{\mathrm{12}}\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot }}$ based on the same model potential as discussed here. Their result, however, was derived solely from the condition that J1539+0239 is bound to the Milky Way while neglecting any other observational constraint that keeps the mass low. In contrast, the mass estimates presented in Table 1 are fully consistent with all the observational data mentioned before.

The best-fit parameters yield a total mass (with 90%-confidence limits computed analogously to the model parameters’ uncertainties) of ${\mathit{M}}_{\mathrm{total}}\mathrm{=}\mathrm{1.}{\mathrm{9}}_{-0.8}^{\mathrm{+}\mathrm{2.4}}\mathrm{\times }{\mathrm{10}}^{\mathrm{12}}\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot }}\mathrm{\approx }{\mathit{m}}_{\mathrm{h}}$. This result is a modest revision of the Przybilla et al. (2010) value, and it confirms their conclusion that a total Galactic mass of about 1 × 10¹² M_⊙ such as found by Xue et al. (2008) possibly underestimates the true value. From Eq. (18) and the fact that Λ = 200 kpc at the best fit, it is obvious that M_total is affected by the restriction Λ ≤ 200 kpc. Although the latter is justified by cosmological simulations on galaxy formation, the specific choice for the upper boundary can of course be subject to discussion. But since even the mass corresponding to Λ = 450 kpc lies within the 90%-confidence interval of M_total, this concern does not seriously change the result. To remove the Λ dependence, consider the quantity M_{R < 50 kpc}, i.e., the total mass within 50 kpc. From the kinematic analysis of a combined sample of field horizontal-branch stars, globular clusters, and satellite galaxies, Sakamoto et al. (2003) have obtained ${\mathit{M}}_{\mathit{R}\hspace{0.17em}\mathit{<}\hspace{0.17em}\mathrm{50}\hspace{0.17em}\mathrm{kpc}}\mathrm{=}\mathrm{5.}{\mathrm{4}}_{-0.4}^{\mathrm{+}\mathrm{0.1}}\mathrm{\times }{\mathrm{10}}^{\mathrm{11}}\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot }}$ agreeing with ${\mathit{M}}_{\mathit{R}\hspace{0.17em}\mathit{<}\hspace{0.17em}\mathrm{50}\hspace{0.17em}\mathrm{kpc}}\mathrm{=}\mathrm{5.}{\mathrm{1}}_{-0.4}^{\mathrm{+}\mathrm{3.3}}\mathrm{\times }{\mathrm{10}}^{\mathrm{11}}\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot }}$ as found in this study.

It is worth noting that J1539+0239 is not at all gravitationally bound to the Milky Way in the best-fitting model. This is possible because of the large uncertainty in its position and velocity determination resulting in a small contribution (χ_∗ = 0.66) to the total ${\mathit{\chi }}_{\mathrm{w}}^{\mathrm{2}}$. More precise kinematic information on the star would set tighter constraints on the Galactic potential, probably reducing the uncertainties in the mass estimates.

We also model the trajectory of the very distant satellite galaxy Leo I based on the full six-dimensional phase space information (position: $\mathit{\alpha }\mathrm{=}{\mathrm{10}}^{\mathrm{h}}{\mathrm{8}}^{\mathrm{m}}\mathrm{28}\stackrel{\mathrm{s}}{.}\mathrm{68}$, δ =  + 12°18′19.7″; distance: 256.7 ± 13.3 kpc; radial velocity: 282.9 ± 0.5 km s^-1 ; proper motion: μ_αcosδ = −0.1140 ± 0.0295 mas yr^-1, μ_δ = −0.1256 ± 0.0293 mas yr^-1, Sohn et al. 2012). Leo I is found to be formally bound to the Galaxy, showing that J1539+0239 is indeed the dominant constraint on the halo mass.

Table 1 also contains some derived quantities that are not directly fitted, such as v_esc, ⊙ or v₀/r_⊙. The escape velocity at the Sun’s position is about 616 km s^-1, thereby lying very close to the 90%-confidence limit 498 km s^-1 < v_esc, ⊙ < 608 km s^-1 of Smith et al. (2007). The ratio of v₀ over r_⊙ is 28.8 km s^-1 kpc^-1 and thus somewhat less than the values found by Reid et al. (2009), 30.3 ± 0.9 km s^-1 kpc^-1, or McMillan & Binney (2010), 29.9−31.6 km s^-1 kpc^-1, from maser sources alone.

As pointed out before, the halo mass distribution is not well constrained by the observations, so different potential shapes of the halo can yield almost equal matches to the data but differ significantly in other properties, as for instance the total mass. To investigate this behavior, two other representations of the dark matter halo are considered in the following sections.

4. Model II

In Model II, the halo is replaced by the truncated, flat rotation curve model used in Wilkinson & Evans (1999) and Sakamoto et al. (2003).

4.1. Components and characteristics

While bulge and disk components have the same shape as in Model I, the halo potential reads as (Wilkinson & Evans 1999) ${\mathrm{\Phi }}_{\mathrm{h}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{{\mathit{M}}_{\mathrm{h}}}{{\mathit{a}}_{\mathrm{h}}}\ln \left(\frac{\sqrt{{\mathit{R}}^{\mathrm{2}}\mathrm{+}{{\mathit{a}}_{\mathrm{h}}}^{\mathrm{2}}}\mathrm{+}{\mathit{a}}_{\mathrm{h}}}{\mathit{R}}\right)\mathrm{·}$(22)The resulting density ${\mathit{\rho }}_{\mathrm{h}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\frac{{\mathit{M}}_{\mathrm{h}}}{\mathrm{4}\mathit{\pi }}\frac{{{\mathit{a}}_{\mathrm{h}}}^{\mathrm{2}}}{{\mathit{R}}^{\mathrm{2}}{\left({\mathit{R}}^{\mathrm{2}}\mathrm{+}{{\mathit{a}}_{\mathrm{h}}}^{\mathrm{2}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}}$(23)is cusped like R^-2 for R ≪ a_h and falls off like R^-5 for R ≫ a_h. In this way, the corresponding rotation curve is flat in the inner regions, and the total halo mass m_h = M_h is finite without invoking a cutoff parameter.

4.2. Results

The properties of the best-fitting Model II are visualized in Figs. 2 and 4 and summarized in Table 2. While the rotation curve reveals the same systematic shortcomings as found in Model I, the match of the remaining constraints is of comparable quality. The unweighted χ² per d.o.f. is 1.7.

Again, there is a tight correlation between the halo parameters M_h and a_h (see Fig. 5). To avoid unphysically large halos, the restriction a_h ≤ 200 kpc has to be imposed based on exactly the same reasoning as the condition Λ ≤ 200 kpc in Model I. Owing to the correlation, this step also sets an upper limit to the parameter M_h, hence to the total halo mass m_h. Unfortunately, the effects on the latter and thus on the total Galactic mass are more severe in this case since its 90%-confidence interval, ${\mathit{M}}_{\mathrm{total}}\mathrm{=}\mathrm{1.}{\mathrm{7}}_{-0.5}^{\mathrm{+}\mathrm{0.2}}\mathrm{\times }{\mathrm{10}}^{\mathrm{12}}\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot }}$, covers a_h values only up to 218 kpc. The mass inside 50 kpc, on the other hand, is almost independent of the chosen upper boundary for a_h. According to Sakamoto et al. (2003), this quantity is much more robust than, e.g., the total mass of the Milky Way. The reason for this is the specific form of the applied potential. It is such that loci of equal masses inside spheres of radii ranging from ~0 to ~50 kpc lie in a very narrow band in the M_h-a_h plane (see Fig. 5). Although the correlation follows contours of equal mass within the central ~20 kpc, it therefore nearly follows contours of equal mass within the central 50 kpc, resulting in small uncertainties for the latter. The value derived here, M_{R < 50 kpc} = 4.6 ± 0.3 × 10¹¹ M_⊙, is slightly lower than in Sakamoto et al. (2003), ${\mathit{M}}_{\mathit{R}\hspace{0.17em}\mathit{<}\hspace{0.17em}\mathrm{50}\hspace{0.17em}\mathrm{kpc}}\mathrm{=}\mathrm{5.}{\mathrm{4}}_{-0.4}^{\mathrm{+}\mathrm{0.1}}\mathrm{\times }{\mathrm{10}}^{\mathrm{11}}\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot }}$.

In general, the gravitational potential of Model II is shallower than in Model I (see Fig. 2) implying systematically lower masses, a lower local escape velocity of about 576 km s^-1 in agreement with Smith et al. (2007), and an unbound orbit for J1539+0239 (χ_∗ = 0.80). Nevertheless, Leo I is still on a bound orbit.

5. Model III

The third halo potential considered here is based on the universal density profile of dark matter halos suggested by Navarro et al. (1997) from cosmological simulations.

5.1. Components and characteristics

Using ${\mathit{\rho }}_{\mathrm{h}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\frac{{\mathit{M}}_{\mathrm{h}}}{\mathrm{4}\mathit{\pi }}\frac{\mathrm{1}}{{\mathrm{(}{\mathit{a}}_{\mathrm{h}}\mathrm{+}{\mathit{R}}^{\mathrm{)}}}^{\mathrm{2}}\mathit{R}}$(24)(Navarro et al. 1997) together with Eq. (2), one obtains a halo potential of ${\mathrm{\Phi }}_{\mathrm{h}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{{\mathit{M}}_{\mathrm{h}}}{\mathit{R}}\ln \left(\mathrm{1}\mathrm{+}\frac{\mathit{R}}{{\mathit{a}}_{\mathrm{h}}}\right)\mathrm{·}$(25)The weighting factor M_h is thereby equivalent to the mass inside a sphere of radius ~5.3 times the scale length a_h. The combination of this halo potential with a Miyamoto & Nagai (1975) disk and bulge component (see Sect. 3.1.1) is hereafter denoted Model III. Owing to the use of Eq. (24) for R → ∞, the total halo mass is formally logarithmically divergent.

5.2. Results

The properties of the best-fitting Model III are visualized in Figs. 2 and 6 and summarized in Table 3. The quality of the fit is similar to Models I and II and yields an unweighted χ² per d.o.f. of 1.7.

In contrast to the two previous models, the rotation curve is rising outside the solar circle, reaching its maximum at about 82 kpc, and falling off beyond this point. This is a consequence of the search for a compromise between the opposing constraints of making J1539+0239 bound, i.e., giving the Milky Way a very high total mass, and, at the same time, limiting the mass in the central ~19 kpc (see Fig. 7) to match the remaining observations. Unlike Model III, Models I and II were designed to give flat rotation curves, avoiding any analogical behavior. As a result, the masses of Model III systematically exceed those of the others, e.g., by more than a factor of 1.5 in terms of M_{R < 50 kpc}. Moreover, because of using Eq. (24) for R → ∞, the total mass M_total is logarithmically infinite. Following Navarro et al. (1997), an alternative measure of the halo mass is given by M₂₀₀, which is the halo mass inside a sphere of radius r₂₀₀ defined to have a mean interior density of 200 times the critical value for closure ρ_crit = 3H²/8πG. Assuming a Hubble constant H of 73 km s^-1 Mpc^-1, values of ${\mathit{r}}_{\mathrm{200}}\mathrm{=}\mathrm{31}{\mathrm{9}}_{-58}^{\mathrm{+}\mathrm{61}}\hspace{0.17em}\mathrm{kpc}$ and ${\mathit{M}}_{\mathrm{200}}\mathrm{=}\mathrm{4.}{\mathrm{0}}_{-1.8}^{\mathrm{+}\mathrm{1.2}}\mathrm{\times }{\mathrm{10}}^{\mathrm{12}}{\mathit{M}}_{\mathrm{\odot }}$ are derived from Model III agreeing with r₂₀₀ = 301 kpc and M₂₀₀ = 3.4 × 10¹²  M_⊙ as used by Abadi et al. (2009). The respective local escape velocity of 812 km s^-1 significantly exceeds the Smith et al. (2007) value but is comparable to the one in Abadi et al. (2009). Although this is high enough to keep Leo I formally bound, it does not hold J1539+0239 within 200 kpc (χ_∗ = 0.20).

Kenyon et al. (2008) considered a Galactic potential that only differs from Model III in the bulge component: the expression R + b_b is used instead of $\sqrt{{\mathit{R}}^{\mathrm{2}}\mathrm{+}{{\mathit{b}}_{\mathrm{b}}}^{\mathrm{2}}}$. Putting special emphasis on the very central region of the Milky Way, the parameters of their model were tweaked to provide a reasonable match to observations from R = 5 pc to R = 100 kpc. Compared to Model III, their resulting model is less massive (M₂₀₀ = 1.6 × 10¹²  M_⊙), which is also reflected by lower values for the local escape velocity (v_esc, ⊙ = 635 km s^-1), the local mass and surface density (ρ_⊙ = 0.046 M_⊙ pc^-3, Σ_1.1 = 46 M_⊙ pc^-2), and the circular motion of the LSR (v₀ = 220 km s^-1).

6. The hypervelocity star HE 0437-5439 revisited

Given the smoothness of the presented potentials, a simple fourth-order Runge-Kutta method with adaptive stepsize control is sufficient to reliably and efficiently solve the equations of motion (as given in Appendix A) numerically. Various trial trajectories, including those of the LSR, were successfully calculated in order to test the self-written ISIS-implementation of the adopted integrator. The hypervelocity star HE 0437-5439 (Edelmann et al. 2005) is analyzed as a first and very interesting application because the spatial origin of this star is still under debate.

When hypervelocity stars were first discovered in 2005 (Brown et al. 2005), dynamical ejection from the supermassive black hole at the GC (Hills 1988) was supposed to be their only origin because these stars move so fast that they are unbound to the Galaxy. This scenario was challenged by the discovery of the third hypervelocity star HE 0437-5439 (Edelmann et al. 2005), a massive B star, because the travel time from the GC to its present position would exceed its lifetime by a factor of 4. Edelmann et al. (2005) notice that the star is much closer to the Large Magellanic Cloud (LMC) than to the Galaxy and might therefore originate in the LMC. This was corroborated by a sophisticated differential abundance analysis of high-quality, high-resolution spectra by Przybilla et al. (2008), who shows that the abundance pattern is inconsistent with that of GC stars but consistent with that of LMC B stars to within error limits. A kinematic investigation for the place of birth of HE 0437-5439 would require precise measurements of the star’s, as well as LMC’s, proper motion. While several studies have been carried out to measure the proper motion of the LMC, no measurements have been available for HE 0437-5439 until recently. Brown et al. (2010) used two epochs of images taken with the Advanced Camera for Surveys (ACS) onboard HST to determine the proper motion of HE 0437-5439 and argue that their measurement rules out a place of birth in the LMC at the 3σ level but is consistent with an origin from the GC. To remedy the time-of-flight versus lifetime problem, HE 0437-5439 needs to be a blue straggler (Edelmann et al. 2005). This would imply that the progenitor was a triple system from which a binary was ejected by the supermassive black hole. Later on its trajectory, the binary system merged to form the blue straggler (Perets 2009). However, the result of Brown et al. (2010) is strongly affected by the treatment of LMC motions, as well as the charge-transfer inefficiency (CTI) correction applied to the data. To demonstrate both issues, the kinematical analysis of HE 0437-5439 is revisited here considering three different cases.

• 1).

  To start with, trajectories for 10⁴ Monte-Carlo realizations, which simultaneously account for statistical and systematic uncertainties in current positions and velocities via Gaussian distribution, are computed for HE 0437-5439 and LMC independently using the input parameters listed in Table 4. The distributions of pericenter distance d_p, its corresponding point in time T_p, and relative velocity v_p derived from the resulting sample of 10⁸ combinations of orbit pairs are insensitive to the choice of the mass model (see run #1a in Table 5), and they confirm the results of Brown et al. (2010) (d_p = 13 kpc, T_p = −23 Myr) based on the Kenyon et al. (2008) Galactic potential and identical input values. Since d_p exceeds the radii of LMC’s central region ${\mathit{R}}_{\mathrm{cen}}^{\mathrm{LMC}}\mathrm{=}\mathrm{3}\hspace{0.17em}\mathrm{kpc}$, as well as the outermost regions ${\mathit{R}}_{\mathrm{out}}^{\mathrm{LMC}}\mathrm{=}\mathrm{10}\hspace{0.17em}\mathrm{kpc}$ (Brown et al. 2010), an origin in the LMC seems unlikely in this context.

  Fig. 7 Visualization of the correlation between Mh and ah in Model III. The meaning of the cross, curves, and shaded regions is the same as in Fig. 3 except that the red dashed line defines loci of constant halo mass inside a sphere of radius 18.8 kpc.

  Table 4

  Kinematic input parameters.

  Fig. 8 Histograms showing the distribution of distances dp, times Tp, and relative velocities vp at periastron of HE 0437-5439 with respect to the LMC (upper panel, run #3a) and the GC (lower panel, run #3b) for Model III. Mean values and standard deviations are given in Table 5. The gray-shaded areas mark regions with ${\mathit{d}}_{\mathrm{p}}\mathrm{\le }{\mathit{R}}_{\mathrm{cen}}^{\mathrm{LMC}}\mathrm{=}\mathrm{3}\hspace{0.17em}\mathrm{kpc}$, ${\mathit{d}}_{\mathrm{p}}\mathrm{\le }{\mathit{R}}_{\mathrm{out}}^{\mathrm{LMC}}\mathrm{=}\mathrm{10}\hspace{0.17em}\mathrm{kpc}$, and 15 Myr ≤ |Tp| ≤ 21 Myr. The last is the lifetime of HE 0437-5439 assuming a single-star nature (Przybilla et al. 2008).

• 2).

  Next, the LMC proper motions of Table 4 are replaced by (μ_αcosδ,μ_δ) = (+1.94 ± 0.29, −0.14 ± 0.36) mas yr^-1 (Kroupa & Bastian 1997) to explore their influence on the outcome. Performing the same Monte-Carlo method as before, all three models give smaller pericenter distances (see run #2 in Table 5) and thus closer encounters of the two objects.

  Table 5

  Results of the kinematic investigation of HE 0437-5439.

• 3).

  Assuming the CTI of ACS to be rising linearly in time, Brown et al. (2010) applied 55% of the epoch-2 CTI correction to their epoch-1 images. According to Massey (2010), however, there was a dramatic increase in the CTI between the two epochs that lead to an overcorrection in the epoch-1 data by Brown et al. implying proper motions of HE 0437-5439 that were larger than stated in Brown et al. (2010). The corresponding effects are roughly estimated here by adding the systematic errors to the mean value and omitting them afterwards, i.e., by using (μ_αcosδ,μ_δ) = ( + 0.86 ± 0.25, + 0.57 ± 0.21) mas yr^-1 as input for the Monte-Carlo calculation. The resulting distributions (see run #3a in Table 5) are still almost model independent and visualized for Model III in the upper panel of Fig. 8. In all three models, about 10% of all orbit pairs yield pericenter passages within the central region, i.e., ${\mathit{d}}_{\mathrm{p}}\mathrm{\le }{\mathit{R}}_{\mathrm{cen}}^{\mathrm{LMC}}$, or 76% in the outermost regions of the LMC (${\mathit{d}}_{\mathrm{p}}\mathrm{\le }{\mathit{R}}_{\mathrm{out}}^{\mathrm{LMC}}$). Out of these, 2% (21%) have shorter flight times than the star’s age of 18 ± 3 Myr (Przybilla et al. 2008). Thus, 0.2% (16%) of the trajectories are consistent with an origin in the inner (outer) LMC without invoking additional requirements, such as a blue straggler nature. The decrease in pericenter distances due to CTI effects (Table 5: run #2 versus run #3a) is stronger than when they are due to a change in LMC proper motions (Table 5: run #1a versus run #2).

For comparison, we now consider the GC as the place of origin. Using the input values of Table 5, HE 0437-5439 passed 95 ± 14 Myr ago within d_p = 14.8 ± 10.1 kpc of the GC (see Model III of run #1b in Table 5). While the travel time is in good agreement with Brown et al. (2010) (T_p = −98 Myr), it is difficult to compare d_p since it is not given in Brown et al. (2010). Nevertheless, the upper lefthand panel of Fig. 3 in Brown et al. (2010) suggests that substantially more than 68% of all trajectories have Galactic plane-crossing locations below 15 kpc. Because the Galactic plane-crossing location is an upper limit for d_p, this fraction increases further when we consider the pericenter distance. However, this contradicts our result that less than 60% of all orbits pass the GC within 15 kpc. As seen from run #3b in Table 5, the situation does not change significantly when accounting for the CTI overcorrection. The respective distributions from the 10⁴ trajectories of HE 0437-5439 are shown in the lower panel of Fig. 8, whereby only ~2% (27%) of the trajectories have pericenter passages within 3 kpc (10 kpc). All of them exceed the stellar lifetime and therefore require HE 0437-5439 to be a blue straggler.

As a consequence, ruling out an origin in the LMC in favor of the GC is disputable and the question of HE 0437-5439’s birthplace remains open.

7. Conclusions

The motions of stars provide important information about the mass distribution in various components of the Galaxy. In particular, they trace the dark matter. Oort’s (1932) measurements of stellar motions were the first to hint at the presence of dark matter in the Milky Way. Soon, the Gaia mission will provide velocity information of a huge number of stars and satellite galaxies with unprecedented precision. Analytical mass models of the Galaxy are utilized to calculate the orbits of stars. Allen & Santillán (1991) derived such a model by making use of observational constraints such as the Galactic rotation curve, the distance to the GC, and the local mass density. During the past twenty years, observational data have greatly improved. Therefore, it was time to revisit the Galactic gravitational potential of Allen & Santillán (1991) and update its parameters by making use of the latest observations. The halo mass is hereby constrained by the most extreme halo star discovered so far (Przybilla et al. 2008). For comparison, two other widely used halo mass distributions – the truncated, flat rotation curve model (Wilkinson & Evans 1999) and a model derived from numerical cosmological simulations (Navarro et al. 1997) – are fitted as well.

All three Milky Way mass models are analytical, simple, and equally capable of reproducing their imposed restrictions. Major discrepancies only become apparent at Galactocentric distances greater than ~15 kpc, which is the region where observational constraints are rare and the halo component dominates. In particular, the depth of the gravitational potential is very sensitive to the form of the dark matter halo seriously affecting for instance predictions of the local escape velocity. The significantly different shapes of the halos allow systematic uncertainties in kinematic investigations to be estimated by comparing the results of orbit computations performed in the three models separately.

Finally, the enigmatic hypervelocity star HE 0437-5439 is re-investigated by inspecting its trajectory in the context of the updated potentials. There are several pros and cons of whether the object originated in the LMC as proposed by the discoverers of the star (Edelmann et al. 2005) or from the GC, which is the suspected place of origin of all hypervelocity stars (Brown et al. 2005). The problem with the latter is that the travel time to the GC exceeds the stellar lifetime and that the chemical composition of HE 0437-5439 differs from what is known about stars in the GC. A birthplace in the LMC, on the other hand, would be consistent with both. However, based on their own measurement of the star’s proper motion, Brown et al. (2010) claim that the LMC – in contrast to the GC – is ruled out as the place of origin. Therefore, we calculated trajectories for the star by making use of the three mass models under study. Using the same input data as Brown et al. (2010), we could confirm their results regarding the LMC, while those for the GC show discrepancies that are independent of the applied mass model. Our trajectories pass the GC within a considerably larger distance, showing that an origin in the GC is much less likely than suggested by Brown et al. (2010). Moreover, we investigated the influence of LMC proper motions and inspected systematic errors in the star’s proper motion measurements stemming from the CTI overcorrection in the ACS images. The latter have not been considered by Brown et al. (2010) but turned out to be crucial. Our calculations show that an origin of HE 0437-5439 in the LMC is at least as likely or unlikely as an origin in the GC. Presently available proper motions are therefore inconclusive, and improved measurements are mandatory for settling this issue.

Acknowledgments

A.I. acknowledges support from a research scholarship by the Elite Network of Bavaria. B.W. and E.T. acknowledge support for a research internship by the German Academic Exchange Service. L.S.’s research internship was made possible by support from the MIT-Germany Program. This research has made use of NASA’s Astrophysics Data System Bibliographic Services. We thank John E. Davis for the development of the slxfig module used to prepare the figures in this paper, Yoshiaki Sofue for providing the H i terminal velocity data in digitized form, and Moritz Böck for helpful discussions. We are very grateful to Norbert Przybilla and Ulrich Heber for valuable discussions in the course of this work and for their suggestions, comments, and contributions to the manuscript.

References

Appendix A: Equations of motion

The motion of test particles in an axisymmetric gravitational potential Φ(r,z) can be conveniently derived from the Lagrangian ℒ of the system: $\mathrm{ℒ}\left(\mathit{r,z,}\mathit{ṙ}\mathit{,}\mathit{\phi ̇}\mathit{,}\mathit{ż}\right)\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{ṙ}}^{\mathrm{2}}\mathrm{+}{\left(\mathit{r}\mathit{\phi ̇}\right)}^{\mathrm{2}}\mathrm{+}{\mathit{ż}}^{\mathrm{2}}\mathrm{)}\mathrm{-}\mathrm{\Phi }\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\hspace{0.17em}\mathit{.}$(A.1)Introducing canonical momenta p_r = ∂ℒ/∂ṙ = ṙ, , and p_z = ∂ℒ/∂ż = ż, the corresponding Lagrange’s equations – written as a system of first-order differential equations – read as $\begin{array}{ccc}& & \mathit{ṙ}\mathrm{=}{\mathit{p}}_{\mathit{r}}\mathit{,}\\ & & \mathit{\phi ̇}\mathrm{=}\frac{{\mathit{p}}_{\mathit{\phi }}}{{\mathit{r}}^{\mathrm{2}}}\mathit{,}\\ & & \mathit{ż}\mathrm{=}{\mathit{p}}_{\mathit{z}}\mathit{,}\\ & & \mathit{ṗ}\mathit{r}\mathrm{=}\frac{{{\mathit{p}}_{\mathit{\phi }}}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{3}}}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{b}}\mathit{r}}{{\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}\mathrm{+}{{\mathit{b}}_{\mathrm{b}}}^{\mathrm{2}}\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{d}}\mathit{r}}{{\left({\mathit{r}}^{\mathrm{2}}\mathrm{+}{\left({\mathit{a}}_{\mathrm{d}}\mathrm{+}\sqrt{{\mathit{z}}^{\mathrm{2}}\mathrm{+}{{\mathit{b}}_{\mathrm{d}}}^{\mathrm{2}}}\right)}^{\mathrm{2}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}}\\ & & \mathrm{-}\frac{\mathit{\partial }{\mathrm{\Phi }}_{\mathrm{h}}\mathrm{\left(}\mathit{R}\mathrm{\right)}}{\mathit{\partial R}}\frac{\mathit{r}}{\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}}}\hspace{0.17em}\mathit{,}\\ & & \mathit{ṗ}\mathit{\phi }\mathrm{=}\mathrm{0}\mathit{,}\\ & & \mathit{ṗ}\mathit{z}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{b}}\mathit{z}}{{\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}\mathrm{+}{{\mathit{b}}_{\mathrm{b}}}^{\mathrm{2}}\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{d}}\mathit{z}\left(\frac{{\mathit{a}}_{\mathrm{d}}}{\sqrt{{\mathit{z}}^{\mathrm{2}}\mathrm{+}{{\mathit{b}}_{\mathrm{d}}}^{\mathrm{2}}}}\mathrm{+}\mathrm{1}\right)}{{\left({\mathit{r}}^{\mathrm{2}}\mathrm{+}{\left({\mathit{a}}_{\mathrm{d}}\mathrm{+}\sqrt{{\mathit{z}}^{\mathrm{2}}\mathrm{+}{{\mathit{b}}_{\mathrm{d}}}^{\mathrm{2}}}\right)}^{\mathrm{2}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}}\\ & & \mathrm{-}\frac{\mathit{\partial }{\mathrm{\Phi }}_{\mathrm{h}}\mathrm{\left(}\mathit{R}\mathrm{\right)}}{\mathit{\partial R}}\frac{\mathit{z}}{\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}}}\mathit{,}\end{array}$whereby $\frac{\mathit{\partial }{\mathrm{\Phi }}_{\mathrm{h}}\mathrm{\left(}\mathit{R}\mathrm{\right)}}{\mathit{\partial R}}\mathrm{=}\begin{array}{c}\\ \frac{\mathit{G}{\mathit{M}}_{\mathrm{h}}}{{\mathit{a}}_{\mathrm{h}}\mathit{R}}\frac{{\left(\frac{\mathit{R}}{{\mathit{a}}_{\mathrm{h}}}\right)}^{\mathit{\gamma }\mathrm{-}\mathrm{1}}}{\mathrm{1}\mathrm{+}{\left(\frac{\mathit{R}}{{\mathit{a}}_{\mathrm{h}}}\right)}^{\mathit{\gamma }\mathrm{-}\mathrm{1}}}& \mathrm{if} \mathit{R}\mathit{<}\mathrm{\Lambda }\mathit{andModel}\mathit{I,}\\ \frac{\mathit{G}{\mathit{M}}_{\mathrm{h}}}{{\mathit{R}}^{\mathrm{2}}}\frac{{\left(\frac{\mathrm{\Lambda }}{{\mathit{a}}_{\mathrm{h}}}\right)}^{\mathit{\gamma }}}{\mathrm{1}\mathrm{+}{\left(\frac{\mathrm{\Lambda }}{{\mathit{a}}_{\mathrm{h}}}\right)}^{\mathit{\gamma }\mathrm{-}\mathrm{1}}}& \mathrm{if} \mathit{R}\mathrm{\ge }\mathrm{\Lambda }\mathit{andModel}\mathit{I,}\\ \frac{\mathit{G}{\mathit{M}}_{\mathrm{h}}}{\sqrt{{\mathit{R}}^{\mathrm{2}}\mathrm{+}{{\mathit{a}}_{\mathrm{h}}}^{\mathrm{2}}}\mathit{R}}& \mathrm{if} \mathit{Model}\mathit{II,}\\ \frac{\mathit{G}{\mathit{M}}_{\mathrm{h}}}{{\mathit{R}}^{\mathrm{2}}}\ln \mathrm{(}\mathrm{1}\mathrm{+}{\frac{\mathit{R}}{{\mathit{a}}_{\mathrm{h}}}}^{\mathrm{)}}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{h}}}{\mathrm{(}{\mathit{a}}_{\mathrm{h}}\mathrm{+}\mathit{R}\mathrm{)}\mathit{R}}& \mathrm{if} \mathit{Model}\mathit{III}\mathit{.}\end{array}$Expressing lengths in kpc, masses in galactic mass units M_gal, and time in Myr, the gravitational constant G has to be replaced by the factor ~1.04598 × 10^-4 thereby converting units from 100 km² s^-2 to kpc² Myr^-2.

All Tables

All Figures