© ESO, 2014

1. Introduction

The existence of a bar in our Galaxy is supported by a variety of studies using data from HI 21cm and CO emission, star counts in the Galactic centre (GC), IR observations from DIRBE (Diffuse InfraRed Background Experiment) on COBE (COsmic Background Explorer) and GLIMPSE (Galactic Legacy Infrared Mid-Plane Survey Extraordinaire) with Spitzer, or microlensing surveys (see Gerhard 2002, for a review). However, previous research has revealed inconsistent results regarding the characteristics of the bar. For example, estimates of its pattern speed range from 40 to 65 km s^-1 kpc^-1 (Gerhard 2011) while the estimates of its orientation with respect to the Sun range from 10 deg (López-Corredoira et al. 2000; Robin et al. 2012) to 45 deg (Hammersley et al. 2000; Benjamin et al. 2005). The presence of a secondary bar in our Galaxy is also currently under debate (Martinez-Valpuesta & Gerhard 2011; Romero-Gómez et al. 2011).

Kalnajs (1991) presented an indirect method to measure the bar properties based on the location of kinematic structures in the solar neighbourhood. He related the velocities of the Hyades and Sirius moving groups to the two types of orbits expected around the outer Lindblad resonance (OLR), and in this way constrained the bar’s pattern speed and its orientation.

Many more substructures in the local velocity distribution were unveiled by the ESA’s astrometric mission Hipparcos (e.g. Dehnen 1998; de Zeeuw et al. 1999; Chereul et al. 1999). Most of these groups were initially thought to be remnants of disrupted clusters (Eggen 1996). However, there is evidence of a large scatter in age and metallicity in some of them (Raboud et al. 1998; Dehnen 1998; Skuljan et al. 1999; Famaey et al. 2005; Bobylev & Bajkova 2007; Antoja et al. 2008). Therefore, it is likely that these substructures formed as a response to the non-axisymmetries of the gravitational potential rather than being groups of stars of a common origin.

Several studies after Kalnajs (1991) have attempted to use these local velocity groups to better constrain the properties of the Galactic bar (e.g. Dehnen 2000, hereafter D00), and also of the spiral structure (e.g. Quillen & Minchev 2005). However, Antoja et al. (2009, 2011) have shown that the groups detected in the solar vicinity can be reproduced by models with different parameters, including bar or/and spiral structure, highlighting that local estimates are subject to degeneracies.

The simulations of Antoja et al. (2011, 2009) as well as those of e.g. Quillen et al. (2011) have shown that the groups’ kinematics change across the disc. Recently, using action-angle modelling McMillan (2013) showed how the local Hyades stream can be due to the effects of different resonances such as the Inner Lindblad Resonance or OLR of a non-axisymmetric pattern in the disc. But these models predict differences in the stream kinematics throughout the disc.

With the advent of data from new surveys such as RAVE (RAdial Velocity Experiment, Steinmetz et al. 2006) the detection of kinematic groups is no longer limited to the solar vicinity. The first example was given by Antoja et al. (2012, hereafter A12) where wavelet transform techniques were used to detect kinematic groups beyond the solar neighbourhood in the RAVE survey. The sampled volume allowed the demonstration that some local groups can be traced at least up to 1 kpc away from the Sun in certain directions and that their velocities change with distance. These discoveries point toward the exciting possibility of using observed velocity distributions in a number of regions of the Galaxy to break degeneracies and eventually constrain the properties of the spiral arms and the bar.

A12 showed that Hercules, a local group of stars moving outwards in the disc and lagging the local standard of rest, has a larger azimuthal velocity inside the solar circle and a smaller one outside. Here we quantify this trend in more detail with the new RAVE DR4, showing that it is consistent with the effects of the bar’s OLR, and we use it to constrain the properties of the Galactic bar.

In Sect. 2 we review the properties of the local Hercules stream and its relation with the effects of the bar’s OLR. Also by extending the modelling work of D00, we derive an analytic expression for the variation of the azimuthal velocity of Hercules as a function of Galactocentric radius for different bar properties. We then use simulations of a barred disc to test this model and the recovery of the simulation’s parameters (Sect. 3). In Sect. 4 we measure the observed azimuthal velocity of Hercules as a function of radius for RAVE stars. We finally compare these measurements with the predictions of the effects of the bar’s OLR and we derive the best fit parameters of the bar (Sect. 5). Section 6 contains a final discussion and conclusions.

2. The Hercules stream

2.1. The local Hercules stream and the OLR

The Hercules moving group (also refereed to as the U-anomaly) was initially identified by Eggen (1958) as a group of 22 stars with velocities similar to the high velocity star ξ Herculis. Also Blaauw (1970) noticed an excess of negative U velocities (with U directed to the Galactic centre) for stars with V ~  − 50 km s^-1, that is a stream of stars with eccentric orbits with a mean outward radial motion. The Hercules stream is also evident in the solar neighbourhood velocity distribution of RAVE stars (Fig. 1). By using photometric data for a sub-sample of Hipparcos stars Raboud et al. (1998) showed that the Hercules metallicity distribution covers the whole range observed in the old disc (− 0.6 dex to + 0.6 dex) and a heterogeneous distribution of ages between 6 to 10 Gyr. Similar conclusions were obtained using colour as a proxy for age (Dehnen 1998), with isochrones for giant stars (Famaey et al. 2005), with ages from Strömgren photometry and spectro-photometric metallicities for F and G dwarfs (Helmi et al. 2006; Bobylev & Bajkova 2007; Antoja et al. 2008), and also with high-resolution abundances for F and G dwarf stars (Bensby et al. 2007). Although there is some discrepancy regarding the lower age limit of the group, which ranges from 1 to 6  Gyr, it is now clear that this group does not originate in a single population or cluster.

The first dynamical models for Hercules were presented in D00 and Fux (2001), and were based on the effects of the bar on the local velocity distribution. D00 proposed that Hercules consists of stars that have been scattered by the OLR. In particular, for certain ranges of pattern speeds and orientations of the bar, a group of unstable orbits (${\mathit{x}}_{\mathrm{1}}^{\mathrm{\ast }}\mathrm{\left(}\mathrm{2}\mathrm{\right)}$ orbits) divides the velocity distribution into two main groups (bi-modality) separated by a valley (Fig. 1, see also Sect. 3). One group is approximately centred on the U–V velocity plane (MAIN mode) and the other one has a slower rotation, mean outward radial motion and is associated to the Hercules moving group (OLR mode).

D00 simulated the velocity distribution at the solar position of a barred potential using two-dimensional (in-plane) test particle orbital integrations with the backwards integration technique. The bar model was a quadrupole potential (his Eq. (3)) rotating with speed Ω_b and orientation angle φ_b with respect to the line Sun-GC. He used a simple underlying potential (his Eq. (2b)) with a power-law rotation curve of the form: ${\mathit{v}}_{\mathrm{c}}\mathrm{=}{\mathit{v}}_{\mathrm{0}}\hspace{0.17em}\mathrm{(}\mathit{R}\mathit{/}{\mathit{R}}_{\mathrm{0}}{\mathrm{)}}^{\mathit{\beta }}\mathit{,}$(1)where R₀ denotes the Sun’s distance from the GC and v₀ the local circular speed. Figure 2 shows the rotation curves of this model for different values of β.

By considering only these axisymmetric power-law potentials (thus neglecting the effect of the quadrupole bar), whose orbital frequencies can be derived analytically and by dismissing terms of O(${\mathit{v}}^{\mathrm{3}}\mathit{/}{\mathit{v}}_{\mathrm{0}}^{\mathrm{3}}$), D00 showed that the stars on unstable resonant orbits exactly on the OLR of the rotating frame form a parabola in velocity space (the valley) described by: $\mathit{V}\mathrm{+}\frac{{\mathit{U}}^{\mathrm{2}}}{\mathrm{2}{\mathit{v}}_{\mathrm{0}}}\cong {\mathit{v}}_{\mathrm{0}}\begin{array}{c}\mathrm{˜}\\ {\mathit{V}}_{\mathrm{OLR}}\end{array}\mathrm{\equiv }\hspace{0.17em}\frac{\mathrm{1}\mathrm{+}\mathit{\beta }}{\mathrm{1}\mathrm{-}\mathit{\beta }}\left[\mathrm{1}\mathrm{-}\frac{{\mathrm{\Omega }}_{\mathrm{b}}\mathit{/}{\mathrm{\Omega }}_{\mathrm{0}}}{\mathrm{1}\mathrm{+}\sqrt{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\beta }\mathrm{)}\mathit{/}\mathrm{2}}}\right]\mathit{,}$(2)where U and V are the velocities with respect to the local standard of rest and Ω₀ is the local circular frequency. This parabola has a maximum at $\mathit{V}\mathrm{=}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{V}}_{\mathrm{OLR}}\end{array}$ occurring at U = 0 (a saddle point). However, in his simulations including the quadrupole bar the saddle point between the two modes appears shifted in U and also in V with respect to the analytic estimate given by Eq. (2). Then he found that the V-velocity of the saddle-point V_OLR could be fitted by: ${\mathit{V}}_{\mathrm{OLR}}\mathrm{\approx }\mathit{a}\hspace{0.17em}\begin{array}{c}\mathrm{˜}\\ {\mathit{V}}_{\mathrm{OLR}}\end{array}\mathrm{-}\mathrm{(}\mathit{b}\mathrm{+}\mathit{c}\hspace{0.17em}\mathit{\beta }\mathrm{)}\hspace{0.17em}{\mathit{v}}_{\mathrm{0}}\mathit{,}$(3)where the values of a,b and c are reproduced in Table 1 and depend on the bar’s orientation φ_b. Using Eqs. (2) and (3) we find: ${\mathit{V}}_{\mathrm{OLR}}\mathrm{\approx }\mathit{a}\hspace{0.17em}{\mathit{v}}_{\mathrm{0}}\hspace{0.17em}\frac{\mathrm{1}\mathrm{+}\mathit{\beta }}{\mathrm{1}\mathrm{-}\mathit{\beta }}\left[\mathrm{1}\mathrm{-}\frac{{\mathrm{\Omega }}_{\mathrm{b}}\mathit{/}{\mathrm{\Omega }}_{\mathrm{0}}}{\mathrm{1}\mathrm{+}\sqrt{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\beta }\mathrm{)}\mathit{/}\mathrm{2}}}\right]\mathrm{-}\mathrm{(}\mathit{b}\mathrm{+}\mathit{c}\hspace{0.17em}\mathit{\beta }\mathrm{)}\hspace{0.17em}{\mathit{v}}_{\mathrm{0}}\mathit{,}$(4)which relates the position of the Hercules saddle point V_OLR to the pattern speed of the bar, its orientation (through parameters a, b and c) and the slope and normalisation of the rotation curve. Using the local observed velocity distribution of Hipparcos stars D00 found this saddle point to be at V_OLR = ( − 31 ± 3) km s^-1.

2.2. Analytic model for the Hercules stream across the Galaxy

Our purpose is now to generalise Eq. (4) to different Galactocentric radii R (i.e. not necessarily the solar neighbourhood). For this, we replace the quantities describing the solar neighbourhood by their respective functional forms, that is v₀ by v₀(R/R₀)^β and Ω₀ by v₀(R/R₀)^β/R (Eq. (1)). In cylindrical Galactocentric coordinates (v_φ = V + v₀), Eq. (4) becomes: $\begin{array}{ccc}{\mathit{v}}_{\mathit{\phi ,}\mathrm{OLR}}\mathrm{\left(}\mathit{R}\mathrm{\right)}& \mathrm{\approx }& \mathit{a}\hspace{0.17em}{\mathit{v}}_{\mathrm{0}}\hspace{0.17em}\mathrm{(}\mathit{R}\mathit{/}{\mathit{R}}_{\mathrm{0}}{\mathrm{)}}^{\mathit{\beta }}\hspace{0.17em}\frac{\mathrm{1}\mathrm{+}\mathit{\beta }}{\mathrm{1}\mathrm{-}\mathit{\beta }}\left[\mathrm{1}\mathrm{-}\frac{{\mathrm{\Omega }}_{\mathrm{b}}\mathit{R}}{{\mathit{v}}_{\mathrm{0}}\mathrm{(}\mathit{R}\mathit{/}{\mathit{R}}_{\mathrm{0}}{\mathrm{)}}^{\mathit{\beta }}}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\sqrt{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\beta }\mathrm{)}\mathit{/}\mathrm{2}}}\right]\\ & & \mathrm{-}\mathrm{(}\mathit{b}\mathrm{+}\mathit{c}\hspace{0.17em}\mathit{\beta }\mathrm{-}\mathrm{1}\mathrm{)}\hspace{0.17em}{\mathit{v}}_{\mathrm{0}}\hspace{0.17em}\mathrm{(}\mathit{R}\mathit{/}{\mathit{R}}_{\mathrm{0}}{\mathrm{)}}^{\mathit{\beta }}\mathit{.}\end{array}$(5)In this way, the position of the saddle point between the Hercules and the MAIN mode is a function v_φ,OLR = f(Ω_b,v₀,β,φ_b,R). Now φ_b is the angle between the considered region (not necessarily the solar neighbourhood) and the bar. As explained in D00, a higher (lower) force of the bar creates more (less) pronounced features in the velocity plane but does not influence significantly v_φ,OLR, and, therefore, it does not appear explicitly in Eqs. (2) and (3).

Figure 3 shows v_φ,OLR as a function of Galactocentric radius R for different bar properties and different rotation curves. For this plot we have set a Galactocentric radius of the Sun of R₀ = 8.05 kpc and a circular velocity at the Sun of v₀ = 238 km s^-1 following recent results by Honma et al. (2012) based on VLBI astrometry of Galactic maser sources. We see that the position of the saddle point decreases with R. This is in agreement with what was reported by A12 for Hercules in the RAVE data. According to the model, the position of the saddle point decreases linearly for β = 0 as: ${\mathit{v}}_{\mathit{\phi ,}\mathrm{OLR}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{\approx }\mathrm{(}\mathit{a}\mathrm{-}\mathit{b}\mathrm{+}\mathrm{1}\mathrm{)}\hspace{0.17em}{\mathit{v}}_{\mathrm{0}}\mathrm{-}\frac{\mathit{a}{\mathrm{\Omega }}_{\mathrm{b}}}{\mathrm{1}\mathrm{+}\sqrt{\mathrm{1}\mathit{/}\mathrm{2}}}\mathit{R}\mathit{.}$(6)We see also that for a given angle (shown by lines with the same colour), a higher pattern speed Ω_b (lower group of curves) produces a smaller v_φ,OLR at a given radius compared to lower Ω_b (upper group). For a fixed pattern speed, larger bar orientations (red curves) lead to a relation between saddle point v_φ,OLR and R that has a different slope compared to smaller angles (blue curves). We also notice a slight dependence on the slope of the rotation curve β: depending on the pattern speed, decreasing (dotted) or increasing (dashed) rotation curves give smaller or larger v_φ,OLR compared to flat (solid) rotations curves. This dependence is due to the resonances moving closer or farther away from the Sun as the rotation curve is changed.

3. Validating the analytical model with simulations

In the previous section we presented an analytical model for the Hercules stream which relied on specific assumptions on the potential, namely the shape of the rotation curve, and on the derived frequencies of the orbits. To test the validity of this model, especially Eq. (5), in this section we use an independent simulation that has been run using a different potential.

3.1. Test particle simulations

We use a simulation similar to that of Antoja et al. (2009), in which the bar’s pattern speed is Ω_b = 47.5 km s^-1 kpc^-1. Our simulation uses the same quadrupole bar as in D00 and is also two-dimensional. However, our axisymmetric potential is given by Allen & Santillan (1991, A91), and composed of a bulge and a flattened disc modelled as Miyamoto-Nagai potentials, and a spherical halo. This axisymmetric model uses a value¹ of R₀ = 8.5 kpc for the solar radius and a local circular speed of v₀ = 220 km s^-1. The resulting circular velocity of the model is shown in Fig. 2 (dotted-dashed line). This curve is fairly different than the power-law models of Eq. (1) by D00, and presents sections with different slopes and normalisations. The inner peak is due to the presence of the bulge. This different underlying model does not have the same orbital frequency dependencies used to derive Eq. (2) and, therefore, allows us to test if the approximations are nonetheless valid for other potentials.

Another important difference of our simulation as compared with D00 is that we use different initial conditions and a different integration scheme. Instead of the backwards integration, we start with 12 × 10⁶ test particles with an initial distribution function as discussed in Hernquist (1993). The density follows an exponential disc and the velocity distribution is adopted as a Gaussian with a radial velocity dispersion decaying exponentially with radius, with value of ~50 km s^-1 at the solar radius. The azimuthal velocity dispersion is related with the radial one through the epicyclic approximation and the asymmetric drift is also taken into account. The initial conditions generated in this way are not fully consistent with the potential and we expect these to change in time until reaching stationarity. To avoid these transient effects we first let our initial conditions evolve in the axisymmetric potential for several   Gyr (see Monari et al. 2014, for a discussion). Afterwards, we introduce the bar abruptly in the potential² and the final distribution is obtained through forward integration of the orbits for 0.4 Gyr (equivalent to ~3 bar’s rotation). We consider the particles in a given volume to study the velocity distribution of a particular position in the disc. This is in contrast to D00, whose results correspond to a single position in configuration space.

Figure 4 shows a sketch of the face-on view of the simulation, with the Sun at X = −8.5 kpc and Y = 0, the bar oriented with φ_b = 20 deg with respect to the Sun, and the Galaxy rotating clockwise. From this simulation we have selected the particles located in 4 different bands with orientations of 20 deg (the assumed Sun’s position), 40 deg, 60 deg and 80 deg with respect to the bar. The width of these bands is Δφ = 4 deg and we add the particles in the symmetric bands at respective angles of 200 deg, 220 deg, 240 deg and 260 deg, which are dynamically equivalent. We now aim to explore the velocity distribution (v_R,v_φ) of these bands as a function of Galactocentric radius. We set v_R positive towards the GC as U. We take bins in radius of a width of Δ_R = 600 pc every 600 pc. The number of stars per bin ranges between 2000 stars for the outermost bins to 10 000 stars for the innermost ones (Fig. 6 top).

In the first column of Fig. 5 we show these velocity distributions for the band with 40 deg orientation (blue band in Fig. 4). These panels reveal a bimodal distribution, with the structure at lower v_φ and negative v_R being the modelled Hercules group. For this band, we take bins from R = 7.6 kpc to R = 10.6 kpc as this is the range for which Hercules can be traced. Simple visual inspection shows that, as predicted by our model, the v_φ velocity of Hercules (or equivalently the velocity of the saddle point v_φ,OLR) decreases as a function of R. In the next section we show how we measure the velocity of the saddle point v_φ,OLR.

3.2. Measuring v_φ,OLR in the simulations

We measure the position of the saddle point v_φ,OLR as illustrated in Fig. 5 for the band at 40 deg:

• 1.

  We rotate the (v_φ, v_R) velocities to align the Hercules structure with the horizontal axis, leading to the new coordinates (${\mathit{v}}_{\mathit{R}}^{\mathit{\theta }}$, ${\mathit{v}}_{\mathit{\phi }}^{\mathit{\theta }}$). Visual inspection shows the rotation angle θ to be between 10 deg and 20 deg depending on the band and radius considered. To simplify the method we use the same angle for all bins and we choose a value of 15 deg for reasons specified in step 3.

• 2.

  We estimate the probability density in this velocity space using the Epanechnikov adaptive kernel density estimator method (Silverman 1986) with an adaptability exponent of 0.1. (Fig. 5 second column).

• 3.

  We integrate over ${\mathit{v}}_{\mathit{R}}^{\mathit{\theta }}$, only inside the range ${\mathit{v}}_{\mathit{R}}^{\mathit{\theta }}\mathrm{=}\mathrm{[}\mathrm{-}\mathrm{130}\mathit{,}\mathrm{-}\mathrm{10}\mathrm{]}\hspace{0.17em}\mathrm{km}\hspace{0.17em}{\mathrm{s}}^{-1}$ (within the black vertical lines in the second column) to avoid contamination from other groups or regions of the velocity plane. The distribution along ${\mathit{v}}_{\mathit{\phi }}^{\mathit{\theta }}$ is shown in the third column. We clearly see the presence of the two peaks separated by a valley. The Hercules peak is indicated with a black vertical line. We see how for small R, Hercules is stronger than the MAIN mode, while as we move outwards in the disc it becomes weaker. Of all rotation angles θ = 10,15,20,25 deg, the angle of θ = 15 deg gives the maximum height of Hercules (black vertical line) for most of the bands and radial bins. This means that for this angle the structures are better aligned with the horizontal axis. In step 6 we estimate the error on the final location of the saddle point derived by assuming this value. Then we locate the position of the minimum ${\mathit{v}}_{\mathit{\phi ,}\mathrm{OLR}}^{\mathit{\theta }}$ (red line).

• 4.

  We estimate the error in ${\mathit{v}}_{\mathit{\phi ,}\mathrm{OLR}}^{\mathit{\theta }}$ by generating 500 bootstrap samples, repeating steps 1 to 3, and computing the standard deviation of the obtained set of ${\mathit{v}}_{\mathit{\phi ,}\mathrm{OLR}}^{\mathit{\theta }}$, which is typically very small (~2 km s^-1, red error bars in Fig. 5 third column). Figure 6 (second panel) shows ${\mathit{v}}_{\mathit{\phi ,}\mathrm{OLR}}^{\mathit{\theta }}$ as a function of R for the 5 radial bins (blue diamonds). These velocities decrease with R. Additionally, we show ${\mathit{v}}_{\mathit{\phi ,}\mathrm{OLR}}^{\mathit{\theta }}$ for the other bands at 20 deg, 60 deg and 80 deg in different colours which depict the same behaviour.

  Fig. 6 Several measurements for the bands at different bar angles of the simulation as a function of radius R. Top panel: number of stars per bin. Second panel: measured ${\mathit{v}}_{\mathit{\phi ,}\mathrm{OLR}}^{\mathit{\theta }}$. Third panel: measured ${\mathit{v}}_{\mathit{R,}\mathrm{OLR}}^{\mathit{\theta }}$. Bottom panel: final determination of vφ,OLR and expected theoretical curves.

• 5.

  To get the ${\mathit{v}}_{\mathit{R}}^{\mathit{\theta }}$ of the saddle point, we derive the distribution along ${\mathit{v}}_{\mathit{R}}^{\mathit{\theta }}$ (fourth column) inside the green rectangles of Fig. 5. These are centred in the valley (${\mathit{v}}_{\mathit{\phi ,}\mathrm{OLR}}^{\mathit{\theta }}$) with a width of 20 km s^-1. We then find the maximum of the curve ${\mathit{v}}_{\mathit{R,}\mathrm{OLR}}^{\mathit{\theta }}$ (green line). We also estimate the error in ${\mathit{v}}_{\mathit{R,}\mathrm{OLR}}^{\mathit{\theta }}$ through the bootstrapping technique and it is typically of 5 − 10 km s^-1 (green error bars). The position of ${\mathit{v}}_{\mathit{R,}\mathrm{OLR}}^{\mathit{\theta }}$ as a function of R is shown in the third panel from the top of Fig. 6 (blue diamonds). In general, this velocity becomes more negative with radius. For larger radii the distribution in ${\mathit{v}}_{\mathit{R}}^{\mathit{\theta }}$ is noisier and has several maxima (bottom panels of Fig. 5). This is because the number of particles for large R decreases and the valley is wider and contains less particles (only 24 particles were inside the green rectangle of the last radial bin). In this case our determination of ${\mathit{v}}_{\mathit{R,}\mathrm{OLR}}^{\mathit{\theta }}$ may not be accurate. For instance, the last bin of the band at 40 deg does not follow the overall trend.

• 6.

  We convert (${\mathit{v}}_{\mathit{R,}\mathrm{OLR}}^{\mathit{\theta }}$, ${\mathit{v}}_{\mathit{\phi ,}\mathrm{OLR}}^{\mathit{\theta }}$) back to (v_R,OLR, v_φ,OLR) by rotating an angle − θ. The position of the saddle point is indicated with red bars in the first and second columns of Fig. 5. The value v_φ,OLR is the observable needed in our modelling. To obtain the errors in v_R,OLR and v_φ,OLR we must consider two contributions. First, the statistical errors which, as explained above, we get from the bootstrapping method (e_bst). Second, the error made by using a fixed value for the rotation angle θ (e_θ). To estimate the latter we repeat steps 1 to 6 using the two extreme angles of θ = 15 ± 5 deg, we compute the maximum difference between the new determinations of v_R,OLR and v_φ,OLR and the ones for θ = 15 deg and assign this difference to the error, which turns out to be ≲ 5 km s^-1. Finally, we add both errors e_bst and e_θ in quadrature.

The measured velocities v_φ,OLR are shown Fig. 6 (bottom) for the four bands. For all bands, the velocity decreases with R. Overlaid on the points are the theoretical curves from Eq. (5) for the input parameters of the simulation (Ω_b = 47.5 km s^-1 kpc^-1, v₀ = 220 km s^-1 and R₀ = 8.5 kpc), along with the four bar orientations of the bands. As the slope of the rotation curve changes with radius for the A91 model (Sect. 3.1), we plot for each band three curves corresponding to β = 0.2,0, − 0.2. We see that our measured v_φ,OLR are consistent with the predictions of Eq. (5), given the errors. We see more discrepancies at large radii where it is more difficult to detect reliably the position of the saddle point for the reasons mentioned above. Note that the discrepancy between the estimated value in the last radial bin for the band at 40 deg and that expected is due to the poor determination of ${\mathit{v}}_{\mathit{R,}\mathrm{OLR}}^{\mathit{\theta }}$ (step 5) and we shall reject this data point in our analysis of Sect. 3.4.

We have also validated the analytical model with other simulations with different values of the pattern speed, which moves the resonances to other positions of the disc, and the bar’s force and obtained similarly satisfactory results. We will now recover the input model parameter from the simulated data.

3.3. Maximisation and parameter space sampling

We compare the determined values of the position of the saddle point v_φ,OLRi with the estimates obtained from Eq. (5) through the chi-square statistic: ${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\sum _{\mathit{i}}(\frac{{{\mathit{v}}_{\mathit{\phi ,}\mathrm{OLR}}}_{\mathit{i}}\mathrm{-}{{{\mathit{v}}_{\mathit{\phi ,}\mathrm{OLR}}}_{\mathit{i}}}^{\mathrm{model}}}{{\mathit{\sigma }}_{\mathit{i}}})2$(7)where the subscript i stands for each data point and σ_i is the error of each point. We assume that the noise associated with the data points can be represented as a Gaussian process and, therefore, we can approximate the likelihood function by $\mathit{prob}\mathrm{\propto }\exp \left(\mathrm{-}\frac{{\mathit{\chi }}^{\mathrm{2}}}{\mathrm{2}}\right)$. In the maximisation of the probability, we consider the pattern speed Ω_b in units of Ω₀, i.e. in practice we fit Ω_b/Ω₀. We take a range of [0,3.4] Ω₀, which corresponds to [0,100] km s^-1 kpc^-1, in steps of 0.0025 or ~0.07 km s^-1 kpc^-1. This range is large enough not to influence the posterior probability density function (pdf). For the slope of the rotation curve β we use the range [ − 0.2,0.2] in steps of 0.01. These limits are the ones considered in D00 for which the fit given in Eq. (3) is valid. The bar’s angle φ_b is explored in the range of [0,80] deg in steps of 0.5 deg. Outside this range, the model of Eq. (4) is not valid as the Hercules structure does not exist or there is just a counterpart at v_R > 0 in the velocity plane (e.g. Fig. 2 in D00). Note that this range is actually larger than the limits considered in D00 to obtain the fit of Eq. (3) (he used up to 50 deg). However, as we showed in Sect. 3.2 (Fig. 6) the extrapolation to larger angles is valid.

3.4. Recovering the parameters of the model

Figure 7 shows the pdf for our toy model in the φ_b-Ω_b plane for the bands at 40 and 60 deg. This pdf has been marginalised over β, since we do not expect to constrain the slope of the rotation curve β to a single value as it varies in our simulations over the distance range considered. Indeed, we find a very flat pdf in the direction of β. In Table 2 we give the maximum of the probability, the mean or expectation of each parameter E(φ_b,Ω_b), together with the standard deviation of the pdf σ_φb and σ_Ωb/Ω0, and the correlation coefficient ρ_φbΩb for each band. In Fig. 7 the 1σ, 2σ and 3σ confidence regions are delimited by the solid, dashed and dotted red lines, respectively. The maximum of the probability is indicated with a orange cross whereas the mean E(φ_b,Ω_b) is shown with a orange asterisk. There is a high correlation between orientation and pattern speed (also noticed in D00) with a correlation coefficient around ρ_φbΩb ~ 0.98, with higher values of Ω_b preferred for larger bar angles.

For the band at 40 deg the maximum and the mean of the pdf are similar. They are also close to the input value of the model (light blue cross), which is in the limit between the 1σ and 2σ confidence regions. However, if we remove the (problematic) last bin in R for this band, the input value lies well inside the 1σ region. For the band at 60 deg, the mean and the maxima of the pdf differ significantly. This is because the probability distribution is flatter and more asymmetric, and in this case the mean can be considered a better estimate and more representative as it takes into account the skewness of the pdf. The input value is close to the mean value and falls well inside the 1σ region. From the values in Table 2 we see that the recovered values can present an offset with respect to the input parameters of around ~10 deg in the orientation but only of 0.04Ω₀ or ~1 km s^-1 kpc^-1 for the pattern speed. Nevertheless, input and recovered values are consistent given the standard deviations. We obtain similar results for the other bands.

The two-dimensional pdf contours are approximately elliptical and can be locally well approximated by a multivariate Gaussian centred on the expected values and with a covariance matrix given by the values of Table 2. This approximation allows us to establish a tighter joint constraint on the set of parameters (orientation and pattern speed). Furthermore, if we have independent constraints on the bar’s orientation φ_b = φ_b1, our conditional best estimate for Ω_b would be: $\mathit{E}\mathrm{(}{\mathrm{\Omega }}_{\mathrm{b}}\mathit{/}{\mathrm{\Omega }}_{\mathrm{0}}\mathrm{|}{\mathit{\phi }}_{\mathrm{b}}\mathrm{=}{{\mathit{\phi }}_{\mathrm{b}}}_{\mathrm{1}}\mathrm{)}\mathrm{=}\mathit{E}\mathrm{(}{\mathrm{\Omega }}_{\mathrm{b}}\mathit{/}{\mathrm{\Omega }}_{\mathrm{0}}\mathrm{)}\mathrm{+}\frac{{\mathit{\rho }}_{{\mathit{\phi }}_{\mathrm{b}}{\mathrm{\Omega }}_{\mathrm{b}}}{\mathit{\sigma }}_{{\mathrm{\Omega }}_{\mathrm{b}}\mathit{/}{\mathrm{\Omega }}_{\mathrm{0}}}}{{\mathit{\sigma }}_{{\mathit{\phi }}_{\mathrm{b}}}}\left({{\mathit{\phi }}_{\mathrm{b}}}_{\mathrm{1}}\mathrm{-}\mathit{E}\mathrm{\left(}{\mathit{\phi }}_{\mathrm{b}}\mathrm{\right)}\right)$(8)with a variance: $\mathit{Var}\mathrm{(}{\mathrm{\Omega }}_{\mathrm{b}}\mathit{/}{\mathrm{\Omega }}_{\mathrm{0}}\mathrm{|}{\mathit{\phi }}_{\mathrm{b}}\mathrm{=}{{\mathit{\phi }}_{\mathrm{b}}}_{\mathrm{1}}\mathrm{)}\mathrm{=}{\mathit{\sigma }}_{{\mathrm{\Omega }}_{\mathrm{b}}\mathit{/}{\mathrm{\Omega }}_{\mathrm{0}}}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{\rho }}_{{\mathit{\phi }}_{\mathrm{b}}{\mathrm{\Omega }}_{\mathrm{b}}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{.}$(9)This linear relation is shown as a green line in Fig. 7. For example, we might put a prior on φ_b to be the exact input value, i.e. 40 or 60 deg. The resulting conditional expected values are indicated in the last two columns of Table 2 (in units of Ω₀ and in   km s^-1 kpc^-1) and we see that we recover with high accuracy (1%) the input pattern speed. This can also be seen in Fig. 7 where the light blue cross almost lies on top of the green line.

If we marginalise the pdf’s of Fig. 7, we obtain the best estimates for each individual parameter independently on the rest of parameters and their corresponding confidence intervals. The maxima of the individual marginalisations are shown with a orange triangle in Fig. 7. Whereas for the band at 40 deg, this yields a peak that is close to the maximum in the two-dimensional pdf, for 60 deg the new peak is completely off. This is because the global pdf is highly degenerate and asymmetric, and the one-dimensional pdf’s do not capture the main correlation between the parameters, giving unsatisfactory results. Therefore our best results are given when the two parameters are simultaneously estimated.

4. Hercules in the RAVE data

We use now the RAVE survey data to measure the position of the saddle point v_φ,OLR as a function of R. Our aim is to establish whether the observed trend is consistent with the analytic model developed in Sect. 2.2 and to constrain the bar properties through best fits to the observations.

4.1. The RAVE data

We use the RAVE DR4 (Kordopatis et al. 2013) and the distance determination method by Zwitter et al. (2010), which leads to a new data set with 315 572 stars³. Alternatively, we may also use the distances from the method by Burnett et al. (2011), presented in Binney et al. (2014) which are in fact the recommended ones for DR4. We shall see later, however, that we obtain very similar results in both cases. The stellar atmospheric parameters of the DR4 are computed using a new pipeline, based on the algorithms of MATISSE and DEGAS, and presented in Kordopatis et al. (2011). Compared to DR3, DR4 is 5 times larger and the spectral degeneracies and the 2MASS photometric information are better taken into consideration, improving the parameter determination (and hence the distance estimation) with respect to previous data releases. We use proper motions from different catalogues, mainly PPMX (Röser et al. 2008) and UCAC2 (Zacharias et al. 2004), choosing from each catalogue the values with the smallest errors.

Following Honma et al. (2012) we use a position of the Sun of (X,Y) = ( − 8.05,0) kpc and a circular velocity at the Sun of v₀ = 238 km s^-1 to compute the positions and cylindrical velocities v_R and v_φ of the stars. We adopt the velocities of the Sun with respect to the local standard of rest of (U_⊙,V_⊙,W_⊙) = (10,12,7) km s^-1 from Schönrich et al. (2010). We examine later on the implications of these adopted values on our results. The value of (v₀ + V_⊙)/R₀ by Honma et al. (2012) is 31.09 ± 0.78 km s^-1 kpc^-1 which is compatible with the value from the reflex of the motion of the Sgr A* 30.2 ± 0.2 km s^-1 kpc^-1 (Reid & Brunthaler 2004).

First, we select 274103 stars in the plane with | Z |  ≤ 1 kpc, as for this range of heights we expect to be able to detect the effects of the bar (Monari et al. 2014). Figure 8 (grey dots) shows the positions of these selected RAVE stars. From the stars within 1 kpc from the plane, we select a band of stars at 6 deg with respect to the line Sun-GC with a width of Δφ = 6 deg (blue dots in Fig. 8). There are in total 71 605 stars in this band, of which 94% are giants. The median relative error in distance for this band is 27%, the median error in transverse velocity is 20 km s^-1, whereas radial velocity errors are smaller than 1.5 km s^-1 for 90% of the stars. If φ_b is the orientation of the bar with respect to the solar neighbourhood, the orientation of this band with respect to the bar is φ_b + 6 deg.

We choose this band because it covers a large range of R while keeping the errors in distance and kinematics small. Ideally, one could use data on different bands, i.e. not restricted to a given φ. However, given the quantity and quality of the current data, this does not improve the fit on the parameters: bands at other angles have less stars and larger kinematic errors or cover a smaller range of radii. First, as one can see in Fig. 8, the RAVE data extends far beyond the blue band selected. However, at these locations the Hercules stream is hardly recognisable. The stars are at least 1.5 − 2. kpc far from the Sun and their distance and kinematic errors are large. One could also take a band for the φ of the solar neighbourhood. However, the data for this band does not cover a large range of R and is so close in angle (only 6 deg) to the band selected that we do not improve the fit by combining the information of the two bands.

As in Sect. 3, we divide the band in bins of R but we now take bins every 0.2 kpc with a width of ΔR = 0.2 kpc. In Fig. 9 (top) we show the number of stars per bin. The range of radii that we probe for this band is [7.8,8.6] kpc (red asterisks in the figure). Outside this range we fail to detect the Hercules structure. This may be due to observational errors, to the fact that Hercules is masked by the other groups and due to the number of stars which decreases substantially. On the other hand, the average height of the stars (bottom panel of Fig. 9) increases significantly outside the mentioned range due to the RAVE fields selection. For large heights above the plane the kinematic structures may be also diluted and, additionally, the behaviour of the orbital frequencies can be different from those in-plane which may invalidate Eq. (5).

4.2. Measuring v_φ,OLR in the RAVE data

We now follow the steps outlined in Sect. 3.2 using the RAVE data. The process is shown in Fig. 10. The first column shows that the Hercules stream is not as clear as in the simulations depicted in Fig. 5. This may be due to several reasons. First, the presence of observational errors dilutes velocity structures. Also the existence of other kinematic groups in the data may mask the presence of Hercules. Second, the bins in the RAVE data are located at relatively large heights from the plane (Fig. 9) which can also wash out the Hercules signal. Third, the test particle simulation of Sect. 3 shows a particularly conspicuous OLR mode but Monari et al. (2014) have shown recently that for longer integration times (i.e. older bars) the distinction between the MAIN and the OLR modes is less clear and more similar to observations. The strength of the bar can also influence in the conspicuousness of Hercules.

Due to the above limitations we introduce a change in our method with respect to Sect. 3. This is because the adaptive kernel density estimator produces a weak signature of the Hercules peak, sometimes seen only as an inflection point. We therefore prefer to use the wavelet transform (WT) instead. This is especially suitable to enhance overdensities and underdensities (Fig. 10 second column), and has been applied extensively for the detection of kinematic groups (Skuljan et al. 1999; Antoja et al. 2008, A12)⁴. We use here a range of scales between 22 − 45 km s^-1 (see A12). The WT detects also other peaks such as Hyades or Sirius apart from the Hercules group in the distribution of ${\mathit{v}}_{\mathit{\phi }}^{\mathit{\theta }}$ (third column in Fig. 10). Figure 11 shows the v_φ,OLR for the different bins which decreases with R as expected if Hercules is caused by the bar’s OLR. The data points are also tabulated in Appendix A.

5. Results: application of the analytic model to the RAVE data

We proceed to obtain the most likely bar properties consistent with the RAVE data. We use the maximisation parameter ranges as explained in Sect. 3.3. Note that we had to assume values for v₀, V_⊙, U_⊙ and R₀ to compute the individual v_φ, v_R and R from the observables. In Sect. 5.1 we keep these parameters fixed, while in Sect. 5.2 we consider also changes in these parameters. In Sect. 5.3 we discuss the effect of the observational errors and possible biases in distance.

5.1. Results for fixed solar parameters

Figure 12 shows the two-dimensional marginalised pdf’s φ_b-Ω_b (left), β-Ω_b (middle) and φ_b-β (right). The first panel presents a well defined peak. By contrast, the other panels show flatter distributions, especially for the slope β of the rotation curve for which we do not obtain any constraint. In Table 2 (Model 1) we give the details of the pdf of φ_b-Ω_b, that is the maximum of the probability, the expectation of each parameter E(φ_b,Ω_b), the standard deviations of the distribution σ_φb and σ_Ωb/Ω0, and the correlation ρ_φbΩb. From Fig. 12 we can observe (as in the simulations of Sect. 3.4) the strong correlation between φ_b and Ω_b with a correlation coefficient of ρ_φbΩb = 0.98. Correlations between other parameters are much smaller: ρ_βΩb = −0.03 and ρ_φbβ = 0.05.

The maximum of the pdf is the orange cross in Fig. 12 and is located at (φ_b,Ω_b/Ω₀) = (45 deg,1.89). The expected values are shown as a orange asterisk. For the choice of the parameters v₀ = 238 km s^-1 and R₀ = 8.05 kpc, the pattern speed of 1.89Ω₀ corresponds to Ω_b = 56.0 km s^-1 kpc^-1. There are no significant differences between the maximum and the mean of the pdf for Ω_b, as they differ only by 1%. For the bar’s orientation φ_b we obtain a broader likelihood distribution than for Ω_b. In the left panel of Fig. 12 we see that the 1σ region (solid red line) covers almost the whole range of φ_b (from ~20 deg to ~80 deg). The maximum of the pdf and its mean differ by 6%. The fit given by the maximum of the two-dimensional pdf is plotted on top of the data points in Fig. 11 (black curves) labelled as “2d max”, for three different values of β.

The fact that we can constrain the value of the pattern speed but not the orientation is mainly due to the small range of radius (~600 pc), in comparison to the toy model (2 − 3 kpc). This is expected from inspection of Fig. 3, where we see that different pattern speeds occupy distinct regions, while curves for different angles or slopes of the rotation curve can be rather close to each other for certain R and may become indistinguishable given the errors. We also performed a test reducing the range of radii for our simulations to establish whether the quality of the constraint obtained depended on the radial range considered. In practise, for each of the four bands (20, 40, 60 and 80 deg) we repeated the fit using only the innermost bins in a 600 pc radius range, and separately, the outermost bins. These tests showed that the effect of reducing the radius range doubles the uncertainty in the constrained φ_b, making it of the order of 20 deg, as for the RAVE data.

Although the errors could also be a cause of our weak constraint in the bar’s orientation for the RAVE data, the precision to which v_φ,OLR is determined is not so different for the simulations and the data. For RAVE most of the errors are between 2 and 3   km s^-1. For the simulations, although the uncertainties in the determination of v_φ,OLR are larger than 3 km s^-1 for 40% of the bins (for the outer bins) we still determine the bar orientation more accurately in the simulations (to 10 deg) than for the RAVE data (~20 deg).

One could marginalise the pdf’s of Fig. 12 to obtain the best estimates for each individual parameter. If we proceed in this way, we get the maxima of the individual marginalisations as shown with a orange triangle in Fig. 12 left. For Ω_b the maximum of the two-dimensional pdf is similar to the one-dimensional maximum, differing only by 4%. However, note that for φ_b the maximum of the one-dimensional pdf is quite different (55%) from that obtained from the two-dimensional φ_b-Ω_b panel. This is analogous to what happened in the case of our simulations in Sect. 3, and is due to the global pdf being degenerate and skewed especially in the φ_b direction. The resulting fit of the one-dimensional pdf’s is the red curve shown in Fig. 11 labelled as “1d max”. This curve fits very poorly our data, showing once more that the one-dimensional pdf’s do not capture the main correlation between the parameters and give misleading results.

Because of the tight correlation between Ω_b and φ_b and the large dispersion in the probability for φ_b, we actually obtain a better fit and a tighter constraint when we use the appropriate combination of parameters. In the same manner as in Sect. 3, under the bivariate normal approximation, using Eqs. (8) and (9) we can establish a linear relation between φ_b and Ω_b that allows us to obtain the best estimate of Ω_b given a particular value of φ_b. We obtain: $\mathit{E}\mathrm{(}{\mathrm{\Omega }}_{\mathrm{b}}\mathit{/}{\mathrm{\Omega }}_{\mathrm{0}}\mathrm{|}{\mathit{\phi }}_{\mathrm{b}}\mathrm{=}{{\mathit{\phi }}_{\mathrm{b}}}_{\mathrm{1}}\mathrm{)}\mathrm{=}\mathrm{1.91}\mathrm{+}\mathrm{0.0044}\left({{\mathit{\phi }}_{\mathrm{b}}}_{\mathrm{1}}\mathrm{\left(}\mathrm{deg}\mathrm{\right)}\mathrm{-}\mathrm{48}\right)$(10)with standard deviation of 0.02. In units of   km s^-1 kpc^-1, when we fix the solar parameters to v₀ = 238 km s^-1 and R₀ = 8.05 kpc, this is: $\mathit{E}\mathrm{\left(}{\mathrm{\Omega }}_{\mathrm{b}}\mathrm{\right|}{\mathit{\phi }}_{\mathrm{b}}\mathrm{=}{{\mathit{\phi }}_{\mathrm{b}}}_{\mathrm{1}}\mathrm{)}\mathrm{=}\mathrm{56.3}\mathrm{+}\mathrm{0.1316}\left({{\mathit{\phi }}_{\mathrm{b}}}_{\mathrm{1}}\mathrm{\left(}\mathrm{deg}\mathrm{\right)}\mathrm{-}\mathrm{48}\right)$(11)with standard deviation of 0.5 km s^-1 kpc^-1. The green line in the left panel of Fig. 12 indicates this linear relation. An example is shown in the last two columns of Table 2 (in units of Ω₀ and in   km s^-1 kpc^-1). For the angle of φ_b = 30 deg we obtain a pattern speed of 54.0 ± 0.5 km s^-1 kpc^-1. This model is shown on top of the data points in Fig. 11 (blue curves) labelled “30 deg” for three different values of β. We see that these curves fit better the data points, compared to the curve for the maxima of the one-dimensional marginalised pdf’s (black curve). In the range of 10 − 45 deg, which as explained in Sect. 1 encompasses the bar orientations independently estimated in the literature, we would obtain a range of pattern speed of 51.3 − 55.9 km s^-1 kpc^-1.

Figure 13 shows the 1σ confidence limits for slices of the three-dimensional probability at different values of β compared to Model 1 (marginalised over β). The curves do not differ significantly, meaning that the dependence on the β parameter is not strong. For instance, when we fix our model to β = 0 (Model 2 in Table 3)⁵, we obtain similar results for the pdf of φ_b and Ω_b when compared to Model 1 (blue and red curves in Fig. 13). For β ≠ 0 the curves are similar only for φ_b around the two-dimensional maxima (orange cross). This is why the fit obtained with the two-dimensional maximum (marginalised over β) fits the trend in Fig. 11 for different values of β. For other angles (away from the 2D maximum), there is though a slight dependence on β. Under the assumption that the parametrisation of the rotation curve of Eq. (1) is valid, this figure demonstrates that φ_b = 30 deg, and β = 0.2 (or in general positive β) are favoured by the data only if the pattern speed is higher. As clearly shown in Fig. 11 with the blue-dashed curve, the linear relation of Eq. (10) is not valid for positive β.

5.2. Varying the solar parameters

As explained before, we had to assume values for v₀, V_⊙, U_⊙ and R₀ to compute the individual v_φ, v_R and R from the observables. The adoption of this specific set of values has an effect on our derived bar’s parameters. In this section, though, we shall see that the effects are in practise little or they lead to only a scaling on the obtained pattern speed.

First, we changed the value of V_⊙ from 12 km s^-1 to 5 km s^-1 (e.g. Dehnen & Binney 1998). To do this, we would have to recompute the velocities of our RAVE data points using this new value and re-run our method to find v_φ,OLR as a function of R. In particular, v_φ (which is one of the required observables) is obtained by adding to the “measured” heliocentric velocities the adopted values of v₀ and V_⊙ and rotating them by an angle that depends on the position of the star in the disc (which in turn depends on the positions in the sky, distances from the Sun and the adopted value of R₀). As a short-cut to this, one can see that for a particular angular band the change on the individual v_φ of a star due to a change in V_⊙ and v₀ will translate into a shift of the measured v_φ,OLR. Therefore, our new determinations v_φ,OLR can be recomputed from the old ones by adding a factor (5 − 12)cos(6 deg) km s^-1, where 6 deg is the angle of the selected band. Model 3 of Table 3 shows the results for this change in V_⊙. We see that this change slightly reduces the expected value of Ω_b to 1.85Ω₀ (Ω_b = 54.7 km s^-1 kpc^-1). Also the conditional value of Ω_b for φ_b = 30 deg is reduced to 52.6 ± 0.5 km s^-1 kpc^-1. On the other hand, for the band considered here (only at 6 deg from the line GC-Sun) the parameter U_⊙ has little influence on the computation of v_φ.

Following the idea of the previous paragraph, we can also turn parameters such as v₀ into free parameters of the model without first having to compute the new velocities v_φ of the data points for many values of v₀, then derive for each case the new determinations of v_φ,OLR and finally do the fitting process. Model 4 in Table 3 is the best fit obtained when v₀ is a free parameter. We explore this parameter in the range of [224,252] km s^-1 (Honma et al. 2012) using steps of 1 km s^-1. This change, however, does not affect the determinations of Ω_b/Ω₀ and φ_b with respect to Model 1. Figure 14 shows the pdf in the Ω_b − v₀ plane, with a correlation coefficient of σ_Ωbv0 = −0.03. For v₀ of 224 or 252 km s^-1 we get similar best fit pattern speeds in the combined φ_b-Ω_b pdf (1.88Ω₀ and 1.90Ω₀, respectively). However, once scaled to the respective Ω₀, the pattern speeds become 52.4 and 59.5 km s^-1 kpc^-1. On the other hand, the pdf for v₀ is very flat as can be seen in Fig. 14.

5.3. Analysis of errors and biases

Here we explore the influence of the observational errors and biases on the analysis of Sect. 5.1 by considering four possible cases. In the first case (Model 5) we consider only stars with velocity errors in the v_R and v_φ directions smaller than 15 km s^-1. This contains 35% of our initial sample and has a total of 26 076 stars. In the second case (Model 6), we use the distances from Binney et al. (2014) obtained with the method by Burnett et al. (2011), instead of the one by Zwitter et al. (2010). Finally, we also explore how a bias in distance would affect our results. We redo the analysis considering the extreme cases of having distances overestimated (thus we reduce the original values) and underestimated (thus we increase the original values) by 30%. These are Models 7 and 8, respectively. The new measured values of v_φ,OLR for these four cases are shown in Fig. 15 with different symbols and colours. Using the same symbols, in Fig. 16 we show the maximum in the φ_b-Ω_b plane, and the respective 1σ confidence limits are marked with a dotted line. In Table 3 we give the details of the two-dimensional pdf’s in the same manner as the previous cases. The data points of each case are tabulated in Appendix A.

The results for these four additional cases are similar to Model 1. For example, for the distance method by Burnett et al. (2011) (Model 6) we find an almost identical two-dimensional maximum and 1σ contour. For the four cases, the two-dimensional maxima are all located inside or very close to the 1σ contour of our standard Model 1. Moreover, the maxima are only shifted along the direction of degeneracy of Model 1 and their confidence regions also follow the same degeneracy. We do see, however, that the bar’s orientation is more sensitive to observational errors and biases. For instance, the sample with the smaller errors (Model 5, blue diamond) has a maximum for a bar orientation that is 4 deg smaller than for Model 1 but the same expectation value. We also find a smaller (larger) bar’s orientation in Models 7 and 8 when we correct the distances supposing that they were overestimated (underestimated), although they are consistent within the errors. These differences are because biases in distance systematically change the slope of the relation between v_φ,OLR and R. We obtain similar results for the best pattern speed, whereas the value for a fixed orientation of 30 deg changes at most by 0.7 km s^-1 kpc^-1.

Finally, we also perform a test to assess whether the observational errors and the resulting blurring of the substructures could produce a bias in the derived models parameters. To this end we convolve the simulations of Sect. 3 with typical RAVE errors. More details are given in Appendix B.

The results of this test indicate that, as expected, the addition of errors produces a blurring of the structures which makes impossible the Hercules detection in certain radial bins, especially those located farthest from the simulated Sun. In the bins where Hercules can still be detected we observe a slight tendency to obtain smaller estimates of v_φ,OLR for the farthest bins and that appears to increase with distance. This seems to be due to contamination from stars located originally farther away and which fall in our sample because of their large distance error more than to kinematic errors only. In our parameter fits, the mean or expectation of both the pattern speed and the bar’s orientation are consistent with the true parameters within 1σ in most of the cases. However, the obtained correlation between the angle and the pattern speed presents a small bias. In particular, if the true bar orientation is assumed, it leads to an overestimation of the pattern speed by an amount between 0.6 and 1.1 km s^-1 kpc^-1 depending on the orientation.

Note that this test is rather simple at least regarding our RAVE error model and selection function and that it is based on a very idealized model. Moreover, the mentioned bias could be different depending on the model parameter region that we are probing. For these reasons one cannot conclude that our obtained bar parameters with RAVE in the case of the linear relation are overestimated by the cited amount. We rather use these results to estimate that the approximate systematic error due to the RAVE data precision in the determination of the pattern speed could be around 1 km s^-1 kpc^-1.

6. Discussion and conclusions

We have derived the pattern speed of the Galactic bar from the analysis of the kinematics of the Hercules stream at different Galactocentric radii, assuming that Hercules is caused by the effects of the bar’s OLR. The crucial observable for this measurement is the azimuthal velocity of the saddle point that separates Hercules from the main part of the velocity distribution.

In particular, starting from the model by D00, we have derived an analytical expression for how the azimuthal velocity of the saddle point changes as a function of position in the Galaxy and its dependence on the properties of a barred potential, namely, the bar’s pattern speed, orientation, and the slope and normalisation of the rotation curve. We then used data from the RAVE survey to measure this velocity as a function of Galactocentric radius. We have found that it decreases with radius in a manner that is consistent with our analytic model. By fitting the measured trend, we have derived the best fit parameters of the Galactic bar. To our knowledge, this is the first time that the information on how a moving group changes as a function of radius is used in deriving the parameters of the non-axisymmetries of the disc.

We tested the reliability of our analysis by comparing the model predictions with the “measurements” of the velocity of the saddle point in a toy model consisting of a test particle simulation. Although the analytical model was derived using the stellar orbital frequencies for simple power-law Galactic potentials, it was found to reproduce well the trends found with a more complex Galactic potential (with three components: halo, bulge and disc). Our method to locate the velocity of the saddle point successfully finds velocities that are consistent with the predictions and we recover the input parameters of our simulation in most cases inside the 1σ confidence region. We emphasise that a much more accurate constraint is obtained when the proper combination of Ω_b and φ_b (which are largely degenerate) and some prior information on φ_b are used.

Our model has provided new constraints for the parameters of the Milky Way bar. The likelihood function of the pattern speed and the bar’s angle is highly degenerate. We find that the combined likelihood is maximum for a bar’s pattern speed of Ω_b = (1.89 ± 0.08) × Ω₀, where the latter is the local circular frequency. Assuming a solar radius of 8.05 kpc and a local circular velocity of 238 km s^-1, this corresponds to a pattern speed of 56 km s^-1 kpc^-1 with a standard deviation of ~2 km s^-1 kpc^-1. Also, because of the high correlation between φ_b and Ω_b, we find that a better description of our best fit results is given by the linear relation E(Ω_b/Ω₀ | φ_b = φ_b1) = 1.91 + 0.0044(φ_b1(deg) − 48) with standard deviation of 0.02. For the angle of φ_b = 30 deg we obtain a pattern speed of 54.0 ± 0.5 km s^-1 kpc^-1, reducing further the uncertainty in this determination. In the range of bar’s orientation of 10 − 45 deg, as other independent studies suggest, we obtain a range of pattern speed of 51.3 − 55.9 km s^-1 kpc^-1. Tests made by adding typical RAVE errors to the the toy model indicate that RAVE-like errors could produce systematic errors in the pattern speed of around 1 km s^-1 kpc^-1 when estimated using this linear relation.

The determination of Ω_b in units of Ω₀ is not very sensitive (typically only by a few centesimal digits) to the assumed Galactic parameters such as the circular velocity curve, the peculiar motion of the Sun, or to different distance determination methods or biases in distance. For instance, using a smaller value for the peculiar velocity of the Sun V_⊙ reduces the best estimate for the pattern speed of the bar by ~2 km s^-1 kpc^-1.

Our result agrees well with the pattern speed derived with other methods compiled in Gerhard (2011). For instance, with a direct determination of the pattern speed using the Tremaine-Weinberg method a value of 2.15 ± 0.54Ω₀ was obtained (Debattista et al. 2002). Several studies (e.g. Englmaier & Gerhard 1999; Fux 1999) that fitted hydrodynamical models to the observed CO lv-diagram give a combined value of 1.89 ± 0.36Ω₀. Other determinations that are also consistent with ours come from the analysis of the kinematics of solar neighbourhood stars. For example, the observed trend with velocity dispersion of the Oort constant C could be due to the bar’s effects with a pattern speed of 1.87 ± 0.04Ω₀ according to Minchev et al. (2007). Also a value of 1.82 ± 0.07Ω₀ could be responsible for the existence of several low-velocity streams in the solar neighbourhood such as Pleiades and Coma Berenices, or Pleiades and Sirius (Minchev et al. 2010). Compared to the value of D00 of 1.85 ± 0.15Ω₀, who fitted the same stream as in the present work but only locally, we have obtained a consistent value but a tighter constraint.

On the other hand, we find that the bar’s orientation cannot be constrained well as the 1σ confidence region extends from 20 to 80 deg. This is mainly due to limitations in the data, in particular the small range of radii that it spans. This conclusion stems from the fact that we could constrain better this parameter for our toy model, which can be studied for a larger range of radii even if the the errors in the determination of the azimuthal velocity of the saddle point are comparable or larger than those of the data. Additional tests reducing the radius range of the simulations to make it similar to RAVE showed that the obtained uncertainties in the bar’s orientation are larger and more comparable to the data.

Alternative interpretations and models to explain the Hercules stream have been proposed in the literature. Fux (2001) in his study with test-particle and N-Body simulations suggested that near the external regions of the OLR the bar generates an overdensity of stars in velocity space that is made of chaotic orbits, which is a different interpretation to the scattering mechanism proposed by D00. Quillen et al. (2011) showed that features similar to Hercules can be associated to the coupling of several spiral structures. Also Antoja et al. (2009) found a group similar to Hercules in test-particle simulations containing only spiral arms, although a bar was required to move this feature to negative U. Our findings of the trend of Hercules with radius provide evidence that this stream may indeed be due to the effects of the bar’s OLR, but whether these other models can account for such a trend remains to be seen.

An analysis with i) a sample covering larger range in radius and regions with different bar’s orientation; ii) with more stars per bin; and iii) with smaller errors in distances and proper motions should allow us to constrain the bar’s orientation and even better the pattern speed. The astrometric data from the ESA’s Gaia mission will provide us with such numerous, extended and precise observations. A clear benefit would also be obtained when observations spanning several bands in azimuth could be used at the same time for tighter constraints.

Online material

Appendix A: v_φ,OLR and errors

In this appendix we tabulate the data points of Figs. 15 and 11 used for our different fits, and their errors (Table A.1).

Appendix B: Effects of the RAVE errors

Here we perform some simple tests to assess whether the observational errors produce any bias in the derived models parameters. To this end we convolve the simulations of Sect. 3 with typical RAVE errors. In particular, we use the median errors of the radial bins that we use in our analysis. That is we take a relative error in distance of 25%, error in proper motion (both in α and δ) of 1.9 mas/yr, and error in radial velocity of 0.9 km s^-1.

More in detail, we first orient the disc of the simulation so that the desired band with different position angles for the bar of 20,40,60 or 80 deg has the same orientation with respect to the solar position (6 deg) as our RAVE band. We select particles in a band with the same angular width as in the observations and we then convolve the positions and velocities with the mentioned errors (assuming they are Gaussian). Afterwards we bin the simulation in the same way as the data. We select 4 bins located in the same (relative) position with respect to the simulated Sun as in the data to have particles with a similar distance distribution (and, therefore, similar errors in distance and in transverse velocity).

An example of the final velocity distribution for the 4 bins for the band of 40 deg is shown in Fig. B.1 (middle). The new bins are smaller in extent and therefore have significantly less particles compared to Fig. 5. In order to do a proper comparison, we include here also the original simulation (i.e. without error convolution) in the same radial bins (left column). In the upper left corner we indicate the median direction of the transverse movement with an arrow. The blurring of the substructures occurs along this direction as this is the one influenced by errors in distance and proper motion, which are significantly larger than errors in line of sight velocity. The velocity distribution is distorted along this direction as explained in A12 and in Fig. 9 of McMillan (2013). For some of the bins (outside the range presented here), the blurring is large enough that the Hercules gap is no longer detected. The final determinations of v_φ,OLR for all bands (black symbols in Fig. B.2) fall close to the expected theoretical lines in most of the bins. Because of this no bias is observed in our final bar’s parameter determination.

We now perform another test where we consider also particles originally located farther away in the disc but that, due to distance errors, end up in the selected band after error convolution. For this we first select particles with a maximum distance of 2.4 kpc from the Sun’s position, we convolve with the RAVE errors, and finally we take the subset of particles in the bands. The limit of 2.4 kpc corresponds to the maximum distance at which a red clump star with absolute magnitude of M_J = −0.87 would be observed by the RAVE survey assuming that the upper magnitude limit of the survey is J ~ 11. With this limit we avoid including particles in our band that were originally very far from the Sun and therefore, that would have never been observed because of the magnitude limit of the survey.

The example for the band of 40 deg is shown in the right column of Fig. B.1. The final determinations of v_φ,OLR for all bands are the colour symbols in Fig. B.2. We observe in this case a slight tendency to obtain smaller estimates of v_φ,OLR for the bins that are far from the simulated Sun (the ones at larger radius) and that seems to increase with distance. However, as the RAVE bins are still quite close to the Sun with a maximum distance of 1.67 kpc, this bias is not very significant.

The recovered parameters for all bands after RAVE error convolution are shown in Table B.1. We see how in most of the cases the correct model parameters both in angle and in pattern speed are recovered within the error bars. However, the recovered correlation between the angle and the pattern speed does present a bias. In particular, the recovered pattern speeds obtained with the linear relation using the correct (true) value of the bar’s orientation are slightly larger than expected. In all cases the difference between the recovered value and the correct one is between 0.6 and 1.1 km s^-1 kpc^-1, increasing with bar’s orientation. This is in all cases equal or less than 2σ. We remark that we do not see this bias in the expectation values of the likelihood.

---

Acknowledgments

We acknowledge funding support from the European Research Council under ERC-StG grant GALACTICA-240271. Funding for RAVE has been provided by: the Australian Astronomical Observatory; the Leibniz-Institut fuer Astrophysik Potsdam (AIP); the Australian National University; the Australian Research Council; the French National Research Agency; the German Research Foundation (SPP 1177 and SFB 881); the European Research Council (ERC-StG 240271 Galactica); the Istituto Nazionale di Astrofisica at Padova; The Johns Hopkins University; the National Science Foundation of the USA (AST-0908326); the W.M. Keck foundation; the Macquarie University; the Netherlands Research School for Astronomy; the Natural Sciences and Engineering Research Council of Canada; the Slovenian Research Agency; the Swiss National Science Foundation; the Science & Technology Facilities Council of the UK; Opticon; Strasbourg Observatory; and the Universities of Groningen, Heidelberg, and Sydney. The RAVE web site is at http://www.rave-survey.org. E.K.G. acknowledges funding from SFB 881 “The Milky Way System” of the German Research Foundation, particularly through subproject A5. We also thank the anonymous referee for the helpful comments and suggestions which significantly improved the first version of this paper.

References

All Tables

All Figures