The disturbed outer Milky Way disc

The outer parts of the Milky Way's disc are significantly out of equilibrium. Using only distances and proper motions of stars from Gaia's Early Data Release 3, in the range |b|<10{\deg}, 130{\deg}<l<230{\deg}, we show that for stars in the disc between around 10 and 14 kpc from the Galactic centre, vertical velocity is strongly dependent on the angular momentum, azimuth, and position above or below the Galactic plane. We further show how this behaviour translates into a bimodality in the velocity distribution of stars in the outer Milky Way disc. We use an N-body model of an impulse-like interaction of the Milky Way disc with a perturber similar to the Sagittarius dwarf to demonstrate that this mechanism can generate a similar disturbance. It has already been shown that this interaction can produce a phase spiral similar to that seen in the Solar neighbourhood. We argue that the details of this substructure in the outer galaxy will be highly sensitive to the timing of the perturbation or the gravitational potential of the Galaxy, and therefore may be key to disentangling the history and structure of the Milky Way.


INTRODUCTION
The European Space Agency's Gaia mission (Gaia Collaboration et al. 2016) has allowed us to see with new clarity the extent to which the Milky Way's disc is out of equilibrium. The data from Gaia includes photometry and spectra, but it is the extraordinarily precise astrometry that sets Gaia furthest apart from any other instrument.
With the releases of Gaia DR2 and EDR3 (Gaia Collaboration et al. 2018a, 2021a, new disturbances have been discovered, and ones that were already known have been mapped in far greater detail. Substructure in velocities parallel to the disc have been mapped locally and across the disc (e.g., Dehnen 1998;Famaey et al. 2005;Antoja et al. 2008;Gaia Collaboration et al. 2018b;Kawata et al. 2018;Trick et al. 2019;Friske & Schönrich 2019;Khanna et al. ★ E-mail: paul@astro.lu.se 2022). Vertical asymmetries have been found in both number counts and velocities, and associated with both the Galactic warp and with bending or breathing modes in the disc (e.g. Widrow et al. 2012;Williams et al. 2013;Schönrich & Dehnen 2018;Poggio et al. 2018;Bennett & Bovy 2019;Carrillo et al. 2019;Khanna et al. 2019;Cheng et al. 2020).
The Gaia phase spiral 1 was discovered in Gaia DR2 data by Antoja et al. (2018) and is an overdensity of stars along a spiral in the − plane, where is the direction perpendicular to the Galactic disc. It is also found as a spiral shaped variation in the average or velocities in the − plane. A number of studies have demonstrated that this can be produced as a perturbation induced by the passage of the Sagittarius (Sgr) dwarf galaxy phase-mixes over time (Binney & Schönrich 2018;Darling & Widrow 2019;Laporte et al. 2019;Bland-Hawthorn & Tepper-García 2021;Gandhi et al. 2022), and this is now the generally favoured formation scenario. Other studies have investigated how this structure persists and varies across the galaxy as a function of chemistry, position and kinematics (e.g., Bland-Hawthorn et al. 2019;Xu et al. 2020;Li 2021, Alinder, McMillan & Bensby in prep.).
In connection with Gaia EDR3, Gaia Collaboration et al. (2021b, henceforth AC21) looked at stars in the Galactic anticentre (170 • < ℓ < 190 • , | | < 10 • ), and showed that in this region the velocity distribution in the outer disc was significantly disturbed. At Galactic radii between 10 and 14 kpc stars predominantly follow a bimodal velocity distribution. One of these modes is at ≈ 10 km s −1 and is mostly populated by stars below the Galactic plane, while the other is at ≈ −5 km s −1 , is mostly populated by stars above the Galactic plane, and lags the rotation of the first group by around 30 km s −1 . We, perhaps inelegantly, refer to these two modes as 'clumps' throughout this study. The strength of the two clumps varies strongly with increasing radius, with the clump at positive becoming increasingly dominant at greater radii.
This behaviour is elegantly, if incompletely, summarised by plotting the distribution in velocity of stars at different angular momenta. AC21 showed that towards the Galactic anticentre this distribution has a clear break from negative to positive at an angular momentum around | | = 2750 km s −1 kpc, corresponding to a guiding centre radius of ∼11.5 kpc.
The AC21 study was limited to the anticentre region because in this region the velocity in the plane of the sky closely corresponds to the Galactocentric azimuthal and vertical velocities. This allowed the Gaia team to gain significant insight into the velocity structure of the outer disc of the Milky Way using a sample of stars without measured line-of-sight velocities 2 . The productive use of proper motion samples alone to study Galactic dynamics has a long history, including the work by Eggen on moving groups (e.g., Eggen 1958); the characterisation of the local velocity distribution by Dehnen & Binney (1998), Dehnen (1998) and recently for the white dwarf population by Mikkola et al. (2022); measurement of velocity asymmetry in the Galactic disc (Antoja et al. 2017); and the determination of the Milky Way's escape velocity by Koppelman & Helmi (2021).
In this study we extend the approach of AC21 to include stars over a wide area in the outer disc. Because we are not limited to stars with measured line-of-sight velocities, we can use a far larger sample of stars than would otherwise be available. This allows us to explore the region with a level of detail that is not otherwise possible. We investigate whether the bimodality of the velocity distribution is only found around the Galactic anticentre, or if it can be found elsewhere in the galaxy. We ask whether we can see any change in it as we move around the galaxy.
The paper is organised as follows: In Section 2 we describe the dataset we study and the approximations we make to allow us to draw conclusions from these data. We also demonstrate the validity 2 Throughout this text we will refer to 'line-of-sight' velocities rather than 'radial' velocities to avoid confusion with the velocity component in the Galactocentric radial direction. of these approximations on mock data. In Section 3 we look at the disturbances in this part of the outer Milky Way disc. In Section 4 we present -body simulations of an impulse-like interaction with the Milky Way disc, and show that it can produce similar behaviour in the outer disc to that which we see. Finally we discuss and summarise our results in Secs. 5 & 6, respectively.

DATA
In this work we investigate the dynamics of stars in the outer galaxy, with Galactic longitudes in the range 130 • < ℓ < 230 • and latitudes | | < 10 • . The range of Galactocentric azimuths covered by these data are illustrated in Fig. 1. The data are taken from Gaia EDR3 and consist of the proper motions * , ; and the photometric magnitudes , BP , RP . We take the distance to the stars to be the photogeometric estimates from Bailer-Jones et al. (2021). These estimates take into account the observed colours and magnitudes of stars, along with their parallaxes (corrected following the recipe from Lindegren et al. 2021) and a prior constructed from a threedimensional model of our Galaxy, to derive distance estimates that can be significantly more precise than those derived without photometry.
Importantly, we do not restrict ourselves to only those stars with measured line-of-sight velocities. Gaia EDR3 contains over 1.4 billion stars with astrometric measurements, of which only 7 million (i.e., less than 0.5 percent) have published Gaia line-ofsight velocities. After Gaia DR3 in June 2022 this number will rise, but only to around 33 million (i.e., around 2.5 percent of the full catalogue). Our choice not to use them, therefore, enables us to work with a very large sample of stars, but at the cost of losing one component of the velocity. As we shall demonstrate, we are able to make approximations that allow us to study the velocity The quoted numbers are numbers of stars in a pixel of width 0.025 mag in colour and height 0.1 mag in absolute magnitude. The arrow shown at the top right corner of the diagram illustrates the assumed reddening vector. Stars will be moved in this direction on the diagram depending on whether the extinction to them has been over-or under-estimated. distribution in the vertical and azimuthal directions with just these data. However these approximations become less reasonable away from = 0 • and ℓ = 180 • .
To ensure that we are using high quality data we require that the astrometric 'Renormalised Unit Weight Error' (RUWE) meets the criterion RUWE < 1.4. We also require that the parallax divided by its error / > 3, so that the photogeometric distance estimation is not completely dominated by the photometric information. With these criteria, we have 20 041 385 stars in our full sample. Bland-Hawthorn et al. (2019) argued that the kinematically colder populations in the Milky Way disc are more strongly affected by the phase spiral than kinematically warmer populations. Further, several authors (e.g., Amôres et al. 2017;Poggio et al. 2018;Romero-Gómez et al. 2019) have shown that the structure of the Galactic warp depends on the age of the tracer population considered. We therefore construct a sub-sample dominated by young blue main-sequence stars as indicated in Fig 2. Extinction is taken as the median extinction at the quoted distance in the dust map of Green et al. (2019), and the conversion from ( − ) to the Gaia bands is approximated as = 2.294 ( − ); ( − ) = 1.309 ( − ), figures we have taken from Sanders & Das (2018).
We define our blue sample with a simple cut in extinctioncorrected absolute magnitude and colour as illustrated in Fig.2. This leaves us with a sample of 1 093 192 stars. The selection of these stars is advantageous for two further reasons: Firstly, because the stars are relatively bright, their Gaia parallax measurements will have relatively small uncertainties; secondly, because their velocity dispersion in the Galactocentric radial direction is typically very small, the uncertainty due to the absence of a measured line-of-sight velocity is reduced (see below).

Velocity estimates
If we have a measured distance and proper motion of a star, then we know two components of the star's velocity, which we can decompose into the direction of increasing Galactic coordinates ℓ and (we refer to these as ℓ and respectively). Without a measured line-of-sight velocity for a star we cannot know its full 3D velocity.
To understand motions in a galaxy it is convenient to consider velocities in terms of their Galactocentric components in cylindrical polar coordinates ( , , ). The relationship between ( ℓ , ) and ( , , ) depends on a star's position in the Galaxy, and towards the Galactic anticentre (ℓ, = 180 • , 0 • ) and centre (ℓ, = 0 • , 0 • ), we have ℓ = ± and = , while proper motions provide no information about . This convenience has led a number of authors to focus studies in these regions (e.g., Kawata et al. 2018;Fernández-Alvar et al. 2021, AC21).
At other Galactic coordinates we do not have this convenient relationship. However, we can make further progress by making a simple approximation: That the typical ≈ 0 at any point in the disc. There is a dispersion around this value, which will produce a corresponding uncertainty in the derived values which increases as we move further from ℓ = 180 • . We know that even the average is not exactly zero across the Galaxy, with Gaia DR2 in particular demonstrating that it varies at the level of ∼10 km s −1 in the extended Solar neighbourhood (Gaia Collaboration et al. 2018b), while Eilers et al. (2020) showed that it varies with an amplitude of ∼ 7 km s −1 in the outer disc. However, this is much smaller than the typical velocity in the disc, so makes a fairly small difference. We quantify the errors introduced by this approximation in Section 2.2.
For a star in any position in the Galaxy with measured ( ℓ , ), there is a line-of-sight velocity, which we call * , that would give it = 0. We can show that this velocity * = ⊥ + ⊥ ( cos + sin ) cos , where ( , , ) are Cartesian components of a star's position in the Galaxy 3 , and ( ⊥ , ⊥ , ⊥ ) are the components of the star's Galactocentric velocity in the , , coordinate system that come only from its velocity in the plane of the sky. We can use this to find our estimates for and , which are * = (− ⊥ sin ℓ + ⊥ cos ℓ) A full derivation of these equations is in Appendix A.
One could consider taking more sophisticated approaches than the approximation that = 0. The value of could be estimated as a function of position using the stars for which there are measurements of the line-of-sight velocity. This would have to be done carefully to avoid adding excess noise, especially far from the Sun where there are few such stars. It would also be taking its estimate from a different (brighter) population than the population being considered. Recently, it has been proposed that a Bayesian neural network can be used to estimate the missing line-of-sight velocities for Gaia stars (Naik & Widmark 2022), but it is not clear    whether this provides a better point estimate than our approach. In this study our simple approach is sufficient, and has the benefit of being straightforwardly reproducible and less susceptible to unexpected systematic errors. We defer the study of improvements we could make for this approximation to future work.
We have to assume a position and velocity for the Sun in the Milky Way to convert our measurements into the Galactocentric frame. We assume the Sun is located a distance = 8.178 kpc from the Galactic centre (Gravity Collaboration et al. 2019) and a height above the Galactic plane of = 0.02 kpc (Bennett & Bovy 2019). We further assume ( , , , , , = (−11.1, −248.5, 7.25) km s −1 for the solar motion (Schönrich, Binney & Dehnen 2010;Reid & Brunthaler 2020).
Throughout this paper we indicate that a quantity has been estimated without a line-of-sight velocity measurement by giving it an asterisk. Our estimate of the angular momentum of a star is * ≡ * , which for most stars in the Milky Way disc is negative.
In the interests of making our results easier to interpret, we augment this throughout this article with a simple estimate of the guiding centre radius * = | * |/(236 km s −1 ).

Demonstration on mock data
To test our assumption that we can proceed using only proper motions, we apply our method to mock data with known properties. We base this on a dynamical model constructed using the software package (Vasiliev 2019). This model consists of two disc components and a stellar halo, and is designed to be in equilibrium within a given axisymmetric gravitational potential (the main model from McMillan 2017). We sample 40 000 000 particles from the model which, since the model is axisymmetric, we can place at any position angle in the Galaxy we are interested in. The samples we construct from this model, which come from 10 • ranges in ℓ between 130 • and 230 • , contain between 2 800 000 and 3 200 000 particles when we apply the constraint that they also feature only stars with | | < 10 • . We perturb the vertical velocities of stars in this model to (very approximately) mimic the structure seen by AC21. We change the vertical velocity of each star by an amount Δ km s −1 that is a function of the component of angular momentum, . For convenience we express this in terms of = | |/( km s −1 kpc): To demonstrate that our method is robust against the variations in in the Galaxy, we further perturb the model with a term that mimics the influence of the spiral arms on the velocity field of the outer Milky Way, guided by the results from Eilers et al. (2020). We assume that the influence of the spiral arms is primarily on , and perturb the velocities of all stars in our model by an amount Δ corresponding to an = 2 spiral perturbation with  These are shown for 10 • ranges in ℓ, ranging from 130 • to 180 • for the upper set of panels and 180 • to 230 • for the lower set of panels. In each case this is also divided into the distribution above and below the plane. Vertical dotted lines are again drawn at * = −2500 & −3000 km s −1 kpc to guide the eye and to ease the comparison between equivalent plots elsewhere in the text. The green circle in the panel for > 0 at 180 • < < 190 • highlights a feature which can be seen in a number of panels. The 'break' seen in the anticentre region at * ≈ 2750 km s −1 kpc is clearly visible in many panels, with noticeable differences above and below the plane.
where we take amplitude max = 7 km s −1 , 0 = 0 • and pitch angle = 12 • . This produces a spiral perturbation that is comparable to that found in the outer disc by Eilers et al. (2020), or the variation found by Cheng et al. (2020).
Finally, we give each star a relative distance uncertainty of 15%, which is slightly larger than the median relative uncertainty of the Bailer-Jones et al. (2021) distance estimates. We then make our velocity estimates * and * for these mock data.
The results are illustrated in Figure 3, which shows the angular momentum-vertical velocity space ( * , * ), coloured by density in * at each * , for bins between ℓ = 130 • and 230 • . The choice to show the density in * rather than in ( * , * ) means that we can focus on the changes in * with * , rather than having the structure of the plot being dominated by the decreasing number of stars as a function of * . These mock observations can be compared to the model distribution, which is shown in Fig. 4. We also show in each panel the typical uncertainty in and (defined by the 16th and 84th percentile of the differences between the true values and the derived values).
It is obvious from first inspection that the structure in this figure is most clear for bins near the anticentre. Similarly, the measurement uncertainties grow as we move further from the anticentre. This is completely expected given that the approximation described in Sec 2.1 is insensitive to the real line-of-sight velocity at ℓ = 180 • and becomes increasingly sensitive to it as we move further away across the sky. Another subtler reason for this is that the influence of distance errors becomes greater as we move further from ℓ = 180 • , because a given Galactocentric radius corresponds to an ever increasing distance from the Sun.
Despite the smearing out of the underlying structure, we are still able to see some sign of it out to 50 • from the anticentre. This informs our understanding when we look at the data -we should expect perturbations of the galaxy to appear significantly smoothed further from the anticentre, not introduced or accentuated by our approximations.
In these tests, the distance errors (and the velocity errors they introduce) are the dominant source of blurring of the features of the distribution out to around ℓ = 150 • or 210 • , at which point the uncertainties due to having no line-of-sight velocity measurements begin to dominate. We demonstrate this in Appendix B, where we isolate the influence of the two effects separately. We also show that the influence of the spiral perturbation is negligible. Reducing distance uncertainties is almost as important as adding line-of-sight , and the feature it highlights is weaker for this kinematically cool population. The 'break' at * ≈ 2750 km s −1 kpc is even more clearly visible than for all stars, and we now see that it extends to ℓ = 130 • , but still shows no sign of being present for ℓ > 210 • . measurements for stars in our efforts to find this structure, even far from the anticentre.

THE OUTER MILKY WAY
First, in Figure 5 we show, for our full sample, angular momentumvertical velocity space ( * , * ), coloured by density in * at each * . This is separated into panels each showing a 10 • range in Galactic longitude ℓ, and also split above and below the Galactic plane. The feature found by AC21 is again visible in these panels, because in many of the panels there is a downturn in * at * ≈ 10 kpc ( * ≈ 2400 km s −1 kpc) which goes down to a minimum at * ≈ 11.5 kpc ( * ≈ 2750 km s −1 kpc) before what appears to be an abrupt break, with * jumping to a positive value. As we look in bins of increasing ℓ, this pattern (a downturn followed by an abrupt break) can first be seen below the plane in the 140 − 150 • bin (the high * part of the 130 − 140 • bin is very poorly populated below the plane). The break appears above the plane for the 150 − 160 • bin, but remains stronger below the plane up to the 160 − 170 • bin. The break is then stronger above the plane for increasing ℓ, and has faded away below the plane by about 190 • and above the plane by about 210 • . This extends our understanding beyond the results of AC21. In that study it was found that the downturn towards lower * was dominated by stars below the plane. We see the same for stars in the area of the sky probed by AC21 (170 • < ℓ < 190 • ), but can now see that this is not consistent across the sky.
The panels of Figure 5 also have another feature worth noting, which is an overdensity near * = −20 km s −1 at * ≈ 10 kpc, particularly clear for the panels at > 0 for 180 • < ℓ < 200 • , and visible elsewhere. We indicate this feature with a circle in one panel of this figure. This feature is faintly visible in the equivalent plot in AC21, but was not commented upon. This feature is investigated further by Alinder, McMillan & Bensby (in prep.).

Young population
Previous studies have shown that the phase spiral is stronger for stars on kinematically colder orbits, i.e., younger, bluer, stars. AC21 also showed that the break at * ≈ 11.5 stands out very clearly in their young population (their Fig. E.8). We therefore show in Figure 6 the ( * , * ) plane for the bright young main-sequence sample described in Sect.2, again coloured by density in * . We note also that these are all bright stars and therefore their Gaia Figure 7. Velocity distributions in the * − * plane for stars that are divided into bins that are 20 • -wide in ℓ, 0.5 kpc in , and contain stars on one side or another of the Galactic plane (as indicated on the left in blue, at the top, and on each panel, respectively). Dotted grid lines are added to ease comparisons between panels. Bimodal structure is visible in many panels, especially below the plane, with the gap being approximately at * = 2750 km s −1 kpc, which is marked on each panel with a red cross at * = 0. parallax uncertainties will be relatively small, allowing for more accurate distance estimates even far from the Sun.
The vertical velocity dispersion of these sources is smaller than for the whole population, which is clear in all of these plots. As seen by AC21, this means that the break at * ≈ 11.5 stands out clearly, and we have extended this result to a far wider range in ℓ. The feature near * = −20 km s −1 at * ≈ 10 kpc is harder to see for this limited population, presumably because it has a low vertical velocity dispersion so does not sample the phase space around * = −20 km s −1 well.
Using these stars as tracers we can see for the first time a hint that the break extends to ℓ = 130 • , and therefore perhaps even beyond. The panel for 130 • < ℓ < 140 • appears to have a break at around * = 12 kpc, but there is no sign of a downtrend before this. There are still no real signs of any break for ℓ > 210, so we tempted to say that it has disappeared. Certainly we can say with some confidence that it is weaker at ℓ = 210 than at ℓ = 150 because, by symmetry, the uncertainties at these two longitudes should be very similar (assuming differences due to differing extinctions are negligible). However, we are very much in the position of having limited data with which to make these judgements, and a more detailed look, and a clearer perspective beyond this range in ℓ, will have to wait for a larger sample of stars with measured line-of-sight velocities.

Velocity distributions
Finally, in Figure 7 we show the velocity distribution in the * - * plane for 20 • slices in ℓ, separated into radial bins of width 0.5 kpc from 10.5 kpc to 13.5 kpc, and further divided into stars above and below the plane.
These plots confirm the results seen by AC21, that the break seen in the * − * plane can equally be thought of as being made by a bimodal velocity distribution, consisting of two clumps of stars in the * − * plane, one tending to move downwards with a slower rotational velocity, one moving upwards with a faster rotational velocity.
To guide the eye, we have put a cross in each panel of these figures at an angular momentum of 2750 km s −1 kpc (for a star at a radius in the middle of the range). The transition we see happens at around this angular momentum, and the gap between the two clumps is at approximately this point in many cases. This confirms that there is a connection between the structure seen here and that seen in Figs 5 & 6. However it also demonstrates that we lose some information by treating the disturbance as a function of , because in some panels the gap in the velocity distribution is at different * . For example, a number of the panels for (12.5 < < 13) kpc have a substantial clump of stars at the slower rotational velocity and lower * that sits almost exactly at * = 2750 km s −1 kpc. Therefore it is clear that the exact position of the break between the clumps is not simply a value of , and the behaviour of the outer Galaxy in this region can not be completely simplified as in Figs. 5 & 6. The break into two clumps can be seen below the plane in all ℓ ranges from 130 • to 210 • , with varying strength and clarity. Above the plane it is only really clear in the range 150 • < ℓ < 190 • .

SIMULATIONS
In an effort to interpret what we are seeing in the outer Milky Way, we turn to a numerical simulation of the interaction of an impulsive mass with a cold stellar disc at a single transit point. This is based on the simulations presented by Bland-Hawthorn & Tepper-García (2021), with the only differences being that we have reduced the mass of the perturbing particle from 2 × 10 10 M to 10 10 M and reduced the number of particles. The reduced number of particles allowed us to explore the parameter space of impacts before settling on the simulation we present here, while still allowing us to resolve the behaviour we sought to examine. We found that a 2 × 10 10 M perturber produced an effect in the outer disc that was significantly larger than that seen in our Milky Way data, while the 10 10 M simulation provided a closer match to the perturbation that we see, as the scale of this perturbation scales approximately linearly with the mass of the perturber. We do not claim to have made an exhaustive search of the possible impacts or timings, to have pinned down the mass of the Sgr dwarf when it last had a major impact on the Milky Way disc, or created an exact analogue of the behaviour seen in Sec 3. That work is reserved for future studies.
The interested reader is directed to Bland-Hawthorn & Tepper-García (2021) for a full description of this simulation but we will summarise the key elements here. The galaxy model consists of a live -body dark-matter halo, stellar bulge and stellar disc with a total mass of ∼1.45 × 10 12 M and a nearly flat rotation curve out to = 20 kpc. It is constructed from an action-based distribution function using such that it is stable in its initial configuration without any relaxation, and its disc is kinematically cold to enhance the signature of perturbations. We populate this model with 10 million particles in the halo, 5 million particles in the disc and 2.25 million particles in the bulge, then evolve it using the -body code (Teyssier 2002). The perturbation is modelled as a nearly impulsive interaction of a point mass which intersects the Galactic plane at ≈ 18 kpc travelling with a velocity of ∼330 km s −1 , and a Galactic azimuth = 0. 4 The disc crossing occurs at a simulation time ≈ 100 Myr, and after = 150 Myr the mass of the perturber begins to exponentially decrease, such that by the time of its second disc crossing it has a mass of ∼ 10 −6 times its original mass. This simplifies the encounter and allows us (like Bland-Hawthorn & Tepper-García 2021) to focus on the effect of a single impact.
Bland-Hawthorn & Tepper-García (2021), following the analytical models of Binney & Schönrich (2018), showed that the impact of a perturber like this is able to produce a phase spiral in the Solar neighbourhood that is qualitatively similar to that found in Gaia data, and that it will have produced similar phase spirals across the disc (mainly differing in their phase). Here we focus instead on the influence on the outer disc and the effects that we have shown and discussed in Sec. 3. Specifically, we focus on a circular sector in the galactic disc of width Δ = 40 • , centred on = 0. This is roughly the same range in azimuth covered by our Milky Way data (see Fig. 1 for a comparison). Figure 8 shows the velocity distribution in this sector for 3 selected snapshots at = 0.35, 0.68 and 1.08 Gyr in our simulation. Like Fig. 7 we show the velocity distribution in radial bins of width 0.5 kpc from 10.5 kpc to 13.5 kpc, and show the distribution separately for stars above and below the plane. These snapshots are selected because they are at times when the velocity distribution in this part of the Galaxy separates into two clumps of the kind seen in the Milky Way data in Fig. 7.
Of all of these velocity distributions, it is the one at 0.35 Gyr    . Velocity distributions of a 20 • segment of the outer disc of the simulated galaxy at three selected times (as indicated, note that the disc crossing of the perturber occurs at = 100 Myr). The panels, equivalent to those in Fig 7, show density in the -plane for radial bins between 10.5 and 13.5 kpc, divided into particles above and below the disc plane. Again, dotted grid lines are added to ease comparisons between panels and with Fig 7. In all these selected times, the velocity distribution is divided into two clumps, but only in the first case are these two clumps displaced from one another in .
(i.e, 250 Myr after the disc crossing of the perturber) that bears the strongest resemblance to the velocity distribution seen in the Milky Way. The two clumps are at different , with the higher | | clump being at positive and the lower | | clump at negative . There is a clear difference between the behaviour above and below the Galactic plane. There is a clearer distinction between two clumps in the velocity space for < 0, and the lower velocity clump has essentially disappeared by the 12.5-13.5 kpc bin for < 0 while this is not true above the plane.
The other two snapshots do show a velocity distribution divided into two clumps, but they do not have a clear difference in their values. We have run this simulation to nearly 4 Gyr, long after the perturber has reduced to a negligible mass and have found that the velocity distribution continues to divide into two clumps, then Z < 0 Figure 9. The distribution in as a function of angular momentum in the outer disc of the simulated galaxy at the same three selected times as Fig. 8. The panels, equivalent to those in Figs. 5 & 6, are again divided into the distribution above and below the plane, and we have again drawn dotted lines * = −2500 & −3000 km s −1 kpc to ease the comparison between these equivalent plots. We see corrugation of the distribution, which appears over shorter and shorter scales as it winds up over the period shown here, but we do not see a clean break of the kind seen in the Milky Way data.
coalesce into a single clump, then divide into two clumps again repeatedly throughout the simulation time. However, the division does becomes weaker, concentrated in the outer bins, and does not show a clear asymmetry either side of = 0 or difference in for the two clumps.
To understand this better, in Fig. 9 we look in the simulations at the distribution in as it changes with , again in the sector of width Δ = 20 • around the assumed Solar azimuth in the Galaxy, and in the same times as the panels of Fig. 8. In each case there are disturbances in the vertical velocities that are of similar scale to the disturbances seen in the Milky Way. We also see, particularly in the two later snapshots, breaks in the distribution which are comparable to those seen in the Milky Way. The panel at = 0.35 Gyr does not show a clear break in the distribution, but, as we have seen in Fig. 8, is the snapshot with a velocity distribution that best matches that seen in the Milky Way.
We would also note that the panels show that the wavelength (in angular momentum) of the disturbance is relatively long shortly after the disc crossing, and gets increasingly short as time passes. This occurs because the disturbance winds up around the disc.
The choice of position to centre this analysis is, to an extent, arbitrary. The initial disc is axisymmetric, and Bland-Hawthorn & Tepper-García (2021) used their simulation to explore the behaviour of the phase spiral at a range of positions in the Galaxy. We have explored other choices of and found that the results are, qualitatively, rather similar. The times at which the velocity distribution forms into two clumps varies somewhat with , as the disturbance wraps its way around the disc. The general pattern remains the same as for the three snapshots illustrated in Figs 8 & 9. For the interested reader we have made videos showing the time evolution of the plot shown in Fig. 8 for choices of at 40 • intervals available at https://www.astro.lu.se/~paul/MyResearch.

Galactic seismology
For over a century, geologists have learned about the interior of the Earth through seismology, studying the waves that pass through it. More recently, helioseismology has revolutionised our understanding of the Sun, and asteroseismology has brought us great insight into the stars of our Galaxy. We have now entered the era of Galactic seismology (or Galactoseismology), with much to learn about our Galaxy from its oscillations. 5 Galactic seismology differs from the other fields in many ways, of which the two most important are perhaps timescale, and influence on the medium through which the waves are passing. The timescale of the oscillations in the Galaxy are measured in Myr, as opposed to oscillations of the Earth or stars, which can be measured in minutes or hours. We therefore can not wait to follow the oscillations of the Galaxy as they go through their phases. The oscillations of the Earth and the stars can be treated as non-invasive -not changing the underlying structure of the medium through which they are travelling. This is not true for the oscillations of the Galaxy. Darling & Widrow (2019) note, following the argument of Hunter & Toomre (1969), that the influence of the perturbation on the unperturbed disc is of the same order as the influence of the unperturbed disc upon the perturbation. They further show that the influence of self gravity slows the growth of the phase spiral. This result that can also be seen by comparing the phase spiral found by Bland-Hawthorn & Tepper-García (2021) with the one seen in the analytic models of Binney & Schönrich (2018) which winds up quickly (see also Hunt et al. 2021).
The bimodality of the velocity distribution in the outer Galaxy discovered by AC21 was, like the phase spiral before it, not predicted. Here we have shown that it can be explained, at least qualitatively, by a model that also produces a phase spiral. These different measurements, or phase space projections in which we conduct our analysis, can be thought of as two separate Galactic seismographs, but they must be intrinsically linked and understanding of our Galaxy will come from models that can consistently explain both.
For now we find ourselves like blindfolded people discovering an elephant, each describing what is right in front of them (a trunk, a tail, a leg, etc.) and not being able to agree on what they have found. 6 Similarly we will probably need to simultaneously explain the phase spiral, the velocity bimodality, the Galactic warp, the ridges in − space (Kawata et al. 2018;Ramos et al. 2018, AC21), the variations in as a function of angular momentum (e.g. Friske & Schönrich 2019) and Galactic (e.g. Eilers et al. 2020), and even the structures of the outer Galaxy (e.g., Xu et al. 2015;Laporte et al. 2022, AC21). This will be complicated by the fact that some of these phenomena, especially those seen in the planar velocities, are affected by the Galactic bar and spiral arms in a way that is likely to be unrelated to the vertical velocity perturbations that we see. Future data from Gaia, or spectroscopic surveys such as WEAVE (Dalton et al. 2016) and 4MOST (de Jong et al. 2019), will no doubt allow us to identify a new tusk or two, but our task will remain that of piecing all these elements into an over-all picture. This study has not described the entire elephant, but has extended our understanding of one part of it. We have shown that the velocity bimodality persists across a range of Galactic angle that goes far beyond those probed by AC21, that it changes in character with varying , and appears to fade away by ≈ 10 • (ℓ ≈ 210 • ) although, as we have discussed, this is at least partly due to the difficulty in measuring it at these longitudes. The simulations we have performed show that, like the phase spiral, bimodal velocity distributions can be induced by the passage of the Sagittarius dwarf. However, we have not carefully matched the timing of the velocity bimodality in our simulation to that of the formation and evolution of the phase spiral, nor to the past history of Sagittarius. That work is reserved for future study.

Density and bending waves
Bland-Hawthorn & Tepper-García (2021) argue that the impact of the Sgr dwarf triggers two different = 2 modes in the inner disc: a density wave, and a corrugated bending wave. They associate the phase spiral with the density wave rolling up and down the bending wave. They note in addition that the outer discs' vertical displacement is better described instead as dominated by an = 1 bending mode.
The density wave can be expected to create gaps (or simply low density regions) in the angular momentum distribution, which translate into gaps in for stars in a given radial bin. The bending wave can mean that the two clumps either side of this gap have different . This is what we appear to see in Fig. 8 for the snapshot at = 0.35 Gyr, but the snapshots at = 0.68 and 1.08 Gyr have gaps, which we associate with the density wave, which do not coincide with a significant vertical velocity difference from the bending wave -the clumps in these two cases have nearly identical vertical velocities. In this picture, the structure seen in the Milky Way can only be produced by a coincidence of a minimum point in the density wave and a point where the gradient of in the bending wave is very large.
It is certainly questionable whether this is exactly what we are seeing in the outer Milky Way. The break in the * − * plane seen in Figs. 5 & 6 does appear to be relatively clean at a given * , but in velocity space (Fig. 7) the divide does not appear to lie along a line of constant * , as it does in the simulation (Fig. 8). Equally, we do not see a clean break in the − plane at = 0.35 Gyr, as we do in the Milky Way. There is also a very clear difference in the Milky Way between what we see above and below the plane. We do not see this here, but this may be because we are averaging over a wider range in .
The timing of the snapshot that seems to best mimic the behaviour in the Milky Way ( = 0.35 Gyr, i.e., 250 Myr after the disc crossing) is also before Bland-Hawthorn & Tepper-García (2021) see a phase spiral forming, though it is consistent with the timing of the phase spiral found using a simpler, analytic, approach by Binney & Schönrich (2018). The later snapshots do have breaks in the − plane, but we can clearly discern in each of them that there is an oscillation with an ever decreasing scale length in , which we can associate with the winding up of the bending wave across the Galactic disc. The equivalent distribution in the Milky Way does not appear to have this oscillation on short scales in , instead having a large scale variation with an abrupt break. We also use a perturber that has a mass which half that of the Bland-Hawthorn & Tepper-García (2021) simulation that reproduced the phase spiral, in order to match the scale of perturbation seen in the outer Milky Way, so cannot claim that we simultaneously model both.
Certainly we have only begun to interpret the velocity bimodality in the outer Milky Way disc. One can hope that this will prove a fruitful area, both for learning about the mass and timing of the most recent impacts on the Milky Way disc, and for measuring the properties of the disc too. If the bimodality can indeed only be explain by a coincidence of the density and bending wave then this sets a strong constraint on the impact's timing relative to orbital frequencies in the disc. We might reasonably expect that in the outer disc, as the disc itself becomes less dominant, the influence of self-gravity becomes less important, allowing us to make a unambiguous estimate of these properties and to use simpler, faster, modelling techniques of the kind employed by Binney & Schönrich (2018).

SUMMARY
In this paper we have shown that the velocity bimodality -or, equivalently, break in the − distribution -found in the Galactic anticentre by AC21 persists across some significant fraction of the Galactic disc. It can be found out to ℓ = 130 • , which corresponds to a Galactocentric angle of ≈ 20 • , but appears to be fading away by ℓ = 210 • ( ≈ 10 • ). Its behaviour is, as seen by AC21, different above and below the Galactic plane, with a stronger break/bimodality usually visible below the plane. We note that selecting young main-sequence stars allows us to see a stronger effect, which is consistent with previous results that show younger, kinematically colder, stars are more affected by disturbances in the vertical velocities of stars (e.g., Bland-Hawthorn et al. 2019).
Using rather idealized -body simulations we have shown that the passage of a Sagittarius-like dwarf galaxy is capable of producing qualitatively similar effects in a Galactic disc. These simulations are also capable of producing a phase spiral akin to the one which has been discovered in the Solar neighbourhood. This study has been conducted entirely without measured lineof-sight velocities for the stars considered. For the foreseeable future we will be in the position that we have measured proper motions for far more stars than line-of-sight velocities. This is due to both the two billion stars observed by Gaia, and the still larger number which would be observed by a successor mission conducting astrometry in the infrared (Hobbs et al. 2021). Exploiting these incomplete 5D phase-space data are absolutely key to building up large samples of stars for study, which allows us to either constrain our models far more accurately or, as in the case of the velocity bimodality, find features of our Galaxy that we have not predicted using our models. Our investigation of the uncertainties inherent in our approach makes the further point that distance uncertainties are as serious a barrier to characterising sharp structure in velocity space as the absence of line-of-sight velocity measurements.
To make further progress we have to introduce our approximation that = 0. From the above, we can see that which, under the approximation that = 0, can be rearranged to give an estimate of , which we refer to as * where * = ⊥ + ⊥ ( cos + sin ) cos (A4) This allows us to find our estimates * = ( ⊥ − * sin cos ) − ( ⊥ − * cos cos ) We can find the error associated with a star having a true Figure B1. Velocity perturbation in applied to our mock data. The perturbation is comparable to that seen in the Milky Way, but somewhat larger, to provide a more stringent test. by placing this in eq. A3 and following the equations through. We find errors = − * and = − * to be = sin − cos cos + sin (A7) and = − tan cos + sin .
(A8) runs from ∼ − 3 at ℓ = 130 • to ∼3 at ℓ = 230 • , with a mild dependence on distance (in the sense that more distant stars are worse affected).

APPENDIX B: DETAILS OF THE MOCK DATA
The dynamical model on which we base the experiments in section 2.2 has two discs and a halo, each of which is defined such that it is in equilibrium within the axisymmetric Milky Way gravitational potential taken from McMillan (2017). We use the Stäckel fudge approach of Binney (2012) to derive the actions = ( , , ≡ ) as a function of position and velocity.
The discs are based on the quasi-isothermal distribution function (Binney & McMillan 2011) as modified in , , with ( , , Ω) being the epicycle frequencies.
• is an estimate of the typical orbital radius, which is given by radius of a circular orbit with angular momentum circ = + + /4. The surface density of this disc is Σ ≈ Σ 0 exp(− / disc ), and the functions˜and˜, which approximately give the velocity dispersions in and (and therefore in the latter case also determine the disc scale height), are of the form 2 ( • ) = 2 ℎ 2 2 ( • ) + 2 min (B2) 2 ( • ) = 2 ,0 exp(−2 • / ) + 2 min ,  Figure B2. Density in at different for mock data with realistic distance uncertainties but with precise line-of-sight velocities and no spiral perturbation. As with Fig. 3, each panel shows the results for different 10 • ranges in ℓ. Error bars indicate the typical uncertainties in and . The error introduced is a significant fraction of the total uncertainty seen in Fig. 3 in all panels, and is the dominant source of uncertainty near the anticentre. No LOS velocity only Figure B3. Density in * at different * for mock data with no distance errors or spiral perturbation, but with no line-of-sight velocity measurement. Each panel shows the results for different 10 • ranges in ℓ. The error introduced by the lack of line-of-sight velocities is never much larger than that introduced by the distance uncertainties (Fig B2), and only becomes the dominant source of uncertainty ∼30 • from the anticentre.  where min = 5 km s −1 is introduced to avoid numerical difficulties at large • . The term ± ( ) is a function which is much larger for negative than the equivalent positive . This gives the model a negative sense of rotation, like the Milky Way. Omitting it would give the distribution function of a disc that has equal numbers of stars rotating in opposite directions.
The thick disc distribution function is taken to be a single quasiisothermal distribution, while the thin disc one is a superposition of quasi-isothermal populations taken to be of ages , weighted proportional to exp( / SFR ) in order to mimic the thin disc's expected decreasing star formation rate. The radial and vertical dispersions of the quasi-isothermal populations both vary to reflect the agevelocity dispersion relation in the form ( ) = 1 +(1− ) 1/ , with 0.33 being the power law associated with the age-velocity dispersion relation, and being the ratio of the velocity dispersion of the oldest population and the youngest population.
The stellar halo is based on the double power-law form introduced by Posti et al. (2015), as modified in where ( ) ≡ + + | |, The interested reader is directed to Vasiliev (2018Vasiliev ( , 2019) for more details. Parameters of the distribution functions of all three components are given in table B1.
The samples drawn from this model are perturbed in the way described in 2.2, which produces the disturbance in illustrated in figure 4 and a disturbance in which is illustrated in Fig B1, and is broadly similar to that found by Eilers et al. (2020).

B1 Different sources of error
Figures B2, B3 & B4 subdivide the effects which blur out the structure in the − plane. Figure B2 shows the effect of the 15 percent distance uncertainty applied to the model on its own (i.e., assuming that line-of-sight velocities are known). This causes blurring in all panels that increases gradually as we move further from the anticentre, as the distance we have to look to see stars at a given (which generally have ≈ ) grows with distance from the anticentre. Note that we give quantities as, for example, rather than * because we are not having to estimate these quantities without line-of-sight velocities. Figure B3 shows the effect of assuming that there is no measured line-of-sight velocity, even if the distance to the star is known precisely. The uncertainty introduced is very small towards the anticentre -far smaller than the uncertainty introduced by the distance uncertainty. However, the uncertainty introduced by this approximation grows with increasing angle from the anticentre faster than the uncertainty associated with the distance, and is larger by ℓ ≈ 150 • . Figure B4 show that the effect of the spiral in addition to the lack of line-of-sight velocity is so small as to be completely negligible. We do not attempt to isolate the effect of the spiral completely because, if the line-of-sight velocity is known, we no longer have to to make the approximation that = 0, so the spiral perturbation has no effect on our results. This paper has been typeset from a T E X/L A T E X file prepared by the author.