The Milky Way's Shell Structure Reveals the Time of a Radial Collision

We utilize the analytical phase-space density model for a caustic structure derived in Sanderson & Helmi (2013),

$$\begin{eqnarray}&&f(r,{v}_{r})\propto \displaystyle \frac{1}{{r}_{s}^{5/2}{{\rm{\Omega }}}_{s}}\sqrt{\displaystyle \frac{\kappa }{{\delta }_{r}^{2}}}\ \exp \left\{\displaystyle \frac{-{[{r}_{s}-r-\kappa {v}_{r}^{2}]}^{2}}{2{\delta }_{r}^{2}}\right\},\end{eqnarray} \tag{ 2 }$$

where r_s is the distance of the shell from the Galactic center, κ describes the curvature of the caustic surface, δ_r is the characteristic width of the shell, and Ω_s is the solid angle spanned by the shell. Sanderson & Helmi (2013) also included a term for the radial velocity of the debris at the surface of the shell in their model, v_s . This term is directly proportional to the total angular momentum of the progenitor galaxy. The total angular momentum of the VRM debris was shown to be very small by Donlon et al. (2019), so we ignore the v_s term in our analysis. It should be noted that we expect our shells to be relatively thick due to the axisymmetric potential of the Milky Way. The curvature of the caustic surface approximately depends on the strength of the underlying potential,

$$\begin{eqnarray}&&\kappa \approx \displaystyle \frac{1}{| 2g({\boldsymbol{r}})| },\end{eqnarray} \tag{ 3 }$$

where $g({\boldsymbol{r}})$ is the gravitational force toward the Galactic center at some position ${\boldsymbol{r}}$ with respect to the center of the Galaxy.

In order to isolate shell stars, we wish to describe the density of stars with some $| {v}_{r}| \lt {v}_{c}$ in a caustic surface, where v_c is chosen to select as large and pure a sample of shell stars as possible. This selection will allow us to determine the corresponding shape of the shells in our radial density histograms. We integrate over our velocity range (−v_c  < v_r  < v_c ) to find our density function:

$$\begin{eqnarray}&&\rho (r)={\int }_{-{v}_{c}}^{{v}_{c}}f(r,{v}_{r}){{dv}}_{r}.\end{eqnarray} \tag{ 4 }$$

We do not require a complete computation of this expression, as we simply wish to determine its approximate form. Assuming that v_c is small, we approximate the integral as a rectangle of width 2v_c and height f(r,0):

$$\begin{eqnarray}&&\rho (r)\approx 2{v}_{c}\cdot f(r,0).\end{eqnarray} \tag{ 5 }$$

Evaluating this expression yields

$$\begin{eqnarray}&&\rho (r)\propto \displaystyle \frac{2{v}_{c}}{{r}_{s}^{5/2}}\sqrt{\displaystyle \frac{\kappa }{{\delta }_{r}^{2}}}\exp \left\{\displaystyle \frac{-{[{r}_{s}-r]}^{2}}{2{\delta }_{r}^{2}}\right\},\end{eqnarray} \tag{ 6 }$$

which is a Gaussian distribution, where the overall amplitude of the density peak depends linearly on our value of v_c . Thus, it is reasonable to approximate the radial density of stars in a caustic structure as a Gaussian for sufficiently small values of v_c .

Figure 2 shows a comparison of this Gaussian approximation and the original model. Note the slight difference in the radial location of the maximum of each model, and that the Gaussian distribution has a smaller maximum than the original model. This means that we have fewer stars available when looking at shells compared to entire caustic surfaces; the trade-off is that shells are substantially easier to isolate.

Zoom In Zoom Out Reset image size

The requirement for the approximation performed in Section 3.1 requires that $\kappa {v}_{r}^{2}\ll | {r}_{s}-r|$. We determine that κ = 4.74 × 10⁻⁴ kpc s² km⁻² after evaluating Equation (3) for the model potential we use for the Milky Way in this work (see Section 4.1) at a distance of 25 kpc from the Galactic center and 25/$\sqrt{2}$ kpc above the plane of the disk. This is roughly the location of the VOD. A reasonable shell width, δ_r , is on the order of 1 kpc (Sanderson & Helmi 2013), which corresponds to a typical separation of $| {r}_{s}-r| \lt$ 0.5 kpc. Solving for the requirement of the approximation, we find that $| {v}_{r}| \,\lt 30$ km s⁻¹. Due to the size of the velocity errors in our data set (${\sigma }_{{v}_{r}}=\pm 21$ km s⁻¹), we require that $| {v}_{r}| \lt 10$ km s⁻¹ in order to reduce contamination from material that is not actually in shells.

Figure 3 shows histograms of galactocentric radius, r, for the observed data, including all of the sample stars in both the VOD and HAC fields with $| {v}_{r}| \lt 10$ km s⁻¹ and $| {L}_{z}| \lt 500\,\mathrm{kpc}$ km s⁻¹. We acknowledge that our radial velocity cut likely removes many stars that are actually in shells from the data sets due to our large radial velocity errors; the entire range of velocity in the selected data is about the same as the 1σ velocity error. Cutting in L_z is less problematic, as our error in angular momentum is less than a quarter of our cut range. The candidate shell star data can be found in Tables 1 and 2.

Zoom In Zoom Out Reset image size

In Figure 3, Gaussian mixture models are fit to the cut candidate shell data using the scikit-learn sklearn.mixture.GaussianMixture implementation (Pedregosa et al. 2011) of the expectation-maximization (EM) algorithm (Dempster et al. 1977). A kernel density estimation (KDE) of the data is shown in Figure 3 to provide a smooth approximation of the underlying density of the sample that is not dependent on binning and is a better representation of what the fitting algorithm "sees." The radial density of each Gaussian is modeled as

$$\begin{eqnarray}&&f(r)=a\exp \left\{\displaystyle \frac{-{\left(r-b\right)}^{2}}{2{c}^{2}}\right\},\end{eqnarray} \tag{ 7 }$$

where the values of a, b, and c are determined by the fitting algorithm.

Since the true number of Gaussian components in our data is not clear, we use the EM algorithm to fit Gaussian mixture models with up to five components to the data. Each quality of each fit was then assessed with its associated Bayesian information criterion (BIC; Schwarz 1978) and Akaike information criterion (AIC; Akaike 1974). These information criteria help prevent over- or underfitting the data. By adding another Gaussian component to a Gaussian mixture model, one will generally find that the model is a better fit to the data. The BIC and AIC weigh the improvement in the likelihood of the fit against a penalty for adding additional parameters. The model where the information criteria are minimized is then the most statistically significant result.

Both the BIC and AIC suffer from the assumption that the number of data points is much larger than the number of parameters in the model. The small number of data points in our sample impacts the validity of the BIC and AIC values, particularly when the number of Gaussians is large. The corrected AIC (AICc; Hurvich & Tsai 1989) is an adjustment of the AIC for small numbers of data points. We find that this value is a better indicator of the quality of a model than the AIC. The BIC is less affected by the small number of data points, and we find that the BIC agrees with the AICc for models with fewer than five Gaussian components.

The EM algorithm preferentially fit Gaussians to individual stars at large distances instead of the interior structures with many stars. Finding a shell with only one star in it is not reasonable, so we omitted data points with r ≥ 30 kpc while fitting the Gaussian mixture models.

Comparing the information criteria of the data for Gaussian mixture models with varying numbers of components, we find that both the VOD and HAC data sets are best modeled by only two Gaussians. This corresponds to four statistically significant shell structures in total. Our Gaussian mixture model fit to the candidate shell star data suggests the existence of shells in the VOD region at r = 13.8 and 24.5 kpc and the HAC region at r = 8.8 and 21.6 kpc. These results show that the shells in the HAC region are closer to the center of the Galaxy than the shells in the VOD region. This supports the idea that the shells are all formed from debris of the same merger event; one expects shells to arise at different distances on opposite sides of the Galaxy, since the material in one shell must have a different energy than the material in the other shells.

In the case of a smooth halo with no substructure, we would expect to see a single Gaussian distribution. If enough cuts are made to the data, the single Gaussian distribution could be quite noisy due to small number statistics. The EM algorithm produces a best-fit single Gaussian at a distance of 15 kpc in the VOD region and 16 kpc in the HAC region. While the information criteria of the data are minimized for a two-component model, we wish to strongly rule out the possibility that the data come from a single Gaussian. To this end, we use two tests: the Anderson–Darling test (Anderson & Darling 1952; Stephens 1974) and Hartigan's dip test (Hartigan & Hartigan 1985). Both of these tests are designed to evaluate the likelihood that a given distribution is derived from a single normal distribution. The Anderson–Darling test yielded p < 0.01 for both regions, meaning that we can confidently reject the null hypothesis that either region is composed of data derived from a single Gaussian. The Hartigan's dip test produced p = 0.16 for the VOD region data and p = 0.15 for the HAC region, providing an ∼85% chance that either region is derived from more than one Gaussian distribution. Between these two statistics, we claim that the shell candidate stars make up substructure in the halo and are not simply a smooth halo background.

Figure 4 shows histograms of the candidate shell stars against the best fits from the EM algorithm. The middle row of Figure 4 shows the data after velocity cuts have been applied to select only shell stars, split into separate distributions of RRL and BHBs. Although BHBs make up a much smaller percentage of the sample in the VOD region than in the HAC region, the distributions for different types of stars in each region appear to be similarly distributed in galactocentric radius. Slight differences in the distributions of the different types of stars may be due to distance errors, which are about the size of a bin width beyond 20 kpc for BHBs and around half a bin width for RRLs. Curiously, the shell at ∼25 kpc in the VOD region appears to be composed almost entirely of RRLs. A comparison of the candidate shell stars in both regions shows that the shells in the VOD region do not lie at the same distance as the shells in the HAC region. We note that BHBs have distance errors of around 10%. At 20 kpc from the Galactic center, a 10% error can move the stars over in the histogram by as much as an entire bin. The slight differences in the distances in the peaks of the RRL and BHB distributions in the HAC region may be due to the errors in distance or simply small number statistics.

Zoom In Zoom Out Reset image size

Some stars in the data sets that were identified as RRLs also satisfied the constraints for BHBs (Section 2). We chose to take the RRL distance data for these stars, as they are probably more accurate than the corresponding BHB distance calculation. However, we want to identify how our analysis changes if these stars are actually BHBs and not RRLs. There were a total of 1994 stars in our data (16.9% of all stars) that were identified as both RRLs and BHBs. Out of these stars, 371 (3.2% of all stars) had a difference between the RRL and BHB distance calculations of greater than 0.1 kpc. For these 371 stars, the mean difference between the calculated distances from the Galactic center is 3.9 kpc, and the standard deviation in the distance differences is 3.4 kpc.

We then recalculated which of these stars are located in shells if we use the BHB distances instead of the RRL distances. In the VOD region, two stars in the 18–20 kpc bin and one star in the 14–16 kpc bin are removed. The VOD region gains one star in the 14–16 kpc bin and two stars in the 12–14 kpc bin. In the HAC region, one star is removed from the 24–26 kpc bin and two stars are removed from the 26–28 kpc bins. The HAC region gains one star in the 20–22 kpc bin and two stars in the 18–20 kpc bin. If we use the BHB distances for the stars in this work that are identified as both types of star, the interior shell of the VOD and the exterior shell of the HAC become more pronounced. Since it seems more likely that stars identified as variable are RRLs rather than BHBs, we choose to adopt the RRL distance values, even though they make the shells slightly less obvious.

Donlon et al. (2019) showed an excess of material at 40 kpc from the Sun in the VOD, and Fardal et al. (2019) showed an "Outer Virgo Overdensity" at ∼75 kpc from the Sun. The data from Vivas et al. (2016) in the bottom right panel of Figure 4 show stars with small v_r at ∼40 kpc from the Galactic center that extend out to larger distances. These structures may be outer shells formed from the VRM in the VOD region, but we are not able to analyze them in this work due to the lack of 6D information for those stars. The L_z and v_r errors are large at these distances, which inhibits our ability to properly locate distant shell substructure in observed data. These distant shells are not incompatible with our current models, as the VRM N-body simulation from Donlon et al. (2019) placed merger material at 40–70 kpc from the Sun in the VOD region, and shells are seen out to 60 kpc from the Galactic center in the simulations described later in this work. We predict that simulations of a progenitor falling in from the virial radius would also form shells at large r values from material in its trailing tidal tail.

When considering stellar velocities, it is important that we recognize that RRLs are pulsating variable stars; light from these stars will have Doppler shifts due to both translational motion and the time-dependent expansion and contraction of the surface of the star. These pulsational velocities can be as large as 120 km s⁻¹ (see Figure 2 of Sesar 2012, Figure 5 of Duffau et al. 2014, and Figure 2 of Vivas et al. 2016). The radial velocities of the RRLs in our catalog that are not from the Liu et al. (2020) catalog are derived from SDSS stack spectra. Since RRLs are typically targeted for SDSS spectra as potential BHB or quasar candidates based on color, SDSS radial velocities are not calculated with pulsating variable stars in mind. This means that our measured velocities, based only on the Doppler shift, could be wildly different than the actual velocities of the stars.

Vivas et al. (2016) created a data set of RRLs in the VOD in a way that cuts out stars with wildly erroneous radial velocity measurements. In this data set, the majority of RRLs had radial velocity uncertainties below 10 km s⁻¹. We matched this data set with Gaia data in order to obtain 3D velocities and then cut it so that $| {v}_{r}| \lt 10$ km s⁻¹ and $| {L}_{z}| \lt 500\,\mathrm{kpc}$ km s⁻¹ in order to isolate shell structure. This is compared with our RRL data in the VOD region in the bottom left panel of Figure 4. In this figure, the Vivas et al. (2016) data imply the existence of shells in the VOD region at the same distances as the two statistically significant shells that were identified in RRLs and BHBs.

One might be concerned that large radial velocity errors would confuse the selection of shell stars using galactocentric radial velocities. However, the velocity errors just from measurement are already quite large, of the order of the random velocity errors due to pulsation of the RRLs. Many stars that are actually in the shell are not selected due to random error or RRL pulsation, reducing our shell signal. Stars that are actually up to 120 km s⁻¹ away from v_r  = 0 could be erroneously included in the sample. However, on average, the radial velocity measured from stack spectra will not be 120 km s⁻¹ off from the actual radial velocity measurement, since the pulsational velocity is only that large for a small portion of the RRL oscillation period, and stack spectra usually include measurements from several different epochs. In practice, the maximum difference between a stack spectrum with multiple measurement epochs and the actual radial velocity is closer to 50 km s⁻¹, and the actual error is usually even smaller (see Figure 3 of Vivas et al. 2016). The right panel of Figure 2 shows that within the same caustic structure, stars that actually have galactocentric radial velocities of v_r  = ±50 km s⁻¹ would only be ∼2 kpc closer than the actual shell structure, which is within the systematic distance error for stars in our data set beyond 20 kpc from the Sun.

The radial velocities of RRLs from the Liu et al. (2020) catalog are calculated by fitting the observed velocity as a function of phase to an empirical pulsating velocity template curve. It is expected that this will reduce the error in radial velocity due to pulsational velocity compared to calculations of radial velocity from stack spectra. In order to test the differences between radial velocities that are calculated using stack spectra and velocities that are calculated by fitting an observed velocity to an empirical velocity curve, we looked at the sample of 369 RRLs that had measured spectral radial velocities in both the Gaia RRL and Liu et al. (2020) catalogs. The radial velocities of both catalogs are generally consistent; the mean difference between the two catalogs' radial velocities is 21 km s⁻¹, and the standard deviation of the differences is 15 km s⁻¹. The mean difference between the catalogs' radial velocities is roughly the same size as the upper limit for the estimated error in the radial velocity from either catalog. This suggests that a sizable fraction of the differences between the two catalogs' radial velocities is due to measurement uncertainties and not necessarily to pulsational velocities at the surface of a star. While radial velocities of RRLs derived from stack spectra are likely less accurate than radial velocities derived from empirical velocity curves, the typical departure of the radial velocity value derived from stack spectra from the "correct" value is in practice much smaller than the 120 km s⁻¹ upper limit for the magnitude of the pulsational velocity. We maintain that radial velocities derived from stack spectra for RRLs are fairly good estimates for the actual radial velocities of those stars, although they should be used with caution.

Especially in the VOD and HAC regions of the sky, RRLs that are not associated with any VRM caustic structure and also have low L_z angular momentum are relatively uncommon. This explains why we are able to recover the shells even though many RRLs could have large radial velocity errors. The similarity between the BHB and RRL data and the positions of shells from RRLs in Vivas et al. (2016) data in the VOD support the conclusion that the RRL radial velocity errors that are derived from SDSS spectra do not significantly impact the outcomes of this work.

It should be noted that using the radial velocity values from the Liu et al. (2020) catalog in place of the radial velocity values from the Gaia/SDSS catalog did not substantially change the histogram of stars (Figure 4) that satisfy the galactocentric radial velocity and angular momentum cuts. This is consistent with our analysis that the radial velocities of RRLs derived from stack spectra that we use in this work are typically not very far off from the "correct" radial velocity measurements, and that the errors in the measured radial velocity values are often just as large as the difference between the stack spectra and velocity curve fit radial velocity measurements. This may not be the case for all data sets, as the quality of the radial velocities derived from stack spectra depends on several factors, such as the number of observations, the quality of the spectra, and the times at which the observations are made.

In order to explore phase mixing in radial mergers, we create a suite of 120 N-body simulations of radial mergers. These simulations are all initialized with a single Plummer profile progenitor dwarf galaxy (Plummer 1911). A third of these simulations use a dynamical mass of 10⁹ M_⊙ and a scale radius of 3 kpc, which is consistent with measurements of typical dwarf galaxies such as the Fornax dSph (Penarrubia et al. 2008). These parameters were previously used as likely parameters for the progenitor of the VRM by Donlon et al. (2019), who claimed that an N-body simulation with this mass and scale radius more accurately recovered the kinematics of the VRM than simulations with a mass of 10⁷ M_⊙ and a scale radius of 0.4 kpc. In order to explore the effect of progenitor mass on radial merger phase mixing, we also perform simulations with a dynamical mass of 10⁸ M_⊙ and a scale radius of 1 kpc, as well as simulations with a dynamical mass of 10¹⁰ M_⊙ and a scale radius of 10 kpc. The progenitors in these simulations are given zero initial velocity, so that the initial energy of the dwarf galaxy is determined solely by its initial distance from the Galactic center, and the merger is guaranteed to be nearly radial.

Since we are starting the dwarf galaxy at rest in an axisymmetric potential, there are only two parameters required to create the full range of substructure for a particular dwarf model. These are the initial inclination angle (i) from the Galactic disk as viewed from the Galactic center and the initial distance from the Galactic center (r₀). The inclination angle is defined so that i = 0° in the disk and i = 90° when the dwarf galaxy is positioned along the Galactic Z-axis. The values we explored for i range from 0° to 90° in increments of 10°, and the values we explored for r₀ are 20, 30, 45, and 60 kpc. The observed VRM debris has an L_z distribution centered on zero (Donlon et al. 2019), so we did not run simulations with nonzero angular momenta. We used a value of Y = 0 for the initial position of the dwarf galaxy, as the azimuthal angle of impact does not change the simulation due to the axisymmetric symmetry of the Galaxy.

The simulations were run on the MilkyWay@home N-body software (Shelton 2018) in a static gravitational potential consisting of a Hernquist bulge, Miyamoto–Nagai disk, and logarithmic halo using potential parameters from Orphan Stream Model 5 in Newberg et al. (2010). Each simulation contained 20,000 bodies and was integrated forward in time for 10 Gyr. For the rest of this work, we will refer to each simulation solely based on its progenitor's mass, initial inclination angle, and initial distance from the Galactic center. When necessary, the results of the simulation are also identified by a time step.

Note that it is likely that the progenitor of the VRM was originally located at or near the Milky Way's virial radius and then fell into the Milky Way along a decaying orbit. As the orbit decayed, each apogalacticon of the dwarf galaxy's orbit would be located closer to the Galactic center than the previous apogalacticon. Eventually, the dwarf galaxy collided with the Milky Way; this occurred after the final apogalacticon on the dwarf's orbit, which determines the final orbital energy of the progenitor dwarf galaxy. After the dwarf galaxy passed through the Milky Way, its constituent stars became bound solely to the Milky Way, and their individual energies were locked in. We need only to constrain the distance of the final apogalacticon pass of the progenitor of the VRM in order to characterize the present-day kinematics of the structure, which corresponds to r₀.

Simulated radial mergers experience phase mixing of the dwarf galaxy debris over large timescales. Figure 5 shows how small variations in the initial energies of stars in a radial merger will cause the material in shells to segregate based on energy. This phase-wrapping causes more shells to form as the age of the radial merger increases. Phase mixing is also caused by slight tangential velocities of stars in each shell due to variations in angular momentum. This causes shells to grow in surface area over time. Eventually, debris from a merger reaches a phase-mixed equilibrium within the Galaxy.

Zoom In Zoom Out Reset image size

Figure 6 shows the simulation of a radial merger with initial inclination angle i = 30°, r₀ = 30 kpc, and a mass of 10⁹ M_⊙. The data for this simulation were cut to mimic our data shown in Figure 4. This simulation shows that over time, the strong peaks demonstrating shells phase mix into a smooth distribution. When a merger has phase mixed completely, shells will no longer be present (that is, the number of caustic surfaces eventually approaches infinity, and it is impossible to detect any one shell in particular).

Zoom In Zoom Out Reset image size

We seek an objective way to describe the amount that a radial merger has phase mixed from its initial state. To this end, we develop a numerical metric for a histogram h, which we call the "causticality" C of that histogram. With n bins in h and N_i counts in the ith bin of h, we define causticality as

$$\begin{eqnarray}&&C=\displaystyle \frac{{\displaystyle \sum }_{i=2}^{n}{\left({N}_{i}-{N}_{i-1}\right)}^{2}}{{\displaystyle \sum }_{i=2}^{n}{\left({N}_{i}+{N}_{i-1}\right)}^{2}}.\end{eqnarray} \tag{ 8 }$$

Causticality consists of a ratio of the sum of the squared differences between each bin and its previous neighbor and the sum of the squared sums of each bin and its previous neighbor. The denominator in Equation (8) imposes that C = C(m · h) for any constant m > 0, which ensures that the causticality only depends on the shape of the input histogram and not the total number of data points in the histogram.

Causticality can take on all values between zero (completely mixed) and unity (completely unmixed). Consider a distribution with particles in only one bin or that has undergone no phase mixing; the measured causticality of this system would be equal to unity. This is also the case for a system where every other bin is empty. On the other end, a uniform distribution has a causticality of zero. Here C is larger for distributions with sharp peaks than distributions with gradual, smooth peaks; Figure 6 shows how phase mixing lowers the measured value of the causticality over time in a simulation of a radial merger.

In practice, the measured causticality of our systems never actually becomes zero. This is because the debris of radial mergers relaxes into a smooth Gaussian distribution, not a uniform distribution. Additionally, the invariance of causticality with respect to the number of particles assumes that a constant bin size is being used. Typically, bin size decreases as more objects are added to a histogram, which can alter the measured value of causticality. Thus, it is important that the bin sizes for the observed and simulated data are identical.

If observable shells are present in our data, we would expect to see a few strong, sharp peaks. As radial mergers age, these few thin shells become wider and more numerous until the overall radial distribution approaches a smooth distribution. We therefore expect the causticality to decrease over time for each simulation of a radial merger until the merger is completely phase mixed, at which point the causticality assumes some small constant value.

Appendix A contains an analysis of the uncertainty in causticality, as well as the effects of bin size on measured causticality values. We determine that 2 kpc bin sizes are best for these data.

This is the first time that causticality has been introduced, and we recognize that it is beneficial to perform a similar analysis of phase mixing using a widely studied statistic. In Appendix B, we use the Kullback–Leibler divergence (KLD; Kullback & Leibler 1951) instead of causticality to evaluate how phase mixed a distribution is. We find that it is more difficult to get an estimate of the merger time of the VRM from measurements of the KLD, but we are still able to determine that the VRM must have occurred within the last 5 Gyr. One advantage of causticality over the KLD is that causticality does not require comparison to some assumed baseline distribution; while the KLD as it is applied in this work only approaches zero when the distribution is uniform, the causticality rapidly approaches zero as a distribution becomes smooth, regardless of its shape.

Figure 7 shows plots of the causticality over time for all 120 radial merger simulations with varying inclination angle, radius of initial infall, and mass. Evaluating the causticality for our observed data, we find that the value of the causticality is 0.20 ± 0.09 for the VOD region and 0.12 ± 0.06 for the HAC region. Any background halo stars or observational errors in these data sets will decrease the measured values of the causticality in the observed data, so these calculated values of the causticality of the observed data are lower limits. Background stars or observational errors will make the observed data appear to be older than they actually are, since the simulated data do not suffer from these issues.

Zoom In Zoom Out Reset image size

Figure 7 shows that simulations with larger initial distances take longer to phase mix, on average, than simulations with smaller initial distances. This means that the approximate time it took the VRM debris to phase mix to its current state depends on what the distance of the final apogalacticon (r₀) of the VRM progenitor's orbit was. In order to get a good estimate of the age of the VRM from this data, we need to determine a likely value of r₀ for the progenitor of the VRM.

Assume that the progenitor dwarf galaxy remains largely intact prior to reaching the final apogalacticon. The individual stars in the dwarf galaxy have a range of energies that are limited by the fact that they are gravitationally bound together. When the dwarf galaxy is tidally disrupted as it passes through the Galactic center, the individual stars have a range of energies near the energy per mass of the original dwarf galaxy and will therefore populate orbits with a range of apogalacticons that are near the final apogalacticon of the dwarf galaxy. For example, Figure 6 shows that the debris of a simulated radial merger with r₀ = 30 kpc primarily populates the halo at distances between 20 and 40 kpc.

If the radial merger event has enough material to dominate the stellar halo, which is true for the VRM, then we expect to see a stellar break in the halo in the range of distances that the radial merger populates. Various works place the stellar breaks at a range of 18–28 kpc (Watkins et al. 2009; Deason et al. 2011; Sesar et al. 2011; Pila-Díez et al. 2015; Xue et al. 2015; Deason et al. 2018). This range in the stellar break radius may be due to measurements in different regions of the sky; certain portions of the Galaxy, such as the VOD and the HAC, will have proportionally more VRM debris than regions without visible overdensities. Measurements of the stellar break radius in the overdense regions will produce values similar to r₀ for the VRM progenitor; measurements in other areas of the sky will produce values of the stellar break that are smaller than r₀ for the VRM progenitor due to the majority of VRM material in that region not being near apogalacticon and therefore populating regions with smaller galactocentric radii. A measurement of the stellar break over the entire halo would then be an average of the actual r₀ value for the VRM progenitor and other regions of the sky and expected to yield a value somewhat smaller than the actual value of r₀ for the VRM progenitor.

These values for the stellar break suggest that the value of r₀ for the VRM must be between 20 and 30 kpc in order to generate a stellar break near those distances. This value of r₀ for the VRM is supported by Villalobos & Helmi (2008), who showed that simulations of dwarf galaxies falling into the Milky Way from distances near the virial radius tend to have a final apogalacticon around r₀ ∼ 20 kpc from the Galactic center before colliding with the Galaxy.

In Figure 7, we find that simulations with r₀ = 20 kpc have similar causticalities to the observed data between 1 and 2 Gyr after the simulation is started. The simulations with r₀ = 30 kpc have causticalities that are similar to the observed data between 3 and 5 Gyr. Since we expect the progenitor of the VRM to have a value of r₀ between 20 and 30 kpc, the time at which the VRM has values of causticality that are similar to the observed data is expected to be within the last 5 Gyr.

This is again consistent with Villalobos & Helmi (2008), who stated that shell structure is present for only ∼2 Gyr after dwarf galaxies collide with the Milky Way in their simulations. Further, we were unable to locate identifiable shell structure beyond 5 Gyr after collision in our simulations with r₀ < 60 kpc. These points suggest that after 5 Gyr, VRM-like mergers have phase mixed beyond what we see in the observed data, and that the VRM cannot have happened that long ago, as we still see shell substructure.

Note that as the mass of the progenitor dwarf galaxy increases, less time is required to phase mix the debris to observed levels. For progenitors with a mass of 10¹⁰ M_⊙, this mixing time is very consistent and does not appear to strongly depend on initial inclination angle or distance. No simulations with this mass had a value of causticality similar to the observed data after 5 Gyr.

For simulations with masses of 10⁹ and 10¹⁰ M_⊙, the time to phase mix increases monotonically with initial distance. However, simulations with a mass of 10⁸ M_⊙ and r₀ = 60 kpc appeared to phase mix more quickly than simulations with identical mass and smaller initial distances. This is due to the distance cut that we made on our simulations when calculating causticality. Since the observed data did not extend beyond 50 kpc, we cut the simulated data to only include objects within 50 kpc. The 10⁸ M_⊙, r₀ = 60 kpc simulations place the majority of their material beyond 50 kpc from the Galactic center, so the measured causticalities only include the small number of objects that remain within 50 kpc. This deflates the measured value of causticality. However, a structure that places almost all of its debris beyond 50 kpc cannot possibly be responsible for shells within 30 kpc of the Galactic center, so simulations with 10⁸ M_⊙ and r₀ = 60 kpc are not viable models for the VRM. Data from the 10⁸ M_⊙ and r₀ = 60 kpc simulations are omitted from Figure 7 to avoid potentially misleading the reader.

Stars in a shell will possess similar energies and will have moved similarly throughout the lifetime of the radial merger. Thus, we can expect a single star in some shell to oscillate with the same period as the rest of the stars in the shell. This means that we can model caustic oscillation in the Galactic potential as a free oscillation of a single model mass at the location of a shell. At any given point in time, we expect to see multiple shells at different radii and positions on the sky for a given radial merger. By assigning model particles with zero initial velocity to each shell and then tracing these particles' motions back in time, we can compare the relative positions of each caustic surface throughout time. At some time in the past, the particles on all of the shells must have been at the same place, bound to the dwarf galaxy. The difference in each shell's distance from one another can then be used to pinpoint the time that a radial merger was still self-bound.

We begin by placing a model particle at each shell found by the BIC fitting algorithm for some radial merger (either simulated or observed data). We must be able to allow motion from one side of the Galaxy to the other, so we arbitrarily select the shells on one side of the Galaxy to be located at positive r and the shells on the other side to be located at negative r. Particles will then oscillate between positive and negative r values. Using Galpy version 1.4.1 (Bovy 2015; http://github.com/jobovy/galpy), we calculate the forces on each model particle at each time step in order to integrate the particles backward in time for 10 Gyr. By integrating each model particle back in time, we expect to find a point in time where each caustic structure is located at the same point in space. Additionally, if each caustic structure is moving in the same direction (has the same signs of dr/dt) at that point, we claim that we have found the point in time where the progenitor of the radial merger was still self-bound. This point in time will be the "merger time" of the radial merger.

In order to measure the overall difference in the location of shells, we introduce δ(t) as a metric to measure the total distance between model particles,

$$\begin{eqnarray}&&\delta (t)=\displaystyle \sum _{i\lt j}^{N}| {{\boldsymbol{r}}}_{i}(t)-{{\boldsymbol{r}}}_{j}(t)| .\end{eqnarray} \tag{ 9 }$$

Here we use vector subtraction in order to ensure that the model particles are spatially close to one another and not simply at similar galactocentric distances. A small δ(t) means the caustic structures are clumped together, while a high δ(t) corresponds to a large separation of the shells. Times with a small δ(t) are more likely to be times when the progenitor was still bound.

As the progenitor of the radial merger collides with the host galaxy, all of the caustic structures are initially falling inward and then transition to moving outward from the host galaxy. We count the number of first derivatives of the model particles that are positive (dr/dt > 0) in order to determine the number of caustic structures moving in the same direction. The passage of the dwarf galaxy through the Galactic center would cause the number of positive first derivatives to go directly from zero to the total number of shells. If a minimum of δ(t) corresponds to a time where the number of model particles with dr/dt > 0 is either zero or the number of shells in the data, then that is a likely time that the dwarf galaxy passed through the Galactic center and was tidally disrupted. This is especially true if the number of model particles with dr/dt > 0 goes directly from zero to the number of shells or vice versa at that time, as this corresponds to the progenitor passing through the Galaxy.

We tested the method outlined in Sections 5.1 and 5.2 by using it to measure merger times for the series of simulations of radial mergers in a Milky Way–like galaxy described in Section 4.1. These simulations had a known evolve time, which allowed us to verify whether or not the method was able to recover the correct time of collision for each different simulation.

We created data sets similar to the ones used for the analysis of the VRM debris by cutting our simulation data in v_r and L_z identically to the way they were cut in the observed data. Additionally, we only looked at data out to r = 50 kpc, in order to ensure that our method could recover the correct infall time despite our limitations in distance in the observed data sets. The majority (>75%) of the simulation data were typically found in this distance range.

Next, we cut the simulated data based on sky position. Our observed data were cut to only two regions of the sky, the VOD and HAC regions. In order to emulate this, we looked at the radial merger simulations in R.A. and decl. and selected the data from two regions that were chosen to be similar in size to the VOD and HAC regions, have roughly the same positions in the sky (within ∼20°), and contain an overdensity of the simulated bodies. The regions had to be chosen by hand, since the simulations did not all place shell structure in the same places on the sky. We had to cut out regions around the shell overdensities because the Galactic potential was not spherically symmetric, so the shells are not at a constant radius; by selecting a limited region of the sky, the shells were at an approximately constant radius over the selected data. We made these cuts in the 40 radial merger simulations with a mass of 10⁹ M_⊙ for a range of evolve times between 1 and 5 Gyr. The correct merger time was calculated for each simulation as when the center of mass of the dwarf galaxy passed through the Galactic center. Only the simulations with a mass of 10⁹ M_⊙ were used, since this mass is similar to the predicted mass of the VRM (see Section 7.1) and still has a wide range of possible phase-mixing times compared to more massive merger events (Figure 7). Using only a third of the simulations also allowed us to cut down on the number of simulations that had to be analyzed manually, while still retaining an idea of the capabilities and limitations of the method.

The next step was to locate the shells in the data. First, we used an EM algorithm to fit a Gaussian mixture model to a histogram of the data and optimized the number of Gaussians using a BIC method, similar to the process outlined in Section 3.2. However, the EM algorithm often preferred to fit Gaussian components to outlier peaks with small star counts if the interior shells were located close to one another. As it was not possible to supervise every fit in order to avoid this behavior, we chose to use a more robust algorithm utilizing a residual sum of squares as the goodness of fit.

This fitting algorithm minimizes a binned residual sum of squares,

$$\begin{eqnarray}&&\mathrm{RSS}=\displaystyle \sum _{i}^{m}\displaystyle \frac{1}{\eta }{\left[{\int }_{{r}_{i}}^{{r}_{i+1}}f(r){dr}-{n}_{i,\mathrm{obs}}\right]}^{2},\end{eqnarray} \tag{ 10 }$$

where m is the number of bins, f(r) is the fit Gaussian mixture model, r_i and r_i+1 describe the bounds of the ith bin, and n_i,obs is the value of the ith bin of the histogram of observed data. The observed data are compared to the integral of the fit model over each bin instead of the value of the model at the center of each bin; this will constrain the fit model to the number of stars in the data, as well as improve the estimate of the width of each shell. We also use a normalization constant

$$\begin{eqnarray}&&\eta =N-k-1,\end{eqnarray} \tag{ 11 }$$

where N is the total number of data points in the histogram and k is the total number of parameters being fit. We used the differential evolution algorithm (Storn & Price 1997) to minimize the residual sum of squares of our data, which provides the optimal fit for our data. Specifically, we use the implementation of the differential evolution algorithm from the scikit-learn python package (Pedregosa et al. 2011).

We employ the BIC in order to determine how many shells are statistically significant. For this case, the BIC is given by

$$\begin{eqnarray}&&\mathrm{BIC}=k\mathrm{ln}(N)+2\mathrm{ln}(\mathrm{RSS}),\end{eqnarray} \tag{ 12 }$$

where N and k are the same as in Equation (11). As before, the BIC prevents overfitting our observed data. We fit up to three Gaussians to both regions individually and compare the corresponding BICs; the fit with the lowest BIC was taken as the most statistically significant result. The positions of the peaks of the best-fit Gaussian mixture model were taken to be the distances of the shells.

We then used the model particle rewinding method on the adjusted simulation data to recover the most likely merger time for each simulation. We placed model particles in the centers of the selected regions at the distances of the shells. These model particles were allowed to move freely in our Milky Way potential, resulting in oscillations in their distance from the Galactic center. Tracing these oscillations throughout this period identified times when the model particles were all close to one another, corresponding to a likely time where the progenitor of the merger was still coherent. The best-fit merger time from the simulation was then selected from the results and compared to the correct merger time.

Figure 8 shows the result of our method on one such simulated radial merger. The galactocentric distances of the model particles are provided in order to show how δ(t) depends on the positions of the model particles. It is clear that the method is able to recover the time at which the oscillating model particles line up, as this is, on average, within a few tenths of a Gyr of the actual merger time of the simulation.

Zoom In Zoom Out Reset image size

In order to examine the dependence of the method on the shape of the Galactic potential, we also rewound the model particles in the Galpy built-in potential MWPotential2014 (Bovy 2015). This uses a spherical power-law potential for the central bulge, a less massive Miyamoto–Nagai disk than the Newberg et al. (2010) model, and a double power-law spherical halo potential. The two model potentials are compared in Figure 9. Calculating merger times with this second potential helps us determine whether or not an "incorrect" model potential impacts the results, as the radial merger simulations were only run in the Newberg et al. (2010) model potential. This helps us understand the limitations of our method, as the Milky Way's actual gravitational potential will not be truly identical to any model potential that we choose.

Zoom In Zoom Out Reset image size

Table 3 and Figure 10 show the results of our method on simulated data using both the "correct" and "incorrect" potentials. It is clear that the method is able to recover the infall times of the simulated radial mergers for a range of initial distances and evolve times. Including simulations that recovered two or more shells in the data, the standard deviation of the differences between the calculated and correct values is 0.52 Gyr for the "correct" Newberg et al. (2010) model potential and 0.53 for the "incorrect" MWPotential2014 model. The standard deviation of the differences decreases to 0.38 Gyr in the Newberg et al. (2010) model and 0.39 Gyr in the MWPotential2014 model if only simulations with three or more shells are considered, and it drops even further to 0.20 Gyr in the Newberg et al. (2010) model and 0.12 in the MWPotential2014 model if only simulations with four shells are considered. There are four shells identified in our observed data, so we estimate the error in infall time to be the larger of the two standard deviations of the differences: σ_t  = ±0.2 Gyr.

Zoom In Zoom Out Reset image size

This value is only the approximate error contributed by the recovery method. Other factors, such as using an incorrect shape of the Galaxy potential, might increase the actual error in our calculated value. However, the two model potentials used in this work did not seem to have a large impact on our method. This leads us to believe that any reasonably shaped model potential would likely produce similar results.

Some combinations of initial angles and distances of the simulated dwarf galaxy progenitors resulted in little to no shell structure. This is why, for example, there are no simulations with i = 0°, r₀ = 30 kpc, and a mass of 10⁹ M_⊙ in Table 3, as this simulation did not produce identifiable shell substructure. The simulations with inclination angle i = 30° appeared similar in R.A. and decl. to the simulations of individual components of the VRM debris in Donlon et al. (2019). It is possible that this inclination angle may be close to the correct value of the inclination angle of the progenitor of the VRM.

We only tested simulated evolve times up to 5 Gyr, as that is near the age of radial mergers at which shells stop being easily located. A more thorough suite of tests would include tests exploring the impacts of different progenitor profiles. Even without this exhaustive testing, we believe that our method is able to determine the correct time of collision of the VRM progenitor and the Milky Way, as the method shows clear success in recovering the correct values for a variety of radial merger simulations.

The concept of the "ancient last major merger" of the Milky Way has occupied the literature on Galactic structure for decades (Kuijken & Gilmore 1989; Gilmore et al. 2002). Gilmore et al. (2002) identified a group of disk stars rotating slower than expected in the local solar neighborhood and claimed that this group was evidence of the last major merger, which had occurred 10–12 Gyr ago. Further, Deason et al. (2013) claimed that a massive merger approximately 10 Gyr ago would explain the break in the density profile of the Milky Way halo at 30 kpc. The theory suggests that in the time since this massive merger event, there has been a period of Galactic quiescence, up until the Sagittarius dSph and Magellanic Cloud mergers that are currently underway.

The last major merger is thought to have formed the Milky Way's thick disk, which has stellar ages older than 10 Gyr. One theory for the creation of a thick disk is that the Milky Way's protodisk was heated by a collision with a large dwarf galaxy merger (Quinn & Goodman 1986; Velazquez & White 1999). This would have caused the existing disk stars to be kicked up on orbits with larger vertical action. Gas in the disk would radiatively cool, so stars formed after this merger event could remain in a cold, younger thin disk (effectively "quenching" the thick disk).

With the release of Gaia DR2, evidence suggested that the last major merger was finally identified. The GSM (Simion et al. 2019) and the Gaia–Enceladus Merger (Helmi et al. 2018) were independently discovered and are interpreted to be the same "last major merger of the Milky Way." The age of this merger event was reported to correspond to the age of the thick disk: between 8 and 11 Gyr. Helmi et al. (2018) found that the youngest age of stars in this merger was 8–10 Gyr, which is consistent with this massive merger indeed having created the thick disk. This agrees with the results of Gallart et al. (2019), who claimed that star formation in the Milky Way transitioned from thick to thin disk star formation around 9.5 Gyr ago, corresponding to the GSM. Gallart et al. (2019) also claimed that both the Gaia Sausage progenitor and the Milky Way's in situ halo had finished the bulk of their star formation by 10 Gyr ago. Belokurov et al. (2019) claimed that star formation in the inner halo ended around the same era, and, assuming that the cause for the quenching of star formation was the merger event itself, this is additional evidence for an ancient merger.

The presumed merger age determines a preferred mass of the ancient merger progenitor. Comparison of the metallicity of the RRLs in the Gaia Sausage to the RRLs in the LMC suggests that the progenitor of the Gaia Sausage had a similar mass as the LMC did 10 Gyr ago (Belokurov et al. 2017; Zinn et al. 2020). Mackereth et al. (2019) analyzed the chemical and orbital properties of the Gaia Sausage in comparison to EAGLE simulations and determined that it likely had a mass of at least 10⁹ M_⊙ and was accreted around 9 Gyr ago.

This large mass suggests that this merger is the main component of the stellar halo; Deason et al. (2019) showed that the total stellar mass of the halo is on the order of 10⁹ M_⊙. Simulations suggested that most of the stars on radial orbits in the halo came from massive satellites approximately 6–10 Gyr ago, and that these stars would dominate the halo between 10 and 30 kpc of the Galactic center (Belokurov et al. 2018). The inner halo is indeed dominated by stars on highly eccentric orbits (Iorio & Belokurov 2019). This is all consistent with the halo being primarily formed from a single merger event.

For all of these reasons, the present view of the Gaia Sausage is that it is an ancient merger event with a stellar mass between 10⁸ and 10⁹ M_⊙ that is responsible for the formation of the thick disk and makes up the majority of the stellar mass of the halo. We will now attempt to reconcile this viewpoint with the merger time of the VRM progenitor that was calculated in this work.

The Gaia Sausage is characterized by a structure in velocity space with high dispersion in v_r and small rotational velocities in the local solar neighborhood (Belokurov et al. 2018). It has been shown that the VRM debris looks very similar to the Gaia Sausage in the local solar neighborhood, so much so that the two structures are likely identical (Donlon et al. 2019). If the VRM debris and Gaia Sausage are identical, then how do we explain the gap between the 2.7 Gyr ago infall time of the progenitor of the VRM and the 8–10 Gyr age of the Gaia Sausage?

In Section 4, we developed an argument that the VRM debris must populate the stellar halo at the distance ranges in which we identify shells. Since the VRM debris dominates the stellar halo, the debris must also be located near and inside the stellar break at ∼20 kpc. To satisfy these conditions, the progenitor of the VRM must have had a final apogalacticon before collision between 20 and 30 kpc. The measured values of the causticality in simulations with r₀ = 20 and 30 kpc in Figure 7 match the measured causticality values for the observed data around 1 and 5 Gyr, respectively. This suggests that the merger time of the VRM is within the last 5 Gyr. Additionally, by 5 Gyr after collision, our method's ability to recover the merger time of a radial merger dropped substantially (Figure 10). Finally, this work recovered a merger time of only 2.7 Gyr for the VRM progenitor. All of this suggests that the progenitor of the VRM collided with the Milky Way more recently than 8–11 Gyr ago.

If the dynamical mass of the VRM is actually closer to 10¹⁰ M_⊙, then even large initial distances constrain its merger time to be within the last 5 Gyr (Figure 7). The GSM is expected to dominate the halo and would then have a dynamical mass of around 10⁹–10¹⁰ M_⊙. This suggests two possible scenarios: if the GSM is indeed ancient, it either cannot dominate the material in the halo or cannot be responsible for the relatively unmixed stellar overdensities in the halo. In this case, the VRM and GSM cannot be identical events, as the VRM debris populates shells in the VOD and HAC. If the VRM and GSM are not the same, then the Milky Way halo is composed of more than one major merger event. On the other hand, if the GSM is responsible for the VOD and HAC and dominates the halo, then it is unlikely that the merger event occurred more than 5 Gyr ago. If this is the case, the VRM and GSM probably describe the same merger event.

Villalobos & Helmi (2008) ran simulations of dwarf galaxy mergers that collided with a model Galactic disk after falling in from near the virial radius. After an orbital decay over a period of ∼1.5 Gyr, the dwarf galaxy finally collides with the Galaxy from a distance of ∼20 kpc from the Galactic center. Thick disks are shown to form approximately 0.25 Gyr after these collisions. The work states that thick disk and satellite disruption reach an equilibrium ∼2 Gyr after the collision occurs, and that major shell structure is only visible for these 2 Gyr. If the thick disk was formed in this way 8–11 Gyr ago, then shells in the resulting merger debris would not be visible at the present day.

In this work, we use the term "merger time" to describe when the progenitor of the VRM's center of mass passes through the Galactic center. This is a collision between the dwarf galaxy and the Milky Way. If the progenitor for the Gaia Sausage crossed the Milky Way's virial radius sometime around 10 Gyr ago ("infall time"), it could have spent the last 8 Gyr in a prolonged orbit about the Galaxy while it shed its initial energy until finally colliding with the Galaxy 2.7 Gyr ago. This would have resulted in the progenitor of the Gaia Sausage making many passes through the Milky Way over its lifetime, which could cause tidal shocks in both structures. We speculate that these tidal shocks from perigalacticon passes of the Gaia Sausage progenitor could potentially kick up thick disk stars from the Milky Way protodisk through large-scale tidal torques, as well as generate star formation through the collapse of gas clouds in the disk. This could explain the findings that the star formation history of the thick disk oscillates with a period of ∼3 Gyr (Gallart et al. 2019), which is a reasonable time between passes for a distant halo object with a highly eccentric orbit.

We did not include dynamical friction in our models of the VRM in this work, because it does not play a role during and after the disruption of the satellite as it passes through the center of the Milky Way. However, dynamical friction would allow the VRM progenitor to shed energy as it orbits and falls inward from the virial radius. A reduction in energy is required for a dwarf galaxy formed outside or near the Milky Way's virial radius to eventually become substantially more bound to the Milky Way and would cause the VRM progenitor to slow down from earlier nonradial passes of the Milky Way before its final radial collision.

An initial tidal shock from a close encounter with the Milky Way could be responsible for the quenching of star formation in the Gaia Sausage 8–11 Gyr ago (the stellar "age" of the progenitor). However, according to Brown et al. (2014), many dwarf galaxies in the Local Group stopped star formation around the same time 12 Gyr ago during reionization of the universe, so alternatively, it could be that the Gaia Sausage had already stopped active star formation even before it became gravitationally bound to the Milky Way. This could also explain why the star formation in both the Milky Way in situ halo and the Gaia Sausage progenitor appear to slow down around 12 Gyr ago (Gallart et al. 2019), even before the proposed accretion time of the progenitor of the Gaia Sausage around 8–11 Gyr ago. If the accretion event were the primary underlying cause of the drop in star formation rates in both structures, as Gallart et al. (2019) claimed, then it is not expected for star formation to begin slowing down before the accretion event occurs. It is possible that the previous metallicity arguments for the age of the Gaia Sausage were not dating the time of the collision between the dwarf progenitor and the Milky Way but instead were measuring other properties of the infall event. Note, however, that our method for determining merger time only measures the time since the collision between the progenitor and the Milky Way and does not constrain when the progenitor became bound to the Milky Way, star formation was quenched in the Milky Way or the progenitor dwarf galaxy, or the thick disk was created.

If the Gaia Sausage did not collide with the Milky Way 8–11 Gyr ago, then what caused the thick disk? Amarante et al. (2019), Beraldo e Silva et al. (2020), and Clarke et al. (2019) all showed that it is possible to form a thick disk in Milky Way analogs in isolation, provided that there is a period of lumpy disk activity early in the Galaxy's life. Additionally, Ma et al. (2017) showed a cosmological simulation of a Milky Way–mass galaxy forming a bifurcated disk without a collision. These systems are able to generate both the kinematic and chemical properties of the thick and thin disks without merger events. The thick disk may form itself without perturbation from a satellite, in which case, there is no need for the Gaia Sausage to be the perturber of the thick disk. Additionally, Rodriguez Wimberly et al. (2019) claimed that it is unlikely that star formation in a thick disk would be quenched by the infall of a large dwarf galaxy. If either of these situations is the case, then the age of the GSM is not required to line up with the time of the quenching of star formation in the thick disk.

The EPO is a recently discovered overdensity with a lower surface brightness than either the VOD or the HAC. It has been suggested that the EPO is associated with the VOD and HAC and that they share a common origin (Li et al. 2016; Donlon et al. 2019). The results of this work are consistent with this idea; Figure 5 shows ejecta in a direction 120° from the two regions containing the majority of the shell substructure, as viewed from the Galactic center. This work suggests a single "trefoil" structure connecting the three overdensities, instead of a model where material oscillates only between the HAC in the south and the VOD.

We note that the VOD appears to visually extend to higher declinations, tracing the Perpendicular Stream (Weiss et al. 2018b). It is possible that the Perpendicular Stream is associated with one or more shells or that it traces a long radial tail of infalling or outgoing material. It may be beneficial to look for shell structure in the upper portion of the Perpendicular Stream in the future in order to further constrain the behavior of the progenitor of the VRM. A study of the southern portion of the HAC would also be beneficial in exploring the full extent of the VRM debris, as well as continuing the search for shells in other regions of the Milky Way.

The phase-space spiral in the disk is a structure in z–v_z phase space that has been characterized as ongoing phase mixing in the local solar neighborhood from a recent perturbation in the disk (Li & Shen 2020). The disk is known to not be in equilibrium (Widrow et al. 2012; Carlin et al. 2013; Williams et al. 2013; Yanny & Gardner 2013; Xu et al. 2015), and the vertical disk disequilibrium presents locally as the phase-space spiral (Li & Shen 2020). The phenomenon that caused the phase-space spiral occurred at least 500 Myr ago, and it is suggested that the structure could only exist in the Milky Way for ∼4 Gyr (Li & Shen 2020). Our merger time for the progenitor of the VRM lies within these constraints, and we suggest that it is possible that the cause of the phase-space spiral is indeed the VRM. Perhaps the VRM "wobbled" the Milky Way's central bulge in a way that would propagate waves in the disk and generate vertical motion in disk stars. Further exploration is required to determine the relationship between the VRM and the phase-space spiral, if any exists.

The "Splash" (Belokurov et al. 2019) is another substructure of the Galactic disk in the local solar region. The Splash is characterized by a large population of metal-rich stars on highly radial orbits in the inner halo. While the formation origin of stars in the Splash is not yet known, it has been hypothesized that the Splash consists of stars thrown out of the Milky Way's (proto)disk from the GSM 9.5 Gyr ago (Belokurov et al. 2019). Since the VRM is a radial merger, it is likely that metal-rich material in the Galactic center would be disrupted from the initial impact of the progenitor of the VRM with the Galactic center, which could explain the high-metallicity material in the Splash. Amarante et al. (2019) showed that the Splash could instead have been created through another disk process, such as a lumpy disk early in the Milky Way's formation history. Either way, there is no need for an ancient infall event to have created this substructure.

Snaith et al. (2014) showed an oscillation in star formation of the Milky Way inner disk from 6 Gyr ago until approximately 2 Gyr ago. There is a strong peak in star formation around 2.5 Gyr ago. It is possible that this oscillation in the star formation rate could be due to tidal shocks from the passes of the VRM progenitor as it fell in from the virial radius, and that the burst of star formation 2.5 Gyr ago was due to the infall of gas from the VRM and perturbations in the disk due to the collision event. Mor et al. (2019) also find a strong peak in star formation around 2.7 Gyr ago that may correspond to the VRM. It is also possible that the increase in star formation around 2.5 Gyr ago could be caused at least in part by an impact of the Sagittarius dSph progenitor with the disk; Ruiz-Lara et al. (2020) attributed a peak in the star formation history of the local solar region around 2.0 Gyr ago to such an event. It should be noted that Mor et al. (2019) do not find these peaks in star formation associated with the passages of the Sgr dwarf.