The GALAH survey and Gaia DR2: Linking ridges, arches and vertical waves in the kinematics of the Milky Way

TL;DR

This study links Gaia DR2 star kinematic patterns with phase mixing of spiral arms, revealing ridges and arches as surfaces of constant energy, and demonstrates their formation through N-body simulations involving satellite perturbations.

Contribution

It shows that phase mixing of spiral arms explains the observed ridges and arches in Milky Way star kinematics, highlighting their energy and angular momentum properties and coupling between planar and vertical motions.

Findings

Ridges are present in multiple kinematic and chemical abundance maps.

Ridges correspond to surfaces of constant orbital energy or angular momentum.

Simulations show ridges can form in both isolated and satellite-perturbed galactic discs.

Abstract

Gaia DR2 has revealed new small-scale and large-scale patterns in the phase-space distribution of stars in the Milky Way. In cylindrical Galactic coordinates $(R,\phi,z)$, ridge-like structures can be seen in the \vphiR{} plane and asymmetric arch-like structures in the \vphivR{} plane. We show that the ridges are also clearly present when the third dimension of the \vphiR{} plane is represented by $\langle z \rangle$, $\langle V_z \rangle$, $\langle V_R \rangle$, $\langle$[Fe/H]$\rangle$ and $\langle[\alpha/{\rm Fe}]\rangle$. The maps suggest that stars along the ridges lie preferentially close to the Galactic midplane ($|z|<0.2$ kpc), and have metallicity and $\alpha$ elemental abundance similar to that of the Sun. We show that phase mixing of disrupting spiral arms can generate both the ridges and the arches. It also generates discrete groupings in orbital energy $-$ the ridges and arches are simply surfaces of constant energy. We identify 8 distinct ridges in the \gaia{} data: six of them have constant energy while two have constant angular momentum. Given that the signature is strongest for stars close to the plane, the presence of ridges in $\langle z \rangle$ and $\langle V_z \rangle$ suggests a coupling between planar and vertical directions. We demonstrate, using N-body simulations that such coupling can be generated both in isolated discs and in discs perturbed by an orbiting satellite like the Sagittarius dwarf galaxy.

Tables

Table 1. Table 1: Data quality cuts on GALAH DR2.

Selection	Comments
	-
0$\le$ Field id	Excludes data without proper selection
	function

Table 2. Table 2: Parameters for the isolated Galaxy (Model P). Column headers are as follows: M t subscript 𝑀 t M_{\rm t} := total mass ( 10 9 superscript 10 9 10^{9} M ⊙ subscript 𝑀 direct-product M_{\odot} ); r s subscript 𝑟 s r_{\rm s} := scalelength (kpc); r tr subscript 𝑟 tr r_{\rm tr} := truncation radius (kpc); N p subscript 𝑁 p N_{\rm p} := number of particles ( 10 6 superscript 10 6 10^{6} ).

	Profile	$M_t$	$r_s$	$r_{\mathrm{tr}}$	$N_p$
Galaxy
DM halo	H	${10}^3$	38.4	250	10
Bulge	H	9	0.7	4	1
Thick disc	MN	20	5.0a	20	2
Thin disc	Exp/Sech	28	3.0b	20	3

Table 3. Table 3: N-body models and properties of the perturber . Column headers are as follows– M tot subscript 𝑀 tot M_{\rm tot} : total mass ( 10 9 superscript 10 9 10^{9} M ⊙ subscript 𝑀 direct-product M_{\odot} ); M tid subscript 𝑀 tid M_{\rm tid} : tidal mass ( 10 9 superscript 10 9 10^{9} M ⊙ subscript 𝑀 direct-product M_{\odot} ); r tr subscript 𝑟 tr r_{\rm tr} : truncation radius of dark matter in kpc; N p subscript 𝑁 p N_{\rm p} : number of particles ( 10 5 superscript 10 5 10^{5} and v 0 subscript 𝑣 0 v_{0} : approximate initial orbital speed in km s -1 . See the notes below the table for more information.

Model	$M_{\mathrm{tot}}$	$M_{\mathrm{tid}}$	$r_{\mathrm{tr}}$	$N_p$	$v_0$
P (unperturbed/isolated galaxy)	0
S (intermed. mass, one transit)	50	30	19	1	360
R (high mass, one transit)	100	60	24	2	372

Equations

$$\varphi =\frac{1}{b}(r-a)+i\frac{\pi }{2},$$

Full text

††pubyear: 2019††pagerange: The GALAH survey and Gaia DR2: Linking ridges, arches and vertical waves in the kinematics of the Milky Way–A

1 Introduction

The second data release of the Gaia astrometric mission (DR2: Gaia Collaboration et al., 2018a) has heralded a new era in the field of Galactic dynamics. The rich dataset provides position, parallax and proper motion ($\alpha,\delta,\omega,\mu_{\alpha},\mu_{\delta}$) for over a billion stars at unprecedented precision (e.g., $\sigma_{\mu}\approx 10$ micro-arcseconds yr*-1* for $G<14$). In addition, a subset of the full dataset includes line-of-sight velocities from the Radial Velocity Spectrometer (Gaia DR2 RVS, Soubiran et al., 2018) for about 7 million stars, thus providing full 6D phase-space information for this sample. The sheer number of objects covered by Gaia DR2, combined with its high precision, allows us to map the Galaxy’s kinematics in a volume more than an order of magnitude larger than that covered by Hipparcos (Perryman et al., 1997).

The revelation of multiple new stellar streams (Malhan et al., 2018; Price-Whelan & Bonaca, 2018), evidence of non-axisymmetry through substructure in velocity (e.g., Gaia Collaboration et al., 2018b; Trick et al., 2019), inter alia, have all helped build a good consensus that the Galactic disc is far from being in dynamic equilibrium. Particularly remarkable was the discovery made by Antoja et al. (2018, A18 hereafter), which revealed a spiral pattern in the $(z,V_{\rm z})$ plane density around the Solar neighbourhood, the so-called phase-spiral111This feature has been variously referred to as the phase-plane spiral (Binney & Schönrich, 2018) and the phase-space spiral (e.g. Khoperskov et al., 2018). Consistent with the traditional use of “phase mixing” rather than “phase-space mixing,” we adopt the more compact language of phase-spiral (e.g. Bland-Hawthorn et al., 2019). The distinctive phase pattern has also been described as the “snail” or “snail shell” (A18). (their Fig. 1). The phase-spiral, seen most strongly when color coded by $V_{\mathrm{\phi}}$, is thought to be a signature of the Galaxy relaxing from a disturbed state, through phase-mixing. Using toy models, both A18 and Binney & Schönrich (2018) suggested that the phase-spiral was evidence of the Galaxy’s interaction with the Sagittarius dwarf galaxy (Sgr), and further constrain the last impact to about 0.5 Gyr ago. Recent N-body simulations (e.g., by Bland-Hawthorn et al., 2019; Laporte et al., 2019) have shown that tidal interaction with Sgr can indeed reproduce the phase-spiral seen in Gaia DR2, and suggest a similar or younger timescale for the interaction. However, Khoperskov et al. (2018) have shown that the phase-spirals can also be generated through an entirely internal mechanism. In their simulations, they show that the buckling of the Galactic bar can generate bending waves in the disc. This is able to create the phase-spiral, and the wave takes about 0.5 Gyr to travel to the outer disc ($\sim$10 kpc). The spirals survive well after the end of the buckling phase, where these bending waves are supported by the disc self-gravity. These results show that it is non-trivial to distinguish between an internal or an external perturbation.

A18 also revealed that the $(R,V_{\phi})$ space has substructure in the form of diagonal ridges (Figure 1a) and that the $(V_{R},V_{\phi})$ space has arches (Figure 1b), some of which are asymmetric about the $V_{R}=0$ line. They suggest that arches are just the projection of ridges in the $(V_{R},V_{\phi})$ velocity space; however, ridges can also be present without any arches. Ramos et al. (2018) identified some of the $(V_{R},V_{\phi})$ arches and traced their median $V_{\phi}$ at different Galactocentric radii $R$, suggesting that the arches and ridges are linked.

We show that physical understanding of the connection between the ridges and arches is still missing. Moreover, any connection between these dynamical excitations and the phase-spiral has yet to be clearly established. All of these phenomena have distinct properties (amplitudes, wavelengths, etc.) and their unification is the topic of a later paper.

Several models using various simulation techniques have been proposed to explain the ridges and arches. Most models explain either the ridges or the arches, but not necessarily both. Resonant scattering by non-axisymmetric features rotating with a fixed pattern speed, e.g., the bar or spiral arms, has been shown to generate arches. Dehnen (2000) showed that one prominent asymmetric arch and two other weak arches can be developed by a bar, which has since been demonstrated by several other simulations (e.g., Antoja et al., 2014; Monari et al., 2017; Hunt & Bovy, 2018; Hattori et al., 2019; Pérez-Villegas et al., 2017). A18 further showed that resonance with a bar can also generate ridges, but only one or possibly two ridges can be seen in the solar neighborhood as compared to the many seen in Gaia DR2. More recently, simulations by Fragkoudi et al. (2019) also showed that the outer Lindblad resonance of the Galactic bar could give rise to one of the prominent ridges in the $(R,V_{\phi})$ plane and a Hercules-like feature in the $(V_{R},V_{\phi})$ plane.

Phase-mixing models have also been used to explain these kinematic features. In such models, test particle simulations are employed. Test particles are set up to mimic a perturbation and are then evolved in a Milky Way-like potential. A18 showed that ridges can be generated using a horizontal phase-mixing model, but did not show if they lead to arches. Moreover, the physical motivation for the model was also not made clear. Minchev et al. (2009) showed that phase wrapping after interaction with a dwarf galaxy can produce multiple arches in the $(V_{R},V_{\phi})$ plane, similar to those seen in the solar neighbourhood (see also Gómez et al., 2012), but they do not explore the occurrence of ridges.

This raises an interesting question: is the impact with a dwarf galaxy necessary to see multiple ridges? Quillen et al. (2018) point out that the arches seen in Gaia DR2 are tilted about the $V_{R}=0$ line, but those generated by the phase-wrapping model of Minchev et al. (2009) are symmetrical. They propose that a model in which the stars that have recently crossed spiral arms at their apocenter or pericenter can explain the asymmetric arches; however, they do not study the ridges.

Hunt et al. (2018) consider a potential with 2D transient spiral arms that wind up over time, and using the backward integration technique of Dehnen (2000), show that this perturbation can give rise to features such as the Hercules stream in the $(V_{R},V_{\phi})$ plane, as well as multiple ridges in the $(R,V_{\phi})$ plane, and multiple asymmetric arches. Transient spiral arms have been shown to develop in self-gravitating disc simulations (Sellwood, 2011). This sets up the motivation to look for ridges and arches in simulations of this kind. Quillen et al. (2011) showed that asymmetric arches can be generated in self-gravitating N-body simulations, but did not study the ridges in $(R,V_{\phi})$. Laporte et al. (2019) studied N-body simulations involving interaction with a dwarf galaxy, and were able to generate ridges in $(R,V_{\phi})$, but only one arch or moving group could be seen in the $(V_{R},V_{\phi})$ plane. An interesting question to ask is whether the source of the ridges and arches is internal or external, and how we could distinguish between the two. Can a phase-mixing model motivated by transient spiral arms explain multiple ridges and multiple asymmetrical arches?

Vertical waves have also been reported in the Gaia DR2 data (e.g., Gaia Collaboration et al., 2018b; Bennett & Bovy, 2019). Already with the limited coverage of Gaia-TGAS, Schönrich & Dehnen (2018) and Huang et al. (2018) found that the vertical velocity in the solar neighborhood varies with angular momentum $J_{\phi}$. They found a large-scale trend of $V_{z}$ increasing monotonically with $J_{\phi}$, which is a signature of the Galactic warp. Superimposed on this large-scale trend, they also found undulations (or corrugations) indicative of a wave-like pattern. Undulations in the profile of $V_{z}$ as a function of Galactocentric radius $R$ were also reported by Kawata et al. (2018). Gómez et al. (2013) and D’Onghia et al. (2016) both show that undulations in the $V_{z}(R)$ profile can be seen in N-body simulations involving interaction with Sgr. However, the variation of $V_{z}$ as a function of angular momentum was not studied. Are these vertical waves linked to ridges and arches? Can these vertical waves be seen in simulations with or without the interaction of Sgr? This is a question we attempt to address.

In this paper, we revisit the $(R,V_{\phi})$ ridges seen in Gaia DR2. First, we dissect and characterize the ridges using radial velocity, vertical height, and vertical velocity. Furthermore, we explore the nature of ridge stars by considering elemental abundances from GALAH and relate this to the nature of the perturbation itself. Next, we simulate phase mixing of spiral arms and show how this model can be used to understand the connection between the $(R,V_{\phi})$ ridges and the $(V_{R},V_{\phi})$ arches. Finally, we carry out N-body simulations of the Galactic disc, both with a Sgr-like perturber and without any perturber and study the phase-space features in these simulations and compare them with those seen in Gaia DR2.

2 Data set and methods

Throughout the paper, we adopt a right-handed coordinate frame in which the Sun is at a distance of $R_{\odot}=8.2$ kpc from the Galactic center (Bland-Hawthorn & Gerhard, 2016), consistent with the new ESO Gravity measurement (Gravity Collaboration et al., 2018), and has Galactocentric coordinates $(X,Y,Z)=(-8.2,0,0.25)$ kpc. The cylindrical coordinate angle $\phi={\rm tan}^{-1}(Y/X)$ increases in the anti-clockwise direction, while the rotation of the Galaxy is clockwise. The heliocentric Cartesian frame is related to Galactocentric by $X_{\rm hc}=X+R_{\odot}$, $Y_{\rm hc}=Y$ and $Z_{\rm hc}=Z$. $X_{\rm hc}$ is negative toward $\ell=180^{\circ}$ and $Y_{\rm hc}$ is positive towards Galactic rotation. For transforming velocities between heliocentric and Galactocentric frames we use $(\dot{X}_{\odot},\dot{Y}_{\odot},\dot{Z}_{\odot})=(U_{\odot},\Omega_{\odot}R_{\odot},W_{\odot})$. Following Schönrich et al. (2010), we adopt $(U,V,W)_{\odot}=(11.1,12.24,7.25)$ km s*-1*, while for the azimuthal component we use the constraint of $\Omega_{\odot}=30.24$ km s*-1kpc-1* which is set by the proper motion of Sgr A*, i.e., the Sun’s angular velocity around the Galactic center (Reid & Brunthaler, 2004). This sets the rotation velocity at the Sun to $V_{\phi,\odot}=-248$ km s*-1*, and thus the circular velocity at the Sun to $V_{c,\odot}=-236$ km s*-1*. We now describe the astrometric and spectroscopic data that we use in this work and the quality cuts that we apply on them.

2.1 Gaia DR2 RVS sample

In this paper we make use of the Gaia DR2 radial velocity sample (Gaia DR2 RVS) which provides full 6D phase space information ($\alpha,\delta,\omega,\mu_{\alpha},\mu_{\delta},V_{\rm los}$). We selected stars with positive parallax and with parallax precision $\sigma_{\omega}/\omega<$ 0.2, which gave a sample of $6376803$ stars. The SQL query used to generate the sample is given in Appendix A. We estimated distance as $1/\omega$, which is reasonably accurate for our selected stars and for the purpose of this paper (Luri et al., 2018).

2.2 GALAH DR2 sample with elemental abundances

The spectroscopic data used here is taken from an internal release of GALAH DR2, which includes, public data (Buder et al., 2018, DR2), and fields observed as part of the K2-HERMES (Wittenmyer et al., 2018) and TESS-HERMES (Sharma et al., 2018) programs. To maintain the survey selection function, we have applied the quality cuts summarised in Table 1, which gives a total of 465870 stars cross-matched with Gaia DR2. This internal release includes non-LTE corrections on [Fe/H] but not on [$\alpha$/Fe]. For the kinematics of this dataset, we make use of the parallax and proper motion ($\omega,\mu_{\alpha},\mu_{\delta}$) from Gaia DR2, but use the highly precise radial velocities from GALAH, which have typical error of 0.1 km s*-1* (Zwitter et al., 2018). Since we are mainly interested in nearby stars, we restrict our GALAH sample only to dwarfs, by applying a surface gravity cut of (\log g$$>3), which results in a final sample of 258289 stars. This avoids any issues related to systematic errors in stellar parameters between dwarfs and giants.

2.3 Phase mixing simulations

To understand the origin of the phase-space substructures like ridges and arches, we perform simulations in which spiral arms phase mix and disrupt over time. The simulations are motivated by the desire to mimic the effect of transient spiral arms. For this we consider an initial distribution of particles confined to four thin spiral arms. The $i^{th}$ arm is setup as an Archimedean spiral, with azimuth:

[TABLE]

where, $0<a<2\pi$ controls the orientation of the spiral, $b=\pi/10$ controls the tightness of the winding, and $r=R/8$ kpc. The radial distribution was assumed to be skew normal with skewness of 10, location parameter of 4 and scale parameter of 6. This is to ensure that there are enough particles in the Solar neighbourhood-like volume. The radial velocity was sampled from $\mathcal{N}(0,20)$ and the azimuthal velocity from $\mathcal{N}(\Theta(R),20)$, where $\Theta(R)$ denotes the circular velocity. For simplicity, the particles were set up in the midplane with zero vertical velocity. A total of 640000 particles were evolved for 650 Myr with galpy (Bovy, 2015) using the MWPotential2014 potential, consisting of an axisymmetric disc, a spherical bulge and a spherical halo. The set up is similar to Antoja et al. (2018), but they start with stars arranged in a single line as compared to four spiral arms used by us. A set of movies showing the evolution of the system is available as Supporting Information in the online version of this paper.

2.4 N-body simulations

Another approach that we adopt, in order to gain insight into the origin of phase-space substructures, is the use of N-body simulations of a multi-component Galaxy. Ideally, we would like to carry on a detailed modelling of every component of the Galaxy, both collisionless (e.g. dark matter and stars) and gas. This is, however, not feasible as it is both computationally expensive and non-trivial to do. We adopt a common, simplifying assumption: we assume that a pure N-body (rather than a full N-body, hydrodynamical) model is sufficient for our purposes. We caution that neglecting the gas components in these type of experiments may not always be appropriate (see e.g. Tepper-García & Bland-Hawthorn, 2018).

We consider here the following two scenarios as the plausible origin of the ridges: i) instabilities internal to the Galaxy; and ii) tidal (external) interactions. In consequence, we focus our attention on one representative model for each of these. On the one hand, we simulate the evolution of an isolated Galaxy starting from some prescribed initial conditions (see below). On the other hand, we simulate how the stars behave in a Galaxy that has been tidally perturbed by the interaction with a smaller system. It has been suggested that Sgr may lie behind many of the kinematic features revealed by Gaia DR2 (e.g. Antoja et al., 2018; Binney & Schönrich, 2018; Laporte et al., 2019; Bland-Hawthorn et al., 2019). It therefore seems natural in our case to simulate the interaction of the Galaxy with a Sgr-like perturber.

Our isolated model Galaxy consists of four collisionless components: a host DM halo; a stellar bulge; a thick stellar disc; and a thin stellar disc. We refer to this model as the ‘isolated’ model (Model P); see Table 2 for details of this model Galaxy. Note that the values for structural parameters (scalelength, scaleheight, etc.) are only consistent with the range of values inferred from observations (Bland-Hawthorn & Gerhard, 2016). We assume the distribution of stars in the thick disc is well approximated by a Miyamoto & Nagai (1975b) (MN) profile. This choice, and its corresponding parameter values, are founded on the work of Kafle et al. (2014), who use precise stellar kinematic information to infer the mass distribution of the Milky Way, assuming that total stellar disc component is well described by a single MN component. They obtain a scalelength and scaleheight of 5 kpc and 0.5 kpc, respectively. Note that Bland-Hawthorn & Gerhard (2016) provide instead mean value of 2 kpc and 0.9 kpc, respectively, which is however predicated by the assumption that the mass distribution in the (thick) stellar disc be well described by an exponential profile. Our adopted values are consistent with the fact that the scalelength of a MN may differ by a up to factor of $\approx$2 compared to the scalelength of an equivalent exponential disc (Flynn et al., 1996) and that the exponential scaleheight can differ, potentially by the same factor of 2, from the MN scaleheight and that the exponential scaleheight can differ, potentially by the same factor of 2,from the MN scaleheight (Smith et al., 2015). For the thin disc we adopt a scalelength and a scaleheight of 3 kpc and 0.3 kpc, respectively. These values agree well with the mean values quoted by Bland-Hawthorn & Gerhard (2016). Our adopted mass for the thin disc is well within the range inferred from observations. The adopted mass for the thick disc is higher than quoted by Bland-Hawthorn & Gerhard (2016). However, it should be noted that our adopted value corresponds to the mass integrated out to 20 kpc. A lower mass is obtained if truncating the disc at a smaller radius. Nevertheless, the total stellar disc mass is consistent with other estimates (e.g., Bovy & Rix, 2013). Overall, our choices of component profiles and the values of their corresponding parameters define a valid model for the Galaxy, comparable to the model successfully adopted by others in numerical studies (e.g., Chequers et al., 2018).

Our interaction models consist essentially of a two-component (DM, stellar spheroid) system orbiting an initially isolated Galaxy along an (unrealistic) hyperbolic orbit. The reason for choosing such an orbit rather than a more realistic orbit for the perturber is that, as we have shown previously (Bland-Hawthorn et al., 2019), each passage of Sgr across the Galactic plane washes out the kinematic signatures of its previous crossing, thus limiting the time span available between crossings. In contrast, by adopting a hyperbolic orbit we ensure that Sgr transits the Galactic plane (disc) once only, thus facilitating the analysis of its effect on the Galactic stars.

We consider perturbers with total masses of 5 or $10\times 10^{10}$ $M_{\odot}$, spanning the mid-to-high range of plausible Sgr masses at infall (e.g Niederste-Ostholt et al., 2010). Both the stellar system and the dark halo are modelled as truncated Hernquist (1990) spheres. Their scale radii are initially set at $0.85$ and 10 kpc, respectively. The stellar system is initially truncated at 2.5 kpc while the truncation radius of the dark halo is listed in Table 3. A simulation with each of these masses was started with the perturber at $(x,y,z)=20.8,0.,45.5)$ kpc on an orbit of eccentricity $e=1.3$ (hyperbolic) and pericentric distance 10 kpc.222The exact initial initial velocities for model R and model S was $(v_{x},v_{y},v_{z})=(-267,0,-260)$ km s*-1* and $(v_{x},v_{y},v_{z})=(-258,0,-251)$ km s*-1*, respectively.

Two key requirements on these type of simulations, imposed by the exquisite detail on the kinematics of stars revealed by the data, are the mass resolution (or particle number) and the limiting spatial resolution. The latter has to be low enough to allow for a correct simulation of the evolution of the dynamically coldest stellar component within $|z|<0.2$ kpc. The former needs to be high enough to allow for a dense enough sampling of the $(R,V_{\phi})$.

We choose values for the particle number and spatial resolution such that we fulfill these requirements while keeping the computational cost of the simulations reasonably low. More specifically, we set the limiting spatial resolution at 30 pc, to be compared with 300 pc, the (initial) scaleheight of the cold stellar disc, which is the smallest length scale in our simulation. The adopted particle number varies from component to component, depending on the total mass of the component and the corresponding particle mass; as per our above discussion, the thick and thin stellar discs have been assigned the absolute highest particle number (see Table 2).

The simulations’ axisymmetric initial conditions, i.e. the particles’ positions and velocities for each component were assigned by the technique of Springel et al. (2005) as implemented in the dice code (Perret et al., 2014). In doing so, all the components are intended to be in dynamical equilibrium with the total potential of the compound system. However, in reality they will in general be slightly out of equilibrium (e.g. Kazantzidis et al., 2004), and even an isolated Galaxy disc will develop some small-scale structure, such as rings and transient spirals. While usually an unwanted numerical artifact in experiments like ours, here we actually need these instabilities in order to investigate their effect on the stellar phase-space kinematics. In real galaxies, such transient features will develop as well, but for reasons yet to be fully understood.

The evolution of the system333The dice and Ramses configuration files used to create our initial conditions and to setup our simulations, respectively, are freely available upon request. in each case is calculated with the adaptive mesh refinement (AMR) gravito-hydrodynamics code Ramses (version 3.0 of the code described by Teyssier, 2002). Simulation data are stored at approximately $\Delta\tau=$10 Myr intervals. A set of movies showing the evolution of the system in each model are provided at Gaia-GALAH-phase-spiral.

3 Results

We begin by studying the $(R,V_{\phi})$ plane using the observed data. Next we compare the observed results with predictions from two type of simulations, phase mixing simulation of disrupting spirals and disc N-body simulations. This is followed by a study of the $(V_{R},V_{\phi})$ plane and subsequently the $(R,V_{R})$ plane. In both cases we compare the observed results with the predictions from the simulations.

3.1 Analysis of the $(R,V_{\phi})$ plane using the observed data

3.1.1 Dissection in kinematics and vertical height with Gaia

Antoja et al. (2018) revealed the diagonal ridge-like structures in the $(R,V_{\phi})$ density distribution. In this section, we further explore this plane using kinematics and vertical height. In Figure 2, we show results using the Gaia DR2 RVS sample. We select stars with ($|\phi-180^{\circ}|<25^{\circ}$) & ($|R-R_{\odot}|/{\rm kpc}<3.5$).

Figure 2a shows the density distribution $\langle V_{R}\rangle$ in the $(R,V_{\phi})$ plane. Multiple diagonal ridges are clearly visible, extending between $6<R/{\rm kpc}<12$. Ridges are prominent at $R=8.2$ kpc, but seem to fade away as we move away from the solar radius, this is because, as we move away from the Sun there is a fall in number density of stars and an increase in uncertainty in $R$ and $V_{\phi}$. Figure 2b shows a map of $\langle V_{R}\rangle$ in the $(R,V_{\phi})$ plane. The ridges are more prominent and are visible even at large distances from the solar neighborhood. They appear at a similar location to that in the density map in Figure 2a). Stars along the ridges are moving either radially outward or inward, with respect to the background distribution, with \langle V_{R}\rangle$$\approx 10 km s*-1*.

Next we explore the properties of the ridges in the vertical direction. Figure 2c shows a map of $\langle V_{z}\rangle$ in the $(R,V_{\phi})$ plane. The ridge structure can again be seen but it is weaker as compared to the $\langle V_{R}\rangle$ map. Three ridges are clearly visible and the $\langle V_{z}\rangle$ associated with the structures is about 2 km s*-1*. Figure 2d shows $(R,V_{\phi})$ mapped by $\langle|z|\rangle$, i.e., the mean of the absolute distance from the mid-plane of the disc. The ridges are primarily composed of stars that lie close to the Galactic plane ($|z|\lesssim 0.2$ kpc), as indicated by the distinctive dark color. It is important to note that if stars at all heights above the plane participated in the ridges, the $|z|$ map in Figure 2d would be completely featureless. This preferential distribution must thus be linked to the nature of the perturber responsible for the ridges. Three ridges are also visible in the map of $\langle z\rangle$, i.e. the mean distance from the plane (Figure 2e).

We have overplotted curves of constant angular momentum (black dotted lines) at $L_{\rm z}=[1350,1600,1800,2080]$ kpc km s*-1* and curves of constant orbital energy $E$ (white dashed lines) at (E-E_{\rm circ}(R_{\odot}))/V_{\rm circ}^{2}(R_{\odot})$$=[-0.112,-0.021,0.097]. The energy was evaluated using the MWPotential2014 potential in galpy (Bovy, 2015). Both curves decrease with $R$ and resemble ridges, and hence either of these physical quantities can be used to label the ridges. The difference between the two is that the constant energy curves are straight lines but the angular momenta ones are not. From these plots, it is difficult to say if the ridges are constant energy or constant angular momentum, we revisit this issue later in subsection 3.2.

We have shown that the ridges are present in maps of density, kinematics and vertical height. We now investigate if the ridges in the maps of different quantities are correlated with each other. For this we select stars in a narrow range in $R$ and then study the one dimensional profiles of various quantities as a function of orbital energy $E$ (Figure 3). We choose $E$ instead of $V_{\phi}$ as ridges are well approximated by curves of constant energy. We use the MWPotential2014 potential in galpy to compute the energy (Bovy, 2015). Instead of directly using $E$, we use the dimensionless form given by E^{\prime}=$$(E-E_{\rm circ}(R_{\odot}))/V_{\rm circ}^{2}(R_{\odot}), where $V_{\rm circ}$ is the circular velocity at a given radius, and $E_{\rm circ}$ is the energy of a star in a circular orbit at $(R,z)=(R_{\odot},0)$. Figure 3(a,b) show the density profiles. At least 8 peaks can be identified and these are marked with vertical dotted lines. The peak at $E^{\prime}=-0.156$ corresponds to the Hercules stream and is shown with a different color. Figure 3(c) shows the profile of mean vertical velocity. A large-scale trend of increase in $V_{z}$ with $E$ can be seen similar to Schönrich & Dehnen (2018) who studied $V_{z}$ as a function of $L$. Note, for a given $R$ and $V_{R}$, $E$ increases monotonically with $L$, and here the range of $R$ is almost constant and $V_{R}$ is small. This large-scale trend of $V_{z}$ is due to the warp.

Besides the large-scale trend, peaks at $E^{\prime}=[-0.021,0.097]$ can also be seen. The location of these peaks matches with peaks seen in density. Figure 3(d) plots median vertical distance z. There is no large-scale trend, but 3 peaks ($E^{\prime}=[-0.21,0.097,0.189]$) are clearly identifiable and they match with the peaks in density. Two of the peaks also match with peaks in $V_{z}$. Figure 3(e) shows the median value of $-|z|$. Almost all density peaks have a corresponding peak in this plot, which is a reflection of the fact that the stars in the density peaks lie close to the Galactic plane. Finally, Figure 3(f) shows the profile of median $V_{R}$. Although all peaks do not match in location with all peaks in density, however, for each undulation in the profile of density there is an undulation in $V_{R}$. This indicates that $V_{R}$ and the density peaks are strongly correlated with each other. Note, a consequence of $V_{R}$ peaks not matching up with density peaks is that stars in a ridge are not symmetrically distributed about $V_{R}$, and arches in the $(V_{R},V_{\phi})$ plane show such a behaviour.

3.1.2 Dissection in elemental abundances with GALAH

We now study the elemental abundance in the $(R,V_{\phi})$ plane. In Figure 4a, $(R,V_{\phi})$ is mapped by [Fe/H]. For the region 200<$$V_{\mathrm{\phi}}/km s*-1*$<$250 the background metallicity is around [Fe/H]$\approx-0.1$, reflecting the local ISM around the solar neighbourhood which is sub-solar (Nieva & Przybilla, 2012). The ridges in this region however, are mainly composed of solar metallicity stars, with typical [Fe/H]$\approx$0.03. In Figure 4b $(R,V_{\phi})$ is mapped by [$\alpha$/Fe]. The ridges stand out as a population with [$\alpha$/Fe]$\approx 0.05$ (close to solar values). This is consistent with the ridges being made of stars that lie predominantly in the plane. Stars close to the plane are younger, and young stars are metal rich and alpha-poor (age-scaleheight and age-metallicity relations, e.g., Mackereth et al., 2017).

Beyond L_{\rm Z}$$=2080 kpc km s*-1* there is a sharp cut-off in metallicity (black dotted curve, Figure 4a). This region is dominated by relatively metal-poor stars with typical [Fe/H] $\approx-0.3$, and is also alpha-enhanced around [$\alpha$/Fe]$\approx 0.1$ (Figure 4b). This suggests that the origin of these stars is different from those along the ridges. L_{\rm Z}$$>2080 kpc km s*-1* corresponds to a guiding radius $R_{\rm G}>9.5$ kpc (assuming a flat rotation curve). These stars thus belong to the outer disc and their low metallicity is consistent with the Galaxy’s negative metallicity gradient with $R$ (Hayden et al., 2014).

Similarly, stars at the bottom of Figure 4(a,b) with V_{\mathrm{\phi}}$$< 150 km s*-1* also show a sharp change in abundances. These stars have large asymmetric drift, are rotating slowly, and have [$\alpha$/Fe] $>$ 0.14 and [Fe/H] $<-0.4$. These properties are consistent with that of the traditional thick disc, which is metal-poor, alpha-enhanced, and kinematically hot (Bensby et al., 2014; Duong et al., 2018).

3.2 Analysis of the $(R,V_{\phi})$ plane using a phase-mixing simulation

We now consider a toy model of phase mixing similar to that used by Antoja et al. (2018) to explain some of the features seen in the $(R,V_{\phi})$ plane. We consider an initial distribution of particles confined to four thin spiral arms; in Antoja et al. (2018) the particles were confined to a single line. The particles are then evolved in time under a multi-component analytic potential. The simulation is designed to mimic phase mixing of perturbations caused by transient spiral arms (for further details see subsection 2.3).

The distribution of stars in the $(X,Y)$ and the $(R,V_{\phi})$ planes are shown in Figure 5 for four different snapshots in time. Also shown are maps of $\langle V_{R}\rangle$ in the $(R,V_{\phi})$ plane; see Figure 5(i-l). As we move forward in time, the spiral pattern decays (Figure 5(a-d)), and the ridges start to form and they increase in number and become more stretched and therefore thinner (Figure 5(e-h)). The ridges can also be seen in maps of $\langle V_{R}\rangle$. The ridges are approximately linear in the $(R,V_{\phi})$ plane and resemble lines of constant angular momentum. The appearance of the ridge structure is a consequence of phase mixing and can be understood in terms of Liouville’s theorem, which states that the full phase-space density (or volume) of a system evolving in a fixed potential is conserved. In the case of our simulation, the phase space is made of $(X,Y,V_{X},V_{Y})$. Initially the density in the $(X,Y)$ space is high while that in $(V_{R},V_{\phi})$ space is low. As the spiral pattern disperses, the density in the $(X,Y)$ plane reduces, but to conserve the phase-space density, the density should increase in other dimensions. The structures in the $(R,V_{\phi})$ are in some sense a reflection of this phenomenon.

Figure 5(m-p) show the distribution of stars in orbital energy and Galactocentric radius $(E,R)$ plane. Discrete energy levels can be seen. Stars in a ridge lie in a narrow energy interval; we explore this further in Figure 6. The figure shows the distribution of stars in the $(E,R)$ and $(J_{\phi},R)$ plane. The top panels show results for Gaia DR2, while the bottom panels are for the phase-mixing simulation where we only show stars belonging to a single spiral. It is clear from Figure 6(c,d) that in phase-mixing simulations, the ridges are curves of constant energy rather than constant angular momentum. According to Figure 6(c,d), near $E^{\prime}=0$, constant energy and constant angular momentum curves are both expected to be flat, it is only at higher values of $|E^{\prime}|$ that one can differentiate between the two cases. For the observed data, ridges between $-0.15<E^{\prime}<0$ (corresponding to the Hercules stream) look flat in $J_{\phi}$, while the rest of the ridges, especially with large values of $|E|^{\prime}$, are flatter in $E^{\prime}$ than in $J_{\phi}$. This suggests that different ridges can originate from different physical processes. The lower most ridge (vertical coordinate of -0.3) in Figure 6b is slanted downwards in $J_{\phi}$ just like in Figure 6d, while the topmost ridges in Figure 6b is slanted outwards just like in Figure 6d. For the lower two ridges with $E^{\prime}<-0.2$, it is clear that the ridges are sharper in energy than in angular momentum, lending further support to the ideas that energy as a quantity is better than angular momentum for characterizing these ridges.

3.3 Analysis of the $(R,V_{\phi})$ plane using disc N-body simulations

We now consider the more realistic N-body simulations of the Galaxy described in subsection 2.4. The first scenario, Model P, is that of an isolated galaxy, i.e., unperturbed by a satellite. The $(X,Y)$ density of four selected snapshots at $\tau=[0.21,0.56,0.97,2.59]$ Gyr are shown in Figure 7(a-d). At $\tau=0.21$ Gyr (Figure 7a), the disc settles into an equilibrium configuration and develops tightly wound spiral arms. Such self-excited instabilities forming spiral arms are a known feature of N-body simulations in disc galaxies (Sellwood, 2012). The corresponding $(R,V_{\phi})$ density map (Figure 7e) is largely uniform and lacks ridge-like substructure as seen in Gaia DR2 (in Figure 2). Similarly, in the velocity maps (Figure 7(i,m)), there are fine-structure blobs in the kinematics with $\langle V_{R}\rangle$ $\approx 10$ km s*-1* and $\langle V_{z}\rangle$ $\approx 3$ km s*-1*, but no ridges can be seen.

By $\tau=0.56$ Gyr, the spiral arms have weakened slightly, they are fewer and thicker (Figure 7b). Interestingly, the $(R,V_{\phi})$ density at this snapshot shows large-scale diagonal stratification, with a span of about 4 kpc (Figure 7f). The $\langle V_{R}\rangle$ map shows multiple thin diagonal ridges with an alternating pattern of radially outward and inward motion (Figure 7(j)).

By the next snapshot at $\tau=0.97$ Gyr, the spiral arms are found to have diffused and weakened (Figure 7c). The corresponding $(R,V_{\phi})$ density map shows several prominent ridges that extend over $5<R/kpc<15$ and have a more linear appearance compared to the previous snapshot (Figure 7g). The ridges are also clearly present in the $\langle V_{R}\rangle$ and $\langle V_{z}\rangle$ maps, where the amplitude of the radial oscillations is again higher than the vertical component.

By the final snapshot, chosen at $\tau=2.59$ Gyr, the spiral arms are found to have significantly decayed. A central bar with half-length of $\sim 2.5$ kpc is visible prominently (Figure 7d). The density, $\langle V_{R}\rangle$, and $\langle V_{z}\rangle$ maps continue to show large scale ridges (Figure 7(h,i,p)). In summary, Figure 7 shows that an unperturbed galaxy can reproduce ridges in the $(R,V_{\phi})$ plane with features similar to that seen in Gaia DR2. The ridges appear as the spiral structure decays, and are maintained as long as this decay is going on. As was already mentioned in subsection 3.2, this is a consequence of Liouville’s theorem which requires that the density in phase-space is always conserved. This suggests that internal instabilities such as transient spiral arms, could be responsible for the ridges seen in Figure 2.

Next, we consider the scenario where the Galaxy is tidally perturbed by an orbiting satellite. Model S simulates the interaction with an intermediate mass Sgr galaxy ($5\times 10^{10}M_{\odot}$), while Model R simulates the interaction with a heavier Sgr galaxy ($10^{11}M_{\odot}$). In both cases, Sgr crosses the disc at around $\tau=0.15$ Gyr, and perturbs the galaxy from its equilibrium state. Previously, in simulations run in Bland-Hawthorn et al. (2019), we noted that disc crossing by Sgr wipes out previous coherent structure and generates new structures in the Galaxy. Evolving the galaxy for $\tau=1.5$ Gyr, allows for enough time to develop, decay, and phase mix the spiral arms as well as the effects of Sgr. For this reason we compare the unperturbed and perturbed scenarios at roughly coeval timestamps of ($\tau=1.5$ Gyr), i.e., allowing for enough time for perturbations to phase mix.

Figure 8 shows the density, $\langle V_{R}\rangle$, and $\langle V_{z}\rangle$ maps in the $(R,V_{\phi})$ plane for the various N-body simulations alongside Gaia DR2 data. We note the presence of ridges in all three simulations (Figure 8(a,b,c)). Ridges are also present in the maps of $\langle V_{R}\rangle$ (Figure 8(f,g)) and $\langle V_{z}\rangle$ (Figure 8(j,k)). The ridges for Gaia DR2 data in (Figure 8d appear to be smeared out at the edges. This is due to observational errors in proper motion and parallaxes that are dominant at larger distances. We would like to point out, that the purpose of the N-body simulations in this paper, is to demonstrate what kinematic signatures can be generated at a given location in the Galaxy. These are not, however, selection function matched snapshots i.e., the particles do not have stellar parameters/magnitudes assigned that we could convolve Gaia DR2-like errors. In any case, our aim her is not to match the exact number of features or their location one-to-one, which would be affected by the smear due to observational errors.

In Figure 3, we saw that for Gaia DR2, ridges are correlated in kinematics and spatial density. In Figure 9 we explore similar correlations for our N-body simulations. We select stars around $(R,z,\phi)=(8.2,0.0,180.0^{\circ})$ and consider the profiles of $\langle V_{z}\rangle$, $\langle V_{R}\rangle$, and $z$ against the dimensionless orbital energy, E^{\prime}=$$(E-E_{\rm circ}(R_{\odot}))/V_{\rm circ}^{2}(R_{\odot}). For all simulations, peaks can be seen in profiles of density, $\langle V_{z}\rangle$, $z$, and $V_{R}$. It is worth noting that the unperturbed model has no tidal interactions, hence, the observed vertical oscillations for the unperturbed model must be due to internal processes.

A number of features seen in Figure 3 for the Gaia DR2 data can also be seen in the simulations. The location of peaks in $z$ match with location of peaks in $v_{z}$. Location of extrema in $v_{z}$ match with location of peaks in density. For the the unperturbed case it is the minima that matches and for the high mass case it is the maxima. For the intermediate mass case we do not see such an association. We note that the matching of the location of peaks in $z$ and $v_{z}$ is not a general feature, because it was only seen at a few special locations within the galaxy.

The $\langle V_{z}\rangle$ profiles (Figure 9(j-l)) show a large scale trend like in the Gaia DR2 data. Such a trend is expected for the presence of a warp. A clear warp was detected in all our simulations. A plot of mean $z$ and $v_{z}$ as a function of $\phi$ showed a sinusoidal pattern with the $v_{z}$ profile being shifted by $90^{\circ}$ with respect to the $z$ profile.

The amplitude of fluctuations for all the plotted quantities (density, $\langle V_{z}\rangle$, $z$, and $V_{R}$), is considerably higher for the high mass Sgr case compared to the other two simulations. A comparison with Figure 3 shows that the amplitude of $\langle V_{z}\rangle$, $z$, and $\langle V_{R}\rangle$ fluctuations for the case of Gaia DR2 is comparable to the case of unperturbed and intermediate mass Sgr simulations, making the case for the high-mass perturber unfavourable.

3.4 Analysis of the $(V_{R},V_{\phi})$ plane: arches

We now study the $(V_{R},V_{\phi})$ plane. Figure 1b shows the distribution of Gaia DR2 stars. Arch-like structures can be seen and they are asymmetrical about the $V_{R}=0$. In Figure 10, we show the distribution of stars in the phase-mixing simulation. Initially, there are no arches, but as time proceeds, arches start to appear, increase in number, and become thinner. In the final snapshot, at 494 Myr, multiple arches are clearly visible, and they also appear to be asymmetrical as in the observed data.

We now study the distribution of stars in the $(V_{R},V_{\phi})$ space using disc N-body simulations. Figure 11 shows the distribution of stars for simulation P. We show snapshots corresponding to time $\tau$ of 250, 500, 1000 and 1500 Myr. For each time, we show distributions at four different locations in azimuth. The simulation starts with a smooth disc and by 250 Myr strong tightly wound spiral arms can be seen, however the velocity distribution is devoid of any substructures at this stage. As the simulation evolves, the velocity distribution becomes irregular and develops substructures. Arches are visible in all snapshots with $\tau>=500$ Myr and they are not symmetric about the $V_{R}=0$.

Figure 12 shows the distribution of stars for simulation S that corresponds to interaction with an intermediate mass satellite. Similarly, Figure 13 shows the distribution for simulation R that corresponds to interaction with a high mass satellite. As compared to the simulation P of the unperturbed galaxy, considerably more substructures and arches can be seen in simulations S and R. The simulation R with high mass shows more arches than simulation S. It is clear that arches can develop even when there is no external perturber, but when a perturber like an orbiting satellite is present, the arches are stronger and more numerous. The fact that the number and strength of arches depends upon the mass of the satellite means that we can use the observed data to put limits on the satellite mass. From our set of simulations, we conclude that simulation P has too few arches and simulation R too many, and it is the simulation S that matches best with the Gaia DR2 data.

3.5 Analysis of the $(R,V_{R})$ plane

Structures have been reported in the $(R,V_{\phi})$ and $(z,V_{z})$ planes, but so far the $(R,V_{R})$ space has not been explored. Figure 14 shows maps of density and $\langle V_{\phi}\rangle$ in the $(R,V_{R})$ phase space. The density map is extremely smooth and shows no substructure in Figure 14(a-d). However, arrow shaped substructures can be seen in the $V_{\phi}$ map in Figure 14(e-h). Phase mixing can explain the substructures seen in this space. Such substructures can also be seen in disc N-body simulations. These substructures provide additional independent constraints on models trying to explain the origin of phase-space substructures in the Galaxy.

4 Discussion and Conclusions

We have explored the ridge-like features in the $(R,V_{\phi})$ plane using position and velocities from Gaia DR2 and elemental abundances from GALAH. We find that ridge-like features are visible not only in the density maps but also in maps of $\langle V_{R}\rangle$ , $\langle|z|\rangle$, $\langle V_{z}\rangle$, [Fe/H] and $[\alpha/{\rm Fe}]$ (Figure 2 and Figure 4). Ridges in the $\langle V_{R}\rangle$ map are more prominent and visible to much larger Galactocentric radii $R$ than in the density map. The $\langle|z|\rangle$ map suggests that the ridges are more prominent for stars close to the mid-plane of the Galaxy. The GALAH data suggest that stars in the ridges are predominantly of higher metallicity than the non-ridge stars ($\sim$solar [Fe/H]) and solar $[\alpha/{\rm Fe}]$ (Figure 4). Since, typically stars close to the plane have values of [Fe/H] and $[\alpha/{\rm Fe}]$ that are close to solar, this explains the trends with elemental abundance. That the ridge stars are predominantly at low $|z|$ could be due to one or all of the following three reasons: a) The ridges are due to transient perturbations (i.e., spiral arms) that are close to the plane and are disrupting and phase mixing with time; b) the ridges are due to interaction of stars with perturbations that are close to the plane; c) stars close to the plane are kinematically cold and it is easier to perturb them.

Our phase-mixing simulation of disrupting spiral arms can explain a wide array of kinematic features in the observed Gaia DR2 data. They simultaneously reproduce the ridges in the $(R,V_{\phi})$ plane (Figure 5(e-h)), the ridges in the $\langle V_{R}\rangle$ maps (Figure 5(i-l)), and the arches in the $(V_{R},V_{\phi})$ plane (Figure 10). They also reproduce the observed asymmetry in the arches. While a bar perturbation has been shown to generate ridges, the number of ridges generated from a bar alone are too few to match the observed data (Antoja et al., 2018; Hunt et al., 2018). Phase mixing generates surfaces of constant energy and this explains the occurrence of both the ridges and the arches.

More realistic N-body simulations of a disc in which spiral arms are naturally generated support the results obtained from phase mixing. In these simulations, the spiral arms grow in strength with time till about 500 Myr, and then start to decay. As the spiral arms decay and get phase mixed, the ridges and arches are found to grow in prominence, a phenomenon that was also seen in the phase-mixing simulation (Figure 7). Our N-body simulations show ridges in the $\langle V_{z}\rangle$ maps as seen in the observed data. Simulations in which the disc is perturbed by the passage of an orbiting satellite also show features similar to the case of an unperturbed disc (Figure 8). However, the ridges are found to be more pronounced, in both the $\langle V_{R}\rangle$ and $\langle V_{z}\rangle$ maps, when the mass of the orbiting satellite is higher. This makes the case of a $10^{11}M_{\odot}$ or higher mass perturber unfavourable but a perturber with $5\times 10^{10}M_{\odot}$ is still consistent with Gaia DR2.

Antoja et al. (2018) tentatively suggest that arches in the $(V_{R},V_{\phi})$ plane are projections of ridges in the $(R,V_{\phi})$ plane. We note that, while ridges do suggest existence of discrete values of $V_{\phi}$ in the solar neighborhood, they do not necessarily suggest the presence of arches. It is impossible to deduce the distribution of $V_{R}$ from the distribution of stars in the $(V_{R},V_{\phi})$ plane. We have shown that phase-mixing simulations of disrupting spiral arms not only generate ridges but also arches. The physical property unifying the two features is the energy. A ridge in $(R,V_{\phi})$ and an arch in $(V_{R},V_{\phi})$ are both curves of constant energy. The phase mixing of disrupting spiral arms generates a regular pattern of peaks in the energy distribution of a sample confined to a narrow range in azimuth (Figure 5 o,p). A curve of constant energy and constant angular momentum both appear as a ridge in the $(R,V_{\phi})$ plane. However, out of the above two, only a constant energy curve will manifest itself as an arch in the $(V_{R},V_{\phi})$ plane. Note, we observe stars in a narrow range of azimuth, only stars of certain discrete values of angular frequency will end up in the chosen azimuth range at a given time. The fact that we see discrete energy levels suggest that the angular frequency is more strongly correlated with energy than with angular momentum.

Two different techniques, our work using phase mixing and work by Hunt et al. (2018) using scattering from a perturbation in the potential, both suggest that transient winding spiral arms can explain the multiple ridges and arches seen in Gaia DR2. Interestingly, transience here is through a process of wrapping up rather than fading away in strength with time as generally thought.

The arches seen in Gaia DR2 are asymmetrical about $V_{R}=0$ in the $(R,V_{\phi})$ plane . Phase mixing is generally thought to produce symmetric arches (Quillen et al., 2018), as was observed by Minchev et al. (2009) in their phase-mixing simulations. This is because, for stars on an arch, the orbital energy is approximately fixed, and since $E\sim V_{R}^{2}+V_{\phi}^{2}$, the arches are symmetric. However, we show that phase-mixing simulations can generate asymmetrical arches, and that the asymmetry is both intrinsic and apparent. The slight intrinsic asymmetry is due to the fact that an arch has a finite width in energy and the $V_{R}$ changes systematically with energy. This occurs in the initial stages when phase mixing is incomplete (Figure 10 b,c). The apparent asymmetry is due to the following reason and is responsible for asymmetry seen at later stages of phase mixing (Figure 10 d). The arch due to a single ridge and a single spiral arm is in general symmetrical in the $(V_{R},V_{\phi})$ plane, but the number density of stars is not symmetrical about $V_{R}=0$. Moreover, the arch is short and does not span the full range of $V_{R}$. When multiple arches from different spiral arms are superimposed they look like a large arch with a strong asymmetry.

We also see asymmetrical arches in N-body simulations in which a disc is evolved in a Milky Way like potential, which includes a live dark matter halo (Figure 11). Asymmetrical arches were reported by Quillen et al. (2011) using similar simulations, but they did not study the effect of an interaction with a satellite. Laporte et al. (2019) studied simulations with an orbiting satellite and reported the presence of ridges but found very few clear arches. We studied simulations both with and without an orbiting satellite. We found that simulations in which the disc is perturbed by an orbiting satellite generates more arches. A high mass satellite generates more arches (Figure 13) than a satellite with lower mass (Figure 12). A $5\times 10^{10}M_{\odot}$ satellite was found to describe the observed data the best. Arches develop within 250 Myr of interaction with a satellite, and are clearly visible even after 1 Gyr. Hunt et al. (2018), using backward integration of test particles in a winding spiral arm potential (Dehnen, 2000), also reach a similar conclusion.

Antoja et al. (2018) used the $V_{\phi}$ separation of consecutive ridges and Minchev et al. (2009) used the $V_{\phi}$ separation of arches to conclude that the perturbation must be older than 1 Gyr and most likely about 2 Gyr. These conclusions are based on the assumption that the ridges are generated by a single perturber. If the ridges and arches are caused by more than one transient spiral arms, then each arm will have its own set of ridges and the separation between the ridges can be smaller as compared to the case of a single perturber for any given age of the perturber. Hence, the $V_{\phi}$ separation cannot be used to reliably date the perturber.

One of the most interesting results of our study is the existence of ridges in the $\langle V_{z}\rangle(R,V_{\phi})$ maps (Figure 2c). At a given $R$ when $\langle V_{z}\rangle$ is plotted as function of angular momentum or energy the ridges show up as undulations with clearly defined peaks and valleys (Figure 3). In addition to undulations, a smooth large scale trend is also seen, the $V_{z}$ increases with $L$ for $L/(V_{\rm circ}(R_{\odot})R_{\odot})>1$. This rise of $V_{z}$ has been associated with the onset of a warp (Poggio et al., 2017; Schönrich & Dehnen, 2018; Poggio et al., 2018). However, the origin of the undulations is not clear. The data shows that the locations of at least two and possibly three $V_{z}$ peaks coincide with the density peaks. This can be interpreted as ridges having a net upward motion. Undulations are also seen in profiles of $z$ with energy. Three peaks are clearly identifiable in $z$ and they match with peaks in $V_{z}$. Such a coupling of peaks between $V_{z}$, $z$ and $V_{\phi}$, is also observed in our N-body simulations, of both the unperturbed and the perturbed disc, but infrequently. We could see such a coupling for only a few locations around the simulated galaxy rather that at all locations.

A 3D phase-mixing simulation with an initial dispersion of $10$ km s*-1*in $V_{z}$ was unable to reproduce the ridges in the \langle V_{z}\rangle$$(R,V_{\phi}) maps. This suggests that the origin of features in $V_{z}$ is dynamical with the self gravity of the disc playing a role. The simulations of both the unperturbed disc and the disc perturbed by a massive satellite show ridges in the $\langle V_{z}\rangle(R,V_{\phi})$ maps (Figure 8). For the two cases of the perturbed disc, the profile of $z$ as a function of orbital energy is also found to show undulations (Figure 9). For the case of the high-mass perturber, the most prominent peak in $\langle V_{z}\rangle$ shows a clear match with the most prominent peak in $\langle z\rangle$. While an interaction with an orbiting satellite can induce coupling between planar and vertical motions, e.g. they are known to generate warps, the case of an unperturbed disc generating such a coupling is intriguing. However, Masset & Tagger (1997) have shown that non-linear coupling between the Galactic spiral waves and the warp waves can lead to outer warps in isolated disc galaxies co-existent with corrugations (undulations) over the inner disc. We propose to investigate this important insight further in the next paper.

Another example of coupling between the vertical and planar motion is the existence of the phase-spiral in the $(z,V_{z})$ and $(z,V_{R})$ planes (Antoja et al., 2018; Bland-Hawthorn et al., 2019). This phase-spiral is seen in density maps, $\langle V_{R}\rangle$ maps and $\langle V_{\phi}\rangle$ maps. Such a coupling can be generated by the impact of a satellite passing through the disc (e.g. Binney & Schönrich, 2018), or due to the buckling of the bar (e.g. Khoperskov et al., 2018). So far there has been no observation or simulation that suggests any link between the arches and ridges, and the $(z,V_{z})$ phase-spiral.

We note that the average $z$ or $V_{z}$ integrated over a single ($z,V_{z}$) phase-space spiral is non-zero and depends upon the orientation of the spiral. So if the orientation of the spiral changes with $L_{z}$, we can expect a change of $\langle V_{z}\rangle$ with $L_{z}$. We find that $L_{z}$ is a more robust quantity to characterize the phase-space spiral compared to $V_{\phi}$. This is because the spiral pattern for a given $L_{z}$ (or orbital energy) is almost invariant with the Galactocentric radius $R$ (also with azimuth $\phi$; Figure 15), within a distance of around 1 kpc around the Sun, but $V_{\phi}$ is not. When the $(z,V_{z})$ plane is studied for different values of $L_{z}$, we find that the spiral pattern is present for a wide range of $L_{z}$ and the orientation of the spiral changes with $L_{z}$ (Figure 16). However, the density distribution along the spiral is not constant and this can override any signatures in $\langle V_{z}\rangle$ generated by the spiral. Therefore, at this stage it is difficult to establish any link between the phase-space spiral and the ridges or the warp.

To conclude, there are many competing and interlocking dynamical processes occurring in the Galaxy. We have a bar, which leaves its imprint on the kinematics through resonances. We have multiple spiral arms, which are thought to be transient and can generate multiple features in kinematics. We have a warp which can couple planar and vertical motions. Other than the above mentioned internal mechanisms to excite the disc, there are also external mechanisms like interaction with orbiting satellites, e.g., Gaia–Enceladus, Sgr, LMC, SMC, and so on (Helmi et al., 2018). These can also couple vertical and planar motions. Clearly, it is important to understand and study the effect of these mechanisms individually. However, in future, we need to devise ways to study the different mechanisms together as the combined effect of mechanisms could be very different.

Acknowledgments

We would like to thank the anonymous referee for their insightful comments.

The GALAH survey is based on observations made at the Australian Astronomical Observatory, under programmes A/2013B/13, A/2014A/25, A/2015A/19, A/2017A/18. We acknowledge the traditional owners of the land on which the AAT stands, the Gamilaraay people, and pay our respects to elders past and present. Parts of this research were conducted by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013.

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.

SK is supported by a Faculty of Science Postgraduate Scholarship at the University of Sydney.

JBH is supported by an ARC Australian Laureate Fellowship (FL140100278) and the ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO-3D) through project number CE170100013. SS is funded by a Dean’s University Fellowship and through JBH’s Laureate Fellowship, which also supports TTG and GDS. MJH is supported by an ASTRO-3D Fellowship. JK is supported by a Discovery Project grant from the Australian Research Council (DP150104667) awarded to JBH.

SB acknowledges funds from the Alexander von Humboldt Foundation in the framework of the Sofja Kovalevskaja Award endowed by the Federal Ministry of Education and Research.

HST-HF2-51425.001 awarded by the Space Telescope Science Institute.

JBH & TTG acknowledge the Sydney Informatics Hub and the University of Sydney’s high performance computing (HPC) cluster Artemis for providing the HPC resources that have contributed to the some of the research results reported within this paper. Parts of this project were undertaken with the assistance of resources and services from the National Computational Infrastructure (NCI), which is supported by the Australian Government.

This research has made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al., 2018). This research has made use of NumPy (Walt et al., 2011), SciPy, and MatPlotLib (Hunter, 2007).

Appendix A Gaia SQL query

Select * from gaiadr2.gaia_source G where G.parallax IS NOT Null AND G.parallax_error/G.parallax < 0.2 AND G.parallax > 0.
where G.radial_velocity IS NOT Null