Mapping Milky Way Halo Substructure Using Stars in the Extended Blue Tail of the Horizontal Branch

When dwarf galaxies fall into the Milky Way, tidal forces pull them apart and their stars become bound to the Milky Way rather than the parent dwarf galaxy. If the dwarf galaxy falls in on a circular orbit, a tidal stream of stars is created, with each star following an orbit similar to the dwarf galaxy progenitor's orbit. Stars that are stripped with lower energy per mass have shorter orbital periods and are pulled ahead, and stars with higher energy per mass lag behind. The entire dwarf galaxy can be dissolved into a long ribbon of stars that encircles the Galaxy. See Newberg & Carlin (2016) for a review of the disruption process and properties of the first couple dozen tidal streams that were discovered.

Dwarf galaxies that fall into a more massive galaxy on radial orbits also disrupt. However, instead of forming tidal streams, they catastrophically disrupt on the first pass through the galaxy center. The stars orbit the host galaxy on radial orbits, piling up at apogalacticon on umbrella-shaped "shells" before falling back through the galaxy center (Hernquist & Quinn 1988; Sanderson & Helmi 2013). Stars from the dwarf galaxy will follow orbits with a range of energies and therefore a range of periods; the galactocentric distance and number of the umbrella-shaped shells will vary with time as stars of different energy reach apogalacticon.

The recent release of Gaia DR2 (Gaia Collaboration et al. 2016, 2018) led to the discovery of the Gaia Sausage merger (Belokurov et al. 2018), the Gaia–Enceladus merger (Helmi et al. 2018), and the Virgo Radial Merger (VRM; Donlon et al. 2019), which were discovered in different Milky Way locations by different techniques. Although the inferred merger timescale differs between different authors, the similarity of the orbits of the stars in each claimed merger event suggests that they originate from the same dwarf galaxy merger event. The relationship might not be quite one-to-one, however; there is a possibility that the Gaia–Enceladus merger as discovered includes debris from both the Gaia Sausage and at least one other halo substructure combined (Belokurov et al. 2020). The Gaia Sausage (Simion et al. 2019) and VRM (Donlon et al. 2019) have both been shown to produce a significant amount of stellar debris in the region of the Virgo Overdensity (VOD; Vivas et al. 2001; Jurić et al. 2008) and the Hercules-Aquila Cloud (HAC; Belokurov et al. 2007), where the stars reach apogalacticon. This radial merger event is thought to have contributed the majority of the stars in the Milky Way inner halo (Deason et al. 2013, 2019; Belokurov et al. 2018) and has changed the way we think about searching for substructure in the outer parts of our galaxy. Recently, Donlon et al. (2020) identified shell structure in the VOD and HAC, where the stellar halo has been known to be overdense, using blue horizontal branch (BHB) and RR Lyrae stars as tracers.

This paper explores the use of the extended blue tail of the horizontal branch, including extreme horizontal branch (EHB) stars, to map structure in the Milky Way galaxy. Greenstein & Sargent (1974) used the terminology extended horizontal branch (EHB) stars, which conflicts with the current common practice of using EHB for extreme horizontal branch stars, which inhabit the extreme blue end of the blue tail. We will refer to stars on the blue extension of the horizontal branch as horizontal branch extension (HBX) stars.

Previous maps of the Milky Way's stellar halo limited the color range of BHB stars to −0.3 < g–r < 0.0 (e.g., Yanny et al. 2000; Sirko et al. 2004). This is because horizontal branch stars that are bluer than g − r = −0.3 become fainter in the g filter and form a "shoulder" where the horizontal branch turns downward toward bluer stars, as the absolute magnitude, M_g , increases (gets fainter). Deason et al. (2011) derived a polynomial fit to get a more accurate fit to the absolute magnitude of horizontal branch stars as a function of color, but they still limited the fit to −0.25 < g–r < 0.0.

The connection between the bluest horizontal branch stars (which are most easily identified as bluer than the shoulder of the horizontal branch in globular clusters) and blue subdwarfs in the field has been studied for many years. Greenstein & Sargent (1974) suggested that the blue end of horizontal branch stars extended from spectral type A to B and all the way to O and explained the origin of the blue (sdB and sdO) halo that they observed. Heber et al. (1986) showed that spectroscopically the field stars and the corresponding EHB stars in globular clusters were identical. However, a solid understanding of the mechanism through which the EHB stars are produced has remained elusive. Heber (2009) gives a review of the field of hot subdwarfs and the struggle to match them to their evolutionary origins.

More recently, the morphology of the horizontal branch has begun to be better understood. Newell & Graham (1976) showed that there were two gaps (and therefore three distinct populations) of high-latitude faint blue stars, which were expected to be part of this extended horizontal branch; the gaps were measured from Strömgren ubvy photometry to be at $\mathrm{log}{T}_{\mathrm{eff}}\sim 4.11$ and $\mathrm{log}{T}_{\mathrm{eff}}\sim 4.34$. These gaps are similar to the Grundahl jump (Grundahl et al. 1998) discovered in M13 in the Strömgren uvbyβ system and the Momany jump (Momany et al. 2002) discovered in NGC 6752 in U versus (U − V). These gaps are described below.

Grundahl et al. (1998) showed that bluer than (u − y)₀ = 0.95 the horizontal branch is 0.4 mag brighter in u than zero-age horizontal branch models would predict. The authors suggested that this would indicate two populations on the horizontal branch. The anomalous bluer horizontal branch stars could be caused by He mixing in the pre−horizontal branch giants, or variations in the CNO abundances. In a subsequent paper, Grundahl et al. (1999) show that this jump is ubiquitous in globular clusters that have a horizontal branch that extends hotter than T_eff ∼ 11,500 K. Because the jump does not depend on properties of the cluster (such as metallicity) and is also seen in the field (as a jump in surface gravity), the later paper suggests that there is no evolutionary difference between the stars redder and bluer than the Grundahl jump, but instead the stellar atmosphere models need to be adjusted to fit the data. Grundahl et al. (1999), Pace et al. (2006), and Brown et al. (2016, and references therein) explore radiative levitation of metals and gravitational settling of helium to explain the Grundahl jump.

Momany et al. (2002) showed that in addition to the Grundahl jump, there is a jump at U − V ∼ −1.0 that was attributed to post−horizontal branch evolution and/or diffusion effects. Momany et al. (2004) showed that this second jump (and possibly a third) was common in clusters with horizontal branches that extended blue enough to populate both sides of the jump. The third gap lies between the EHB stars and the blue hook stars, also known as hot He flashers (D'Cruz et al. 1996). Brown et al. (2016) analyze photometry for 53 globular clusters and find that these jumps are remarkably consistent, at the same temperature and color, between the clusters. Because these clusters have a range of stellar populations, the consistency of the jumps suggests that they are a result of changes in the stellar atmospheres of core helium fusing stars at different temperatures. Although extreme differences in He abundance were found to shift the Grundahl jump in two cases, the Momany jump remained stable for all clusters with a sufficiently hot BHB. Brown et al. (2017) showed that the stars hotter than the Momany jump had lower Fe in their atmospheres, which could also be related to atmospheric affects such as radiative acceleration.

We will show that the Grundahl jump is close to (g − r) = −0.3, the bluest color for which horizontal branch stars have been used to map the Milky Way halo. The Momany jump, located at ∼20,000 K, can be identified within the sample of stars we select, and it nominally separates stars that we will call "BHB tail" from EHB stars. Since there are very few blue hook stars compared to other stars on the HBX, their presence or absence in our data set makes little difference to our result. Note that for our purposes it is not important to understand why stars land at one temperature or another on the horizontal branch or HBX; it is only important that these stars are present in particular populations of stars and that the distance can be roughly estimated from the colors. We show that distances to HBX stars can be found to about 20%, which is considerably better than distance measurements to turnoff stars, which have also been used to trace density substructure in the halo (Newberg et al. 2002).

We were inspired to search for EHB stars as halo tracers because it was suggested by Scibelli et al. (2014) that the very blue stars that match the Sloan Digital Sky Survey (SDSS) Data Release 14 (DR14) A0p template might be EHB stars. After exploring the relationship between SDSS spectral type and color, we determined that most of the A0p stars identified by the SDSS pipeline were probably BHB tail stars. However, we recovered a larger number of EHB and BHB tail stars by including all of the blue stars that were not white dwarfs, main-sequence stars, or subdwarf stars in the sample. The analysis we present does not restrict the stellar types of the horizontal branch stars determined from template matching, but it does use the template matching to remove stars that are not horizontal branch stars.

We show how field EHB and BHB tail stars can be identified, provide a function with which to estimate the absolute magnitude (and therefore the distance) to each star, and then use these stars to explore the substructure of the Milky Way halo. We find a large number of HBX stars in the HAC that have similar orbits to those of stars in the VRM, and we observe an overdensity of halo stars at low Galactic latitude that stretches from the HAC to the VOD.

2.1. Sample Selection

2.2. Overdensity

Figure 1 shows the sky positions of our presumed HBX stars for different apparent magnitude ranges. Two-thirds of the total number of stars in the data selection are included in the 16 < g₀ < 19 magnitude range shown in the figure. The footprint of the SDSS is evident in the diagram. The majority of the stars with spectra are in the north Galactic cap. Especially for the brighter magnitudes, there is a clear increase in the star counts toward the Galactic plane. Note the apparent overdensity that is highlighted by the black box in the 18 < g₀ < 19 panel of Figure 1. This overdensity is located in the same direction as the northern end of the HAC region (R.A., decl.) = (250°, 30°) and is not present in other magnitude ranges.

Zoom In Zoom Out Reset image size

Figure 2 shows the line-of-sight, Galactic standard of rest velocity, V_GSR, of the stars in the black box of Figure 1. In this figure, the distribution of V_GSR for the stars is compared to a background Gaussian centered at 0 km s⁻¹ with a dispersion of 120 km s⁻¹, which is characteristic of generic halo stars and is normalized so that the area under the Gaussian integrates to the 57 stars in the selection box. However, the velocity distribution of stars does not look very Gaussian. It looks like there are at least one or two peaks above any Gaussian background: one is located at −90 km s⁻¹, and the other is at +30 km s⁻¹. If the V_GSR peaks are substructure, then they do not belong to the background halo, and we should not include them in the normalization of the underlying Gaussian distribution. In order for the Gaussian background plus the excess counts above the background in each of these peaks to sum to the 57 stars observed, the background Gaussian amplitude must be decreased by a factor of 0.77 (dashed line in the left panel of Figure 2). For this case, there are 12 stars where we expect 4.3 stars in the V_GSR = −90 km s⁻¹ peak (3.7σ) and 11 stars where we expect 5.5 stars in the V_GSR = +30 km s⁻¹ peak (2.3σ). This calculation is meant to give a rough number of stars expected in each "overdensity," and the exact statistics are not important. For example, in the more significant peak, we could hope to find a substructure with eight stars. The HBX stars in these overdensities are in the color range −0.53 < (g–r)₀ < −0.34. The middle and right panels of Figure 2 show the relationship between V_GSR and sky position. While the +30 km s⁻¹ stars are found spread across the selected range of right ascensions, the stars near V_GSR = −90 km s⁻¹ are possibly concentrated near (R.A., decl.) = (248°, 33°). However, a tight group of eight stars is not apparent.

Zoom In Zoom Out Reset image size

We first ask whether either of the V_GSR = −90 km s⁻¹ or +30 km s⁻¹ velocity peaks is consistent with the previously published velocity of the HAC. Belokurov et al. (2007, hereafter B07) measured a V_GSR of +180 km s⁻¹ for the northern end of the HAC, using data from a preliminary version of SDSS DR5. This velocity is not reproduced with newer data releases (Martin 2016), as we will now show. We selected stars in the HAC region, 20° < l < 75° and 20° < b < 55°, with r₀ < 19 and 0.0 < (g − r)₀ < 1.0.

Then, following the method outlined in B07, candidate giant stars were identified by selecting only stars between two M92 ridgelines (Clem 2005), shifted to 10 and 20 kpc. Though the majority of the field stars are expected to be main-sequence stars, the selection between the ridgelines favors stars at the distance of the HAC. The top panel of Figure 3 shows the original results from B07, using the same data that were used for that paper. The fact that 222 candidate HAC RGB stars have a different velocity distribution gives the impression that the selection of the HAC RGB stars was successful and that the velocity peak indicates the line-of-sight velocity of the HAC. However, we show below that the peak was just an artifact of early data pipelines.

Zoom In Zoom Out Reset image size

The middle panel of Figure 3 shows the result of the same selection, but with a newer data release: SDSS DR14. In DR14, there are many more stars present than were originally available in DR5. The newer data selection includes 5858 photometrically selected candidate giant stars, out of roughly 26,000 stars in the field. In the newer data, the V_GSR histograms for field and giant stars no longer differ from each other. The velocity of the peak at 120 km s⁻¹ is consistent with the assumption that most of the stars in both the field and candidate RGB stars are main-sequence thick-disk stars. A comparison of the velocities of individual stars in the B07 data set to the velocities of the same stars in SDSS DR14 shows that some of the values in the original data set were erroneous and resulted in the shifted location of the peak in the top panel of Figure 3. These erroneous velocities were corrected in later data analysis pipelines.

We further corroborate this point by utilizing SDSS DR14 surface gravity information, which was also not available in DR5. Imposing a surface gravity cut of 0.0 < log g_ADOP < 3.5 removed most of the main-sequence thick-disk stars from the sample; the sample was reduced from 5858 to 1519 stars. The bottom panel of Figure 3 shows the line-of-sight velocity histogram for the cleaner data set. The relatively smooth Gaussian distribution is consistent with our expectation for halo stars. Evidently, even the nearly pure sample of giant stars does not give us an indication of a group motion of the stars in the HAC.

In Donlon et al. (2019, 2020), the VRM was shown to be connected to the HAC. Donlon et al. (2020) showed that this merger generated caustic structures with stars moving back and forth through the center of the Milky Way along radial orbits. Stars with lower energy are more tightly bound to the Milky Way and oscillate with higher frequency. When the stars reach apogalacticon, they have a radial velocity of 0 km s⁻¹ with respect to the Galactic center. Closer to the Galactic center but in the same caustic structure, some of the stars will be moving toward the Galactic center and some will be moving away. The range of different velocities of stars in the same part of the Galaxy produces a complex line-of-sight velocity distribution as seen from the solar position, especially when stretched over several distance ranges. Therefore, the stars in the HAC would not be expected to have similar velocities; the velocity of each star in the overdensity is dependent on its own individual radial orbit. Viewing the velocity structure in Figure 2 through this lens, it seems plausible that we are seeing both positive and negative line-of-sight velocity peaks from a radial merger event. Similar velocity structures have been identified in the VOD and are indicative of radial mergers.

4.1. Identifying Known HBX Stars

In order to compute 6D information, we need to calibrate distances to the stars in our sample using the properties of known HBX stars in globular clusters. We used photometric data from An et al. (2008); their DAOPHOT reductions of SDSS imaging data for globular clusters allowed for more accurate crowded field photometry than was obtained with the regular SDSS pipeline.

Spectra in the vicinity of six clusters that had a blue horizontal branch that turned fainter on the blue end and were within the footprint of SDSS were collected from Smolinski et al. (2011) and Morrison et al. (2016). They used the Sloan Extension for Galactic Understanding and Exploration (SEGUE; Yanny et al. 2009) data and determined cluster membership using radial velocities and several other known properties of the globular clusters.

Figure 4 displays the color–magnitude diagrams (CMDs) for M15, M53, M13, M3, M2, and M92, with the photometric data as small gray points (10% of the data points for stars with (g − r)₀ > 0). The remaining clusters presented in An et al. (2008) either did not have a well-populated BHB shoulder or had no verified spectroscopic data. Stars with early A-type spectral classifications are shown as blue squares, one B9 star is shown as a green diamond, and the remaining stars with spectra are shown as red circles. There are many spectra of horizontal branch stars and a few spectra of HBX stars in these clusters. The spectral types of the stars on the HBX and in the blue end of the horizontal branch are all early A (A0/A0p) or very late B (B9). Also apparent is the large number of early A-type stars found in the turnoff region that are redder than our cutoff. Visual examination of these turnoff star spectra led us to conclude that these much fainter stars were misclassified in the SDSS pipeline; the spectra were noisy and did not resemble the A0 and A0p templates.

Zoom In Zoom Out Reset image size

In each cluster, our color cut of (g − r)₀ < −0.3 includes only HBX stars, some possible blue stragglers, and white dwarfs. White dwarfs are easily removed from our data set because of their characteristic spectra. The spectral types for the stars in or near the HBX in these clusters are consistent with the spectral types of candidate HBX stars found in our data set, reinforcing the idea that we have in fact identified field analogs of HBX stars. Also note that there are very few blue straggler stars that are bluer than our color cutoff.

To further connect our data set to stars of known stellar type, we identify the horizontal branch and HBX stars in M13, along with the Grundahl and Momany gaps, and compare the color distribution of actual HBX stars with the color distribution of our field stars. M13 was chosen because its HBX extends the farthest on the blue end and is the only one of the six clusters known to host EHB stars (Grundahl et al. 1999). Figure 5 shows the u₀ versus (u − g)₀ CMD. A histogram of the (u − g)₀ colors of the stars that were selected by eye to be horizontal branch stars (black points in the top panel of Figure 5) is shown by the black line in the middle panel of Figure 5. The red line in the middle panel shows the (u − g)₀ colors for the subset of the M13 horizontal branch stars with (g − r)₀ < −0.3, which is the selection criterion for our field HBX stars.

Zoom In Zoom Out Reset image size

To better interpret this figure, we will estimate the locations of the Grundahl and Momany jumps in SDSS filters. Momany et al. (2002) put the Grundahl jump at U − V ∼ −0.30 and the Momany jump at U −V ∼ −1.0. Caloi (1999) place the Grundahl jump at B −V ∼ 0.0. We estimate the Momany jump to be at B −V ∼ −0.15 from Figure 1 of Momany et al. (2002). Therefore, the Momany jump is expected to be at U −B ∼ −0.85 and the Grundahl jump is expected to be at U −B = −0.3. Using the Jester et al. (2005) conversion from U −B to (u −g) for stars with R_c −I_c < 1.15 and U −B < 0, the Momany jump should be at (u −g)₀ ∼ 0.05 and the Grundahl jump should be at (u −g)₀ ∼ 0.76. These approximate positions for the jumps are indicated as vertical dashed lines in Figure 5.

The apparent horizontal branch jumps in the top panel of Figure 5 and the corresponding gaps in the black histogram in the middle panel are slightly blueward of the estimates from our transformation. The difference is within the errors of the transformation, but it is unsettled as to whether the difference is caused by the transformation, an actual shift of the gaps of this cluster, or small irregularities in the An et al. (2008) photometry. The correspondence between the gaps and the jumps indicates that the M13 horizontal branch stars with (u −g)₀ ≲ 0 are EHB stars, while those with 0 ≲ (u −g)₀ ≲ 0.7 are what we are calling BHB tail stars between the Momany and Grundahl jumps, and those redder than (u −g)₀ > 0.7 are the normal BHB stars that are typically used to trace Galactic structure. The position of the jumps is not exact. In the bottom panel of Figure 5, there is a gap in the (u −g)₀ distribution of field HBX stars at exactly the position at which we expected to find the Momany jump, (u −g)₀ = 0.05. In addition, there are very few stars in our sample that are redder than the Grundahl jump.

For the remainder of this paper, we will assume that all of the sample stars bluer than (u −g)₀ = 0.05 are EHB stars, and sample stars redder than (u −g)₀ = 0.05 will be referred to as BHB tail stars. The two populations together are identified as HBX. Note that our sample primarily consists of EHB stars; 1088 of 1579 stars (69%) are bluer than (u −g)₀ = 0.05.

4.2. Distance Calibration

In order to estimate the distances to HBX stars, we derive the absolute magnitudes for known HBX stars as a function of (g −r)₀ color. Figure 6 shows a fourth-order polynomial fit to M13's HBX. M13 was selected because its BHB tail spans the largest range of (g −r)₀ of the six globular clusters studied. The distance to the globular cluster is 8.4 kpc (Harris 2010). When fitting the data, the mean and standard deviation were used to reject data outside of 2σ iteratively until the mean was a constant. This analysis was carried out in color bins of size 0.05. A polynomial fit, found using an ODR package (Virtanen et al. 2020), yielded

$$\begin{eqnarray*}\begin{array}{rcl}{M}_{g} & = & 17.36+172.43(g-r)+612.74{\left(g-r\right)}^{2}\\ & & +\ 855.10{\left(g-r\right)}^{3}+421.98{\left(g-r\right)}^{4}.\end{array}\end{eqnarray*}$$

This fit to the HBX in M13 is shown in the left panel of Figure 6. The average deviation for the actual stars from the line of best fit was 0.45, which gives us an estimate of the error in the derived absolute magnitudes.

Zoom In Zoom Out Reset image size

To further show that this fit approximates well the absolute magnitude of the HBX, fits were made to six other globular clusters (right panel of Figure 6). These clusters have BHB tails that overlap in (g −r)₀ with M13, allowing for the absolute magnitudes to be compared.

This polynomial fit was then used to estimate the distance to M53 and M92, using BHB tail stars with a color in the range (g −r)₀ < −0.25, the upper bound to which the HBX was fit. These are the only clusters (other than M13) with horizontal branch spectra that are near our color cut of (g −r)₀ < −0.3. In M53 there were two spectra in this cut with calculated distances of 17.2 and 19.0 kpc, for an average distance of 18.1 kpc and a standard deviation of 1.2 kpc. In Harris (2010), the distance for M53 is given as 17.9 kpc, which is within the errors of the distance measured from BHB tail stars. The distances measured for the two BHB tail stars differ from the published distance to the cluster by 4% and 6%, which is well within the ∼20% error expected if the absolute magnitude is measured to an accuracy of 0.45 mag. M92 had three stars that fall within this color cut. From the apparent magnitudes of these stars, we derive a distance to M92 of 6.7, 7.9, and 7.8 kpc for an average of 7.5 kpc with a standard deviation of 0.9 kpc. The distance to M92 is 8.3 kpc (Harris 2010), which is different from each of the three distance measurements by 19%, 5%, and 6%, respectively. For both clusters that had spectra in the color range within which the fit is valid, the distances calculated were within the expected error range of approximately 20%.

The only globular cluster in our sample in which a spectrum of an EHB star was obtained is M13. The B9 spectral type is the correct temperature for EHB stars, but the spectral classification should have been sdB (Greenstein & Sargent 1974). Although it is a bit circular, our distance relation provides an estimate of the distance of 7.6 kpc to that star. This reproduces the distance to M13 within 10%.

5.1. 6D Phase Space

Using radial velocities and angular positions from SDSS, proper motions from Gaia DR2, and distance estimates based on the calculated absolute magnitude of each candidate HBX star, full 6D phase-space information was calculated using the transforms outlined in Johnson & Soderblom (1987). We then transform the result to a right-handed, Galactocentric convention where the Sun is at X = −8 kpc and Y is in the direction of the solar motion. V_X , V_Y , and V_Z are Galactocentric velocities that are positive in the direction of positive X, Y, and Z, respectively. This is the convention that will be used throughout this paper. SDSS DR14 was combined with Gaia DR2 by matching stars within 1'' of each other in position using TOPCAT (Taylor 2005).

Figure 7 shows 3D velocities and positions for field HBXs in the identified overdensity, which is located in the sky position of the northern portion of the HAC. The diamonds represent the group of stars at −90 km s⁻¹, and green squares represent the group at +30 km s⁻¹. Since the stars in the negative-velocity group appear to be separated in velocity, the diamonds are colored to identify the higher-velocity stars (magenta) and the stars with velocities in the usual range expected for halo stars (blue). The bottom panel of Figure 7 shows energy, E, and angular momentum in the Z direction, L_z , which are integrals of motion for an axisymmetric potential. These parameters were calculated using a potential composed of a Miyamoto–Nagai disk, a Hernquist bulge, and a logarithmic halo with the parameters from Orphan Stream Model 5 in Newberg et al. (2010). The energies were renormalized by subtracting 60,000 km² s⁻² to match the energy scale in Donlon et al. (2019).

Zoom In Zoom Out Reset image size

We start by looking for halo substructure in the negative-velocity group, where we see an excess of about eight stars over a reasonable Gaussian background distribution of velocities. Looking first at the magenta diamonds in Figure 7, we see that two of them have nearly zero angular momentum. These are the leftmost two magenta points in the top panels of Figure 7, and the two lower-velocity stars in the velocity panels. The other three magenta points have wildly discrepant energy and angular momentum, which is why they do not appear on the L_z versus E panel. The measured angular momentum for these stars ranges between −4400 and −3000 kpc km s⁻¹ (negative L_z is corotating with the disk in our right-handed coordinate system) with errors of about 1000 kpc km s⁻¹, and the energy ranges from 700 to 44,000 km² s⁻² with errors of about 20,000 km² s⁻². Note that the zero-point of potential energy (and therefore the total energy) is set arbitrarily, so positive energies do not indicate stars that are unbound from the Milky Way; in the next section we will fit an orbit with an apogalacticon of 28.7 kpc to these stars. These three stars are much farther from the Galactic center and have much larger velocities than the other stars in the negative-velocity group. These three stars could represent a stream-like halo substructure.

We now focus our attention on the blue points in the negative-velocity group. Of the seven stars in this group, there are five with very similar positions, velocities, and energies. These stars have nearly zero angular momentum about the Z-axis and are located at (X, Y, Z) = (−4, 5, 6) kpc. A tight moving group like this could also represent a stream-like halo substructure. However, the distribution of angular momentum and energy for the stars in the entire HAC region suggests another possibility. From an N-body simulation of the VRM, Donlon et al. (2019) found that there is a large energy dispersion (−100,000 to −60,000 km² s⁻²) and a small range of angular momentum (roughly within ±500 kpc km s⁻¹ for the VRM). Since the VRM is believed to populate the VOD, HAC, and Eridanus-Phoenix Overdensity, it seems possible that a large fraction of the stars in the overdensity might be part of this known radial merger. This interpretation is supported by the large number of stars in the original overdensity and its coincident location with the northern part of the HAC. There are only five stars in the negative-velocity range that are tightly clustered in phase space, while the overdensity of stars in the right panel of Figure 1 contains many more stars. The HAC could be a component of a large radial merger event, which is the primary component of the stellar halo. If this is the case, we would expect the majority of the stars in the HAC to belong to this merger debris and that there would be very little smooth halo background.

This might also explain why our velocity distribution in Figure 2 does not look like it has an underlying smooth Gaussian distribution. Even if we could find a substructure of approximately eight stars that explains the negative-velocity peak, it would not explain what looks like a factor of two (or higher) stars in the HAC region than in adjacent low-latitude regions. We would need all of the stars in both peaks to be in one or more substructures to explain the observed overdensity. If there is more than one substructure, they would have to be clustered in one sky area.

We will henceforth refer to the high-velocity (magenta) group of three stars within the negative-V_GSR group as the high-energy group, which we will show is on a stream-like orbit. We will refer to the lower-velocity (blue) group of five stars in the negative-V_GSR group as the low-energy group, which we will associate with the VRM.

5.2. Stream-like and Shell-like Orbits in the Negative-velocity Group

To further explore the two negative-velocity moving groups that were identified in the previous section, we calculated orbit fits. We again used a potential composed of a Miyamoto–Nagai disk, a Hernquist bulge, and a logarithmic halo with the parameters from Orphan Stream Model 5 in Newberg et al. (2010).

We fit an orbit to the high-energy group of three stars using the procedure described in Willett et al. (2009). This method minimized the sum of the squares of the differences between the model orbit and the observed sky positions, distances from the Sun, and line-of-sight velocities of the observed stream stars. Figure 8 shows the resulting orbit, which has an eccentricity of e = 0.67, a perigalacticon of 5.6 kpc, and an apogalacticon of 28.7 kpc. The phase-space parameters for the fit orbit are (l, b, d, V_X , V_Y , V_Z ) = (2752, 3773, 25.66 kpc, 0.01 km s⁻¹, 94.51 km s⁻¹, 81.88 km s⁻¹).

Zoom In Zoom Out Reset image size

Also shown in Figure 8 is a comparison between our orbit fit to the high-energy stars and the orbit fit by Donlon et al. (2019) to the Cocytus stream (Grillmair 2009). Donlon et al. (2019) identify a group of stars in the vicinity of the VOD that were near the position of the Cocytus stream in that region of the sky and made the tentative connection. Their orbit has an eccentricity of e = 0.61, which is similar to that of our orbit fit, with a perigalacticon of 11 kpc and an apogalacticon of 44 kpc, which are somewhat different from our orbit. Both our orbit for the high-energy stars and the orbit fit in Donlon et al. (2019) are very different from the orbit that Grillmair (2009) fit to his photometric data. Grillmair (2009) identified a faint overdensity of stars that stretches from Virgo ((R.A., decl.) = (186°, −3°)) to Hercules ((R.A., decl.) = (259°, 20°)) at a distance of about 11 kpc. The stars identified by Donlon et al. (2019) are very well aligned with one end of the photometrically identified Cocytus stream, but their velocities make it impossible for the stream to stay at 11 kpc all the way to Hercules. The three high-energy stars we identify in this paper are possibly at the right sky position and distance to be associated with the other end of the Cocytus stream, but they are not on the same orbit as the stars in Donlon et al. (2019) and are not along a stream that stays at a distance of ∼11 kpc all the way to Virgo.

Grillmair (2009, 2014) finds six streams (Hermus, Hyllus, Cocytus, Lethe, Styx, and Acheron) that pass through the region of the HAC where we identified an overdensity of HBX stars. Of these, all but Styx, which was found a distance of 45 kpc from the Sun, are within 20 kpc of the Sun. Cocytus and Lethe were found at 11 and 13 kpc, respectively. None of these relatively nearby streams have well-determined orbit fits. Given the complexity of the halo in this region and the possibility that the identified streams might result from the combination of actual streams, we have deferred a study of this and related stream-like substructures to a future paper. However, we do show that our BHB tail stars have the potential to trace halo substructure.

Figure 9 shows the orbit fit to the low-energy group of five stars that have near-zero L_z . Since the phase-space coordinates for these stars are similar, the orbit parameters were determined by the average X, Y, Z, V_X , V_Y , and V_Z coordinates for the stars: (l, b, d, V_X , V_Y , V_Z ) = (5643, 4213, 8.50 kpc, −13.18, −57.58, −72.43 km s⁻¹). This orbit is highly radial (eccentricity e = 0.82) and passes very close to the Galactic center (perigalacticon of 1.0 kpc, apogalacticon of 10.7 kpc). The orbit is bent by gravity from the disk, so that consecutive apogalacticons are on the north side of the disk. The highly radial orbit is consistent with our initial guess that this moving group might be associated with the VRM. In the next section, we will use N-body simulations along this orbit to make more detailed comparisons with the VRM.

Zoom In Zoom Out Reset image size

5.3. Connecting HAC BHB Tail Stars to the VRM

An N-body simulation of a dwarf galaxy was evolved along the orbit of the low-energy HBX group with negative line-of-sight, Galactic standard of rest velocity, V_GSR. The simulated dwarf galaxy was created using parameters from Donlon et al. (2019), with the exception that the scale radius was reduced to 0.2 kpc and the evolution time was increased to match the merger time estimate from Donlon et al. (2020). The progenitor was modeled as a 10,000-body Plummer sphere with a mass of 10⁷ M_⊙, a scale radius of 0.2 kpc, and an evolution time of 2.4 Gyr. The small mass was selected by Donlon et al. (2019) to trace this group and is not representative of the dwarf galaxy progenitor, which is likely orders of magnitude more massive.

The left panel of Figure 10 shows 10% of the total bodies simulated. These results can be compared to N-body simulations with identical parameters (right panel of Figure 10) evolved along the Donlon et al. (2019) orbit for Group D, which represents a low-energy portion of the VRM in the VOD. Group D was selected for comparison since its stars have energies similar to those of our HBX stars. The simulation along the orbit of the HBX stars was run by starting a single point at the position of the HBXs, running it backward through the potential for 2.4 Gyr, placing a simulated dwarf galaxy at that position, and then running the simulation along the forward orbit for 2.4 Gyr. The simulation of Group D was run by starting a single point at the position of the Group D stars in the VOD, running it backward through the potential for 2.4 Gyr, placing an identical simulated dwarf galaxy at that position, and then running the simulation along the forward orbit for 2.4 Gyr. In R.A. and decl., the simulations for both Group D and our HBXs place the stars in similar locations. In particular, there are many stars in both simulations around (250°, −30°). The simulation in the left panel additionally puts stars in the region of the sky where we see an excess of HBXs at around (R.A., decl.) = (250°, 30°), and the right panel puts stars in the region of the VOD ((R.A., decl.) = (190°, 0°)).

Zoom In Zoom Out Reset image size

The connection between the low-energy group and Group D in the VRM is apparent from Figure 11, which shows these same N-body simulations in the X- Y plane. In this figure, the orbits of both simulations are plotted over the simulation data. The forward orbits for both are shown as blue solid lines, and the backward orbits are shown as red dashed lines. Going forward in time on the left panel by 0.3 Gyr, which is approximately the galaxy crossing time for stars on this orbit, places stars directly where the orbit in the right panel begins, at the location of the VOD. Similarly, going backward in time from the right panel by 0.3 Gyr places stars at the location where we found the negative-V_GSR group. The time difference between the two simulations results from the fact that each simulation was started from a different position on the orbit but evolved for the same amount of time. This suggests that both the negative-V_GSR group and Group D from Donlon et al. (2019) are on the same orbit; Group D is just farther along the orbit.

Zoom In Zoom Out Reset image size

We have argued that one piece of the HBX substructure in the HAC, and likely much of the HBX overdensity in the HAC, is part of the VRM. But the VRM is expected to be a significant halo substructure that populates all parts of the halo. In this section we will analyze the HBX sample as a whole. We will identify stars that are likely disk HBX stars and likely halo HBX stars and discuss the distribution of eccentricities in the stellar halo.

Figure 12 shows the 1597 candidate HBX stars, spread over the entire SDSS footprint. We calculated an orbit for each star individually using the Galpy Python package (Bovy 2015). These orbits were integrated both forward and backward in time for 0.5 Gyr. Then, we used the Galpy built-in numerical orbital eccentricity method to estimate the orbital eccentricity of each star's orbit. The bottom left panel (solid black line) shows the eccentricities for all of the stars in the sample; 68% of the stars have an eccentricity above 0.5. Our expectation is that disk stars will have an eccentricity less than 0.5 and VRM stars will have an eccentricity greater than 0.5, so we first separate the stars into these two samples. In the top left panel of Figure 12, we show the angular momentum in the Z direction as a function of energy. High-eccentricity stars are shown in magenta. The low-eccentricity stars are separated into two groups: stars that hug the low-L_z edge of the distribution for their energy (negative L_z in our right-handed coordinate system is corotating with disk stars) are assumed to be disk stars (black points), and the rest are assumed to be low-eccentricity halo stars (cyan points). There are 352 stars in the disk sample. If one looks at the region of positive L_z that is located at the symmetrical position on the other side of L_z = 0 from the disk stars, there are very few (of order 20) stars with low eccentricity in the halo. We show the distribution of eccentricities of disk stars in the bottom left panel of Figure 12. Note that the disk stars dominate the low-eccentricity stars; a very small fraction of the stars in the sample are low-eccentricity halo stars.

Zoom In Zoom Out Reset image size

In the right panel of Figure 12 we show the sky positions of the stars in the disk and halo samples. The disk stars (black points) are found predominantly at low Galactic latitudes. Note, however, that there are also a large number of high-eccentricity halo stars at low Galactic latitudes between the HAC region, (R.A., decl.) = (250°, 30°), and the VOD, (R.A., decl.) = (190°, 0°). In the top left panel one sees that the low-eccentricity halo stars are spread across the allowed region of the L_z versus E diagram, but there is a large concentration of high-eccentricity halo stars near L_z = 0 and with energies less than E = −60,000 km² s⁻². These high-eccentricity stars populating the region from the HAC to VOD are likely members of the VRM and are a large fraction of the HBX sample. From Figure 1 it is evident that the low-latitude halo structure extending south from the HAC is composed of brighter stars and is therefore closer to the Sun. These stars were not originally identified as an overdensity because an increase of stars toward the plane, particularly toward the Galactic center, is generally expected to be contamination from disk stars. Note that many VRM halo stars on high-eccentricity orbits will populate other regions of the sky. For example, the N-body simulation in the right panel of Figure 10 puts stars near (R.A., decl.) = (190°, 10°).

Figure 12 shows that we can identify likely HBX stars in the disk and the halo; the halo population appears to be predominantly related to the VRM, but a portion of the stars are unlikely to be associated with the VRM owing to their eccentricity, energy, or location in the sky. The samples identified here can be used to trace old stellar populations in the Milky Way. In future work we plan to combine these stars with other sparse tracers of these populations, such as BHB and RR Lyrae stars, to map the structure of the old stellar population in our Galaxy.

With the HBX sample now divided into Milky Way components, we can explore the morphology of the horizontal branch as a function of component. The bottom panel of Figure 5 shows the (u −g)₀ distribution for the disk, high-eccentricity halo, and low-eccentricity halo. We notice that there are very few EHB stars (bluer than (u −g)₀ = 0.05) in the low-eccentricity halo. A total of 74 of 154 low-eccentricity halo stars (48%) have EHB colors. This distribution is reasonably comparable to the distribution in M13, which is thought to have a comparatively blue horizontal branch compared to other clusters. On the other hand, the disk has a very high fraction of EHBs; 283 of the 351 presumed disk stars (81%) are EHBs. Clearly, this indicates a different population. The stars in the high-eccentricity halo, which could be associated with the VRM, also have a high fraction of presumed EHB stars, though not as high a fraction as the disk. A total of 731 of the 1074 stars in the high-eccentricity halo (68%) have EHB colors.

The EHB fraction in each component might tell us something interesting about stellar populations or help us understand the origins of EHB stars themselves. EHB stars (and field sdB stars) have measured masses of about 0.5 M_⊙ (Heber 1986); to reach the helium core burning stage with this mass, they are thought to have experienced significant mass loss, possibly through the transfer of mass during binary evolution (Han et al. 2002, 2003). The large fraction of EHB stars in the disk might tell us that the disk has a high fraction of close binaries.

Whatever the formation mechanism is for HBX and EHB stars, the fact that they follow a reasonably well-defined color–absolute magnitude relation and are associated with specific populations allows us to use them to trace halo substructure. Since the EHB stars strongly favor some stellar populations over others, they could help to separate different halo populations. If the VRM is the only component of the halo that has a large number of EHBs, then the distribution of EHB stars in the halo would trace the VRM density, for example.

In this paper we outline a method for identifying HBX stars in the field, derive a formula to estimate their absolute magnitudes from color, and use identified HBX stars to explore the 6D phase-space distribution of stars in the Milky Way. We find many likely members of the VRM, which we believe is predominantly the same merger remnant as the Gaia Sausage and Gaia–Enceladus structure. We create a set of candidate HBX stars that trace the disk, the halo, and the VRM.

We selected 1579 likely HBX stars from SDSS DR14 by choosing point sources with (g −r)₀ < −0.3, and then we eliminated white dwarfs, blue main-sequence stars, and subdwarf O stars identified by spectral template matching. This requires the availability of both photometric and spectroscopic data. We connect the HBX stars to the horizontal branch of M13 and show that the stars blueward of (u −g)₀ ≲ 0.05 are likely EHB stars.

A formula to calculate approximate HBX absolute magnitudes was derived by fitting the absolute magnitude of the HBX stars in M13 as a function of (g −r)₀ color. Other clusters approximately follow this same polynomial fit. We recover distances to globular clusters to an accuracy of about 20% by this method.

We identify an overdensity of presumed HBX stars with 18 < g₀ < 19 in the northern HAC region. The velocity distribution in this overdensity is not consistent with a smooth halo distribution and therefore indicates substructure. We identify three high-energy stars that appear to be members of a stream-like halo substructure. We traced orbits for a low-energy comoving group of five candidate HBX stars in the direction of the HAC and show that these stars are less than 0.3 Gyr behind, but are on similar orbits to, the stars in the low-energy component of the VRM. An N-body simulation evolved along the orbit of these candidate HBX stars showed similar placements of debris across the Milky Way to a simulation of VRM stars of similar energy.

A large fraction of our sample (1074 stars) are high-eccentricity halo stars. A large fraction of these stars are low-L_z stars concentrated toward the Galactic plane between the HAC and VOD and are likely members of the VRM. A small fraction of the stars (153 stars) are low-eccentricity halo stars. We also find a significant number of HBX stars (352 stars) with disk-like kinematics, which are preferentially located near the Galactic plane. The morphology of the horizontal branch in these three samples is very different. The disk HBXs have the highest fraction of EHB stars (81%), the high-eccentricity halo has the next-highest fraction (68%), and the low-eccentricity halo has the lowest EHB fraction (32%). These field HBX stars can be used to trace Milky Way substructure.

We thank the anonymous referee for encouraging us to explore the morphology of the horizontal branch. We thank V. Belokurov for providing us with the HAC data used in the Belokurov et al. (2007) paper. This work was supported by NSF grant AST 19-08653; contributions made by the Marvin Clan, Babette Josephs, and Manit Limlamai; and the 2015 Crowd Funding Campaign to Support Milky Way Research. T.D. was partially supported by a NASA/NY Space Grant fellowship.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. 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, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, 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.