Radial Wave in the Galactic Disk: New Clues to Discriminate Different Perturbations

We now turn to test particle simulations involving internal and external perturbations separately, to examine the difference in the resulting phase shift $\Delta L_{Z}$ variation pattern of the $L_{Z}-\langle V_{R}\rangle$ wave. Compared to $N$-body simulations, test particle simulations neglect self-gravity, but do allow us to run more particles and test the effect of varying specific parameters on the orbit evolution.

We sample eight million disk particles from the quasiisothermal df (Carlberg & Sellwood, 1985; Binney, 2010; Binney & McMillan, 2011) implemented in agama (Vasiliev, 2019). These particles are distributed using a radial scale length of 3 kpc within the Milky Way potential model, MWPotential2014 (Bovy, 2015). The central velocity dispersion $\mathrm{(\sigma_{R}|_{R=0},\sigma_{Z}|_{R=0})=(90,110)\,km\,s^{-1}}$ gives a disk that is dynamically colder than the actual Milky Way disk. The scale lengths of the $\sigma_{R}$ and $\sigma_{Z}$ profiles are 6 and 7 kpc, respectively. The $L_{Z}-\langle V_{R}\rangle$ wave signal is more prominent, which simplifies the task of identifying extrema points for the $\Delta L_{Z}$ measurement. To investigate the impact of dynamical hotness on the $L_{Z}-\langle V_{R}\rangle$ wave morphology, we define those particles with $J_{Z}\mathrm{<8\,km\,s^{-1}\,kpc}$ as the dynamically “cold” population, and those with $J_{Z}\mathrm{>23\,km\,s^{-1}\,kpc}$ as the dynamically “hot” population. Under such division, the particle numbers of the cold and hot populations are roughly the same. Note that this criterion of orbital hotness is different from that used in the observation. For the observation, the criterion is determined empirically to enable sufficient extrema points and coherent $L_{Z}-\langle V_{R}\rangle$ wave patterns for phase shift ($\Delta L_{Z}$) measurement. This is mainly due to the differences in the distribution function between the observation and simulation, and also the sample selection bias and velocity measurement error in the observational data, which is not present in the simulation results. Since our main goal is not to quantitatively reproduce the observational $\Delta L_{Z}$ variation trend but to qualitatively understand the physical origin of the $\Delta L_{Z}$ feature, we choose different $J_{Z}$ bounds for the sample division in the observation and simulation.

The external perturber (Milky Way satellite) is a $2.5\times 10^{10}M_{\odot}$ Plummer sphere with a half-mass radius of 3 kpc. Its position $\vec{x}=(4,8,18)\,\mathrm{kpc}$ and velocity $\vec{v}=(-339,-44,76)\,\mathrm{km}\,\mathrm{s}^{-1}$ at the first pericenter passage are from the E1 model of de la Vega et al. (2015). Backward orbit integration for 1 Gyr in the MWPotential2014 potential using galpy’s ChandrasekharDynamicalFrictionForce routine to account for dynamical friction provides the initial condition for our simulation. At 2 Gyr, we reduce the mass and half-mass radius of the satellite to $1\times 10^{10}M_{\odot}$ and 1 kpc, respectively, to mimic the mass loss after the first pericenter passage. The second pericenter passage occurs at 3.1 Gyr when the satellite crosses the disk plane at $R=15\,\mathrm{kpc}$ with $\mathrm{V_{Z}\approx 300\,km\,s^{-1}}$. We adopt this simplified setup of the satellite mass loss for better simulation efficiency. A qualitatively similar prescription has been commonly used in previous works to model the satellite perturbation on the disc (e.g., Xu et al., 2020; Bennett & Bovy, 2021; Gandhi et al., 2022). Our tests, assuming a gradual mass loss history, produce qualitatively similar $L_{z}-\langle V_{R}\rangle$ wave features and $\Delta L_{Z}$ variation trend. The total integration time is set to 4 Gyr to cover these two pericenter passages. The full satellite orbit trajectory color coded by time is illustrated in Figure 2.

Using test particle simulations, we also investigate the influence of internal perturbations with the combination of a steadily rotating bar and two-armed transient spirals. The spiral arm potential reaches its maximal amplitude at 1200 Myr. To differentiate the role of different perturbers in generating the $\Delta L_{Z}$ variation trend, we complement this simulation with a test where the transient spirals act as the sole perturber.

In this test, the spiral potential has a higher amplitude and reaches the maximal value at the start. The two-armed transient spirals follow the default model described in Section 2.2 of Hunt et al. (2018) which follows the potential form given by Cox & Gómez (2002). The winding and decay of the spiral arms are configured by the CorotatingRotationWrapperPotential and GaussianAmplitudeWrapperPotential routines of galpy. Under this setup, spiral arms have a lifetime of about 100 Myr and approximately corotate with stars in the investigated radial range. The morphological parameters of the SpiralArmsPotential model are $\mathrm{N=2,R_{s}=0.3,C_{n}=1,H=0.125,\theta_{sp}=12^{\circ}}$ ($N,\,R_{s},\,C_{n},\,H,\,\mathrm{\theta_{sp}}$ represent the number of spiral arms, radial scale length, sinusoidal potential profile, scale height, and pitch angle, respectively). The bar is modeled as the 3D generalization of the Dehnen bar potential (Dehnen, 2000; Monari et al., 2016) with the same pattern speed (40 $\mathrm{km\,s^{-1}\,{kpc}^{-1}}$, invariable with time) and bar length (4.5 kpc) as the “fiducial” model presented in Trick et al. (2021).

The integration time is 2 Gyr for the case with a bar and 0.4 Gyr without a bar (spirals only). Particle orbits are integrated using the galpy code (Bovy, 2015). In the first 1 Gyr, the bar is the main driver of the phase space structure formation. Therefore, it can be used to illustrate the $\Delta L_{Z}$ variation due to the effect of the bar. A fully grown bar generates noticeable $L_{Z}-\langle V_{R}\rangle$ wave signals (extrema points for the phase shift measurement) only at the 2:1 and 1:1 OLR of the bar, which could make the number of measurable extrema points too limited to ascertain a clear $\Delta L_{Z}$ variation trend (See Figure 11 of Chiba et al., 2021). On the other hand, generating more extrema points in the $L_{Z}-\langle V_{R}\rangle$ wave requires higher-order Fourier modes from the bar potential (Monari et al., 2019) or a decrease in the pattern speed of the bar (Chiba et al., 2021), which is beyond the scope of this work. We exclude stars with an initial radius smaller than 1.5 kpc to save integration time. This will have little impact on the simulation results since the number of discarded particles capable of entering the investigated $L_{Z}$ range is negligible.

With the above setup, we divide each simulation snapshot into eight equally spaced azimuthal ranges and extract extrema points from the “cold” and “hot” waves. We focus on the $L_{Z}-\langle V_{R}\rangle$ wave pattern in the range of $L_{Z}\in[600,3000]\thinspace\mathrm{km\thinspace s^{-1}\thinspace kpc}$, divided into 165 equally spaced bins. The Gaussian kernel size used to smooth the $L_{Z}-\langle V_{R}\rangle$ wave for $\Delta L_{Z}$ measurement is $\sigma=4$. We compare the $\Delta L_{Z}$ variation trends in different azimuths to conclude a universal variation pattern. If there is none, we try to unveil the $\Delta L_{Z}$ variation feature existing in most azimuthal ranges, which may also be valuable for discriminative purposes. Due to the subtle and discrete nature of the measurable extrema points, conclusions drawn from the analysis on the $\Delta L_{Z}$ variation trend are reliable in the qualitative sense, but any quantitative conclusion should be treated with caution. Providing a quantitative match of the $\Delta L_{Z}$ variation curve between simulation results and observation data is difficult since the phase shift $\Delta L_{Z}$ is influenced by measurement settings, distribution function, lack of self-gravity, and other factors.

In the case of the transient spirals, the phase shift amplitude $\Delta L_{Z}$ exhibits no universal trends with $L_{Z}$ among different azimuthal ranges. As Figure 6 demonstrates, within the majority of azimuthal ranges, the phase shift amplitude $\Delta L_{Z}$ displays irregular oscillation with $L_{Z}$, with an amplitude between 30 and 80 $\mathrm{km\,s^{-1}\,kpc}$. We cannot observe any similar $\Delta L_{Z}$ variation pattern shared by both snapshots, either. Interestingly, the decreasing trend with $L_{Z}$ resembling the trend produced by external satellite perturbation exists for very few azimuthal ranges of the two snapshots. In such cases, the slope of the $\Delta L_{Z}$ curve is shallower than those generated from the satellite impact, which still allows us to distinguish from the $\Delta L_{Z}$ variation trend of external perturbation.

The disk evolution shown in Figure 7 is when the perturbations of the bar and transient spirals are combined. The central depression in density is the result of excluding the innermost particles from orbit integration. Before the transient spiral starts growing ($T\mathrm{<1000\,Myr}$), the disk dynamics are mainly influenced by the bar, as shown by the butterfly pattern in the $\langle V_{R}\rangle$ map. When the spiral potential reaches maximal amplitude at 1200 Myr, the two-armed spiral pattern becomes most prominent in both the density and $\langle V_{R}\rangle$ maps. The amplitude of $\langle V_{R}\rangle$ decreases as the spirals wrap up and decay. The bar resumes domination as the time approaches the end of our simulation.

As shown in the leftmost column of Figure 8, when only the bar potential is active, the $\Delta L_{Z}$ variation trends in different azimuthal ranges do not share common features. The $\Delta L_{Z}$ amplitude exhibits smaller dispersion than the transient spiral case and increases to a larger value at the $L_{Z}$ range below the corotation resonance. However, the wave shape consistency between the two subpopulations becomes worse within the lower $L_{Z}$ range, making this feature susceptible to systematic error.

After the spiral potential becomes active, although the overall $\Delta L_{Z}$ pattern illustrated in the three right columns of Figure 8 varies drastically with time, a common feature in the phase shift variation pattern emerges across most of the azimuthal ranges: a characteristic peak at $L_{Z}\sim 2000\,\mathrm{km\,s^{-1}\,kpc}$ with variable prominence. Interestingly, the location of the $\Delta L_{Z}$ peak coincides with the location of the 2:1 OLR of the bar (marked by a vertical blue dashed line) in most azimuthal ranges. In a minority of the azimuthal ranges, the $\Delta L_{Z}$ peak is less prominent or absent, which makes the $\Delta L_{Z}$ variation trend qualitatively similar to those produced by the transient spirals. Despite the caveats mentioned above, the sufficiently prominent $\Delta L_{Z}$ peak could provide the most likely range of the 2:1 OLR to constrain the bar pattern speed $\Omega_{b}$.

Due to the fundamental difference between internal and external perturbations, our toy model of radial phase mixing cannot account for the $\Delta L_{Z}$ variation trends generated by the bar and spirals, which may be caused by the dependence of resonance lines on dynamical hotness ($J_{Z}$). As illustrated in Figure 10 of Trick et al. (2021), the location of bar resonance (the so-called “ARL”, axisymmetric resonance line) shifts toward lower $L_{Z}$ for stars of higher $J_{Z}$. In the region of the bar’s 2:1 OLR, the majority of stars are forced toward higher $L_{Z}$, and the changes of $L_{Z}$ ($\delta L_{Z}[L_{Z,0}]$) are smaller than corotation resonance. Within corotation resonance, radial ($L_{Z}$) migration can proceed in both directions, which could partially cancel out the effect of the ARL shift. This is not the case for Lindblad resonances. Therefore, the bar’s 2:1 OLR can generate a distinct $\Delta L_{Z}$ peak in the presence of the spirals. No prominent $\Delta L_{Z}$ peaks emerge under transient spiral perturbation since most stars approximately corotate with the spirals and are influenced by corotation resonance. However, the explanation above cannot account for the absence of the $\Delta L_{Z}$ peak feature in the bar-only simulation. More comprehensive theoretical studies are needed to fully comprehend the complex interplay between spiral and bar perturbations in the future.

To briefly summarize, the $L_{Z}-\langle V_{R}\rangle$ wave pattern displays different dependence on dynamical hotness depending on the nature of dominant perturbations. Utilizing the amplitude of phase shift ($\Delta L_{Z}$) between the wave patterns of different dynamical hotness, we could quantify this dependence and associate its variation trend versus $L_{Z}$ with particular perturbation types. The external satellite generates a general trend of decreasing $\Delta L_{Z}$ in the lower $L_{Z}$ range in almost all azimuths during the whole time. Such universality does not exist in the case of the transient spirals or the bar. Under the joint perturbation of the bar and the transient spiral, $\Delta L_{Z}$ displays a characteristic peak around 2:1 OLR, which qualitatively reproduces the observed $\Delta L_{Z}$ feature. Controlled tests demonstrate that the internal perturbation (the combination of the bar and the spirals) is more favorable than the external one to reproduce the observed phase shift variation pattern.

Despite capturing the essential physical processes of radial phase mixing, our toy model presented in Section 3.2.2 is not exempt from limitations. Firstly, it outlines the backbone of the phase space substructure but falls short of providing a detailed phase space density distribution. As the system becomes fully phase mixed in the lower $L_{Z}$ range, the $\langle V_{R}\rangle$ amplitude should be close to zero, which is not seen in the toy model. Neglecting the azimuthal dependence of the external perturbation, our toy model does not reflect the azimuthal phase-mixing process as the model of Antoja et al. (2022) does. Employing the epicycle approximation, our model can only account for the $\Omega_{R}$ difference between cold and hot orbits, prohibiting us from explaining the larger $\Delta L_{Z}$ at the high $L_{Z}$ end.

It is also imperative to be cautious while interpreting the $\Delta L_{Z}$ variation trends generated by test particle simulations, given the absence of self-gravity effects. Future investigations employing tailored $N$-body simulations (e.g., Bland-Hawthorn & Tepper-García, 2021; Hunt et al., 2021) are required to validate our main results.

We should note that the $\Delta L_{Z}$ variation trend can only differentiate different perturbations in the probabilistic sense because not all azimuthal ranges exhibit the common $\Delta L_{Z}$ variation features unique to specific perturbations. Under the satellite impact, the decreasing $\Delta L_{Z}$ trend predicted by our toy model is less pronounced in some of the azimuthal ranges. In the case of bar and spiral arms, the 2:1 OLR of the bar does not generate the characteristic $\Delta L_{Z}$ peak in all azimuthal ranges. Therefore, broader azimuthal coverage of the Galactic disk in future surveys is essential for robust diagnostics of the perturbation history.

Friske & Schönrich (2019) measured the phase shift by fitting the observed data with the sinusoidal function and comparing the best-fit phase angles with each other. The unit of their phase shift is the degree, not the angular momentum used ($\mathrm{km\,s^{-1}\,kpc}$) in our work. While their approach is advantageous in utilizing complete wave information, it is susceptible to systematic errors due to its model-dependent nature. During fitting, they set the wavelength of the sine function to fixed values, essentially making their approach akin to peak finding in principle. Our method does not assume any waveform, but the results of peak finding could vary with the width of $L_{Z}$ bins. Both approaches have the issue of defining the phase shift when the wave shapes deviate significantly from the sinusoidal shape. New mathematical tools that can robustly measure the phase shift in such nonsinusoidal waves are indispensable for future research.

The comparison between the observational data and our test particle simulation results suggests that the $L_{Z}-\langle V_{R}\rangle$ wave is more likely to originate from internal perturbations rather than external satellites. In the case of external perturbation, extrema points with the most prominent $\Delta L_{Z}$ are located at both ends of the investigated $L_{Z}$ range, which disagrees with the key observational feature. In contrast, the characteristic $\Delta L_{Z}$ peak generated by the 2:1 OLR of the bar persists regardless of the evolution phase of the transient spiral arms. Therefore, if the observed phase shift variation trend does represent the trends within most of the azimuthal ranges, the combination of the bar and the transient spiral is more likely to be the dominant driving force. It is important to note that this conclusion is not definitive, as other scenarios, such as a decelerating bar (Chiba et al., 2021) have not been explored in this work. Therefore, the $\Delta L_{Z}$ variation trend is not a sufficient but necessary criterion to discriminate different perturbing mechanisms. Nevertheless, our study provides a novel approach to differentiating the roles of internal and external perturbations in sculpting the observed phase space substructures.

According to our simulation results, the combination of the bar and the transient spiral arms is more likely to be the dominant perturbation, where the peak of the phase shift amplitude $\Delta L_{Z}$ is associated with 2:1 OLR. This $\Delta L_{Z}$ diagnostic could pave a new way to constrain the pattern speed of the Galactic bar. Utilizing test particle simulations with bar as the sole perturber, Trick et al. (2021) (also in Mühlbauer & Dehnen, 2003) proposed that the $\langle V_{R}\rangle$ map in the $L_{Z}-J_{R}$ action space displaying the sign change is the signature of 2:1 OLR. However, owing to the interference of the transient spirals, the sign change in the $L_{Z}-\langle V_{R}\rangle$ wave pattern could also occur in $L_{Z}$ ranges other than specific types of resonances, as demonstrated by Hunt et al. (2019) in the $R-V_{\phi}$ projection. Our work could help break this degeneracy since it provides a viable way to detect the bar resonance signature in the presence of spirals. Among the four candidates of the 2:1 OLR of the bar proposed by Trick et al. (2021), the “Hat” and “Sirius” moving groups are the two closest to the $\Delta L_{Z}$ peak of the observational data, i.e. the 2:1 OLR of the bar. This gives a crude pattern speed estimate of $\sim 33-41$ $\mathrm{km\,s^{-1}\,kpc^{-1}}$, which favors the long/slow bar model. The above constraint is consistent with previous works using the solar neighborhood kinematics to indirectly infer the bar pattern speed (Binney, 2020; Chiba & Schönrich, 2021; Trick, 2022). The uncertainty of our estimation is greater than conventional methods because we can only give the most probable $L_{Z}$ range of the 2:1 OLR location from discrete extrema points.

Wang et al. (2020) classified the $R-V_{\phi}$ ridges into two types depending on whether they exhibit significant variation with stellar age in the $\langle V_{R}\rangle$ map. Coincidentally, the ridge displaying the most prominent variation with age roughly follows the constant $L_{Z}$ line at $\mathrm{2180\thinspace km\thinspace s^{-1}\thinspace kpc}$ where $\Delta L_{Z}$ exhibits a peak feature. Given the tendency of older stars to be dynamically hotter, this prominent variation could be interpreted through the 2:1 OLR of the Galactic bar, from the dynamical perspective.