Exploring the Perturbed Milky Way Disk and the Substructures of the Outer Disk

Y. Xu; C. Liu; H. Tian; H. J. Newberg; C. F. P. Laporte; B. Zhang; H. F. Wang; X. Fu; J. Li; L. C. Deng

doi:10.3847/1538-4357/abc2cb

1. Introduction

There is a great deal of evidence that the Milky Way disk is in a state of disequilibrium. Near the Sun, there are vertical oscillations both in the density and bulk velocities of disk stars (Widrow et al. 2012; Carlin et al. 2013; Williams et al. 2013; Carrillo et al. 2018). There are kinematic substructures and a radial velocity gradient (Siebert et al. 2011; Tian et al. 2015, 2017; Gaia Collaboration et al. 2018; Wang et al. 2018). Ring-like overdensities appear alternately on the north and south sides of the disk (Newberg et al. 2002; Rocha-Pinto et al. 2004; Xu et al. 2015; Price-Whelan et al. 2015). From an analysis of Gaia DR2 data, López-Corredoira & Sylos Labini (2019) provided evidence that the kinematic distribution of disk stars in the Milky Way deviates from a "simple stationary configuration in rotational equilibrium."

To explain the perturbations, both internal processes (perturbation and resonance of the bar and the spiral arms) and external ones (accretion and passage of satellites) have been proposed (Quillen et al. 2011; Siebert et al. 2012; Gómez et al. 2013, 2016; Debattista 2014; Faure et al. 2014; Kawata et al. 2014; Bovy et al. 2015; Grand et al. 2015; Monari et al. 2015; Tian et al. 2017; Widrow et al. 2014; Khoperskov et al. 2019).

Some studies have considered the coupled effects between the vertical perturbation and planar perturbation. Recent simulations show that the rotating bar and spiral can produce vertical perturbation (Monari et al. 2015) and the passage of a satellite can also produce radial and azimuth perturbation (D'Onghia et al. 2016).

Some of the substructures of the outer disk have been well studied, such as the famous "Monoceros overdensity" (Newberg et al. 2002; Crane et al. 2003; Ibata et al. 2003; Yanny et al. 2003; Conn et al. 2005; Martin et al. 2006; Ivezić et al. 2008; Li et al. 2012; Meisner et al. 2012). The ring-like Monoceros overdensityz was originally thought to be tidal debris (Martin et al. 2004; Bellazzini et al. 2006; Peñarrubia et al. 2005) or part of a warp or flare (Momany et al. 2004; López-Corredoira & Molgó 2014; Wang et al. 2018), but more recent simulations suggest the ring could be induced by the gravitational interaction of the Sagittarius dwarf galaxy with the disk (Purcell et al. 2011; Laporte et al. 2018).

Other studies show that the outer disk overdensities could be related. Deason et al. (2018) and de Boer et al. (2018) suggested that the Anti-Center Stream (ACS) and Eastern Banded Structure (EBS; Grillmair 2006, 2009) are parts of Monoceros overdensity, even though they appear to be distinct structures in photometric data. From star ages, Laporte et al. (2020) showed that thin structures like the ACS and EBS are tidal structures excited by a satellite interaction with the disk, while the Monoceros ring is a distinct structure formed through the gradual flaring of the disk over a more extended disk-satellite interaction. Price-Whelan et al. (2015) suggested that the Triangulum Andromeda Overdensity (TriAnd) is a disk component, since the stellar population is more similar to the disk than to tidal debris from a dwarf galaxy. Bergemann et al. (2018) show that TriAnd and A13 have chemical abundance patterns similar to those of the thin disk. Laporte et al. (2020) showed that the Monoceros ring, ACS, and EBS have [Fe/H]∼-0.7 with no vertical gradient in [Mg/H] versus [Fe/H], which is consistent with chemistry of thin disk and not consistent with the "knee" typically seen in [Mg/H] versus [Fe/H] for dwarf spheroided galaxies (Hayden et al. 2015). Li et al. (2017) suggested that all of the substructure in the outer disk: the Monoceros overdensity, TriAnd1 (Sheffield et al. 2014), A13 (Sharma et al. 2010), and TriAnd2; which appear alternately north and south of the galactic plane with increasing distance from the Galactic center, are vertical oscillations of the disk. Xu et al. (2015) showed that the outer disk overdensities which appear at 2 kpc (north near structure) and 5 kpc (south middle structure) from the Sun can be fit with a star count toy model for a bending wave in the disk. Simulations show that the passage of the Sagittarius dwarf spheroidal galaxy (Sgr dSph) through the Milky Way halo could produce oscillations of the disk with ring-like structures similar to those that are observed (Purcell et al. 2011; Gómez et al. 2017; Laporte et al. 2018).

Antoja et al. (2018) found a Z − V_Z phase space spiral, believed to result from nonequilibrium phase mixing from the passage of a satellite through the disk some several hundred million years ago. Binney & Schönrich (2018) and Bland-Hawthorn et al. (2019) made test particle simulations that can produce a phase space spiral that is quite similar to that found in Gaia DR2 and the Galactic Archaeology with HERMES (GALAH) data. Laporte et al. (2019) showed four density wraps on the Z − V_Z phase space of Gaia DR2 data within 1 kpc from the Sun; their simulation of the Milky Way disk interacting with Sgr dSph can reproduce many features, which are revealed in the Gaia data. Khanna et al. (2019) showed the ridge-like structures in the R − V_ϕ phase space correspond to constant energy or constant angular momentum, and show that phase mixing can produce ridges of constant energy. Chequers et al. (2018) suggested that the cumulative effect of the interaction between the disk and the satellites in a clumpy halo can produce long-lived bending waves.

All of the observations mentioned above were made in a volume of the disk that is limited to 3–4 kpc from the Sun. Large Sky Area Multi-Object Fibre Spectroscopic Telescope (LAMOST) K giants can provide far more distant data, up to at least 20 kpc from the Galactic center to compare with the results of simulations. Through studying the LAMOST K giants, we found a strong connection between the disk kinematic substructures and the phase space spiral. Our work provides new insight into the origin of the outer disk substructures.

In Sections 2–4, we describe our data selection, distance calibration, and calculation of the 3D velocities for the K giant stars. In Section 5, we analyze the V_ϕ kinematic maps projected onto the R – Z plane; we explore the velocity distribution in phase space and the connection between the substructures of K giants and the phase space spiral. Sections 6–8 show the V_R, V_Z, and metallicity distribution in the R – Z plane. In Section 9, the observational kinematic features are compared with the results of the simulations. Sections 10 and 11 respectively present the discussion and conclusion.

2. Sample Selection

K giant stars from LAMOST DR5 were selected using the selection criteria provided in Liu et al. (2014). The stellar parameters for the K giants were obtained from the Stellar LAbel Machine (SLAM) procedure (Zhang et al. 2020). We then crossmatched the LAMOST DR5 K giants and Gaia DR2 stars using Tool for Operations on Catalogues And Tables (TOPCAT). The high-precision line-of-sight velocities measured by LAMOST DR5 (Cui et al. 2012; Deng et al. 2012; Zhao et al. 2012), and the proper motions from Gaia DR2 (Perryman et al. 2001; Prusti 2018; Brown 2018), are combined to obtain 3D velocities.

We remove K giants with true metallicity [M/H] < −1.2 from the sample to eliminate halo stars. We also remove duplication and red clump stars that were identified in Ting et al. (2018). Red clump stars are removed because they have different absolute magnitudes from other K giants with the same color. We require a high signal-to-noise in the g band (snrg > 15) to assure the accuracy of derived stellar parameters and radial velocity errors of less than 10 km s⁻¹, and proper motion errors less than 0.2 mas yr⁻¹ are required to assure the precision of the velocity calculation. Stars with M_K > 0 are eliminated because they have systematic distance errors larger than 10% compared to distances from Gaia DR2 (Bailer-Jones et al. 2018), possibly due to contamination from dwarf K stars. Using all of these selection criteria, we construct a sample of 429,500 K giant stars.

3. Distance Calibration

Our results rely on creating accurate maps of the kinematic distribution of K giants as a function of position in the Milky Way. To achieve this, we need accurate distances for the stars in our sample, which can be observed out to 20 kpc from the Galactic center. We derive the distances to the K giants by applying the Bayesian estimation method from Carlin et al. (2015), using the stellar parameters from the SLAM pipeline; the distance estimates from Bayesian estimation with parameters from SLAM were much closer to Gaia distances than those produced using the standard LAMOST pipeline stellar parameters.

In Carlin et al. (2015), they calibrated the distances derived from LAMOST data with those from the Hipparcos sample to determine that the distance estimates have an uncertainty of about 20%. However, Hipparcos data could only validate distances within several hundred parsecs from the Sun. The distances derived from Gaia DR2 parallaxes can validate the distances within 4 kpc from the Sun. So, we recalibrate the distance of the LAMOST DR5 data with the best estimated distances provided by Gaia DR2 parallaxes as found by Bailer-Jones et al. (2018). Stars from Gaia DR2 within 3 kpc of the Sun and with relative distance errors of less than 10% are matched to our data set of LAMOST K giant stars. The distance error is defined by 1σ of the Gaussian distribution of estimated distance probability (Bailer-Jones et al. 2018). But keep in mind that the distance obtained by Bailer-Jones et al. (2018) is not model-free; the Bailer-Jones et al. (2018) distance is still influenced by the adopted length scale model.

Figure 1 shows the distribution of the relative difference in distance for the two catalogs, Δ_d, as a function of the Gaia DR2 distance for each star. The relative difference in distance is defined as Δ_d = (d_Carlin−d_BJ)/d_BJ. d_Carlin is the distance estimated by procedure outlined in Carlin et al. (2015), applied to the LAMOST DR5 data. d_BJ is distance obtained from Gaia DR2 parallaxes, as determined by Bailer-Jones et al. (2018). The red dots indicate the median value of Δ_d in each distance bin of width 100 pc. The median value of median(Δ_d) is about 0.006, and the variation of the median(Δ_d) of the calibration sample with distance is quite small. From this analysis, we determine that the systematic distance error can be ignored, since it is considerably smaller than the random error.

We further checked our distances against a cleaner sample of Gaia parallaxes. Schönrich et al. (2019) found that the Gaia DR2 parallaxes have a bias of 0.054 mas yr⁻¹ and produced a sample of corrected Gaia DR2 parallaxes (Schönrich et al. 2019) that were free from this bias. This sample uses the entirely safe criteria from Section 8 of Schönrich et al. (2019) to select stars from the Gaia sample, thus assuring that the estimated Gaia distances are accurate. We crossmatch the LAMOST K giants sample with the selected sample. The right panel of Figure 1 shows the distribution of the relative difference of distances, Δ_d, versus the estimated distance from Schönrich et al. (2019). The median value of median(Δ_d) is about 0.004, confirming that the LAMOST distances are also consistent with the estimated distance of Schönrich et al. (2019), as expected. However, there appears to be a trend with distance; at larger distances, the Bayesian distance estimates are lower than the corrected Gaia measurements, with a maximum deviation of 5%.

4. Cylindrical Velocities: V_r, V_ϕ, and V_z

The 3D velocities are calculated with using line-of-sight velocities from LAMOST DR5, proper motions from Gaia DR2, and LAMOST distances calculated using the Bayesian technique from Carlin et al. (2015). A known systematic error of 4.4 km s⁻¹ in radial velocity is subtracted from the LAMOST catalog values (Wang et al. 2018). The 3D velocities in the galactic Cartesian coordinates (U, V, W) are calculated following the method used in Johnson & Soderblom (1987). Positive U, V, and W are oriented toward the Galactic center, the direction of the galactic rotation, and the north Galactic pole, respectively. We adopt a distance from Sun to the Galactic center of R_⊙ = 8.34 kpc (Reid et al. 2014), the height of Sun above the galactic plane of Z₀ = 27 pc (Chen et al. 2001), the peculiar velocity of the Sun with respect of the local standard of rest of (U_⊙, V_⊙, W_⊙) =(9.58, 10.52, 7.01) km s⁻¹ (Tian et al. 2015), and a circular velocity at the solar radius of V_c = 238 km s⁻¹ (Schönrich et al. 2010). Then the 3D Cartesian galactocentric velocities can be obtained from U = U + U_⊙, V = V + V_⊙ + V_c, W = W + W_⊙. The cylindrical velocities V_R, V_ϕ, and V_z are obtained following equations A13, A14, and A15 of Williams et al. (2013).

5. Median V_ϕ Distribution

In this section, we study both the spatial (as a function of R and Z) and phase space (as a function of Z and V_Z) distribution of V_ϕ. We also connect the bulk motion of disk stars, as a function of spatial position, with the phase space spiral. The spatial distribution of V_ϕ is illustrated in Section 5.1. The phase space distribution of V_ϕ is illustrated in Section 5.2. The relationship between the bulk motion and phase space spiral is illustrated in Section 5.3.

5.1. Median V_ϕ as a Function of (R, Z)

Figure 2 shows the spatial distribution of our selected LAMOST K giant stars. They extend to more than 20 kpc from the Galactic center in the radial direction and extend several kiloparsecs from the galactic plane in the Z direction. The right panel of Figure 2 shows that the LAMOST K giants concentrate toward the anticenter direction (−20° < ϕ < 10°, where ϕ is the galactocentric azimuthal angle from the anticenter direction).

**Figure 2.** The spatial distribution of LAMOST K giant stars. The left panel shows galactocentric Y vs. X. The middle panel shows galactocentric Z vs. X. The right panel shows that the stars are heavily concentrated toward the Galactic anticenter. Here, ${phi}=\arctan (Y/X)$ , which is the azimuthal angle from the anticenter directions, as seen from the Galactic center. Positive ϕ are found in the third quadrant, and negative ϕ are found in the second quadrant.
Download figure:
Standard image High-resolution image

In this subsection, we create and analyze two-dimensional projections of the bulk velocities of these stars in the R − Z space. This is a convenient projection because K giants concentrate in a narrow range around the anticenter direction. Figure 3 shows the number of K giant stars in our sample as a function of (R, Z). All of the bins contain at least 10 stars; within R < 20 kpc, most bins have at least 50 stars.

**Figure 3.** Heat map of the number of K giants in each bin of the R, Z map using logarithmic scaling. Note that most of the stars are located closer to the Sun. The lines and square are for comparison with the structures in Figures 4–7.
Download figure:
Standard image High-resolution image

Heat maps of the median V_R, V_ϕ, and V_Z as function of (R, Z) are now constructed to study the kinematic distribution of disk stars as a function of spatial position. The bin size in our maps is (ΔR, ΔZ) = (0, 5, 0.5) kpc.

5.1.1. Disk Substructure in V_ϕ

The first row of Figure 4 ¹¹ shows the median(V_ϕ) and σ(V_ϕ) (the standard deviation of V_ϕ) for our sample, as a function of (R, Z). From this figure, the disk and halo populations can be clearly distinguished; the disk has a high V_ϕ and a low σ(V_ϕ).

The right panels of Figures 4 also illustrate the distribution of σ(V_R), σ(V_Z), and σ([M/H]), respectively. In these panels the disk-like stars also show the same characteristics; the disk-like stars are kinematically colder and have a smaller metallicity scatter compared to those of halo-like stars. The top right panel of Figure 4 shows that in σ(V_ϕ) there is a smooth transition from the kinematically cold area with disk-like stars to the kinematically hot area with halo-like stars, especially before R < 13 kpc. Beyond R = 13 kpc, the transition is sharper.

If one looks at stars that are kinematically colder and have smaller metallicity scatter across the entire R – Z plane, one finds that the distribution of median(V_ϕ) and σ(V_ϕ) of disk-like stars is horn-shaped. The scale height of disk stars grows quickly with R in the outer disk, characteristic of a disk flare. The flare is consistent with the one seen in blue stragglers in the Canada-France Imaging Survey, crossmatched to Gaia DR2 and Sloan Extension for Galactic Understanding and Exploration (SEGUE) and LAMOST (Thomas et al. 2019).

In the solar neighborhood, the V_ϕ distribution is quite similar to the known disk kinematics, which are "dominated by rotation with a smooth vertical gradient" (Bond et al. 2010).

Along the midplane, the median V_ϕ decreases with increasing R, which is consistent with the result of Tian et al. (2016). In the range of R < 11 kpc, the median V_ϕ for bins in the midplane of the disk reaches to 230 km s⁻¹. The median V_ϕ in the disk midplane decreases to around 210 km s⁻¹ when R = 15 kpc.

Beyond R = 11 kpc, the most significant structure in disk kinematics is that the high V_ϕ stars are split into three branches, which we call the main branch, north branch, and south branch. The main branch bends southward; the ridgeline of the main branch is distributed roughly along the line from (R, Z) = (12, −0.5) kpc to (R, Z) = (14, −1) kpc. The north branch is found along the slope from (R, Z) = (13, 1.5) kpc to (17, 5) kpc. The south branch is found along the line from (R, Z) = (11, −1) kpc to (R, Z) = (16, −4) kpc. The south branch has higher median V_ϕ than that of north branch. Both of the north and south branches stretch along the boundary of the flare.

It is apparent from the σ(V_ϕ) panel of Figure 4, in the range of R < 13 kpc, that the ridge of the minimum standard deviation in V_ϕ is located slightly above the midplane, roughly at Z = 0.25 kpc. When R > 13 kpc, the ridgeline of minimum velocity dispersion skews toward the south, following the main branch. This is consistent with a vertical displacement of the midplane of the disk. We will discuss this in more detail in Section 5.1.2.

Also apparent in the σ(V_ϕ) panel of Figure 4 is an area with quite small σ(V_ϕ) in the range of 16.5 < R < 20 kpc, 1.5 < Z < 5 kpc. This location is consistent with famous substructure called the Monoceros overdensity. The location of the main branch, north branch and south branch kinematic structures, identified from the median(V_ϕ) panel of Figure 4, as well as the "Monoceros area" identified from the σ(V_ϕ) panel of Figure 4, are labeled by lines and a square in each panel the figure.

Estimates of the kinematic characteristics of the north branch, south branch, and Monoceros area from Figure 4 are summarized in Table 1. To measure these properties, we selected stars most likely associated with the substructures, in the regions of parameter space delineated in Figure 4. The measured properties are the median over the the selected bins in V_ϕ, V_R, and V_Z.

Table 1. Characteristics of Kinematic Substructures of Disk Stars

Kinematic Features	Median(V_ϕ)	σ(V_ϕ)	Median(V_R)	σ(V_R)	Median(V_Z)	σ(V_Z)
	(km s⁻¹)	(km s⁻¹)	(km s⁻¹)	(km s⁻¹)	(km s⁻¹)	(km s⁻¹)
(i) North branch	213	32.5	−1.6	35	−4.3	33.9
(ii) South branch	213	40.7	0.9	42	−2.5	39.3
(iii) Monoceros area	195.5	26.4	−12.8	30	−7	31.

Download table as: ASCII Typeset image

5.1.2. The Oscillating Disk Traced by Kinematic Features

Using Sloan Digital Sky Survey (SDSS) main-sequence star counts, Xu et al. (2015) found that the midplane of the disk stars appears to oscillate vertically across the galactic plane, as a function of distance from the Galactic center. The oscillating disk asymmetries appear at D = 2, 5, 10, and 15 kpc, where D is the distance from the Sun. The two nearer oscillation asymmetries can be approximated by model with an oscillating disk, in which the disk is offset up by 70 pc at a distance of 2 kpc from the Sun toward the anticenter, and down by 170 pc at a distance of 5 kpc from the Sun toward the anticenter.

In this work, the midplane is traced by the locus of the minimum standard deviation in V_ϕ and the minimum standard deviation in V_Z. We slice the data into bins of different galactocentric radius, R, and find the Z location of the minimum standard deviation for each bin. We use a bin size of ΔR = 0.2 kpc for R < 12 kpc and 0.5 kpc for R > 12 kpc, due to fewer stars and larger statistical noise in the measurement of the minimum at large radius.

In Figure 5, the red triangles and the blue triangles show the locations of the minimum standard deviation of V_ϕ and V_Z, respectively. The Z position of the minimum standard deviation in V_ϕ and V_Z follow each other very closely. The midplane shifts north in the range of R = 10–13 kpc with amplitude of up to 300 pc, and then again shifts south. The black curve is the "oscillation" derived from fitting disk main-sequence star counts in Xu et al. (2015). The trend found in the oscillation model obtained by fitting SDSS main-sequence stars is consistent with that of the midplane oscillation identified kinematically using LAMOST K giants. However, the oscillation in K giant stars suggests a larger oscillation amplitude. Wang et al. (2018) identified the south middle structure from LAMOST K giant star counts. They clarified that the south middle structure defined by Xu et al. (2015) is an extended structure from (R, Z) = (10, −0.5) kpc to (15, −3) kpc, and found that the asymmetry across the plane could be erased by shifting the midplane south by about 300 pc at R = 14 kpc. This suggested a midplane offset and is indicated in Figure 5 with a green plus sign.

We traced the midplane using the location of minimum standard deviation of V_ϕ and V_Z for R < 14.5 kpc. After 14.5 kpc, the locus of the minimum standard deviation in V_ϕ is dominated by the substructure with high Z. For example, the Monoceros area has a small standard deviation in V_ϕ with a median Z of about 3 kpc, which is more consistent with a substructure kicked out of the disk than it is with an oscillation.

5.2. The V_ϕ Phase Space Spiral and the Outer Disk Substructures

Antoja et al. (2018) discovered an impressive snail-like spiral in Z versus V_Z phase space, providing direct evidence of phase mixing produced by a strong perturbation of the disk. The Sgr dSph was identified as the most likely galactic intruder to produce this perturbation, due to the match between the very recent impact time and the phase mixing timescale.

In this section, we will show the variation of the V_ϕ phase space spiral with galactocentric radius, and show the connection between the V_ϕ phase spiral and the outer disk kinematic features.

5.2.1. Comparison of the V_ϕ Phase Space Spiral in LAMOST K Giants with That of Gaia DR2 Data in the Solar Neighborhood

Before we construct the phase spiral maps for different ranges of galactocentric radius, it is necessary to compare the V_ϕ phase spiral seen in LAMOST K giants with that of Gaia DR2 data in the solar neighborhood.

The upper panel of Figure 6 shows the median V_ϕ distribution in the Z − V_Z phase space for 620,000 stars crossmatched between LAMOST DR5 and Gaia DR2, that have galactocentric distances in the 8.24 < R < 8.44 kpc range. We refer to this as the "total sample" in this subsection. The upper panel of Figure 6 perfectly reproduces the two and half circles of phase spiral discovered by Antoja et al. (2018). To produce this panel we used distances from Bailer-Jones et al. (2018).

The upper right panel of Figure 6 shows the 12,870 K giants selected from the total sample with the additional selection criteria of Section 2. Using these stars we can only barely see the phase spiral. The K giants in the disk might not be quite as good as main-sequence stars for illustrating the phase spiral, because they are a kinematically hotter disk population with a larger scale height and therefore less influenced by the Sgr dSph impact. Also, there are many fewer stars in the K giant sample. K giant stars comprise only 1/20th of the total sample within 8.24 < R < 8.44 kpc; the subset of K giants that meet the criteria outlined in Section 2 are only 1/50th of the total sample in this radial slice.

To determine the main reason for the difference in clarity between the upper panels of Figure 6, we tried randomly extracting 1/50th of the stars from the total sample and building the Z − V_Z phase space map. The histograms in the second row of Figure 6 show the distribution in Z of the total sample and K giant stars. The K giants are more spread out in the Z direction, and less concentrated right at the solar position. This is because they are intrinsically bright and can be observed to farther distances than many of the stars in the total sample. To determine whether the wider range of distances is important, we subsampled 12,870 stars from the total sample in two ways. In the left panel of the third row, the stars are randomly sampled. In the right panel of the third row, the stars are extracted from the total sample with a distribution in Z that matches the K giant stars' Z distribution. The two panels in the third row are also produced with distances from Bailer-Jones et al. (2018).

Both of the phase spirals of random samples in the third row of Figure 6 are similar, and only a little clearer than that of the K giant sample in the upper right panel. Extracting randomly from the total sample or extracting based on a broader distribution in Z did not make much difference, but the lower star number did make a big difference. From this experiment, we illustrate the importance of sample size in tracing the phase space spiral. Since K giants are much more sparsely distributed in space than main-sequence stars, we cannot gather as many K giants as main-sequence stars in a thin slice. Because the phase space spiral changes with galactocentric distance, there is a trade-off between using a thicker slice that will include more stars, and a thinner slice that will sample a narrower range of spiral properties; the phase spiral will be blurred in thicker slices of R. We choose 1 kpc bins in R as a trade-off between gathering more stars and cutting the slice as thinly as possible.

The lower panels of Figure 6 show the K giant phase spiral for 8 < R < 9 kpc. The sample size is about 1/12th of the total sample in the upper left panel. The lower left panel is constructed with distances from Carlin et al. (2015); it shows a big central "bulge"(centroid), which is the blurred inner circle of the phase spiral and also shows the wide strong, tail of the phase space spiral. In Figure 1 of Antoja et al. (2018), the width of the inner circle is about ΔZ = 0.1–0.2 kpc and the outer circle width is about 0.2–0.3 kpc. Because the inner portion of the spiral has a finer structure, the inner phase spiral cannot be distinguished from the LAMOST K giant sample in solar neighborhood when the sample size is too small. The lower right panel shows the same data as the lower left panel, but constructed instead with the Gaia DR2 distances (Bailer-Jones et al. 2018), which is more precise than the LAMOST distances (Carlin et al. 2015) in the solar neighborhood. In the lower right panel, we see that the outer portion of the spiral more closely matches the position of the outer spiral in the upper left plot. But the inner portion is not sampled, as can be understood from a comparison of the histograms in the second row of the figure; the central region of the spiral is heavily populated with a large number of intrinsically fainter stars.

Given that the LAMOST K giants are much sparser tracers but sample a much larger volume, we can explore the phase spiral over a larger portion of the Milky Way, but with lower resolution than was seen in local Gaia DR2 data. We will use LAMOST K giants to determine whether there is a phase space spiral in volumes within 20 kpc of the Galactic center and trace the outline of any observed spiral, but we will not be able to observe the phase space spiral's detail.

5.2.2. The Median(V_ϕ) Phase Spiral as a Function of Distance from the Galactic Center

We are now in a position to look for variations in the phase space spiral as a function of position in the Galaxy. First, we separate the disk stars into 1 kpc bins along the projected galactic radius, R, in the range of 6 < R < 20 kpc. In Figure 7, the stars in each bin are plotted in Z –V_Z phase space, color coded by median(V_ϕ).

The phase spiral appears in each bin from 7–15 kpc; after that, the spiral disappears.

At first glance, the phase space spiral in Figure 7 appears to be sketched by the high median(V_ϕ) stars. In the gap of high median(V_ϕ) spiral, there are relatively low V_ϕ stars. From 7–15 kpc, each phase spiral map is similar to the last one, with a small but significant change. The R = 14–15 kpc bin is the last bin in which we can see a full circle of the phase spiral. There is still a high V_ϕ spiral standing out from the lower V_ϕ background, though it is blurry and not smooth. After R = 15 kpc, there are only a few high V_ϕ segments in the R = 15–16 kpc bins.

In the caption of Figure 7, the sample size of each bin is labeled. It is noticeable that sample size is only 5,000 in the R = 14–15 kpc bin. In the upper right panel of Figure 6, the phase spiral can only barely be seen when the sample size decreases to around 12,000. However, the phase spiral can still be seen in the panels of Figure 7 in the bin R = 14–15 kpc. This is because the phase spiral is in the range of − 0.5 < Z < 0.5 kpc in the R = 8 kpc bin, and the interval between the different phase spiral circles is narrow so the phase spiral is quite easily blurred. In the R > 12 kpc bin, the phase spiral is spread over a wider distance from the galactic plane, − 3 < Z < 3 kpc, with a different distribution. The spiral is still visible because the K giants are more spread out in Z when R is larger.

From 7–15 kpc, the winding of the phase spirals becomes looser with increasing R. Moreover, in the 7 < R < 15 kpc range, the phase space gets narrower (compression) in the V_Z direction and wider (stretching) in the Z direction with each increasing R bin. For example, in the bins from 10 < R < 11 kpc to 13 < R < 14 kpc, the high V_ϕ spiral spreads from the 2 kpc extent to 3.5 kpc in the Z dimension. This is consistent with the simulations of Laporte et al. (2019) and Bland-Hawthorn et al. (2019), who explained this behavior as due to the weaker potential with increasing R. As R increases, the stars travel further in the Z direction and move more slowly.

5.2.3. Comparing the V_ϕ Phase Spiral in LAMOST K Giants with Previous Studies

With Figure 7, we showed that the LAMOST K giants can be used to trace the global features of the phase spiral, which change with increasing R. The global features of the phase spirals are consistent with the theoretical predictions, but the detailed features of the phase spiral cannot be fully explored with LAMOST K giants data. Here we try to corroborate the characteristics of the V_ϕ phase spiral in LAMOST K giants by comparing them with Gaia DR2 data within R < 10 kpc.

Laporte et al. (2019) constructed the V_ϕ phase spiral at R = 6, 8, and 10 kpc with the Gaia DR2 stars in Figure 14 in their paper. Comparing their phase spiral with that of the R = 8–11 kpc bin of LAMOST K giants, the shape of the outline is similar and we also can see the outer circle tail from our data, but the inner circle data is also blurred into big centroid, and it is clear that the phase of the phase spiral in the Gaia DR2 data changes with galactocentric radius. The number of phase spiral circles also changes with galactic radius.

Wang et al. (2019) studied the shape of the phase spiral as function of R using stars in common between the LAMOST and Gaia samples, in the range of R < 12 kpc. They also see different phases of spiral shapes with different R. Before R = 9.34 kpc, there are 2.5 spirals and after that radius there is only one circle in Wang et al. (2019). From Figure 7 of this work, we also can see the phase of spirals changing with R, though the phase spiral is blurry. In the first three bins, with 7 < R < 10 kpc, the spiral warp leads toward (Z, V_Z) = (0 kpc, 0 kms⁻¹), where we see a centroid with a strong and wide tail protruding in the clockwise direction in the first and second quadrants. In the 10 < R < 11 kpc bin, the spiral is composed of a big centroid and a long, narrow "tail" in the second and third quadrants. The centroid fades with increasing R, and the centroid deviates from (Z, V_Z) = (0 kpc, 0 km s⁻¹) to the location of Z < 0 kpc.

From Figure 14 of Laporte et al. (2019), the tail of the phase spiral splits when R > 8 kpc. For example, the lower middle panel of Figure 14 of Laporte et al. (2019) shows a branching of spiral in the range of − 1 < Z < 1 kpc, 30 < V_Z < 60 km s⁻¹ for stars in Gaia DR2 data that are within 1 kpc of the Sun. We see similar branching in the 8 < R < 9 kpc bin of Figure 7; at the position of (Z, V_Z) = (0.7 kpc, 30 km s⁻¹), the spiral separates into several branches. We also see that the tail of phase spiral bifurcating in the third quadrant in the bin 10 < R < 11 kpc. We barely see that the tail of phase spiral bifurcating in the fourth quadrant in the 11 < R < 14 kpc bin. This means that either Z – V_Z space is not the correct projection space in which to count the wraps (Laporte et al. 2019), or the simple harmonic oscillation becomes chaotic, or it is not the simple harmonic oscillation but has higher modes.

5.2.4. Cautions

The phase spirals observed in LAMOST K giants are blurry and not smooth. The smaller sample size and larger distance error compared with the Gaia DR2 data explains part of the blurring, but there is another reason for this blurring effect. When we construct the phase spiral, the adopted bin in R – Z space is ΔR, ΔZ = (1 kpc, 0.125 kpc) when R < 12 kpc, ΔR, ΔZ = (1 kpc, 0.15 kpc) when 12 < R < 15 kpc, ΔR, ΔZ = (1 kpc, 0.2 kpc) when 15 < R < 16 kpc, ΔR, ΔZ = (1 kpc, 0.25 kpc) when 16 < R < 18 kpc, and ΔR, ΔZ = (1 kpc, 0.3 kpc) when 18 < R < 19 kpc. The ΔZ bin size is larger with increasing R in order to get enough stars in the sample, since the number of stars per radial increment decreases with increasing R. When we construct the phase spiral, there is an underlying assumption that there is no selection effect inside the (R, Z) bin and the sampled V_ϕ distributions are not biased. This assumption deviates more from the real situation when the bin size is larger. In the R < 10 kpc bin, the phase spiral from LAMOST K giants reproduces most of the features of phase spiral from Gaia DR2, but the properties of the phase spiral beyond this range still need to be confirmed by a larger and more precise data set.

5.3. Connection between the V_ϕ Phase Spiral and Kinematic Substructures

In order to study the relevance between the V_ϕ kinematic feature in R − Z and phase space, we rotate the panels from Figure 7 clockwise and arrange them horizontally in order of increasing R (Figure 8). Then, we compare Figure 8 with the median(V_ϕ) panel of Figure 4.

(i) In Figure 8, the centroid of the phase spiral becomes less and less pronounced with increasing R. This is partially due to the smaller number of K giants with increasing R, and partially because of the lower restoring force which causes the phase spiral to be more spatially extended with larger R (Laporte et al. 2019).

The phase spiral is centered at (Z, V_Z) = (0 kpc, 0 km s⁻¹) when R = 8–9 kpc and (Z, V_Z) = (−0.5 kpc, 10 km s⁻¹) when R = 14 to 15 kpc. In Figure 5, the midplane of the disk as identified by the position of the minimum standard deviation of V_ϕ and V_Z. The mean values of the location of minimum V_ϕ standard deviation from Figure 5 in each 1 kpc R bin are also shown as black triangles in Figure 8. We found that the black triangles trace the centroid of the phase spiral. We can connect this with the observation that the median(V_ϕ) of the main branch decreases with R and bends south after R > 13 kpc in the median(V_ϕ) panel of Figure 4. That is consistent with the prediction that the phase spiral traces the bending wave (Laporte et al. 2019).

(ii) The north branch and south branch of Figure 4 have higher V_ϕ than the adjacent disk-like stars. The location of north branch and south branch are just at the location where the high median(V_ϕ) phase space spiral passes through this projection.

For example, in the R = 13–14 kpc bin the segment of the phase spiral within −50 < V_Z < 50, Z < 0 kpc has a small Z range; after integrating along V_Z it corresponds to a pronounced high V_ϕ locus at (R, Z) = (13.5, −2) kpc. If we imagine compressing each R bin in the V_Z direction, we see that the kinematic substructure of the south branch in the median(V_ϕ) panel of Figure 4, which stretches from (R, Z) = (12, −1.5) kpc to (15, −4) kpc, is equivalent to the location of the high median(V_ϕ) spiral at each radius.

Similarly, the integration of the segment of the spiral with −50 < V_Z < 0 km s⁻¹ and Z > 0 kpc, in bins in the range of 12 < R < 16 kpc, corresponds to the north branch in the median(V_ϕ) panel of Figure 4. From this comparison, we know the north branch and south branch kinematic substructures are a projection of the phase space spiral into the R − Z plane.

The peak lines of north branch and south branch in the median(V_ϕ) panel of Figure 4 are overplotted on Figure 8 by crosses to show the corresponding relationship.

(iii) There are different theories for how a flare is produced, such as secular evolution (Narayan & Jog 2002; Minchev et al. 2012, 2012, 2015) or the cumulative effect of interactions with passing dwarf galaxies (Kazantzidis et al. 2008; Villalobos & Helmi 2008; Laporte et al. 2018).

In Figure 4, the disk-like stars flare with larger R. In Figure 8, the edge of the flare aligns with the location of outer circle of phase space spiral, especially after 11 < R < 15 kpc. Since the phase spiral is most likely induced by a recent passage of the Sgr dSph through the disk, it is reasonable to infer that the last passage of the Sgr dSph contributed stars to the flare.

6. Median V_R Distribution

The second row of Figure 4 shows the distribution of the median and standard deviation of V_R as a function of position in the (R − Z) plane. In the first row of Figure 4, both the maps of median(V_ϕ) and that of σ(V_ϕ) show the distinction between disk-like stars and halo-like stars. There is no significant distinction between disk-like stars and halo-like stars in the V_R map. So in the median(V_R) panel of Figure 4, the boundary within which the bins of median(V_ϕ) panel of Figure 4 have a median(V_ϕ) >160 km s⁻¹ is labeled to guide the eyes. The region of disk-like stars defined by median(V_ϕ) lies between these dark green lines.

Observing the disk-like stars, we see the median(V_R) is positive when R < 8 kpc. The median(V_R) is slightly lower than 0 km s⁻¹ at R = 9 kpc. The median(V_R) is positive in the range of 10 < R < 13.5 kpc, ∣Z∣ < 3 kpc. Then the median(V_R) is negative again. To show it more clearly, the radial variations of V_R within ∣Z∣ < 1 kpc and ∣Z∣ < 3 kpc are plotted in Figure 9. The median V_R of the two positive V_R areas in the range of R < 8 kpc, ∣Z∣ < 3 and 10 < R < 13.5 kpc, ∣Z∣ < 3 kpc is about 5 km s⁻¹.

**Figure 9.** Radial variation of the median(V_R). The black squares show radial variation of median(V_R) within ∣Z∣ < 1 kpc. The red squares show the radial variation of median(V_R) within ∣Z∣ < 3 kpc at R < 15 kpc. From Figure 4, the median(V_R) is seriously influenced by stars of the Monoceros area after R > 15 kpc, especially at higher Z, so the red squares are limited to R < 15 kpc. The error bars show the standard deviation of the values of median(V_R) for subsamples obtained by 1000 times bootstrap.
Download figure:
Standard image High-resolution image

These V_R features are consistent previous studies. Using The Radial Velocity Experiment (RAVE) data, Siebert et al. (2011) were the first to find a negative V_R gradient in the range of 6 < R < 9 kpc; the negative gradient is consistent with that of the median(V_R) of disk-like LAMOST K giants; it is positive in the range of R < 8 kpc and negative in the range of 8 < R < 9 kpc near the galactic plane.

Tian et al. (2017) and Liu et al. (2017) observed the same negative V_R gradient when R < 9 kpc and another positive gradient after R > 9 kpc using RAVE-Tycho-Gaia Astrometric Solution (TGAS) data, LAMOST red clump data, and LAMOST RGB data around the anticenter direction. The positive gradient after R > 9 kpc is consistent with the positive V_R area in the range of 10 < R < 13.5 kpc and ∣Z∣ < 3 with LAMOST K giant tracers.

The Gaia Collaboration et al. (2018) showed the X − Y distribution of disk radial velocities in the range of 4 < R < 12 kpc using giant stars from Gaia DR2. Figure 10 of the Gaia Collaboration et al. (2018) shows that there are two positive V_R regions at 4 < R < 8 kpc and R > 10 kpc, which is also consistent with our result.

**Figure 10.** The V_ϕ distribution from the test particle simulation. Each panel shows the median(V_ϕ) distribution in Z − V_z space for one time and one range of galactocentric radius for one wedge of the test particle simulation with 315° − 20° < ϕ < 315° + 20°. The top row shows data within 8 < R < 9 kpc, the middle row shows 11 < R < 12 kpc, and the bottom row shows 14 < R < 15 kpc. From left to right the simulation times shown are t = 0.12, 0.2, 0.28, 0.36, and 0.4 Gyr. This figure shows that the oscillation propagates outwards with increasing R. The phase spiral gradually moves from smaller R to larger R after the impact time.
Download figure:
Standard image High-resolution image

López-Corredoira et al. (2019) showed that the V_R velocity of Apache Point Observatory Galactic Evolution Experiment (APOGEE) stars increased toward the anticenter in the R = 9–13 kpc range, reaching a maximum of 6 km s⁻¹. V_R decreased at larger R, crossing from positive V_R to negative V_R around R = 15 kpc, which is the position of the outer spiral arm.

The V_R gradients have variously been explained as a perturbation due to spiral arms (Gaia Collaboration et al. 2018; López-Corredoira et al. 2019) or a resonance with the bar (Tian et al. 2017). Cheng et al. (2019) reexamined the V_R distribution in the X − Y plane with LAMOST OB stars, over a larger range in Y. They found that the two positive stripes are not aligned with the spiral arm, and therefore suggest that a satellite impact such as the Sgr dSph would better explain the V_R ridges observed in the X − Y plane.

We were not able to clearly trace the V_R phase spiral, which is not as significant as the V_ϕ phase space spiral, using LAMOST K giants, so it is not included here; a larger amount data is needed to trace the V_R phase spiral.

7. Median V_Z Distribution

The third row of Figure 4 shows the median and standard deviation of V_Z as a function of (R − Z). The stars around the anticenter direction are near the node of the warp, where V_Z is expected to be positive. The midplane stars show positive median V_Z of about 5 km s⁻¹, which is consistent with the presence of a warp as detected from Gaia DR2 kinematics (Poggio et al. 2018; Wang et al. 2020).

From the third row of Figure 4, we can see that the area of disk-like stars with positive median(V_Z) bends to the south after 12 kpc, following the main branch as defined in Figure 4. The fact that both the locus of high median V_ϕ and the locus of positive V_Z bend southward suggests that the disk midplane bends toward the south outside of 12 kpc from the Galactic center.

8. Distribution of Metallicities Abundance

The fourth row of Figure 4 shows the metallicity distribution, [M/H], as a function of (R, Z) in our sample. [M/H] is obtained from LAMOST spectra using the SLAM pipeline (Zhang et al. 2020). The median([M/H]) panel in Figure 4 shows that the value of [M/H] is the maximum in the galactic plane and decreases with increasing ∣Z∣. Also, the flared area in the (R − Z) map, where we find stars with the kinematics of disk stars, also has disk-like metallicity; the median([M/H]) value for stars in the north branch, south branch, and Monoceros area is about −0.5 dex. Not only do these stars have a high [M/H], but they also have a small spread in metallicity compared to the halo-like stars in our sample, as illustrated in the σ[M/H] panel of Figure 4.

9. Comparison with Simulations

To better understand the observational data, we compare our results using LAMOST K giants with the results of both a test particle simulation and an N-body simulation (Laporte et al. 2018) of the Sgr dSph galaxy gravitationally interacting with a Milky Way disk. Test particle simulations have relatively high efficiency and less computational cost; because the model is highly simplified, it is easy to observe the direct influence of the intruder to the disk. The N-body simulation (Laporte et al. 2018) is more realistic, because it includes self-gravity, which cannot be ignored in many situations (Darling & Widrow 2019b). In addition, the N-body simulation includes the effects of the Sgr dSph passing through the disk more than once. In this section, we describe this test particle simulation and N-body simulation, and then compare observational data with the results.

9.1. Description of the Test Particle Simulation

We reproduced Binney & Schönrich's (2018) toy model with galpy (Bovy 2014). In the test particle simulation, the intruder passes perpendicularly through the disk only once, and the gravitational effect of the dwarf galaxy is calculated as an impulse to each of the bodies in the disk.

The galpy potential MKPotential2014, which is a convenient approximation for the Milky Way potential, is adopted. The potential MKPotential2014 includes three parts: a bulge model with spherical potentials that are derived from power-law density models, a disk model that follows the Miyamoto–Nagai potential, and a halo model that follows the Navarro–Frenk–White potential. The parameters and properties of MWPotential2014 are summarized in Table 1 of Bovy (2014). The virial mass of the Galaxy is 0.8 × 10¹² M_⊙. The mass of the disk is 6.8 × 10¹⁰ M_⊙. The scale length of the disk is 2.6 kpc.

The galpy distribution function (df.quasiisothermaldf) is adopted in this work. It is an approximately isothermal distribution function based on action angle variables. MKPotential2014 and df.quasiisothermaldf are basically self-consistent; the simulated disk is stable for several hundred million years of integration.

The dwarf galaxy is modeled as a point mass, simulated with the galpy Kepler potential. The influence of the Milky Way on the passing dwarf galaxy is not considered. The point mass is 2 × 10¹⁰ M_⊙, as in Binney & Schönrich (2018). It passes through the disk from the north side of the disk, starting from Z = 10 kpc. The speed of the point mass remains constant at 300 km s⁻¹. The influence of the dwarf galaxy disappears after 66 Myr, when it arrives on the south side of the disk at z = −10 kpc.

The time at the end of the impact is defined as t = 0.0 Gyr in our simulation. The impact starts at −0.066 Gyr and it ends at 0.0 Gyr. The orbit integration, however, starts at −0.5 Gyr (434 Myr before the impact) and ends at 0.5 Gyr (500 Myr after the impact).

The goal of this section is to search across the whole disk in the simulated result to look for similar kinematic distributions to those found in observational data in both R − Z space and Z − V_Z space during the orbit integration time. We observe each snapshot at each time grid in the directions ϕ = 0°, 180°, 45°,135°, 225°, and 315°, with a bin size of ±20°. Finally, the stars in the wedge with 315° – 20° < ϕ < 315° + 20° are selected to show the detail of kinematic features, because the phase space spiral several hundred million years after the impact in this direction is most similar to the observations. We will compare the observational data with simulation results in Section 9.3. If the simulation was an exact replica of the encounter that produced the phase space spiral, then the ϕ direction in the simulation that matches the data would tell us something about the place or time of the impact that caused the spiral. But since this is only a toy model, the correspondence is only expected to be approximate.

A detailed description of the results of the test particle simulation is available in the Appendix.

9.2. Description of Laporte et al.'s (2018) N-body Simulation

The N-body simulation of Laporte et al. (2018) is a high-precision simulation that considers the entire orbit of the Sgr dSph after falling into the virial radius (R₂₀₀) of the Milky Way, including all disk passages, focusing on the reaction and evolution of the disk to the Sgr dSph. The results of the Laporte et al. (2018) simulation reproduce the vertical density and kinematic oscillation observed in the solar neighborhood and also the ring-like structure in the outer disk.

Laporte et al. (2018) built four kinds of N-body simulation models (L1, L2, H1, and H2) with different initial mass and different radial profiles for the Sgr dSph. The L2 model can best reproduce the spatial distribution of stars in the outer disk and flaring (Thomas et al. 2019) as well as the amplitude of density residual in the solar neighborhood. In Laporte et al. (2019), the L2 model is adopted to explain the observed disk oscillation and phase spiral. We will use the simulation results for the L2 model in this work to compare with the observational data.

The L2 model has a progenitor mass of M200 = 6 × 10¹⁰ M_⊙ and is twice as concentrated as the mean of the mass-concentration relation (Gao et al. 2008), with c = 28. In the L2 model, the Sgr dSph travels around or passes through the Galaxy five times from 5 Gyr ago to the present day. Each passing can produce a disk oscillation and ring-like structures. It can also erase the imprint of the previous passing to a large degree, depending on a trade-off between the rest mass of the Sgr dSph and the galactocentric distance of the impact position. Each passing, the Sgr dSph losses mass, reducing the influence of the satellite on the galactic disk. At the same time, the Sgr dSph gets nearer to the Galactic center, increasing the influence. These two effects compete with each other. From the simulation, the disk passage occurring at present day has a trivial effect on the disk, because the current mass of main body of Sgr dSph is only 10⁹ M_⊙. While the last passage happened more than 0.5 Gyr ago, is important because it reset the velocity distribution that we observe today in both physical space and phase space.

In Laporte et al. (2019), the N-body simulation results are compared with Gaia data in the range of 6 < R < 10 kpc mainly in X − Y space and Z − V_Z space. We will compare with the LAMOST K giants in larger volume of the Galaxy, 8 < R < 20 kpc in R − Z space and Z − V_Z phase space.

The time window of the last Sgr dSph impact is about 0.46–0.8 Gyr ago (Laporte et al. 2019), so we observe snapshots at 0.46, 0.63, and 0.8 Gyr in the directions ϕ = 0°, 180°, 45°, 135°, 225°, and 315°, with a bin size of ±20°. For the N-body simulation, the V_ϕ phase spiral is most similar to that of observational data in the wedge with 225° −20° < ϕ < 225° +20°.

The definition of coordinates used by Laporte et al. (2019) is the same as the test particle simulation. The position of the Sun is (X, Y) = (−8, 0) kpc. The azimuth angle ϕ is 0 in the direction of the negative X-axis. ϕ increases in the direction of disk star rotation.

9.3. Comparison of the Observational Data to the Results of Simulations

In Section 9.3.1, we compare the observed phase spiral with both the test particle and the N-body simulations (Laporte et al. 2018); each simulation provides different insights into the dynamical interaction. In Section 9.3.2, other observational kinematic features in the data are compared with the test particle simulation only.

9.3.1. Comparison of the Observed and Simulated Phase Spirals

Both simulations can qualitatively reproduce the V_ϕ phase space spiral. For the test particle simulation, the wedge with 315° − 20° < ϕ < 315° + 20° is selected as the most similar to the data. For the N-body simulation, the wedge with 225° − 20° < ϕ < 225° + 20° is selected as the most similar.

For the test particle simulation, in the direction of 315° − 20° < ϕ < 315° + 20°, Figure 10 shows the median(V_ϕ) distribution in Z − V_Z space in three different galactic radii and five different times between 120 and 400 Myr after the impact.

Figure 10 shows that the V_ϕ phase spiral appears first at smaller galactic radii, and moves to larger radii with time. In this simulation, the phase spiral starts to appear in the bin with 8 < R < 9 kpc when t = 120 Myr after the impact. In the panel with (R, t) = (8 kpc, 120 Myr), the phase of the modeled spiral is consistent with that of the observed spiral; the phase space spiral is centered at (Z, V_Z) = (0 kpc, 0 km s⁻¹) and trails in the counterclockwise direction. At t = 200 Myr after the impact, the phase space spiral starts to appear at larger R, where R = 14–15 kpc.

Figure 10 also shows that as the phase spiral appears at larger radii, it disappears at smaller radii. The phase spiral has already disappeared from the 8 < R < 12 kpc region by t = 360 Myr after the impact. The phase spiral is just fading from the 14 < R < 15 kpc region when t = 400 Myr.

The fact that the phase spiral moves to successively larger values of R is consistent with a perturbation that is propagating outwards. In Figure 10, the low median(V_ϕ) spiral appears at R = 11–12 kpc at t = 200 Myr. Then the low median(V_ϕ) spiral appears at R = 14–15 kpc at t = 280 Myr. This is consistent with Figure A1 which shows that in the direction of 315° − 20° < ϕ < 315° + 20°, the low median(V_ϕ) ring propagates outwards from R = 10–24 kpc as the time evolves from t = 180–400 Myr after the impact. Similarly, the high median(V_ϕ) phase spiral appears at radius R = 11–12 kpc at t = 280 Myr, and then moves outward to R = 14–15 kpc when t = 360 Myr.

Figure 11 shows the phase spiral distribution over the full range of galactocentric radius at a point t = 0.21 Gyr after the encounter with the dwarf galaxy. Note that the phase space spiral can be seen over a wide range of radius (9–15 kpc).

In the N-body simulation result, the phase space spiral matched the observational phase space spiral in the solar neighborhood 0.4–0.8 Gyr after the Sgr dSph impact (see Figure 10 of Laporte et al. 2019). In this time period, the distance range in which the phase space spiral appears changes from 8 < R < 15 kpc to 8 < R < 20 kpc, depending on azimuth angle. Figure 12 shows an example of phase spiral of each R bin of the azimuth slice with 225° − 20° < ϕ < 225° + 20° in the N-body simulation result when t = 0.8 Gyr.

**Figure 12.** The result of the N-body simulation (Laporte et al. 2019). The V_ϕ distribution in Z − V_Z phase space of stars in the range of 225° − 20° < ϕ < 225° − 20° at a time of 0.8 Gyr.
Download figure:
Standard image High-resolution image

The phase space spirals produced by both the test particle simulation and the N-body simulation qualitatively match the observational phase spirals. Similar to the behavior of the observational phase space spiral, the phase space spirals in simulations also are more extended in the direction of Z and more contracted in the direction of V_Z as R increases. The winding of phase space spiral gets looser with increasing R. In Figure 11 and 12, there are two different phase space spiral shapes, one before and one after R = 12 kpc. This finding is also similar to the observational phase spirals.

The most significant difference between the results of the test particle simulation and the N-body simulation is the persistence timescale of phase space spirals. In the test particle simulation, the phase space spiral only exists for a short time; it survived only 400 Myr within the range of R < 15 kpc. But in the Laporte et al. (2019) N-body simulation (Figure 10 of their paper), the phase space spiral survived more than 800 Myr within the same distance range. This is consistent with the conclusion of Darling & Widrow (2019a), who compared the responses of model disks with and without self-gravity to the same perturbation; they found the bending wave can last 1 Gyr in the simulated disk with self-gravity (as is present in the N-body simulation) while the bending wave damps out within 500 Myr in the simulated disk without self-gravity (similar to the test particle simulation).

At galactocentric distances larger than 16 kpc, the complete phase space spiral is not seen in the result of N-body simulations in the direction of 225° − 20° < ϕ < 225° + 20°; only "hatched chunks" are visible in Z − V_Z space. This is because the phase space spiral has been refreshed by the last impact of Sgr dSph in the inner disk, while in the outer disk the phase space shows the result of perturbations from multiple impacts (Laporte et al. 2019). The hatched chunks are not seen in the results of the test particle simulation outside of 16 kpc. This is because there is only one time impact in the test particle simulation.

9.3.2. Comparison of the Other Observed Kinematic Features with the Test Particle Simulation

(i) The test particle simulation reproduces vertical disk oscillations several hundred million years after the impact. For example, Figure 13 shows the median Z distribution in the X − Y plane at time t = 0.21 Gyr after the impact in the test particle simulation. One sees an oscillating disk around ϕ = 270° (the direction of the negative Y axis). The disk bends to the south at R = 13 kpc, then bends to the north at R = 17 kpc, with an oscillation amplitude of about 300 pc.

(ii) The test particle simulations can reproduce V_R ripples similar to those found by Cheng et al. (2019) or the V_R ripple in median(V_R) panel of Figure 4 of this work. Cheng et al. (2019) studied the kinematic distribution of O, B stars in the X − Y plane. They found a ripple pattern: V_R is −8 km s⁻¹ at R = 9 kpc, 0 km s⁻¹ at R = 12 kpc, and −10 km s⁻¹ when R > 13 kpc. This ripple is similar to the V_R distribution in the R − Z space of Figure 4 and the V_R distribution in R of Figure 9 of our paper. In Figure 4 and Figure 9, we see that the median(V_R) of disk-like stars is positive when R < 8 kpc, dips lower at R = 9 kpc, and is positive again when R > 10 kpc.

In the test particle simulation wedge with 315° − 20° < ϕ < 315° + 20°, the V_R ripple pattern appears after t = 0.14 Gyr. The upper right panel of Figure 14 shows an example V_R distribution for 315° − 20° < ϕ < 315° + 20° in R − Z space when t = 0.24 Gyr. The V_R ripple is apparent in this panel. The median(V_R) is about 10 km s⁻¹ when 2.5 < R < 4 kpc, about −15 km s⁻¹ when 5 < R < 7.5 kpc, and larger than 10 km s⁻¹ when 7.5 < R < 14 kpc. This ripple is illustrated more clearly in Figure 15, which shows the V_R variation of stars within ∣Z∣ < 1 kpc as a function of galactocentric radius. Comparing Figure 15 with Figure 9, the two positive V_R peaks in the test particle simulation result are quite similar to those in the observations. The full sequence of snapshots of the V_R distribution in R − Z space is shown in Figure A6 of the Appendix.

**Figure 15.** Ripple in median(V_R) as a function of R within ∣Z∣ < 1 kpc and azimuthal angle 315° −20° < ϕ < 315° + 20° for the test particle simulation when t = 0.24 Gyr. This variation is comparable to the ripple in the median(V_R) in the data, as shown in Figure 9.
Download figure:
Standard image High-resolution image

(iii) The test particle simulation can reproduce the high median V_ϕ substructure in R − Z space, similar to the north branch and south branch defined in Section 5.1. If we study the phase space spiral in the simulation results, we also find a correspondence between the projection of the phase space spiral in the R − Z space and the high median V_ϕ structures.

Figure 16 shows an example of the median V_ϕ distribution in the R − Z plane (the upper panel) and the V_ϕ phase space spiral in the Z − V_Z plane for each R bin (the lower panel) of snapshot of the test particle simulation with 315° − 20° < ϕ < 315° + 20° and t = 0.21 Gyr. The panels of the phase space spiral map are rotated clockwise and lined up in the direction of increasing R. The upper panel of Figure 16 shows the high median V_ϕ substructure in the range of 10 < R < 13 kpc has Z around 1.5 kpc, which is quite similar to the north branch of Figure 4. There is a less well populated south branch with Z ∼ −2 kpc in that same distance range. The lower panel of Figure 16 shows that the high median V_ϕ branch is equivalent to the projection of the phase space spirals in the 10 < R < 13 kpc bins.

(iv) The test particle simulation exhibits a transient flare that is excited by the Sgr dSph impact. The upper panel of Figure 16 shows an example of a flare excited by a dwarf galaxy impact in the test particle simulation. Comparing the upper and lower panels, we see that the boundary of the flare is defined by the outer ring of the phase space spiral, just as in Figure 8. Figures A4, A6, and A7 show that the flares grow as a function of R, with the oscillation propagating outward. From Figure A4, we can see that there is no flare before the Sgr dSph impact.

(v) The test particle simulations produce Monoceros-like substructures. For the test particle simulation, Figure 14 shows kinematically cold substucture at (R, Z) = (15, 2) kpc with median V_ϕ = 160 km s⁻¹, σ V_ϕ = 15 km s⁻¹, and median V_R = 20 km s⁻¹. Figure A4 shows that this substructure is puffer and has moved to higher radius at t = 0.46 Gyr.

We use the test particle simulation to explore the particle dynamics qualitatively, but the test particle simulation does not include self-gravity, which is important in reducing substructure and maintaining waves over long periods of time (Darling & Widrow 2019a). These concerns are mitigated by the fact that we are only looking at times within a few hundred million years of the impact, but the test particle results in this section should be verified in the future with full N-body simulations.

9.3.3. Summary of Simulation Results

Comparing observations with simulations, the Sgr dSph impact can qualitatively explain most of the observed disk substructure phenomena. From the test particle and N-body simulations, the influenced stars are rotated and dragged into rings, and show a phase space spiral several hundred Myr after the impact. The phase space spiral first appears at small R, and then gradually moves to larger R because the azimuthal frequency decreases with R. Spirals wind up quickly after a few orbital periods.

The passage of the Sgr dSph through the disk can produce the observed substructure. It can produce vertical and radial waves (Gómez et al. 2013; D'Onghia et al. 2016; Laporte et al. 2019); V_R, V_ϕ ripples (D'Onghia et al. 2016); (Laporte et al. 2018) phase space spirals (Binney & Schönrich 2018; Bland-Hawthorn et al. 2019; Laporte et al. 2019); and stream-like rings like the Monoceros ring. In summary, the Sgr dSph passing through the disk is likely the reason for the disk substructures, as was reported in Laporte et al. (2018). We show here in great detail how that substructure, seen in many different phase space projections, is related.

10. Discussion

10.1. Comparison with Previous Studies

Six-dimensional phase space information from Gaia DR2 allowed the local sample of stars to be studied in many projections. Antoja et al. (2018) studied the V_ϕ and V_r distributions in Z − V_Z phase space and discovered the phase space spiral. Laporte et al. (2019) found the number counts in the Z − V_Z phase space spiral. And they found that the vertical density oscillation (Widrow et al. 2012; Yanny & Gardner 2013) is consistent with the 1D projection of the number counts for the spiral in the Z direction. Schönrich & Dehnen (2018) studied the W versus L_Z distribution with Gaia-TAGS data, and found that a feature of the warp is that W is positive in the direction of the anticenter. Also, they found that there is an oscillation in the increase of W with L_Z. They found that the W − V_ϕ oscillation is just the projection of a spiral in Z − V_Z phase space. Cheng et al. (2019) study LAMOST K giants in the X − Y plane. From our observational data and previous study, we know the dominant features of the velocity field in R − Z, X − Y, Z − V_Z, and R − V_ϕ phase space; the density distribution in the solar neighborhood along Z and the density distribution in Z − V_Z phase space are different projections of the same kinematic feature imprinted by the same perturbation.

From our observational data and previous study, the bar is ruled out as the main reason for phase space spirals. This is because (i) The observational phase spiral in LAMOST K giants is apparent 7–15 kpc from the Galactic center, while the phase spiral produced by bar buckling is apparent only 4–10 kpc from the Galactic center, even 4 Gyr after the buckling (Khoperskov et al. 2019). (ii) Stars with different ages coexist in the same phase space spiral at the same R (Tian et al. 2018; Laporte et al. 2019). (iii) From the simulations, the bar cannot produce the observed size of the bending wave (Monari et al. 2015). (iv) A pattern speed of 60 km s⁻¹ kpc⁻¹ is required to fit the observed radial variation of the median(V_R) in the test particle simulation of Liu et al. (2018), which is large compared with recent measurements of 34–47 km s⁻¹ kpc⁻¹ for the pattern speed (Bland-Hawthorn & Gerhard 2016).

A spiral structure is also ruled out as the main formation mechanism for the phase space spiral. If the heavy spiral arms produce ripples in the disk, the ripple should follow those spiral arms (Siebert et al. 2012; Debattista 2014; Faure et al. 2014; Monari et al. 2016; Gaia Collaboration et al. 2018). However, in the larger volume of LAMOST data we see that the ripples do not follow the spiral arms (Cheng et al. 2019).

In the Milky Way disk, there are still kinematic substructures which cannot be explained in detail as a projection of the high V_ϕ phase spiral, and which need more precise data and more research to explore. For example, the "south middle opposite" found by Wang et al. (2019) is located in an apparent gap in the phase spiral. In addition, Carrillo et al. (2019) found that there is a strong azimuthal gradient in V_R, which is not predicted by the phase space spiral.

10.2. The Monoceros Substructure

The Monoceros overdensity is characterized by a velocity streaming feature. The stars in the Monoceros area have disk-like metallicity ([M/H] ≈ −0.5 dex).

Morganson et al. (2016) estimated that the mass of the stars on the north side of the disk in the Monoceros overdensity is 4 × 10⁶ M_⊙, and the mass of stars in "south Monoceros," which is identified as a distinct structure from the south middle structure of Xu et al. (2015), is 4 − 4.8 × 10⁷ M_⊙. The total stellar mass of the Monoceros overdensity is less than 5.2 × 10⁷ M_⊙, even when the south Monoceros structure is included. However, the simulation of Grebel (2005) demonstrated that a dwarf galaxy stellar mass should be at least 10⁸ M_⊙ to produce [Fe/H] as high as −0.5 dex, which is an index of the total metallicity. This metallicity estimate for LAMOST K giants is additional evidence pointing to the identification of the Monoceros overdensity as being part of the disk rather than an accreted satellite.

Laporte et al. (2020) also found that the Monoceros ring is an extension of the outer disk, and not accreted tidal debris, from a more detailed chemical composition of the stars. They detect the iron abundance and the α-abundance of the Monoceros ring and ACS with APOGEE data, and found that both of the structures have abundances in the range of −0.8 < [Fe/H] < −0.3, 0 < [Mg/Fe] < 0.15. Stars with these abundances are at the metal-poor end of the thin disk branch in the map of [Fe/H] versus [Mg/Fe] in Figure 6 of Laporte et al. (2019).

J. Li et al. (2019, in preparation) showed the same result using LAMOST K giant and M giant stars to trace the Galactic anticenter substructure (GASS), which includes Monoceros ring stars. They found that the mean and dispersion of the GASS iron abundance are −0.56 and 0.5 dex, respectively. The α abundance is concentrated in the range of 0–0.2 dex. The lower panel of Figure 7 in J. Li et al. (2019, in preparation) also shows a evidence that the dwarf galaxies of the Milky Way and GASS do not occupy the same region of the [α/Fe] versus [Fe/H] plane. In addition, the high metallicity of the Monoceros ring is inconsistent with a dwarf galaxy of the mass of Fornax; if one uses the mass–metallicity relation for the dwarf spheroidal from Kirby et al. (2013), the implied stellar mass of a galaxy of metallicity [M/H] ∼ −0.5 would be 10¹⁰ M_⊙.

The metallicity suggests that the substructures are very likely connected to the disk rather than with debris from the accretion event.

11. Conclusion

We validate the Gaia proper motions further than R = 20 kpc, where R is the distance from the Galactic center. The distance and radial velocities of the LAMOST K giants are validated to similar galactocentric distance. The combination of Gaia proper motions, LAMOST radial velocities, and distances calculated from LAMOST spectra allow us to study the velocity field within the range of 6 < R < 20 kpc. Because LAMOST K giants are mostly observed within ϕ = ±20° of the anticenter direction, it is most convenient to project velocities in R − Z space, while slicing the data in R. We study the velocity distribution in the R − Z plane, and compare it with projections in Z − V_Z phase space. From the data we learn:

(1) There are plenty of features in the map of the V_ϕ distribution in the R − Z plane. The main characteristics are (i) the midplane as determined by the location of the minimum standard deviation V_ϕ is oscillating in the Z direction, (ii) the outer disk bends southward, (iii) the oscillation in the centroid of the phase spiral traces the large scale bulk motion of the disk (see Figure 8), and (iv) the high V_ϕ locus is trifurcated after R = 13 kpc.

(2) The V_ϕ spiral patterns in the Z − V_Z plane extend from R = 7–15 kpc. Spirals are more tightly wound at smaller R and more loosely wound at larger R. In addition, spirals more elongated along Z and are more squashed along V_Z with increasing R. The simulations predict the same trends due to decreasing of self-gravity and longer dynamical timescale with increasing R (Bland-Hawthorn et al. 2019; Laporte et al. 2019). In the observational data, after R = 15 kpc, there are no complete spirals and the Z − V_Z phase space is full of hatched chunks. Similar to the prediction of Laporte et al. (2019), the imprint of multiple passages of the dwarf galaxy through the disk are present in the stars of the outer disk.

(3) The most significant bulk motion is consistent with the 1D projection of the V_ϕ phase space spiral integrated along V_Z in each radial bin.

(4) The observed flaring of the disk is substantial. The flare in the outer disk is even thicker than the thick disk in the solar neighborhood. From Figure 8, the boundary of the flare is defined by the extension along Z of the phase space spiral. Since the Sgr dSph was found to induce the phase space spiral, this suggests that the recent impact of the Sgr dSph also contributes to the flare.

We cannot rule out other possible mechanisms (such as Narayan & Jog (2002), Minchev et al. (2012)) to produce a flare except the explanation of Sagittarius passage.

(5) The test particle simulation and N-body simulation can reproduce the above observational characteristics qualitatively, including the phase spiral and bulk motions.

In summary, we observed the one-to-one relationship between the phase spiral and the most significant bulk motion of the disk stars. Because the most likely origin of the phase spiral is the last impact of the Sgr dSph, the bulk motion of the disk stars is likely the result of the same impact.

This work is supported by National Key R&D Program of China No. 2019YFA0405500. C.L. thanks the National Natural Science Foundation of China (NSFC) with grant No. 11835057. H.F.W. is supported by the LAMOST Fellow project, National Key Basic R&D Program of China via 2019YFA0405500 and funded by China Postdoctoral Science Foundation via grants 2019M653504 and 2020T130563, Yunnan province postdoctoral Directed culture Foundation, and the Cultivation Project for LAMOST Scientific Payoff and Research Achievement of CAMS-CAS and the National Natural Science Foundation of China grant No. 12003027. H.J.N. is supported by the US National Science Foundation grant AST 19-08653. X.F. acknowledges the support of China Postdoctoral Science Foundation No. 2020M670023, the National Natural Science Foundation of China under grant No. 11973001, and the National Key R&D Program of China No. 2019YFA0405504. This work was supported in part by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan.

Guoshoujing Telescope (the LAMOST) is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the project has been provided by the National Development and Reform Commission. LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences.

This 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.

This work benefited from the International Space Science Institute (ISSI/ISSI-BJ) in Bern and Beijing, thanks to the funding of the team Chemical abundances in the ISM: the litmus test of stellar IMF variations in galaxies across cosmic time (PI D. Romano and Z-Y. Zhang). Thanks to Yang Chengqun who ran the Carlin's code to get the distance of LAMOST K giants.

Appendix A: Detailed Results of the Test Particle Simulation

Figures A1, A2, and A3 show the median(V_ϕ), median(V_R), and median(V_Z) distributions, respectively, of the simulation results in the X − Y plane. The Sun is located at (X, Y) = (−8, 0) kpc. The disk stars rotate clockwise. The azimuth angle (ϕ) is defined so that ϕ = 0 in the direction from the Galactic center to the Sun and ϕ increases in the clockwise direction. The toy model is simplistic, but allows us to study the effect of the impact itself, without complications from other physical effects.

**Figure A1.** The median(V_ϕ) distribution on the X − Y plane in a series of snapshots from the test particle simulation. The time after the impact is labeled on top of each panel. The wedge in the lower right panel shows the direction of 315° −20° < ϕ < 315° + 20°. Arcs at 10, 20, and 30 kpc from the Galactic center are shown in the wedge.
Download figure:
Standard image High-resolution image

**Figure A2.** The V_r distribution on the X − Y plane in the test particle simulation. The time after the impact is labeled on top of each panel. The wedge in the lower right panel shows the direction of 315° −20° < ϕ < 315° + 20°. Arcs at 10, 20, and 30 kpc from the Galactic center are shown in the wedge.
Download figure:
Standard image High-resolution image

The motion of stars on the X − Y plane are described by median(V_ϕ) and median(V_R). From Figures A1 and A2, we can see that the disk feels the disturbance right after the impact begins. The stars are accelerated toward the intruder; stars moving toward the intruder develop a faster speed and stars moving away from the intruder are slowed. Inside the radius of the position of impact, the median(V_R) > 0. Outside that radius, the median(V_R) < 0. After the impact, the high median(V_ϕ) accelerated stars and the low median(V_ϕ) of decelerated stars are dragged into rings. The high median(V_ϕ) stars overtake the low median(V_ϕ) stars at about 0.02 Gyr. Then the high median(V_ϕ) stars are separated into two parts with median(V_R) > 0 and median(V_R) < 0. The part with median(V_R) < 0 moves inward and rolls up, and the part with median(V_R) > 0 becomes a ring at about 0.1 Gyr. The low median(V_ϕ) forms a ring in between two high median(V_ϕ) rings at about 0.18 Gyr. The rings generally propagate outwards.

The distribution of median(V_Z) on the X − Y plane of Figure A3 shows us that the stars around the point of impact have positive median(V_Z) when the intruder is on the north side of the disk. The median(V_Z) of stars around the point of impact is negative, while the intruder is on the south side of the disk. The stars that gain an extra V_Z deviate from their original planar orbit to oscillate above and below the disk. The influenced stars are then dragged to rings as the disk rotates, and the rings propagate outwards.

Figures A1–A3 show that the entire disk of stars is oscillating after several hundred million years. We observed snapshots in the direction of ϕ = 0°, 180°, 45°,135°, 225°, and 315°. The stars in the wedge of 315° − 20° < ϕ < 315° + 20° were selected to show the detailed of kinematic features in R − Z space (Figures A4–A7). This azimuthal angle was selected because the phase space spiral in this direction, several hundred million years after the impact, is most similar to the observed one. From the snapshots, we see the evolution of kinematic substructures induced by the impact of the passing dwarf galaxy. The position of the wedge is shown in the last panel of Figures A1–A3.

**Figure A4.** The median(V_ϕ) distribution of data with azimuthal angle 315° − 20° < ϕ < 315° + 20° in the R − Z plane in the test particle simulation. The time after impact is labeled on top of each panel.
Download figure:
Standard image High-resolution image

**Figure A5.** The standard deviation of the median(V_ϕ) distribution of data in the range of 315° − 20° < ϕ < 315° + 20° in the R − Z plane in the test particle simulation. The time after impact is labeled on top of each panel.
Download figure:
Standard image High-resolution image

**Figure A6.** The median(V_r) distribution of data in the range of 315° − 20° < ϕ < 315° + 20° on the R − Z plane in the test particle simulation. The time after impact is labeled on top of each panel.
Download figure:
Standard image High-resolution image

**Figure A7.** The median(V_z) distribution of data in the range of 315° − 20° < ϕ < 315° + 20° on the R − Z plane in the test particle simulation. The time after impact is labeled on top of each panel.
Download figure:
Standard image High-resolution image

Exploring the Perturbed Milky Way Disk and the Substructures of the Outer Disk

Article metrics

Permissions

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Sample Selection

3. Distance Calibration

4. Cylindrical Velocities: V_r, V_ϕ, and V_z