EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 664 (August 2022) A&A, 664 (2022) A148 Full HTML
Open Access
Issue
A&A
Volume 664, August 2022
Article Number A148
Number of page(s) 15
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/202142251
Published online 23 August 2022

© Z. Prudil et al. 2022

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe-to-Open model. Subscribe to A&A to support open access publication.

1. Introduction

Our Galaxy, the Milky Way (MW), harbors stellar populations with a range of dynamical properties that reflect the underlying mass distribution. In particular, stars on the edge (and beyond) of the MW escape velocity distribution provide an important perspective on both the small and large dynamical scales of the MW environment.

Focusing on the small scale, stars in a galaxy can be accelerated by dynamical interactions with their close surroundings. Principally, two main production classes for high-velocity stars have been introduced. These “runaway stars” are mainly young (O-B type) stars found far away from a star-forming galaxy disk, deep in the stellar halo (Blaauw 1961). The mechanism responsible for their ejection involves two-body (supernova explosion in a binary, e.g., Blaauw 1961; Portegies Zwart 2000; Evans et al. 2020) and multi-body systems (close encounters in stellar systems such as young massive star clusters, e.g., Poveda et al. 1967; Leonard 1991; Gvaramadze et al. 2009). Typically, the runaway stars from the aforementioned production channels can attain velocities from several dozens up to low-hundreds of km s−1 (Portegies Zwart 2000; Perets & Šubr 2012), with very few exceeding the Galactic escape velocity (so called hyper-runaway stars, e.g., Perets & Šubr 2012; Brown et al. 2015; Evans et al. 2020).

The second type of extreme velocity stars is hypervelocity stars, which are the fastest moving stars in a galaxy. The ejection mechanism of hypervelocity stars is linked to a supermassive black hole (SMBH) in the Galactic center (Hills 1988), where during a close encounter, a binary system can be disrupted. Then one of the components can be ejected with an extreme velocity (up to thousands of km s−1, e.g., Hills 1988; Yu & Tremaine 2003; Bromley et al. 2006) exceeding the galaxy escape velocity. The first candidate hypervelocity stars were found by Brown et al. (2005) and since then, hundreds of more candidates have been discovered (e.g., Edelmann et al. 2005; Kollmeier et al. 2009; Brown et al. 2014, 2015; Vickers et al. 2015). Especially with the advent of the Gaia astrometric mission, the number of hypervelocity candidates has increased (e.g., Hattori et al. 2018a; Marchetti et al. 2019; Marchetti 2021; Li et al. 2021), with some of the objects reaching velocities beyond > 1500 km s−1 (Koposov et al. 2020). For a full review on high-velocity and particularly hypervelocity stars, we refer to Brown (2015).

The outlined ejection mechanisms are not intrinsic only to the MW but to other star-forming galaxies as well. For example, star-forming dwarf galaxies like the Large and Small Magellanic Clouds (LMC and SMC) can be a source of stars with extreme velocities that were once ejected (runaway stars) and are now passing through the MW halo (Boubert et al. 2017; Hattori et al. 2018b; Erkal et al. 2019a). The infall of globular clusters and dwarf galaxies toward the Galactic center can lead to the ejection of stars with extreme velocities (Abadi et al. 2009; Capuzzo-Dolcetta & Fragione 2015; Fragione & Capuzzo-Dolcetta 2016) during their interaction with the SMBH. Theoretical predictions also suggest a sub-population of hypervelocity stars ejected from the Andromeda galaxy by its SMBH (Sherwin et al. 2008).

Not all stars accelerated by the mechanisms described above surpass the Galaxy’s escape velocity and become unbound. Some objects populate the high-velocity tail of the Galaxy’s velocity distribution. The escape velocity boundary at the high-velocity tail decreases with increasing Galactocentric radius and is set by the Galaxy’s potential and thus mass distribution. Therefore, weakly bound halo stars can be utilized to estimate the escape velocity at a given radius and, in turn, for a selected gravitational potential; thus, they can provide a constraint on the MW mass. The mass of our Galaxy is one of the large-scale fundamental parameters that still remain uncertain, and current estimates center around 1012 M⊙ (for an illustration of the MW mass estimates see Fig. 7 in Callingham et al. 2019). The extensive dynamical studies of the MW satellite galaxies (e.g., Watkins et al. 2010; Lépine et al. 2011; Koch et al. 2012; Callingham et al. 2019), globular clusters (e.g., Eadie & Harris 2016; Sohn et al. 2018), stellar streams (e.g., Pearson et al. 2017; Malhan & Ibata 2018; Erkal et al. 2019b), and halo stars (e.g., Shen et al. 2022) yield numerous mass estimates for the MW.

The study by Leonard & Tremaine (1990) put forward an approach to modeling the velocity distribution based on a power law. This technique has been often utilized to estimate the Galactic escape velocity and subsequently, the mass of the MW (e.g., Smith et al. 2007; Piffl et al. 2014). Particularly with the availability of astrometric data from the Gaia mission (Gaia Collaboration 2016), many studies have been focusing on estimating the Galactic escape velocity in the Solar neighborhood (i.e., Monari et al. 2018; Deason et al. 2019; Grand et al. 2019; Koppelman & Helmi 2021; Necib & Lin 2022b). In the studies mentioned above, the estimates of the Galactic escape velocity fall in the range from 480 to 640 km s−1 where a subtle difference in the analysis (prior assumptions of model parameters, minimization technique) can play a role in the resulting escape velocity. The treatment of substructures, such as the Gaia-Sausage-Enceladus (Belokurov et al. 2018; Helmi et al. 2018) and stellar streams (e.g., Ibata et al. 2019) around the Solar neighborhood can also lead to discrepancies in estimates of the Galactic escape velocity (Grand et al. 2019; Necib & Lin 2022a). These differences can lead to underestimating the escape velocity by approximately 10% (e.g., Smith et al. 2007; Grand et al. 2019); thus, caution is needed in interpreting the derived escape velocity values.

In the following, we explore both the large- and small-scale dynamical aspects of the MW environment by utilizing RR Lyrae variables. RR Lyrae pulsators are old (above 10 Gyr, see, e.g., Savino et al. 2020) helium-burning giants, located on the horizontal branch inside the instability strip (IS), exhibiting periodic variations on scales from a few hours up to one day. RR Lyrae stars can be divided into three main groups based on the pulsation mode, where fundamental-mode pulsators are denoted as RRab, first-overtone variables are referred to as RRc stars, and double-mode pulsators are marked as RRd. The periodicity in brightness changes is driven by the interplay between pressure and temperature and can be utilized together with metallicity to estimate the stars’ luminosity (PLZ relations) and, subsequently, their distance (Bono et al. 2019). All these properties make them invaluable distance indicators within the Local Group, and probes of the Galactic structure and dynamics (e.g, Koposov et al. 2019; Price-Whelan et al. 2019; Prudil et al. 2020).

RR Lyrae stars with high velocities have been identified in the Galactic bulge by the Bulge Radial Velocity Assay for RR Lyrae stars (BRAVA-RR, Kunder et al. 2015, 2020; Prudil et al. 2019) and one of them was subsequently followed up by a spectroscopic study (Hansen et al. 2016). All of these appear to be bound and kinematically associated with the MW halo (although currently passing through the Galactic bulge), nonetheless, they provide important insight into the high-velocity distribution of the old population in the Galactic center. Unlike the majority of known unbound stars in the MW, RR Lyrae variables are typically metal-poor stars ([Fe/H] ≈ −1.55 dex, Hansen et al. 2011; Crestani et al. 2021) associated with the old MW population (> 10 Gyr). It is important to mention that some old, metal-poor stars with extreme velocities have been detected in the Solar neighborhood (e.g., Hattori et al. 2018b; Du et al. 2018; Huang et al. 2021), with ages above 1 Gyr, where one of the objects is a known RR Lyrae star (Gaia-T-ES14 in Hattori et al. 2018b).

We present the second paper of our series focused on the old population pulsators as a tool to study the MW’s formation history. The present study aims at identifying high-velocity RR Lyrae stars in the MW halo that appear to be unbound and at modeling the velocity distribution at the Solar radius and beyond. The paper is organized as follows: Section 2 describes the data set, criteria imposed on the spatial and kinematical parameters, and outlines the setup for the dynamical analysis. In Sect. 3, we examine the properties of unbound RR Lyrae candidates. In Sect. 4, we attempt to model the high-velocity tail the of the RR Lyrae distribution and estimate the Galactic escape velocity at different radii. Section 5 focuses on estimating the mass of the MW using the calculated escape velocities. Finally, Sect. 6 summarizes our findings.

2. Data set and analysis

As an initial data set, we utilized our sample of 4247 RR Lyrae from the first paper of our series (Prudil et al. 2021), for which we have the complete 6D information on their position and motion (equatorial coordinates, α, δ, proper motions, μα*, μδ, distances, d, systemic velocities, vsys) together with metallicities derived using the calibration by Crestani et al. (2021) and we used also in Fabrizio et al. (2021). Our data set draws from three sources; proper motions from the third early data release of the Gaia astrometric mission (EDR3, Gaia Collaboration 2021a), distances derived based on the i-band mean magnitudes obtained by the Panoramic Survey Telescope and Rapid Response System (PanSTARRS-1, Chambers et al. 2016; Sesar et al. 2017) and systemic velocities corrected for radius variation, determined from the low-resolution spectra collected by the Sloan Digital Sky Survey in its fifteenth data release (SDSS DR15, York et al. 2000; Aguado et al. 2019). Details on the data set, including estimations of the uncertainties, can be found in Prudil et al. (2021).

To increase our chances of finding RR Lyrae variables with high velocities, we expand our initial data set with the low-resolution spectra from the LAMOST1 Experiment for Galactic Understanding and Exploration (LEGUE, Deng et al. 2012; Liu et al. 2014) survey – specifically from its sixth data release (DR6, version 2). Similarly to the SDSS, LAMOST uses a multi-object fiber-fed spectrograph providing low-resolution spectra (R ∼ 1800) with a wavelength coverage between 3700 Å and 9100 Å (Zhao et al. 2012). Each spectrum contains flux, inverse variance of the flux, and vacuum wavelengths in the heliocentric reference frame. Any object observed by the LAMOST survey has at least one co-added spectrum composed of at least one and up to six individual exposures (typically two or three), where the exposures range from a few minutes up to 45 min. Due to the overlapping fields of view for better coverage, some objects have been observed several times with different plates and fibers.

Using our RR Lyrae catalog based on the variability sample from Gaia data release 2 (DR2, Clementini et al. 2019) and the Catalina sky survey (CSS, Drake et al. 2009, 2013a,b; Drake et al. 2014; Abbas et al. 2014) we found matches for 3003 objects with a total of 4754 spectra, which also have i-band mean magnitudes in the Sesar et al. (2017) catalog of RR Lyrae stars. To estimate systemic velocities and distances for the matched stars, we proceeded similarly as in Prudil et al. (2021, see Appendices A and B).

Using the iSpec package (Blanco-Cuaresma et al. 2014; Blanco-Cuaresma 2019) we synthesized a spectrum with a characteristic set of stellar properties for a halo type RR Lyrae star (see, e.g., For et al. 2011; Sneden et al. 2017; Crestani et al. 2021). We used the following stellar parameters: [Fe/H] = −1.5 dex, Teff = 6600 K, log g = 2.25 dex, and microturbulence velocity ξturb = 3.5 km s−1. For the spectral synthesis, we used a python wrapper implemented in the iSpec module for the MOOG radiative transfer code (February 2017 version, Sneden 1973), a line list from VALD2, the solar reference scale of Asplund et al. (2009), and ATLAS9 model atmospheres (Castelli & Kurucz 2003).

In the spectral synthesis, we focused on the four most prominent features in the low-resolution spectra of RR Lyrae stars, the Balmer lines (Hα, Hβ, Hγ, Hδ) and their close neighborhood (defined as the region around each Balmer profile ±100 Å).

Wavelengths for each spectrum were transformed from the vacuum to air wavelengths using the same prescription as for the SDSS (Ciddor 1996).

Individual line-of-sight velocities were obtained using iSpec by cross-correlating each Balmer line with its synthetic counterpart. Uncertainties on the respective line-of-sight measurements were calculated through a Monte-Carlo error simulation, in which we varied fluxes using their uncertainties (under the assumption that they follow a Gaussian distribution).

The systemic velocities were determined by the minimalization process that utilizes spectroscopic products (line-of-sight velocities), photometry (pulsation type and ephemerides), and line-of-sight velocity templates (Braga et al. 2021).

We slightly modified our approach of estimating the systemic velocities due to little to no time overlap between the spectroscopic and photometric observations. Time offsets between photometric and spectroscopic observation can cause an issue since RR Lyrae stars experience changes in their pulsation periods as they evolve off the zero-age horizontal branch. These variations can lead to shifts in the time of mean brightness, M0, which potentially results in a phase offset and subsequently affects the determination of systemic velocities. The accurate estimate of the pulsation phase during the spectroscopic observation plays a crucial role. Thus, for the case of the LAMOST spectra, we increased the boundaries for our uniform (𝒰) prior on M0 from 0.1 (as for SDSS spectra see Appendix B in Prudil et al. 2021) to 0.2 to better sample the uncertainties on the posterior:

p ( θ n ) = U ( − 0.2 < Δ M 0 , n < 0.2 ) . $$ \begin{aligned} p(\boldsymbol{\theta }_{n}) = \mathcal{U} (-0.2 < \Delta M_{0,n} < 0.2). \end{aligned} $$(1)

The systemic velocity, vsys, of a given star was estimated as a weighted average between systemic velocities determined from single Balmer profiles.

In the case of distances, we used i-band mean magnitudes from Sesar et al. (2017), together with the information on metallicity from Crestani et al. (2021) for the crossmatched LAMOST portion of our data set. To account for reddening, we used recalibrated extinction maps from Schlafly & Finkbeiner (2011). All of the aforementioned quantities are a fit with the PLZ relation from Sesar et al. (2017, see Table 1) and subsequently lead to distances for individual objects from our LAMOST data set.

From the 3003 initial LAMOST RR Lyrae cross-match (within 1 arcsec), we removed 242 objects due to their problematic spectra (low signal-to-noise ratio, unphysical behavior in the spectrum, etc.). Additionally, 821 variables we previously analyzed in the SDSS data set; for those stars we used the values provided by the SDSS data products3. The remaining 1940 objects were added to our data set of 4247 RR Lyrae stars from the SDSS sample, increasing our original sample to 6187 RR Lyrae variables with 7D information on their spatial, dynamical, and metallicity distribution.

2.1. Astrometric cuts

To ensure the purity of our data sample, we employed a cut on the significance of proper motions similar to the one used in Prudil et al. (2021). Unlike Prudil et al. (2021), in our study, we required higher proper motion significance and we introduced four additional metrics on the quality of the astrometric solution based on Fabricius et al. (2021), Lindegren et al. (2021), and Gaia Collaboration (2021b):

∑ V 2 / tr ( Σ pm ∗ ) > 4.0 , $$ \begin{aligned}&\sqrt{\sum \boldsymbol{V}^{2}/\mathrm{tr}({\boldsymbol{\Sigma }^{*}_{\rm pm}}) } > 4.0,\end{aligned} $$(2)

ipd _ gof _ harmonic _ amplitude < 0.4 , $$ \begin{aligned}&\mathtt{ipd\_gof\_harmonic\_amplitude } < 0.4,\end{aligned} $$(3)

ipd _ frac _ multi _ peak < 2 , $$ \begin{aligned}&\mathtt{ipd\_frac\_multi\_peak } < 2,\end{aligned} $$(4)

astrometric _ excess _ noise _ sig ≤ 2.0 , $$ \begin{aligned}&\mathtt{astrometric\_excess\_noise\_sig } \le 2.0,\end{aligned} $$(5)

visibility _ periods _ used ≥ 13 . $$ \begin{aligned}&\mathtt{visibility\_periods\_used } \ge 13. \end{aligned} $$(6)

Here, V and Σ pm * $ {\boldsymbol{\Sigma}^{\ast}_{\mathrm{pm}}} $ represent the star’s proper motion vector and its diagonalized proper motion covariance matrix scaled by the re-normalized unit weight error (RUWE factor), respectively. Scaling the covariance matrix by the RUWE factor enabled us to weed out stars with unreliable astrometric solutions (mainly stars with RUWE > 1.4). The reference ipd_gof_harmonic_amplitude denotes the amplitude of changes in the goodness-of-fit in image parameter determination, and together with ipd_frac_multi_peak, it indicates the amount of image asymmetry, which implies possible binarity. The metric on how well the astrometric model matches the observations is expressed in the astrometric_excess_noise_sig flag, which can indicate whether the object is astrometrically well-behaved. The last criterion on visibility_periods_used ensured that our selected variables were well-observed sources. The criteria in Eqs. (2)–(6) reduced our initial sample from 6187 RR Lyrae stars down to 48474, which we refer to henceforth as the study sample.

2.2. Estimating of dynamical properties and analysis setup

To find and examine the possible extreme velocity candidates among our study sample, we utilized the python package for Galactic dynamics, galpy v1.65 (Bovy 2015). In our setup, we assumed the distance between Sun and the Galactic center to be equal to R⊙ = 8.2 kpc (Bland-Hawthorn & Gerhard 2016) and we placed the Solar system at Z⊙ = 20.8 pc (Bennett & Bovy 2019), thus slightly above the Galactic plane. For the Solar peculiar motion we used (U⊙, V⊙, W⊙)=(11.1, 12.24, 7.25) km s−1 from Schönrich et al. (2010) and we set the rotation velocity to 230 km s−1 (Eilers et al. 2019).

With the purpose of examining the orbits and possible ejection mechanism of high-velocity RR Lyrae stars, we took advantage of the composite MW potential implemented in the galpy package (MWPotential2014). The MWPotential2014 is a combination of three gravitational potentials representing the Galactic disk (Miyamoto & Nagai 1975):

Φ d ( R , z ) = − G M d R 2 + ( a d + z 2 + b d 2 ) 2 , $$ \begin{aligned} \Phi _{\rm d}(R, z) = -\frac{\mathrm{G} M_{\rm d}}{\sqrt{R^2 + \left(a_{\rm d} + \sqrt{z^2 + b_{\rm d}^2}\right)^2}}, \end{aligned} $$(7)

where G is the gravitational constant, Md, represents the mass of the disk and ad and bd stand for scale length and height, respectively.

The Galactic bulge in the MWPotential2014 is implemented as a power-law density profile with an exponential cut-off potential:

ρ r = G M b ( r 1 r ) α exp ( − ( r r c ) 2 ) , $$ \begin{aligned} \rho _{\rm r} = {G} M_{\rm b} \left(\frac{r_1}{r}\right)^{\alpha } {\exp }\left(-\left(\frac{r}{r_{\rm c}}\right)^{2}\right), \end{aligned} $$(8)

where rc represents the cut-off radius, r1 stands for a reference radius, and α serves as the inner power-law index. Finally, Mb represents the mass of the Galactic bulge.

The spherical MW halo is defined as the Navarro-Frenk-White spherical dark matter halo (NFW, Navarro et al. 1997):

ρ r = G M h 4 π a 3 1 ( r / a ) ( 1 + r / a ) 2 , $$ \begin{aligned} \rho _{\rm r} = \frac{{G} M_{\rm h}}{4\pi a^3} \frac{1}{(r/a) (1+ r/a)^2}, \end{aligned} $$(9)

where Mh defines the mass of the Galactic halo, and the scale radius is represented by the variable a. For a detailed description of the physical properties of the potential on the individual MW substructures; see Table 1 in Bovy (2015). The original implementation of the MWPotential2014 adopts a virial mass of the MW of Mvir = 0.8 × 1012 M⊙, which lies on the low-end of the virial mass estimates for the MW (see Fig. 7 in Callingham et al. 2019). For a comparison on how our results depend on the assumed MW mass, we explored two additional setups of the MWPotential2014 with slightly larger virial mass; we refer to them in a similarly way to that of the initial potential as MWPotential2014M (25% mass increase) and MWPotential2014H (50% mass increase).

Using galpy, we estimated the Galactocentric rectangular coordinates (x, y, and z) and velocities (vx, vy, and vz). We used a multivariate Gaussian distribution where individual values for astrometric properties, where distances, and systemic velocities are treated as a mean vector. For the full covariance matrix of the individual quantities, we utilized uncertainties and correlations provided in the Gaia EDR3 together with our estimates on distance and systemic velocity (where off-diagonal elements for distances and systemic velocities were set to zero).

We estimated the full covariances between individual coordinates and velocities by running a Monte Carlo (MC) error sampling (for 100 000 iterations). We also assessed the chance for a given star to be unbound (Punb) at the given escape velocity, vesc, the profile of the MWPotential2014 and the star’s Galactocentric radius. In this case, we use a definition of escape velocity from the galpy module, where it is defined as a velocity necessary to escape to infinity6:

v esc ( R ⊙ ) = 2 | Φ ( R ⊙ ) − Φ ( INF ) | . $$ \begin{aligned} v_{\rm esc}(R_{\odot }) = \sqrt{2\left| \Phi (R_{\odot }) - \Phi (\mathrm{INF})\right|}. \end{aligned} $$(10)

Here, the potential at infinity is Φ(INF)=0 and R⊙ represents the Solar radius. For each star, we estimated whether it exceeds vesc at its Galactocentric radii in each MC iteration. For a given star to be considered “unbound” in the assumed MW potential, we required it to exceed vesc in at least half of the cases of the MC simulations:

P unb > 0.5 . $$ \begin{aligned} P_{\rm unb} > 0.5. \end{aligned} $$(11)

Each RR Lyrae star that satisfied the above condition was then further examined. We refer to them henceforth as high-velocity RR Lyrae stars (HV-RRL). We note that the variable Punb does not represent the probability that each star is bound or unbound considering our model.

3. High-velocity variables

Based on our setup from the previous section (Sect. 2), we found a total of nine RR Lyrae variables that fulfilled the condition in Eq. (11). Their complete list, with calculated properties and associated covariances, as well as Gaia astrometric properties, can be found in Tables B.1 and B.2. Besides the quality cuts in Sect. 2.1, all nine candidates have RUWE < 1.4, and none of them had been flagged as a duplicated_source in Gaia EDR3. In the following subsections, we discuss their observed and calculated properties (displayed in Fig. 1) and their RR Lyrae classification.

thumbnail Fig. 1.

Observed and estimated properties for our study sample of RR Lyrae stars (light blue dots) and the candidates for high-velocity RR Lyrae stars (HV-RRL, with individual markers color-coded on the basis of the Punb). Top-left panel: depicts the color-magnitude diagram based on the Gaia passbands, the bottom left panel: spatial distribution in Galactic coordinates of our found HV-RRL candidates. The Galactic rest-frame velocity, vGC, with respect to the Galactocentric distance, RGC, is depicted in the bottom right panel. For each point fulfilling the condition in Eq. (2), we plotted also its 1σ covariance in the form of an error ellipse using the same color-coding as for the data point. Bottom left panel: the Galactic plane is denoted by the solid black line and the red dashed line outlines the Sagittarius stream. Individual blue lines in the bottom right panel represent escape velocity curves calculated using a different gravitational potential for the MW with the solid line representing the MWPotential2014, the dashed line denoting the MWPotential2014M, and the dotted line marking the MWPotential2014H.

3.1. Position in the color-magnitude diagram and variable classification

In the top left panel of Fig. 1 we display the color–magnitude diagram (CMD) for the study sample of RR Lyrae variables, via the Gaia photometry dereddened using the reddening maps from Schlafly & Finkbeiner (2011) and the extinction coefficients from Casagrande & VandenBerg (2018). In the CMD we see two populations separated by color, one at (GBP − GRP)0 ≈ 0.4 mag (characterizing first-overtone RR Lyrae variables) and a second one at (GBP − GRP)0 ≈ 0.6 mag (representing fundamental-mode RR Lyrae stars).

We see that two stars (HV-RRL-03, HV-RRL-04), shown in the top left panel of Fig. 1, with the highest value of Punb (Punb = 1.0, in all three assumed MW potentials) are located beyond the assumed red edge of the IS for the distribution of the RR Lyrae stars. Both of these stars have been identified as RR Lyrae stars in the Catalina sky survey (CSS, Drake et al. 2009) in their variable stars catalog Drake et al. (2014). They were also found in the catalog of RR Lyrae stars based on the Panoramic Survey Telescope and Rapid Response System (PanSTARRS-1, PS1 catalog, Chambers et al. 2016; Sesar et al. 2017). Unlike the CSS variable catalog, we reclassified HV-RRL-04 from a first-overtone pulsator to a fundamental-mode RR Lyrae star, based on the re-analysis of its CSS photometric data. We established a different pulsation period in comparison to the CSS catalog (for details see Prudil et al. 2021). The new classification and pulsation period are in line with the classification by the PanSTARRS-1 survey. It also agrees with the dereddened colors illustrated in the CMD, where fundamental pulsators, on average, lie closer to the red edge of the IS. This is in contrast with the first-overtone pulsators, which are generally located closer to the blue edge of the IS. To verify our classification for HV-RRL-03 and HV-RRL-04, we used a recently published catalog of variable stars in the Zwicky Transient Facility (ZTF, Masci et al. 2019; Bellm et al. 2019; Chen et al. 2020), a time-domain survey of the entire northern sky in the g and r passbands, respectively. The ZTF variable stars catalog contains photometry for both of the aforementioned stars and it was utilized in our effort to confirm or refute our initial classification.

The HV-RRL-03 (in ZTF: SourceID = 764936) is classified as a BY Dra type variable7. We find that both ours and the ZTF period analyses yield a similar period of P = 0.61771 day and low amplitudes of light changes: Amp(g)=0.233 mag, Amp(r)=0.203 mag, and AmpVCSS = 0.17 mag. Under the assumption that HV-RRL-03 is an RR Lyrae star, its low amplitude would be in agreement with its position at the red edge of the IS. If we assume its classification as a BY Dra variable, we would expect it to have a redder color (lower effective temperature, Teff, approximately between 2500 and 5200 K); while using Gaia colors and a prescription for Teff from Mucciarelli & Bellazzini (2020), we obtained temperatures between 5400 and 5800 K for a range of metallicities (expected for RR Lyrae stars from 0 dex up to −3 dex Crestani et al. 2021) and for both giants and dwarfs. This temperature range is also supported by the color-temperature relation for (V − Ks)0 from Ramírez & Meléndez (2005) and González Hernández & Bonifacio (2009). From the latter, we obtained a dereddened color for HV-RRL-03 (V − Ks)0 = 1.41 mag8 using V-band measurements from the CSS survey in combination with Ks-band measurements from the Two Micron All Sky Survey (2MASS, Cutri et al. 2003). We conclude that the Teff of HV-RRL-03 lies around ≈5900 K. This range of effective temperatures is more in line with HV-RRL-03 being an RR Lyrae star, although from photometry its classification remains uncertain. On the other hand, in Gaia astrometry, the parallax value of ϖ = 0.45 ± 0.11 mas suggests that HV-RRL-03 lies at d = 2182 − 488 + 625 $ d = 2182^{+625}_{-488} $ pc (Bailer-Jones et al. 2021), which is thus closer than our estimate based on the period-metallicity-luminosity relation (d = 19 ± 1 kpc). Such a broad difference in the distance would significantly change the calculated vGC and RGC, and subsequently mark HV-RRL-03 as bound.

The star HV-RRL-04 (in ZTF: SourceID = 190951) is classified in the ZTF catalog as an EA variable9 with a period P = 1.2928597 days10, which we subsequently confirmed by re-analyzing the ZTF photometry (using the Period04 software, Lenz & Breger 2004). We corroborated the pulsation period and variability type from the ZTF survey based on the dispersion of the points in the phased light curves. Both stars, HV-RRL-03 and HV-RRL-04, are marked with an asterisk in Tables B.1 and B.2 since their classification as RR Lyrae stars is highly uncertain (in the case of HV-RRL-03) or incorrect (in the case of HV-RRL-04); however, for the remainder of this paper, we treat them as RR Lyrae stars.

The remaining seven RR Lyrae variables are located inside the IS and their position in the CMD matches their assumed RR Lyrae sub-type. Three RR Lyrae stars from our sample (HV-RRL-02, HV-RRL-05, HV-RRL-09) have been identified by the ZTF survey as variable stars; in all three objects, both pulsation types and periods agree between our sample and the ZTF survey. One RR Lyrae variable, HV-RRL-08, has been identified as a variable star candidate by the ZTF survey with a period of 0.33348 days, while in our analysis we found a pulsation period of 0.500229 days. In the frequency space, both pulsation periods lie just 1 cd−1 11 apart, implying that one or the other is a daily alias of the true pulsation period. A careful examination of the light curves phased using both pulsation periods, together with a high amplitude of for HV-RRL-08 (AmpCSS = 0.980 mag), supports its classification as a fundamental-mode RR Lyrae star with a pulsation period of 0.500229 days. The complete list of pulsation parameters, metallicities and classification into RR Lyrae sub-types for the nine RR Lyrae variables can be found in Table 1. In addition, the phased light curves of the nine high-velocity RR Lyrae candidates are shown in Fig. A.1.

Table 1.

Pulsation and chemical properties for high-velocity RR Lyrae candidates found in our analysis.

3.2. Kinematical and spatial properties of the high-velocity sample

The lower panels of Fig. 1 depict the kinematical properties for our study sample of RR Lyrae variables. Pulsators identified potentially as unbound with the selected potential (see Sect. 2) are color-coded by their Punb. The solid, dashed, and dotted lines denote escape velocity curves for three different gravitational potentials of the MW (see Sect. 2). We see that besides the two HV-RRL stars with the highest Punb – at the same time also probably misclassified as RR Lyrae stars (Punb = 1.00, for HV-RRL-03 and HV-RRL-04) – the rest of our high-velocity sample lies beyond 40 kpc in Galactocentric distance. The large distance of these seven RR Lyrae stars leads to potential issues of whether they can genuinely be considered unbound. Mainly due to their considerable error in vGC and RGC depicted with error ellipses in the right-hand panel of Fig. 1, and because none of them appears to be unbound due to its sizeable systemic velocity, their unbound status is uncertain. These seven pulsators might as well be ≈1−2σ outliers, especially since their unbound status is based heavily on their large tangential velocities and not on vsys. Notably, the uncertainties in tangential velocities at such distances are dominated by uncertainties in proper motions, highlighting their probably spurious unbound origin.

Moreover, we varied the value for the rotation velocity (previously set to 230 km s−1, Eilers et al. 2019) between 220 km s−1 and 240 km s−1 to assess the influence of rotation velocity on Punb. We found that some of our HV-RRL do not fulfill the condition in Eq. (11) when the rotation velocity is changed, particularly those with Punb in the range between 0.5 and 0.6. In addition, looking at the number of identified HV-RRL candidates out of the entire sample (nine out of 4847) casts a doubt on high-velocity nature specially looking at the yield of high-velocity stars in samples. For instance, Li et al. (2021) reported 591 high-velocity star candidates out of ≈10 million stars or the case of Marchetti et al. (2019), where 125 high-velocity stars were reported out of ≈7 million stars. Based on these numbers, generally speaking, we would not have expected even a single RR Lyrae star from our sample to be unbound. In addition, from the statistical point of view, based on a large data set of 4847 objects, it is likely to have three or more σ outliers from the general kinematical distribution, therefore, resulting in high Punb.

Regarding the possible origin of our HV-RRL candidates we integrated their orbits backward for 6 Gyr. During the integration, we observe each crossing of zGC = 0.0 kpc and mark the crossing radius:

r cr = x GC 2 + y GC 2 . $$ \begin{aligned} r_{\rm cr} = \sqrt{x^{2}_{\rm GC} + { y}^{2}_{\rm GC}}. \end{aligned} $$(12)

This approach helped us to identify objects that could have been ejected from the Galactic disk (hyper runaway stars) or the Galactic center (hypervelocity stars). We estimated the minimum crossing radius (rmin, calculated using the MC simulation of orbital integration) for individual HV-RRL and compared it to an arbitrary threshold marking the edge of MW disk (usually assumed 25 kpc, see, e.g., Marchetti et al. 2019; Li et al. 2021; Marchetti 2021). Our analysis has shown that none of the found candidates for high-velocity RR Lyrae stars appears to have originated in the two central MW stellar substructures. This is expected, since our sample of HV-RRL candidates is located at large distances and has highly tangential orbits.

The bottom right-hand panel of Fig. 1 depicts the spatial distribution of our HV-RRL candidates in Galactic coordinates. Here, we see that seven out of nine HV-RRL candidates are located along the Sagittarius stellar stream in the coordinate space. Three stars in the HV-RRL candidates data set (HV-RRL-02, HV-RRL-07, HV-RRL-08) had been previously linked with the Sagittarius stream based on the spatial study of the Sagittarius stream by Hernitschek et al. (2017). Based on this overlap, we decide to test whether some of the HV-RRL candidates can be associated with the Sagittarius dwarf galaxy (from here on referred to as Sgr) based on their orbit using the galpy package.

We implemented the aforementioned satellite as a Hernquist profile (Hernquist 1990) with a mass and scale radius from Table 1 in Garrow et al. (2020, and references therein). We incorporated the dynamical friction affecting the orbit of the Sgr within the MW by utilizing the Chandrasekhar prescription (Chandrasekhar 1943) implemented in galpy. To determine whether any of the HV-RRL can be associated with Sgr, we integrated their orbits backward in time for 3 Gyr utilizing the MC error simulation (1000 iterations). In each step, we varied the orbital properties of Sgr12 and the given HV-RRL, by their uncertainties by assuming they follow a Gaussian distribution. We also utilized the full covariance matrix as provided in the Gaia EDR3 data products. In addition, in each alteration, we randomly selected one of the three potentials implemented for the MW (see Sect. 2.2 for details). To assess the possible association of individual HV-RRL with the Sgr, in each MC realization, we estimated whether a given HV-RRL becomes associated with the Sgr by measuring the total distance between HV-RRLs and the Sgr. When the total distance decreased below two times the Sgr’s scale radius we considered given HV-RRL linked with the Sgr. The resulting association probabilities were estimated based on the number of instances when we linked Sgr and a given HV-RRL divided by the number of iterations. We note that each HV-RRL is represented as a mass-less particle which we follow over the past 3 Gyr.

Considering only stars that exhibited a probability of being bound in more than ten % of the instances, we found two stars that may be weakly linked with the Sgr. The stars HV-RRL-02 and HV-RRL-08 both show probabilities of 0.13, of being linked to Sgr in the past 3 Gyr, marking a possibility that they have been ejected recently in one of its pericentric approaches (see Fig. 9 in Vasiliev et al. 2021). From the metallicity point of view (listed in Table 1), two HV-RRL stars fit reasonably well to the metallicity distribution of the Sgr stream and Sgr (Hayes et al. 2020). Both of them had been previously marked as members of the Sgr stellar stream based on their spatial properties (Hernitschek et al. 2017). Therefore, two stars in our sample could be tentatively linked with the Sgr in our simplified simulation. A more sophisticated approach taking into account the disruption of Sgr, its mass-loss, and the LMC in the MW halo could constrain more realistically the link between these two RR Lyrae stars and the Sgr.

4. High-velocity tail

In this section, we focus on those RR Lyrae variables that possess a high vGC value, but remain bound to the MW. In the bottom right panel of Fig. 1, we can see that a few hundred stars from our study sample have a value of vGC >  250 km s−1 and populate the high-velocity tail of the halo RR Lyrae stars. The vGC distribution of stars that occupy a region close to the escape velocity can be described by a power-law distribution function, f, as derived by Leonard & Tremaine (1990):

f ( v GC | v esc , k ) ∝ ( v GC − v esc ) k , v cut < v GC < v esc , $$ \begin{aligned} f\left(v_{\rm GC}| v_{\rm esc}, k\right) \propto \left(v_{\rm GC} - v_{\rm esc}\right)^{k}, v_{\rm cut} < v_{\rm GC} < v_{\rm esc}, \end{aligned} $$(13)

where vesc is the escape velocity and the free parameter k is the exponent. The vcut represents the lowest threshold for vGC (usually assumed between 250 km s−1 to 300 km s−1), where the velocity distribution can still be described by Eq. (13).

Equation (13) assumes a velocity distribution that is populated up to the escape velocity and that does not contain any disk contaminants and unbound objects. The first assumption is often broken since the distribution is not always smooth due to substructures and the velocity distribution is often truncated below the true escape velocity (Smith et al. 2007; Grand et al. 2019; Koppelman & Helmi 2021). Due to this effect, the impact of underestimating the escape velocity is around 7% to 10% below the actual escape velocity (Grand et al. 2019; Koppelman & Helmi 2021).

The last assumption on the purity of the sample is somewhat easier to ensure, since unbound objects have been indentified in Sect. 3. The condition on excluding disk objects can be achieved using the Toomre diagram and probabilistic procedures described in Bensby et al. (2003, 2014) or by the Toomre selection (based on the velocity of the local standard of rest, used, e.g., in Nissen & Schuster 2010; Bonaca et al. 2017; Koppelman & Helmi 2021). In our work, we used the Toomre selection to select only halo RR Lyrae stars that fulfilled following condition:

| v − v LSR | ≥ 230 km s − 1 , $$ \begin{aligned} \left| \boldsymbol{v} - \boldsymbol{v}_{\rm LSR}\right| \ge 230\,\mathrm{km\,s}^{-1}, \end{aligned} $$(14)

where vLSR represents a velocity vector with respect to the local standard of rest vLSR = (0, 230, 0) km s−1 and v = (vx, vy, vz) is a velocity vector.

4.1. Modeling approach

We adopted the formalism prescribed in Koppelman & Helmi (2021), which takes into account the uncertainties in vGC (similarly to the study by Necib & Lin 2022a) by convolving the distribution function, f, with a Gaussian function that is represented by the measured vGC, its uncertainities, σvGC, and the true Galactocentric velocity, v GC ′ $ v^\prime _{\rm GC} $:

C ( v GC , σ v GC , θ M ) = ∫ v cut v esc f ( v GC ′ | θ M ) g ( v GC ′ | v GC , σ v GC ) d v , $$ \begin{aligned} C\left(v_{\rm GC}, \sigma _{v_{\rm GC}}, \boldsymbol{\theta }_{\rm M} \right) = \int _{v_{\rm cut}}^{v_{\rm esc}} f\left(v^\prime _{\rm GC}| \boldsymbol{\theta }_{\rm M} \right) g(v^\prime _{\rm GC} | v_{\rm GC}, \sigma _{v_{\rm GC}}) \mathrm{d}v, \end{aligned} $$(15)

where θM represents the model parameters (k, vesc, vcut), and g( v GC ′ $ v^\prime _{\rm GC} $|vGC, σvGC) is defined as:

g ( v GC ′ | v GC , σ v GC ) = 1 2 π σ v GC exp ( − ( v GC − v GC ′ ) 2 2 σ v GC 2 ) · $$ \begin{aligned} g(v^\prime _{\rm GC} | v_{\rm GC}, \sigma _{v_{\rm GC}}) = \frac{1}{\sqrt{2\pi \sigma _{v_{\rm GC}}}} \exp \left(-\frac{\left(v_{\rm GC}-v^\prime _{\rm GC} \right)^{2}}{2\sigma _{v_{\rm GC}}^{2}}\right)\cdot \end{aligned} $$(16)

The likelihood of this approach is described as:

P ( v GC , σ v GC | θ M ) = ∏ i = 1 N C ( v GC i , σ v GC i , θ M ) , $$ \begin{aligned} P\left(v_{\rm GC}, \sigma _{v_{\rm GC}} | \boldsymbol{\theta }_{\rm M} \right) = \prod _{i=1}^{N} C\left(v^{i}_{\rm GC}, \sigma ^{i}_{v_{\rm GC}}, \boldsymbol{\theta }_{\rm M} \right), \end{aligned} $$(17)

with N representing the number of stars fulfilling the condition vGC >  vcut. In addition, since the vesc varies with the Galactic radius (see, e.g., Monari et al. 2018; Koppelman & Helmi 2021) and due to the broad range over which we determine vesc, we modified our approach to take into account changes in escape velocity over radii. We parametrized vesc following the prescription by Williams et al. (2017) and Deason et al. (2019):

v esc = v esc 0 ( R / R ⊙ ) − γ / 2 , $$ \begin{aligned} v_{\rm esc} = v_{\rm esc}^{0} (R/R_{\odot })^{-\gamma /2}, \end{aligned} $$(18)

where v esc 0 $ v_{\mathrm{esc}}^{0} $ stand for escape velocity at the Solar radius. Parameter γ represents the gravitational potential slope and enters θM = (k, γ, vesc, vcut) as a model parameter. Applying the Bayesian formalism to the likelihood function above yields:

P ( θ M | v GC , σ v GC ) ∝ P ( k ) P ( γ ) P ( v esc ) P ( v GC , σ v GC | θ M ) , $$ \begin{aligned} P\left(\boldsymbol{\theta }_{\rm M} | v_{\rm GC}, \sigma _{v_{\rm GC}} \right) \propto P(k)P(\gamma )P(v_{\rm esc}) P\left(v_{\rm GC}, \sigma _{v_{\rm GC}} | \boldsymbol{\theta }_{\rm M} \right), \end{aligned} $$(19)

where P(k), P(γ), and P(vesc) represent the priors on k, γ, and vesc.

The prior assumption on the two quantities mentioned above significantly influences the resulting posterior distribution. It has been the subject of debate since the introduction of the distribution function describing the high-velocity tail. In their seminal work, Leonard & Tremaine (1990) preferred a lower range for the exponent k (between 0.5 and 2.5), which was later supported by the study of Deason et al. (2019) based on the Auriga simulation suite (Grand et al. 2017). Other studies, such as those of Smith et al. (2007), Piffl et al. (2014) or Monari et al. (2018), used a higher range for k (between 2.3 up to 4.7). Some recent studies, such as Koppelman & Helmi (2021) and Necib & Lin (2022a,b), assumed an agnostic approach on the value of k by using a large range for k (1 to 6 in Koppelman & Helmi 2021) and (1 to 15 in Necib & Lin 2022a).

In our analysis, we selected a similarly broad uniform prior on k as in studies by Koppelman & Helmi (2021) and Necib & Lin (2022a) with the following properties:

P ( k ) = 0.1 < k < 10 . $$ \begin{aligned} P(k) = 0.1 < k < 10. \end{aligned} $$(20)

The consensus on the prior for the escape velocity in the literature has seemingly converged to:

P ( v esc ) ∝ 1 / v esc , $$ \begin{aligned} P(v_{\rm esc}) \propto 1/v_{\rm esc}, \end{aligned} $$(21)

and we chose to utilize it in our study as well. We also included a range of possible values for vesc in the following form:

v cut < v esc < 1000 km s − 1 , $$ \begin{aligned} v_{\rm cut} < v_{\rm esc} < 1000\,\mathrm{km\,s}^{-1}, \end{aligned} $$(22)

where in our analysis we used vcut = 250 km s−1. For the parameter γ, we adopted simple uniform prior:

P ( γ ) = 0 < γ < 1 . $$ \begin{aligned} P(\gamma ) = 0 < \gamma < 1. \end{aligned} $$(23)

4.2. Results for the Solar neighborhood

In the following, we determined k and vesc for the region around the Solar position (4 <  RGC <  12 kpc) and with vGC above 250 km s−1 by maximizing the log-likelihood in Eq. (19). From our study sample, 234 stars13 fulfilled these two conditions. We note that we removed all nine HV-RRL from the subsequent analysis since they would tamper with the results, particularly when analyzing radii outside this local sample.

We utilized the Markov chain Monte Carlo (MCMC) sampler, implemented in emcee (Foreman-Mackey et al. 2013) to find the posterior distribution of k and vesc using 50 walkers, an initial 3000 steps as a burn-in, and restarted the sampler for an additional 2000 steps. Using the emcee implementation we estimated the maximum autocorrelation length and thinned the chains by τ = 19, which yielded 12 600 posterior samples of the k and vesc posterior distributions.

In Fig. 2, where the posterior distributions of k and vesc are displayed as corner plots, the degeneracy between both parameters is obvious. Since the distributions for both parameters have skewed distributions (non-Gaussians), we used their medians and percentiles (15.9 and 84.1) to represent the resulting values of the fit. We found v esc = 512 − 38 + 92 $ v_{\mathrm{esc}} = 512^{+92}_{-38} $ km s−1 and k = 3 . 1 − 0.8 + 1.9 $ k = 3.1^{+1.9}_{-0.8} $. The values for k and vesc agree very well with the values derived by Koppelman & Helmi (2021) and Necib & Lin (2022b). On the other hand, the comparison with studies by, e.g., Deason et al. (2019) and Monari et al. (2018) shows that our estimate on vesc lies slightly on the lower end, but still within the 16th and 84th of our percentile range.

thumbnail Fig. 2.

Corner plot for the fitted parameters k, vesc and γ, based on Eq. (19). Color contours denote the confidence intervals (CIs). The numbers above the histograms represent medians of distributions and their subscripts and superscripts stand for the 16th and 84th percentiles of the distribution.

4.3. Beyond the 12 kpc radius

The vast majority of the studies measuring the escape velocity of the MW focused on the Solar radius (from 4 kpc up to 12 kpc from the Galactic center) mainly due to the limitation on distance information (the lack thereof, or due to large uncertainties). Our data set does not have such a limitation since our derived distances are not dependent on parallaxes but, rather, on the period-metallicity-luminosity relation in the near-infrared i-passband. In addition, uncertainties on distances in our sample vary around 5% to 6%. Therefore, equipped with a sufficient number of targets, we can explore the high-velocity tail beyond the 12 kpc radius.

Applying the same criteria on vGC as in the previous subsection (using only variables with vGC >  250 km s−1) and using two ranges for RGC:

12 < R GC < 20 kpc $$ \begin{aligned} 12 < R_{\rm GC} < 20\,\mathrm{kpc} \end{aligned} $$(24)

and

20 < R GC < 28 kpc , $$ \begin{aligned} 20 < R_{\rm GC} < 28\,\mathrm{kpc}, \end{aligned} $$(25)

we obtained 102 and 97 variables, respectively, that were used to determine k, vesc, and γ beyond the Solar neighborhood. We proceeded in the same way as in Sect. 4.2, with the same set of priors and the same setup for the MCMC sampling. Figure A.2 shows results of our analysis for two regions beyond 12 kpc. We see that with a decreasing number of stars that populate the high-velocity tail, the degeneracy between the parameters k and vesc starts to dominate the posterior distribution.

In order to counter the degeneracy in both parameters, we explored the possibility of utilizing the posterior distribution of k for the Solar neighborhood in regions further away from the Sun’s position. This appears to be an appropriate path since k only weakly varies with distance (see Fig. 7 in Koppelman & Helmi 2021) based on the Aurigaia catalogues (Grand et al. 2018) generated from the Auriga suite of cosmological magneto-hydrodynamical zoom simulations (Grand et al. 2017). Instead of fixing a value for k, as in Monari et al. (2018), we used the entire posterior distribution of k determined in the previous subsection prescribed by a kernel-density estimate (KDE)14 and evaluated the probability density function at each sampled value of k.

In Fig. A.3, we show the posterior distributions for both RGC ranges obtained using the same MCMC setup as in Sect. 4.2. We see that the approximation for k using a KDE improved the resulting posterior for the vesc and in what follows, we use these estimates of escape velocity for regions beyond 12 kpc15.

5. Mass of the MW based on the RR Lyrae velocity tail

The measured escape velocity at a given radius is connected with the gravitational potential of a galaxy, Φ, and can be used to determine that galaxy’s properties, for instance, its mass. In Sect. 3, we defined unbound stars as those exceeding the MW escape velocity estimated at infinite radius (r = ∞). In what follows, we assume that the limiting radius for a star to be considered unbound is two times the virial radius r200, which is a radius where the average galaxy density is 200 times the critical density. This modification leads to a change in Eq. (10) to:

v esc ( R ⊙ ) = 2 | Φ ( R ⊙ ) − Φ ( 2 r 200 ) | . $$ \begin{aligned} v_{\rm esc}(R_{\odot }) = \sqrt{2\left| \Phi (R_{\odot }) - \Phi (2\,r_{200})\right|.} \end{aligned} $$(26)

To estimate the mass of the MW, we used MWPotential2014 (as introduced in Sect. 2.2) and we define the NFW potential by the virial mass, M200, and concentration, c. The escape velocity mainly provides insights into the mass distribution beyond the Solar radius, thus offering a second constraint for the inner mass distribution. We utilize the circular velocity, which in our case is equal to v⊙ = 230 km s−1 (Eilers et al. 2019). We tested different values for v⊙ (e.g., 220 km s−1 and 240 km s−1) and their effect on M200 and concluded that they have a small effect on the final M200, which is still within one σ of the estimated value.

To explore the range of possible parameters for our MW potential, we again employed the Bayesian formalism, as in the previous section (see Sect. 4):

P ( θ n | v obs , σ obs + , σ obs − ) ∝ P ( θ n ) P ( v obs , σ obs + , σ obs − | θ n ) . $$ \begin{aligned} P\left(\boldsymbol{\theta }_{n} | v_{\rm obs}, \sigma _{\rm obs}^{+}, \sigma _{\rm obs}^{-} \right) \propto P(\boldsymbol{\theta }_{n}) P\left(v_{\rm obs}, \sigma _{\rm obs}^{+}, \sigma _{\rm obs}^{-} | \boldsymbol{\theta }_{n} \right). \end{aligned} $$(27)

In this case, we used Eq. (27) to estimate the mass and concentration based on the calculated escape velocities (at 8, 16, and 24 kpc, which represent the centers of our radius ranges) and the circular velocity at the position of Sun. In Eq. (27), vobs stands for observed velocities (in our case escape and circular velocities) and vmod represent the model velocities derived based on a set of parameters θn = (M200, c) for the MW potential. The upper and lower uncertainties ( σ obs + $ \sigma^{+}_{\mathrm{obs}} $ and σ obs − $ \sigma^{-}_{\mathrm{obs}} $) are determined as the 15.9th and 84.1th percentiles of the posterior distributions for the escape velocities. For the circular velocity, we used a symmetric uncertainty of 10 km s−1 similarly to Necib & Lin (2022b).

For the priors, P(θn), we used the following uniform (𝒰) priors on mass and concentration:

P ( θ n ) = U ( 11.4 < log 10 ( M 200 ) < 12.5 ) ∩ U ( 1 < c < 30 ) . $$ \begin{aligned} P(\boldsymbol{\theta }_{n})\,=\;&\mathcal{U} (11.4 < \mathrm{log}_{10}(M_{200}) < 12.5)\,\cap \nonumber \\&\mathcal{U} (1 < c < 30). \end{aligned} $$(28)

Here, θn stands for model parameters. Since the posterior distributions for our calculated escape velocities are highly asymmetric (as a result of the degeneracy between k and vesc, and our set priors), we use the following prescription for log-likelihood that includes a treatment for asymmetric errors (see Sect. 6.3.1 in Barlow 2020, for more details):

P ( v obs , i , σ obs , i + , σ obs , i − | θ n ) ∝ exp ( − 1 2 ( ( v obs , i − v mod , i ( θ n ) ) ( σ obs , i + + σ obs , i − ) 2 σ obs , i + σ obs , i − + ( σ obs , i + − σ obs , i − ) ( v obs , i − v mod , i ( θ n ) ) ) 2 ) · $$ \begin{aligned}&P\left(v_{\mathrm{obs},i}, \sigma _{\mathrm{obs},i}^{+}, \sigma _{\mathrm{obs},i}^{-} | \boldsymbol{\theta }_{n} \right)\nonumber \\&\propto \exp {\left(-\frac{1}{2} \left(\frac{\left(v_{\mathrm{obs},i} - v_{\mathrm{mod},i}\left(\boldsymbol{\theta }_{n} \right) \right) \left(\sigma ^{+}_{\mathrm{obs},i} + \sigma ^{-}_{\mathrm{obs},i} \right)}{2\sigma ^{+}_{\mathrm{obs},i}\sigma ^{-}_{\mathrm{obs},i} + \left(\sigma ^{+}_{\mathrm{obs},i} - \sigma ^{-}_{\mathrm{obs},i} \right) \left(v_{\mathrm{obs},i} - v_{\mathrm{mod},i}\left(\boldsymbol{\theta }_{n} \right) \right)} \right)^2\right)}\cdot \end{aligned} $$(29)

We note that if the upper and lower uncertainties were equal, the likelihood in Eq. (29) would turn into a regular likelihood. The posterior probability of our model is then:

P ( θ n | v obs , σ obs + , σ obs − ) ∝ P ( M 200 ) P ( c ) ∏ i = 1 4 P ( v obs , i , σ obs , i + , σ obs , i − | θ n ) . $$ \begin{aligned}&P\left(\boldsymbol{\theta }_{n} | v_{\rm obs}, \sigma _{\rm obs}^{+}, \sigma _{\rm obs}^{-} \right)\nonumber \\&\propto P(M_{200}) P(c) \prod _{i=1}^{4} P\left(v_{\mathrm{obs},i}, \sigma _{\mathrm{obs},i}^{+}, \sigma _{\mathrm{obs},i}^{-} | \boldsymbol{\theta }_{n} \right). \end{aligned} $$(30)

Using the emcee module, we maximized the log-likelihood in Eq. (30). We used 150 walkers with the initial 2000 steps as a burn-in and then restarting the sampler we ran the walkers for 3000 additional steps. We thinned our chains using the autocorrelation τ = 16, and obtained 30 000 samples for the posterior distributions of concentration and mass. Figure 3 displays the resulting distributions. We find c = 11 . 3 − 2.0 + 2.3 $ c = 11.3^{+2.3}_{-2.0} $ and a mass of log 10 ( M 200 ) = 11 . 92 − 0.09 + 0.13 $ \mathrm{log}_{10}(M_{200}) = 11.92^{+0.13}_{-0.09} $, which translates to M 200 = 0 . 83 − 0.16 + 0.29 × 10 12 M ⊙ $ M_{200} = 0.83^{+0.29}_{-0.16} \times 10^{12}\,M_{\odot} $. Since our analysis covers Galactic radii up to 28 kpc we also estimated the mass within 20 kpc and 30 kpc. We found masses for the MW equal to M MW ( r < 20 kpc ) = 1 . 9 − 0.1 + 0.2 × 10 11 M ⊙ $ M_{\mathrm{MW}} \left(r < 20\,\mathrm{kpc} \right) = 1.9^{+0.2}_{-0.1} \times 10^{11}\,M_{\odot} $ and M MW ( r < 30 kpc ) = 2 . 6 − 0.2 + 0.4 × 10 11 M ⊙ $ M_{\mathrm{MW}} \left(r < 30\,\mathrm{kpc} \right) = 2.6^{+0.4}_{-0.2} \times 10^{11}\,M_{\odot} $.

thumbnail Fig. 3.

Posterior distribution for sampled parameters c and log10(M200) with color-coding representing the confidence intervals of our distribution.

Here, it is important to emphasize that this estimate is based on three estimates of the escape velocity and a single estimate of the circular velocity at the Solar position. If we consider only one single measurement of the escape velocity in the Solar neighborhood, we would obtain a similar result for the MW mass ( M 200 = 0 . 87 − 0.33 + 0.82 × 10 12 M ⊙ $ M_{200} = 0.87^{+0.82}_{-0.33} \times 10^{12}\,M_{\odot} $), only with significantly larger uncertainties. Therefore, including other regions beyond 12 kpc is valuable to increase the precision of the resulting mass estimate.

Based on results from simulations (see, e.g., Smith et al. 2007; Grand et al. 2019; Koppelman & Helmi 2021), we can suspect that our calculated escape velocities are underestimated. If we then artificially boost our escape velocity estimates by 10%, (Koppelman & Helmi 2021), we obtain a mass estimate that seems to be more in line with other mass tracers; log 10 ( M 200 ) = 12 . 10 − 0.08 + 0.12 $ \mathrm{log}_{10}(M_{200}) = 12.10^{+0.12}_{-0.08} $, which translates to M 200 = 1 . 26 − 0.22 + 0.40 × 10 12 M ⊙ $ M_{200} = 1.26^{+0.40}_{-0.22} \times 10^{12}\,M_{\odot} $. For the regions where we are actually able to measure the escape velocity (up to 30 kpc), we obtain the following masses for the MW: M MW ( r < 20 kpc ) = 2 . 1 − 0.1 + 0.2 × 10 11 M ⊙ $ M_{\mathrm{MW}} \left(r < 20\,\mathrm{kpc} \right) = 2.1^{+0.2}_{-0.1} \times 10^{11}\,M_{\odot} $ and M MW ( r < 30 kpc ) = 3 . 0 − 0.2 + 0.4 × 10 11 M ⊙ $ M_{\mathrm{MW}} \left(r < 30\,\mathrm{kpc} \right) = 3.0^{+0.4}_{-0.2} \times 10^{11}\,M_{\odot} $. Considering only the boosted vesc and circular velocity for the Sun’s position, we obtain a similar result for M 200 = 1 . 45 − 0.46 + 0.82 × 10 12 M ⊙ $ M_{200} = 1.45^{+0.82}_{-0.46} \times 10^{12}\,M_{\odot} $, which again leads to a greater uncertainty on the mass.

6. Conclusions

In this work, we utilize our data set of RR Lyrae variables with complete 7D information on their chemo-dynamical properties, described in our previous work (Prudil et al. 2021; Crestani et al. 2021; Fabrizio et al. 2021). We expanded our sample by using low-resolution spectra from the LAMOST survey, which yields a total of 6187 RR Lyrae stars with complete information on their kinematics and chemistry. After making cuts based on the quality of the astrometric solution, our study sample ended up consisting of 4847 variables.

We estimated the kinematical and spatial properties like the Galactocentric velocity and radius for each object. We compared the calculated velocity with the escape velocity determined based on a simple axisymmetric MW potential. We found nine RR Lyrae stars with Punb higher than 0.5. Two of them were likely misclassified as RR Lyrae stars, while the remaining seven objects appear to be genuine RR Lyrae stars, based on their light curve shape and Gaia colors. These variables appear to be on tangential orbits, which indicates that they do not originate in the Galactic disk and bulge.

Here it is important to emphasize that these seven HV-RRL candidates may be ≈1−2σ outliers in the vGC and RGC space, since their “unbound” status is grounded in their high tangential velocities where uncertainties in proper motions dominate the error budget at large distances. Moreover, from the pure probabilistic perspective, one would not expect to find nine unbound candidates out of 4847 of stars, in comparison with studies like Marchetti et al. (2019) and Li et al. (2021), we would expect to find less than a single unbound RR Lyrae stars. In addition, the unbound classification of some of the HV-RRL candidates can be influenced by modification of the rotation velocity, furthermore casting doubt on their unbound nature.

Seven out of nine high velocity RR Lyrae candidates follow the trail of the Sagittarius stellar stream. Three of these (HV-RRL-02, HV-RRL-07, and HV-RRL-08) were previously identified as members of the stream by Hernitschek et al. (2017) based on their spatial properties. We can tentatively link HV-RRL-02 and HV-RRL-08 to the Sagittarius dwarf spheroidal galaxy using their orbits. Since they both had already been associated with the Sagittarius stream by Hernitschek et al. (2017), they could have been ejected quite recently.

Some high-velocity candidates could be reclassified as bound if we consider a more massive MW halo. We explored this possibility by determining the MW escape velocity based on the high-velocity tail of our RR Lyrae sample from which both disk contamination and our candidate unbound RR Lyrae have been removed. Using a power-law distribution to describe the high-velocity tail (Leonard & Tremaine 1990), we expanded our escape velocity analysis beyond the Solar circle (4 <  RGC <  12 kpc) in two regions reaching up to 28 kpc. Thus, we obtained v esc = 512 − 37 + 94 $ v_{\mathrm{esc}} = 512^{+94}_{-37} $ km s−1 for the Solar neighborhood, which is lower than some of the more recent studies (see, e.g., Piffl et al. 2014; Monari et al. 2018; Deason et al. 2019), but more in line with some of the current ones such as Koppelman & Helmi (2021) and Necib & Lin (2022b). It is important to emphasize that the large positive error on our escape velocity estimate is still within the 16th and 84th percentiles with previous measurements. The measurement beyond 12 kpc yields v esc = 436 − 22 + 44 $ v_{\mathrm{esc}} = 436^{+44}_{-22} $ km s−1 and v esc = 393 − 26 + 53 $ v_{\mathrm{esc}} = 393^{+53}_{-26} $ km s−1 at 16 kpc and 24 kpc, respectively.

Together with the circular velocity at the Sun’s position, all three escape velocity estimates entered into our analysis of the MW mass, where we worked with a MW potential that included the NFW halo, for which we varied the mass and concentration. The resulting mass was equal to log 10 ( M 200 ) = 11 . 92 − 0.09 + 0.13 $ \mathrm{log}_{10}(M_{200}) = 11.92^{+0.13}_{-0.09} $ ( M 200 = 0 . 83 − 0.16 + 0.29 × 10 12 M ⊙ $ M_{200} = 0.83^{+0.29}_{-0.16} \times 10^{12}\,M_{\odot} $), which falls at the low end of MW mass estimates (e.g., Watkins et al. 2019; Callingham et al. 2019; Shen et al. 2022), but agrees reasonably well with some of the recent estimates (Eadie & Harris 2016; Eadie et al. 2018; Eilers et al. 2019; Necib & Lin 2022b). Based on escape velocity estimates from simulations, we expect the measured escape velocity to be underestimated by approximately 10%. Using this approximate correction for our escape velocities, we find the mass of the MW to be log 10 ( M 200 ) = 12 . 10 − 0.08 + 0.12 $ \mathrm{log}_{10}(M_{200}) = 12.10^{+0.12}_{-0.08} $ ( M 200 = 1 . 26 − 0.22 + 0.40 × 10 12 M ⊙ $ M_{200} = 1.26^{+0.40}_{-0.22} \times 10^{12}\,M_{\odot} $), which is in agreement with the previous mass estimates.

Mass estimates for the regions where we measure the escape velocity (up to 28 kpc) are centered around M MW ( r < 20 kpc ) = 1 . 9 − 0.1 + 0.2 × 10 11 M ⊙ $ M_{\mathrm{MW}} \left(r < 20\,\mathrm{kpc} \right) = 1.9^{+0.2}_{-0.1} \times 10^{11}\,M_{\odot} $ and M MW ( r < 30 kpc ) = 2 . 6 − 0.2 + 0.4 × 10 11 M ⊙ $ M_{\mathrm{MW}} \left(r < 30\,\mathrm{kpc} \right) = 2.6^{+0.4}_{-0.2} \times 10^{11}\,M_{\odot} $. Considering the possible 10% bias in the escape velocity, estimates for the MW mass shift to M MW ( r < 20 kpc ) = 2 . 1 − 0.1 + 0.2 × 10 11 M ⊙ $ M_{\mathrm{MW}} \left(r < 20\,\mathrm{kpc} \right) = 2.1^{+0.2}_{-0.1} \times 10^{11}\,M_{\odot} $ and M MW ( r < 30 kpc ) = 3 . 0 − 0.2 + 0.4 × 10 11 M ⊙ $ M_{\mathrm{MW}} \left(r < 30\,\mathrm{kpc} \right) = 3.0^{+0.4}_{-0.2} \times 10^{11}\,M_{\odot} $. Particularly values for MMW(r< 20 kpc) allow for comparison with other tracers where, e.g., Malhan & Ibata (2018) found MMW(r< 20 kpc) = 2.5 ± 0.2 × 1011 M⊙ (using the GD-1 stream) and Bovy et al. (2016) estimated mass of MMW(r< 20 kpc) = 1.1 ± 0.1 × 1011 M⊙ based on the Pal 5 and GD-1 stellar streams. Our value falls between both aforementioned mass estimates and agrees well with the MW mass calculated by Watkins et al. (2019) M MW ( r < 20 kpc ) = 2 . 1 − 0.3 + 0.4 × 10 11 M ⊙ $ M_{\mathrm{MW}} \left(r < 20\,\mathrm{kpc} \right) = 2.1^{+0.4}_{-0.3} \times 10^{11}\,M_{\odot} $.

An increase in the mass of the MW halo could potentially change the classification of the unbound stars HV-RRL-02 and HV-RRL-08 to bound stars. On the other hand, stars such as HV-RRL-05, HV-RRL-06, HV-RRL-07, and HV-RRL-09 would require an even higher MW mass to be considered bound (up to M200 = 1.4 × 1012 M⊙ and higher, similarly to the mass estimates from Hattori et al. 2018b; Watkins et al. 2019). HV-RRL-03 and HV-RRL-04 are likely misclassified RR Lyrae variables and, thus, they are likely to be bound. In that case, it is only HV-RRL-01 that remains as an object for which the ejection mechanism is unclear; however, at d = 53 kpc, it stands out as unbound RR Lyrae star with vGC = 655 km s−1.

Lastly with the ongoing and upcoming large spectroscopic surveys, such as 4MOST, WEAVE, and SDSS-V (de Jong et al. 2014; Dalton et al. 2014; Kollmeier et al. 2017), large populations of variable stars in the MW halo will be targeted and these could serve as tracers of the high-velocity tail and subsequently improve current estimates, especially across larger Galactic radii.

1

The LAMOST acronym stands for the Large Sky Area Multi-Object Fiber Spectroscopic Telescope, which refers to the Guoshoujing telescope.

3

For comparison, the difference between systemic velocities determined on the SDSS and LAMOST spectra for overlapping data set centers at −9 km s−1 (lower than usual one σ uncertainty of a given systemic velocity) with a dispersion of 28 km s−1.

4

We note that our astrometric cuts yield similar results as using a classifier of spurious astrometry (Rybizki et al. 2022) in 98.7% cases.

6

Infinity (INF), in this case, is described as the radius at 1012.

7

BY Draconis-type stars belong to the rotation class of dwarf variables with a spectral classification ranging from K to M with an occasional chromospheric activity in the form of starspots.

8

We note that this value lies close to the edge of applicability as defined both in Ramírez & Meléndez (2005) and González Hernández & Bonifacio (2009).

9

This label refers to a class of evolved detached binaries, β Persei type (Algol-type).

10

In our analysis we found P = 0.649351 day based on the CSS data.

11

Counts-per-day.

12

We used uncertainties on proper motions, distances and systemic velocities from McConnachie & Venn (2020) and Vasiliev & Belokurov (2020).

13

We note that we used a comparable sample size to studies such as Monari et al. (2018) and Deason et al. (2019).

14

We utilized the scipy (Virtanen et al. 2020) implementation of 1D kernel-density estimate using Gaussian kernels.

15

We note that for regions beyond the Solar radius Eq. (18) was modified to represent parametrization of vesc at larger radii, by replacing the R⊙ with values of 16 kpc and 24 kpc that better represented the two regions (see Eqs. (24) and (25)).

Acknowledgments

Z.P. acknowledges fruitful discussion with Helmer Koppelman and Andrea Kunder that significantly improved the paper. Z.P., A.J.K.H., B.L., and E.K.G. acknowledge support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID 138713538 – SFB 881 (“The Milky Way System”, subprojects A03, A05, A11). Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the US Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High Performance Computing at the University of Utah. The SDSS website is www.sdss.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics – Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU)/University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional/MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. The CSS survey is funded by the National Aeronautics and Space Administration under Grant No. NNG05GF22G issued through the Science Mission Directorate Near-Earth Objects Observations Program. The CRTS survey is supported by the US National Science Foundation under grants AST-0909182 and AST-1313422. Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope LAMOST) is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the project has been provided by the National Development and Reform Commission. LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences. This research made use of the following Python packages: Astropy (Astropy Collaboration 2013, 2018), dustmaps (Green 2018), emcee (Foreman-Mackey et al. 2013), galpy (Bovy 2015), IPython (Pérez & Granger 2007), Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), scikit-learn (Pedregosa et al. 2011), and SciPy (Virtanen et al. 2020).

References

  1. Abadi, M. G., Navarro, J. F., & Steinmetz, M. 2009, ApJ, 691, L63 [Google Scholar]
  2. Abbas, M. A., Grebel, E. K., Martin, N. F., et al. 2014, MNRAS, 441, 1230 [Google Scholar]
  3. Aguado, D. S., Ahumada, R., Almeida, A., et al. 2019, ApJS, 240, 23 [Google Scholar]
  4. Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
  5. Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
  7. Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147 [Google Scholar]
  8. Barlow, R. J. 2020, in Proceedings of the 2018 Asia–Europe–Pacific School of High-Energy Physics, 5 [Google Scholar]
  9. Bellm, E. C., Kulkarni, S. R., Graham, M. J., et al. 2019, PASP, 131, 018002 [Google Scholar]
  10. Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E., & Deason, A. J. 2018, MNRAS, 478, 611 [Google Scholar]
  11. Bennett, M., & Bovy, J. 2019, MNRAS, 482, 1417 [NASA ADS] [CrossRef] [Google Scholar]
  12. Bensby, T., Feltzing, S., & Lundström, I. 2003, A&A, 410, 527 [CrossRef] [EDP Sciences] [Google Scholar]
  13. Bensby, T., Feltzing, S., & Oey, M. S. 2014, A&A, 562, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Blaauw, A. 1961, Bull. Astron. Inst. Neth., 15, 265 [NASA ADS] [Google Scholar]
  15. Blanco-Cuaresma, S. 2019, MNRAS, 486, 2075 [Google Scholar]
  16. Blanco-Cuaresma, S., Soubiran, C., Heiter, U., & Jofré, P. 2014, A&A, 569, A111 [CrossRef] [EDP Sciences] [Google Scholar]
  17. Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529 [Google Scholar]
  18. Bonaca, A., Conroy, C., Wetzel, A., Hopkins, P. F., & Kereš, D. 2017, ApJ, 845, 101 [NASA ADS] [CrossRef] [Google Scholar]
  19. Bono, G., Iannicola, G., Braga, V. F., et al. 2019, ApJ, 870, 115 [NASA ADS] [CrossRef] [Google Scholar]
  20. Boubert, D., Erkal, D., Evans, N. W., & Izzard, R. G. 2017, MNRAS, 469, 2151 [NASA ADS] [CrossRef] [Google Scholar]
  21. Bovy, J. 2015, ApJS, 216, 29 [NASA ADS] [CrossRef] [Google Scholar]
  22. Bovy, J., Bahmanyar, A., Fritz, T. K., & Kallivayalil, N. 2016, ApJ, 833, 31 [NASA ADS] [CrossRef] [Google Scholar]
  23. Braga, V. F., Crestani, J., Fabrizio, M., et al. 2021, ApJ, 919, 85 [NASA ADS] [CrossRef] [Google Scholar]
  24. Bromley, B. C., Kenyon, S. J., Geller, M. J., et al. 2006, ApJ, 653, 1194 [NASA ADS] [CrossRef] [Google Scholar]
  25. Brown, W. R. 2015, ARA&A, 53, 15 [Google Scholar]
  26. Brown, W. R., Geller, M. J., Kenyon, S. J., & Kurtz, M. J. 2005, ApJ, 622, L33 [Google Scholar]
  27. Brown, W. R., Geller, M. J., & Kenyon, S. J. 2014, ApJ, 787, 89 [Google Scholar]
  28. Brown, W. R., Anderson, J., Gnedin, O. Y., et al. 2015, ApJ, 804, 49 [NASA ADS] [CrossRef] [Google Scholar]
  29. Callingham, T. M., Cautun, M., Deason, A. J., et al. 2019, MNRAS, 484, 5453 [Google Scholar]
  30. Capuzzo-Dolcetta, R., & Fragione, G. 2015, MNRAS, 454, 2677 [NASA ADS] [CrossRef] [Google Scholar]
  31. Casagrande, L., & VandenBerg, D. A. 2018, MNRAS, 479, L102 [NASA ADS] [CrossRef] [Google Scholar]
  32. Castelli, F., & Kurucz, R. 2003, ArXiv e-prints [arXiv:astro-ph/0405087] [Google Scholar]
  33. Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, ArXiv e-prints [arXiv:1612.05560] [Google Scholar]
  34. Chandrasekhar, S. 1943, ApJ, 97, 255 [Google Scholar]
  35. Chen, X., Wang, S., Deng, L., et al. 2020, ApJS, 249, 18 [NASA ADS] [CrossRef] [Google Scholar]
  36. Ciddor, P. E. 1996, Appl. Opt., 35, 1566 [Google Scholar]
  37. Clementini, G., Ripepi, V., Molinaro, R., et al. 2019, A&A, 622, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Crestani, J., Fabrizio, M., Braga, V. F., et al. 2021, ApJ, 908, 20 [Google Scholar]
  39. Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog: II/246 [Google Scholar]
  40. Dalton, G., Trager, S., Abrams, D. C., et al. 2014, Proc. SPIE, 9147, 91470L [Google Scholar]
  41. Deason, A. J., Fattahi, A., Belokurov, V., et al. 2019, MNRAS, 485, 3514 [Google Scholar]
  42. de Jong, R. S., Barden, S., Bellido-Tirado, O., et al. 2014, Proc. SPIE, 9147, 91470M [NASA ADS] [CrossRef] [Google Scholar]
  43. Deng, L.-C., Newberg, H. J., Liu, C., et al. 2012, Res. Astron. Astrophys., 12, 735 [Google Scholar]
  44. Drake, A. J., Djorgovski, S. G., Mahabal, A., et al. 2009, ApJ, 696, 870 [Google Scholar]
  45. Drake, A. J., Catelan, M., Djorgovski, S. G., et al. 2013a, ApJ, 763, 32 [NASA ADS] [CrossRef] [Google Scholar]
  46. Drake, A. J., Catelan, M., Djorgovski, S. G., et al. 2013b, ApJ, 765, 154 [Google Scholar]
  47. Drake, A. J., Graham, M. J., Djorgovski, S. G., et al. 2014, ApJS, 213, 9 [Google Scholar]
  48. Du, C., Li, H., Newberg, H. J., et al. 2018, ApJ, 869, L31 [NASA ADS] [CrossRef] [Google Scholar]
  49. Eadie, G. M., & Harris, W. E. 2016, ApJ, 829, 108 [NASA ADS] [CrossRef] [Google Scholar]
  50. Eadie, G., Keller, B., & Harris, W. E. 2018, ApJ, 865, 72 [NASA ADS] [CrossRef] [Google Scholar]
  51. Edelmann, H., Napiwotzki, R., Heber, U., Christlieb, N., & Reimers, D. 2005, ApJ, 634, L181 [Google Scholar]
  52. Eilers, A.-C., Hogg, D. W., Rix, H.-W., & Ness, M. K. 2019, ApJ, 871, 120 [Google Scholar]
  53. Erkal, D., Boubert, D., Gualandris, A., Evans, N. W., & Antonini, F. 2019a, MNRAS, 483, 2007 [NASA ADS] [CrossRef] [Google Scholar]
  54. Erkal, D., Belokurov, V., Laporte, C. F. P., et al. 2019b, MNRAS, 487, 2685 [Google Scholar]
  55. Evans, F. A., Renzo, M., & Rossi, E. M. 2020, MNRAS, 497, 5344 [Google Scholar]
  56. Fabricius, C., Luri, X., Arenou, F., et al. 2021, A&A, 649, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Fabrizio, M., Braga, V. F., Crestani, J., et al. 2021, ApJ, 919, 118 [NASA ADS] [CrossRef] [Google Scholar]
  58. For, B.-Q., Sneden, C., & Preston, G. W. 2011, ApJS, 197, 29 [Google Scholar]
  59. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
  60. Fragione, G., & Capuzzo-Dolcetta, R. 2016, MNRAS, 458, 2596 [NASA ADS] [CrossRef] [Google Scholar]
  61. Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Gaia Collaboration (Brown, A. G. A., et al.) 2021a, A&A, 649, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Gaia Collaboration (Smart, R. L., et al.) 2021b, A&A, 649, A6 [EDP Sciences] [Google Scholar]
  64. Garrow, T., Webb, J. J., & Bovy, J. 2020, MNRAS, 499, 804 [NASA ADS] [CrossRef] [Google Scholar]
  65. González Hernández, J. I., & Bonifacio, P. 2009, A&A, 497, 497 [Google Scholar]
  66. Grand, R. J. J., Gómez, F. A., Marinacci, F., et al. 2017, MNRAS, 467, 179 [NASA ADS] [Google Scholar]
  67. Grand, R. J. J., Helly, J., Fattahi, A., et al. 2018, MNRAS, 481, 1726 [Google Scholar]
  68. Grand, R. J. J., Deason, A. J., White, S. D. M., et al. 2019, MNRAS, 487, L72 [Google Scholar]
  69. Green, G. 2018, J. Open Sour. Softw., 3, 695 [NASA ADS] [CrossRef] [Google Scholar]
  70. Gvaramadze, V. V., Gualandris, A., & Portegies Zwart, S. 2009, MNRAS, 396, 570 [Google Scholar]
  71. Hansen, C. J., Nordström, B., Bonifacio, P., et al. 2011, A&A, 527, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Hansen, C. J., Rich, R. M., Koch, A., et al. 2016, A&A, 590, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  73. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [NASA ADS] [CrossRef] [Google Scholar]
  74. Hattori, K., Valluri, M., & Castro, N. 2018a, ApJ, 869, 33 [CrossRef] [Google Scholar]
  75. Hattori, K., Valluri, M., Bell, E. F., & Roederer, I. U. 2018b, ApJ, 866, 121 [CrossRef] [Google Scholar]
  76. Hayes, C. R., Majewski, S. R., Hasselquist, S., et al. 2020, ApJ, 889, 63 [Google Scholar]
  77. Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Nature, 563, 85 [Google Scholar]
  78. Hernitschek, N., Sesar, B., Rix, H.-W., et al. 2017, ApJ, 850, 96 [NASA ADS] [CrossRef] [Google Scholar]
  79. Hernquist, L. 1990, ApJ, 356, 359 [Google Scholar]
  80. Hills, J. G. 1988, Nature, 331, 687 [Google Scholar]
  81. Huang, Y., Li, Q., Zhang, H., et al. 2021, ApJ, 907, L42 [NASA ADS] [CrossRef] [Google Scholar]
  82. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
  83. Ibata, R. A., Malhan, K., & Martin, N. F. 2019, ApJ, 872, 152 [Google Scholar]
  84. Koch, A., Lépine, S., & Çalışkan, Ş. 2012, Eur. Phys. J. Web Conf., 19, 03002 [CrossRef] [EDP Sciences] [Google Scholar]
  85. Kollmeier, J. A., Gould, A., Knapp, G., & Beers, T. C. 2009, ApJ, 697, 1543 [NASA ADS] [CrossRef] [Google Scholar]
  86. Kollmeier, J. A., Zasowski, G., Rix, H. W., et al. 2017, ArXiv e-prints [arXiv:1711.03234] [Google Scholar]
  87. Koposov, S. E., Belokurov, V., Li, T. S., et al. 2019, MNRAS, 485, 4726 [Google Scholar]
  88. Koposov, S. E., Boubert, D., Li, T. S., et al. 2020, MNRAS, 491, 2465 [Google Scholar]
  89. Koppelman, H. H., & Helmi, A. 2021, A&A, 649, A136 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  90. Kunder, A., Rich, R. M., Hawkins, K., et al. 2015, ApJ, 808, L12 [NASA ADS] [CrossRef] [Google Scholar]
  91. Kunder, A., Pérez-Villegas, A., Rich, R. M., et al. 2020, AJ, 159, 270 [NASA ADS] [CrossRef] [Google Scholar]
  92. Lenz, P., & Breger, M. 2004, in The A-Star Puzzle, eds. J. Zverko, J. Ziznovsky, S. J. Adelman, & W. W. Weiss, IAU Symp., 224, 786 [Google Scholar]
  93. Leonard, P. J. T. 1991, AJ, 101, 562 [NASA ADS] [CrossRef] [Google Scholar]
  94. Leonard, P. J. T., & Tremaine, S. 1990, ApJ, 353, 486 [NASA ADS] [CrossRef] [Google Scholar]
  95. Lépine, S., Koch, A., Rich, R. M., & Kuijken, K. 2011, ApJ, 741, 100 [CrossRef] [Google Scholar]
  96. Li, Y.-B., Luo, A. L., Lu, Y.-J., et al. 2021, ApJS, 252, 3 [NASA ADS] [CrossRef] [Google Scholar]
  97. Lindegren, L., Klioner, S. A., Hernández, J., et al. 2021, A&A, 649, A2 [EDP Sciences] [Google Scholar]
  98. Liu, X. W., Yuan, H. B., Huo, Z. Y., et al. 2014, in Setting the Scene for Gaia and LAMOST, eds. S. Feltzing, G. Zhao, N. A. Walton, & P. Whitelock, 298, 310 [NASA ADS] [Google Scholar]
  99. Malhan, K., & Ibata, R. A. 2018, MNRAS, 477, 4063 [Google Scholar]
  100. Marchetti, T. 2021, MNRAS, 503, 1374 [NASA ADS] [CrossRef] [Google Scholar]
  101. Marchetti, T., Rossi, E. M., & Brown, A. G. A. 2019, MNRAS, 490, 157 [Google Scholar]
  102. Masci, F. J., Laher, R. R., Rusholme, B., et al. 2019, PASP, 131, 018003 [Google Scholar]
  103. McConnachie, A. W., & Venn, K. A. 2020, Res. Notes Am. Astron. Soc., 4, 229 [Google Scholar]
  104. Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533 [NASA ADS] [Google Scholar]
  105. Monari, G., Famaey, B., Carrillo, I., et al. 2018, A&A, 616, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  106. Mucciarelli, A., & Bellazzini, M. 2020, Res. Notes Am. Astron. Soc., 4, 52 [Google Scholar]
  107. Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [Google Scholar]
  108. Necib, L., & Lin, T. 2022a, ApJ, 926, 188 [NASA ADS] [CrossRef] [Google Scholar]
  109. Necib, L., & Lin, T. 2022b, ApJ, 926, 189 [Google Scholar]
  110. Nissen, P. E., & Schuster, W. J. 2010, A&A, 511, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  111. Pearson, S., Price-Whelan, A. M., & Johnston, K. V. 2017, Nat. Astron., 1, 633 [NASA ADS] [CrossRef] [Google Scholar]
  112. Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, J. Mach. Learn. Res., 12, 2825 [Google Scholar]
  113. Perets, H. B., & Šubr, L. 2012, ApJ, 751, 133 [Google Scholar]
  114. Pérez, F., & Granger, B. E. 2007, Comput. Sci. Eng., 9, 21 [Google Scholar]
  115. Piffl, T., Scannapieco, C., Binney, J., et al. 2014, A&A, 562, A91 [CrossRef] [EDP Sciences] [Google Scholar]
  116. Portegies Zwart, S. F. 2000, ApJ, 544, 437 [NASA ADS] [CrossRef] [Google Scholar]
  117. Poveda, A., Ruiz, J., & Allen, C. 1967, Boletin de los Observatorios Tonantzintla y Tacubaya, 4, 86 [Google Scholar]
  118. Price-Whelan, A. M., Mateu, C., Iorio, G., et al. 2019, AJ, 158, 223 [Google Scholar]
  119. Prudil, Z., Dékány, I., Grebel, E. K., et al. 2019, MNRAS, 487, 3270 [NASA ADS] [CrossRef] [Google Scholar]
  120. Prudil, Z., Dékány, I., Grebel, E. K., & Kunder, A. 2020, MNRAS, 492, 3408 [Google Scholar]
  121. Prudil, Z., Hanke, M., Lemasle, B., et al. 2021, A&A, 648, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  122. Ramírez, I., & Meléndez, J. 2005, ApJ, 626, 465 [Google Scholar]
  123. Rybizki, J., Green, G. M., Rix, H.-W., et al. 2022, MNRAS, 510, 2597 [NASA ADS] [CrossRef] [Google Scholar]
  124. Savino, A., Koch, A., Prudil, Z., Kunder, A., & Smolec, R. 2020, A&A, 641, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  125. Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103 [Google Scholar]
  126. Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829 [CrossRef] [Google Scholar]
  127. Sesar, B., Hernitschek, N., Mitrović, S., et al. 2017, AJ, 153, 204 [NASA ADS] [CrossRef] [Google Scholar]
  128. Shen, J., Eadie, G. M., Murray, N., et al. 2022, ApJ, 925, 1 [Google Scholar]
  129. Sherwin, B. D., Loeb, A., & O’Leary, R. M. 2008, MNRAS, 386, 1179 [NASA ADS] [CrossRef] [Google Scholar]
  130. Smith, M. C., Ruchti, G. R., Helmi, A., et al. 2007, MNRAS, 379, 755 [NASA ADS] [CrossRef] [Google Scholar]
  131. Sneden, C. A. 1973, PhD Thesis, The University of Texas at Austin, USA [Google Scholar]
  132. Sneden, C., Preston, G. W., Chadid, M., & Adamów, M. 2017, ApJ, 848, 68 [Google Scholar]
  133. Sohn, S. T., Watkins, L. L., Fardal, M. A., et al. 2018, ApJ, 862, 52 [NASA ADS] [CrossRef] [Google Scholar]
  134. Vasiliev, E., & Belokurov, V. 2020, MNRAS, 497, 4162 [Google Scholar]
  135. Vasiliev, E., Belokurov, V., & Erkal, D. 2021, MNRAS, 501, 2279 [NASA ADS] [CrossRef] [Google Scholar]
  136. Vickers, J. J., Smith, M. C., & Grebel, E. K. 2015, AJ, 150, 77 [NASA ADS] [CrossRef] [Google Scholar]
  137. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Meth., 17, 261 [Google Scholar]
  138. Watkins, L. L., Evans, N. W., & An, J. H. 2010, MNRAS, 406, 264 [Google Scholar]
  139. Watkins, L. L., van der Marel, R. P., Sohn, S. T., & Evans, N. W. 2019, ApJ, 873, 118 [NASA ADS] [CrossRef] [Google Scholar]
  140. Williams, A. A., Belokurov, V., Casey, A. R., & Evans, N. W. 2017, MNRAS, 468, 2359 [Google Scholar]
  141. York, D. G., Adelman, J., Anderson, J. E., Jr., et al. 2000, AJ, 120, 1579 [NASA ADS] [CrossRef] [Google Scholar]
  142. Yu, Q., & Tremaine, S. 2003, ApJ, 599, 1129 [Google Scholar]
  143. Zhao, G., Zhao, Y.-H., Chu, Y.-Q., Jing, Y.-P., & Deng, L.-C. 2012, Res. Astron. Astrophys., 12, 723 [CrossRef] [Google Scholar]

Appendix A: Additional figures

thumbnail Fig. A.1.

Phased light curves for the high velocity RR Lyrae stars in our sample from the CSS survey. The color-coding of the individual phased light curve represents the Punb condition of them being unbound.

thumbnail Fig. A.2.

Posterior distribution of k, vesc and γ, for 12 <  RGC <  20 kpc (left panel) and for 20 <  RGC <  28 kpc (right panel). Details are similar to those in Fig. 2.

thumbnail Fig. A.3.

The corner plots for the fitted parameters to the velocity distribution. Both panels display the posterior distribution of k, vesc and γ, with the left panel showing results for 12 <  RGC <  20 kpc and the right plot showing results for 20 <  RGC <  28 kpc. To derive these distributions, we used a kernel-density estimate (Sect. 4.2) as a prior for the k parameter.

Appendix B: Additional tables

Table B.1.

List of likely unbound RR Lyrae variables in our sample with Punb >  0.5. The first two columns represent our adopted object ID and Gaia EDR3 ID. Columns 3, 4, 5, and 6 list the heliocentric distances and systemic velocities together with their assumed uncertainties, as estimated in Prudil et al. (2021). Column 8 contains estimated Punb of the RR Lyrae stars being unbound fulfilling the condition in Eq. (13). The calculated Galactocentric velocities and their covariances are listed in columns 9, 10, 11, 12, and 13 starting with the angular velocity vy followed by the quadrature of the radial and vertical velocity, v x 2 + v z 2 $ \sqrt{v_{x}^2 + v_{z}^2} $ both accompanied by their uncertainties and correlations. The last five columns (14, 15, 16, 17, and 18) represent the total Galactocentric velocity vGC and distance RGC followed be their covariances.

Table B.2.

Table of unbound RR Lyrae candidates with Punb >  0.5 and their astrometric and photometric properties from Gaia catalog. The first column represents our object ID with an asterisk marking variables uncertain classification as RR Lyrae stars. Columns 2, 3, 4, 5, and 6 lists Gaia EDR3 ID, equatorial and Galactic coordinates. The proper motions in right ascension and declination, together with their errors, are listed in columns 7, 8, 9, and 10. Columns 11 and 12 contain parallax and its uncertainty. The last three columns (13, 14 and 15) represent re-normalized unit weight error (RUWE) and Gaia photometry in G-band and color GBP − GRP.

All Tables

Table 1.

Pulsation and chemical properties for high-velocity RR Lyrae candidates found in our analysis.

In the text
Table B.1.

List of likely unbound RR Lyrae variables in our sample with Punb >  0.5. The first two columns represent our adopted object ID and Gaia EDR3 ID. Columns 3, 4, 5, and 6 list the heliocentric distances and systemic velocities together with their assumed uncertainties, as estimated in Prudil et al. (2021). Column 8 contains estimated Punb of the RR Lyrae stars being unbound fulfilling the condition in Eq. (13). The calculated Galactocentric velocities and their covariances are listed in columns 9, 10, 11, 12, and 13 starting with the angular velocity vy followed by the quadrature of the radial and vertical velocity, v x 2 + v z 2 $ \sqrt{v_{x}^2 + v_{z}^2} $ both accompanied by their uncertainties and correlations. The last five columns (14, 15, 16, 17, and 18) represent the total Galactocentric velocity vGC and distance RGC followed be their covariances.

In the text
Table B.2.

Table of unbound RR Lyrae candidates with Punb >  0.5 and their astrometric and photometric properties from Gaia catalog. The first column represents our object ID with an asterisk marking variables uncertain classification as RR Lyrae stars. Columns 2, 3, 4, 5, and 6 lists Gaia EDR3 ID, equatorial and Galactic coordinates. The proper motions in right ascension and declination, together with their errors, are listed in columns 7, 8, 9, and 10. Columns 11 and 12 contain parallax and its uncertainty. The last three columns (13, 14 and 15) represent re-normalized unit weight error (RUWE) and Gaia photometry in G-band and color GBP − GRP.

In the text

All Figures

thumbnail Fig. 1.

Observed and estimated properties for our study sample of RR Lyrae stars (light blue dots) and the candidates for high-velocity RR Lyrae stars (HV-RRL, with individual markers color-coded on the basis of the Punb). Top-left panel: depicts the color-magnitude diagram based on the Gaia passbands, the bottom left panel: spatial distribution in Galactic coordinates of our found HV-RRL candidates. The Galactic rest-frame velocity, vGC, with respect to the Galactocentric distance, RGC, is depicted in the bottom right panel. For each point fulfilling the condition in Eq. (2), we plotted also its 1σ covariance in the form of an error ellipse using the same color-coding as for the data point. Bottom left panel: the Galactic plane is denoted by the solid black line and the red dashed line outlines the Sagittarius stream. Individual blue lines in the bottom right panel represent escape velocity curves calculated using a different gravitational potential for the MW with the solid line representing the MWPotential2014, the dashed line denoting the MWPotential2014M, and the dotted line marking the MWPotential2014H.
In the text
thumbnail Fig. 2.

Corner plot for the fitted parameters k, vesc and γ, based on Eq. (19). Color contours denote the confidence intervals (CIs). The numbers above the histograms represent medians of distributions and their subscripts and superscripts stand for the 16th and 84th percentiles of the distribution.
In the text
thumbnail Fig. 3.

Posterior distribution for sampled parameters c and log10(M200) with color-coding representing the confidence intervals of our distribution.
In the text
thumbnail Fig. A.1.

Phased light curves for the high velocity RR Lyrae stars in our sample from the CSS survey. The color-coding of the individual phased light curve represents the Punb condition of them being unbound.
In the text
thumbnail Fig. A.2.

Posterior distribution of k, vesc and γ, for 12 <  RGC <  20 kpc (left panel) and for 20 <  RGC <  28 kpc (right panel). Details are similar to those in Fig. 2.
In the text
thumbnail Fig. A.3.

The corner plots for the fitted parameters to the velocity distribution. Both panels display the posterior distribution of k, vesc and γ, with the left panel showing results for 12 <  RGC <  20 kpc and the right plot showing results for 20 <  RGC <  28 kpc. To derive these distributions, we used a kernel-density estimate (Sect. 4.2) as a prior for the k parameter.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.