Age dissection of the vertical breathing motions in Gaia DR2: evidence for spiral driving

Gaia DR2 has revealed breathing motions in the Milky Way, with stars on both sides of the Galactic mid-plane moving coherently towards or away from it. The generating mechanism of these breathing motions is thought to be spiral density waves. Here we test this hypothesis. Using a self-consistent, high-resolution simulation with star formation, and which hosts prominent spirals, we first study the signatures of breathing motions excited by spirals. In the model, the breathing motions induced by the spiral structure have an increasing amplitude with distance from the mid-plane, pointing to an internal cause for them. We then show that, at fixed height, the breathing motion amplitude decreases with age. Next, we investigate the signature of the breathing motions in the Gaia DR2 dataset. We demonstrate that, at the location with a consistently large breathing motion, the corresponding amplitude increases monotonically with distance from the mid-plane, in agreement with the model. Furthermore, we show that at the same location, the breathing motion amplitude decreases with age, again similar to what we find in the model. This strengthens the case that the observed breathing motions are driven by spiral density waves.

The Milky Way is a barred galaxy (Weinberg 1992) which also hosts spiral structure (e.g., Gerhard 2002, also see Vallée (2005Vallée ( , 2008). The presence of large-scale, non-zero mean vertical motions of Solar Neighbourhood stars has been reported in various surveys, for example, using main-sequence stars in the Sloan Extension for Galactic Understanding and Exploration (SEGUE) survey (Widrow et al. 2012), F-type stars in the Large Sky Area Multi-Object Fibre Spectroscopic Telescope (LAMOST) survey (Carlin et al. 2013), and red clump stars in the Radial Velocity Experiment (RAVE) data (Williams et al. 2013). These vertical motions display two distinct types: one in which stars on either side of the mid-plane move coherently in the same direction along the perpendicular to the disc (bending motion) and the other in which stars move coherently in the opposite direction, thereby compressing towards the mid-plane or expanding away from it (breathing motion). In addition, Widrow et al. (2012) reported evidence for wave-like north-south asymmetries in the number counts of stars in the Solar Neighbourhood, which were further explored by Yanny & Gardner (2013).
The Gaia mission (Gaia Collaboration et al. 2016Collaboration et al. , 2018a has revolutionised the study of Galactic archaeology by measuring the stellar kinematics of stars in the Solar Neighbourhood with an unprecedented precision. The stellar kinematic study from the second Gaia Data Release (hereafter Gaia DR2) revealed rich kinematic substructure in the phase-space distribution of Solar Neighbourhood stars (Antoja et al. 2018), as well as the presence of large-scale non-zero vertical motions (∼ 10 km s −1 in magnitude), and associated bending and breathing motions (Gaia Collaboration et al. 2018b). Using the RAVE and the Tycho-Gaia astrometric solution (TGAS) catalogue, Carrillo et al. (2018) showed the presence of both bending and breathing motions in the Solar Neighbourhood. Bennett & Bovy (2019) revisited the north-south asymmetry in the Solar Neighbourhood stars using the position-velocity measurements from Gaia DR2, and confirmed the presence of periodic over-and under-densities of stars in the vertical distribution, regardless of their colours. For an axisymmetric potential, the bulk radial and vertical motions should be zero (e.g. Binney & Tremaine 2008). Therefore, the question arises what mechanism(s) can induce these non-zero bulk vertical motions in our Galaxy.
Theoretical efforts to understand the cause of vertical motions have explored both external and internal mechanisms. An interaction with a satellite or a dark matter subhalo can excite large-scale coherent vertical bending motions (e.g., see Hunter & Toomre 1969;Araki 1985;Mathur 1990;Weinberg 1991;Gómez et al. 2013;Widrow et al. 2014;D'Onghia et al. 2016;Chequers et al. 2018). However bending motions can be produced by purely internal mechanisms, such as the bar (Khoperskov et al. 2019) or the warp (Khachaturyants et al. submitted). On the other hand, internal causes such as spiral density waves (Faure et al. 2014;Debattista 2014;Monari et al. 2016) or the Galactic bar (Monari et al. 2015) as well as external mechanisms, such as a fly-by encounter with a satellite or passing of a dark matter sub-halo (Widrow et al. 2014), have been proposed for exciting the breathing motions. Earlier studies by Siebert et al. (2011Siebert et al. ( , 2012 showed that spirals and bars can drive a large-scale gradient in the bulk radial velocity in the Milky Way. Strong spirals drive large-scale vertical breathing motions (| | ∼ 5 − 20 km s −1 ) as first shown by Faure et al. (2014). Their semi-analytic models assumed a thin spiral structure with a small vertical extent. Their analytic treatment found that the breathing motions vary with distance from the mid-plane as tanh( / 0 )sech 2 ( / 0 ), where 0 is the scale-height of the spiral. This peaks at ∼ 0.7 0 . At larger heights, this decreases as exp(−2 / 0 ) (see Eqs 10 & 12 in Faure et al. 2014). The self-consistent simulations of Debattista (2014), instead, produced vertically extended spirals, which resulted in vertical motions that increase with distance from the mid-plane. The relative sense of the bulk motions, whether compressing or expanding changes across the corotation resonance (hereafter CR) of the spiral. Inside the CR, the vertical motions are compressive behind the peak of the spiral and expanding ahead of the spiral's peak. Outside the CR, the compressive and expanding breathing motions reverse their sense of occurrence (Faure et al. 2014;Debattista 2014).
In this paper, we continue to explore the possibility that the breathing motions observed in Gaia DR2 are caused by the Milky Way's spiral structure. Using a high-resolution, self-consistent, star-forming simulation of a Milky Way-like galaxy, we study the characteristics of the breathing motions driven by spiral structure. Previous studies of breathing motions either used perturbation theory (Monari et al. 2016) and test particle simulations (Faure et al. 2014) or used selfconsistent simulations but lacked the high mass resolution necessary to sub-sample the stellar populations (Debattista 2014). The novelty of our paper, aside from comparing directly with the Gaia DR2 data, lies in the fact that our high-resolution model allows us to explore how different populations participate in the breathing motions.
The paper is organised as follows: Section 2 provides the details of the simulation, including quantifying the spiral structure present at our chosen snapshot, which, for simplicity, we refer to as 'the model'. Section 3 presents the breathing motions of the model. Section 4 provides the details of a novel technique to detect/identify breathing motions from the variation of the bulk vertical motions. Section 5 compares the vertical kinematic signatures of Milky Way stars with the trends obtained in the model. Section 6 discusses the unique fingerprint of spiral-driven vertical breathing motions, while Section 7 summarises our main conclusions.

STAR-FORMING SIMULATION
We use an -body+smooth particle hydrodynamics (SPH) simulation with gas and star formation to motivate some of the analysis of the Gaia DR2 data. This is a higher mass resolution version of the models described in our previous papers (e.g. Roškar et al. 2008;Loebman et al. 2016). This model starts with a gas corona in pressure equilibrium with a co-spatial Navarro-Frenk-White (Navarro et al. 1996) dark matter halo which constitutes 90% of the mass. The dark matter halo has a virial radius 200 200 kpc and a virial mass 200 = 10 12 M . Gas velocities are given a spin of = 0.065 (Bullock et al. 2001), with specific angular momentum proportional to radius. Both the gas corona and the dark matter halo consist of 5 × 10 6 particles; gas particles have softening = 50 pc, while that of dark matter particles is = 100 pc. All stars form out of cooling gas, inheriting the softening of the parent gas particle. We evolve the simulation for 13 Gyr with (Wadsley et al. 2004(Wadsley et al. , 2017. As gas cools it settles into a disc; once the gas density exceeds 0.1 cm −3 and the temperature drops below 15,000 K, star formation and supernova blast wave feedback commence (Stinson et al. 2006). Supernovae feedback couples 40% of the 10 51 erg per supernova to the interstellar medium as thermal energy. Gas mixing uses turbulent diffusion as described by Shen et al. (2010).
We use a base timestep of Δ = 5 Myr with timesteps refined such that = Δ /2 < √︁ / , where we set the refinement parameter = 0.175. We set the opening angle of the tree-code gravity calculation to = 0.7. Gas particle timesteps also satisfy the condition = ℎ/[(1 + ) + ], where = 0.4, ℎ is the SPH smoothing length set over the nearest 32 particles, and are the linear and quadratic viscosity coefficients and is described in Wadsley et al. (2004). Star particles represent single stellar populations with a Miller-Scalo initial mass function.

Spiral structure in the model
We use a snapshot from this model at = 12 Gyr. We choose this snapshot strictly because it allows a good division of its stellar populations, not because of any special properties of the model at this time. While the simulation we use has not been published before, the general evolution of the spiral structure in this class of models is explored in Roškar et al. (2012).
In order to facilitate comparison with the Milky Way, we quantify the spiral amplitude in the model to show that the spirals in the model are not unreasonably strong. We also quantify the strength of the spiral as a function of stellar age. We divide the stars into three broad age bins: a young stellar population (with ages between 0-4 Gyr), an intermediate-age stellar population (with ages between 4-8 Gyr), and an old stellar population (with ages between 8-12 Gyr). The top panels of Fig. 1 show the surface density of stars in the face-on, ( , ), view for the three age populations. The bottom panels show the residual surface density,Σ( , ), calculated as: where Σ avg ( ) is the average density at radius , and the averaging is carried out over the whole azimuthal range.
The density plots of Fig. 1 show that the youngest population exhibits a very prominent two-armed spiral structure, which is also present in the intermediate-age population, but is indistinct in the old population. However the residual density plots show that the = 2 spiral is present even in the old population, albeit weaker. This therefore is a genuine density wave which propagates in all stellar populations. The estimated pitch angle, , of the spiral in this snapshot is ∼ 33 • . In comparison, the Milky Way's spirals are tighter, having a pitch angle ∼ 10 • (Siebert et al. 2012).
To quantify the strength of the spirals, we measure the radial profiles of the Fourier moments of the surface density, defined as: where is the coefficient of the ℎ Fourier moment of the density distribution, is the mass of the ℎ particle and is its cylindrical angle. The peak of the Fourier = 2 amplitude, ( 2 / 0 ) max 0.19 for our model. In comparison, Rix & Zaritsky (1995) found that, for their sample of galaxies with = 2 spirals, the peak value of 2 / 0 varies in the range 0.15 to 0.6, with ∼ 28% of their sample having ( 2 / 0 ) max ≤ 0.2. The spiral structure in the model therefore is sufficiently realistic to provide a basis for comparing the kinematic signature of the spiral in the model with that in the Milky Way. The amplitudes for the = 2 and = 4 Fourier harmonics for all three stellar populations at = 12 Gyr are shown in Fig. 2. The spiral has a significant = 2 amplitude at 4 ≤ / kpc ≤ 8. While the amplitude decreases with stellar age, all stellar populations take part in the spiral structure, as opposed to material spirals present in only the young stellar populations (for detailed discussion, see e.g. Binney & Tremaine 2008;Dobbs & Baba 2014). The = 4 amplitude is likewise strongest in the young population, and is almost zero for the old population. Thus the spiral in the old population varies more gently in azimuth while that in the young population is more sharply delineated, as can be seen in Fig. 1.
We also calculate the residual density distribution Σ( , ) within a radial annulus of width Δ as where Σ avg ( ; Δ ) is the average density, and the averaging is carried out over the whole azimuthal range within the radial annulus of radial width Δ , at a fixed height . Fig. 3   maps in the ( , ) plane at different annuli of Δ = 1 kpc. At each annulus, the peak density in the mid-plane ( = 0) is ahead (relative to the sense of rotation) of the peak density at larger heights from the mid-plane. Also, the density distribution at fixed and is azimuthally asymmetric relative to the peak density, in agreement with the findings of Debattista (2014). The spiral extends to the full vertical extent, rather than being limited to close to the midplane. This behaviour was also found in pure -body simulations by Debattista (2014), and will be important for our subsequent analysis.

BREATHING MOTIONS INDUCED BY SPIRALS
We start by studying the breathing motions induced by the spirals in the star-forming model in order to explore how the breathing motions vary between different age populations. The self-consistent simulation of Debattista (2014) showed that the amplitudes of the breathing motion increases with height from the mid-plane; here we extend that earlier work but exploring how the amplitudes change as a function of both height and stellar age.
We define the breathing velocity, breath , as (Debattista 2014;Widrow et al. 2014;Gaia Collaboration et al. 2018b): where ( , , Δ ) is the mean vertical velocity at position ( , ) in the galactocentric cartesian coordinate system, averaged over a horizontal layer of thickness Δ (for details see Gaia Collaboration et al. 2018b). We choose Δ = 400 pc, and calculate breath in three such vertical layers: | | = [0, 400] pc, | | = [400, 800] pc and | | = [800, 1200] pc. In this definition of breath , the average in a particular layer is carried out before the difference between the two sides is taken. In the observational data this helps reduce the effect of any selection function differences between the two sides of the disc. This also removes the trivial effect of bending waves, which would otherwise give rise to a spurious signature by displacing stars to one side of the mid-plane. Eq. 4 implies that when breath > 0 then stars are coherently moving away from the mid-plane (expanding breathing motion), while breath < 0 when stars are moving towards the midplane (compressing breathing motion). For the sake of comparison, we also consider the bending velocity: The star-forming model, which is not externally perturbed, lacks substantial coherent bending waves, so the model's bending velocities are not shown here.

Bulk vertical motions of the different age populations
We examine the star-forming model's vertical motions as a function of stellar age. The distribution of the bulk vertical motions in the meridional, ( , ), plane, calculated at four different azimuthal wedges, for the three stellar populations at = 12 Gyr, are shown in   mid-plane. There are not enough stellar particles in the young population at large heights to measure a reliable mean vertical velocity for them there. The intermediate-age population exhibits a large | | at larger heights; | | gets somewhat weaker for the old stellar population. All three populations are coherently moving away from, or towards, the disc mid-plane indicating breathing motions predominate in all stellar populations. There are also bending motions present, but these are weaker compared to the breathing motions. These bending motions are studied in detail elsewhere (Khachaturyants et al. submitted). All three age populations exhibit prominent bulk vertical breathing motions in the ( , ) plane. However, the azimuthal wedges used in Fig. 4 are large, therefore the coincidence of strong breathing motion with the spiral structure is not readily apparent in this figure. This is further demonstrated in Section 4.
We then select stars in the radial annulus 5 ≤ / kpc ≤ 7, calculate in the ( , ) plane for the different stellar populations, and plot these in annulus. The vertical motions switch from compressive to expanding at 0 • , which Fig. 1 shows is the location of the peak density in the youngest population. The amplitude of the vertical motions decreases with age, with the old population having velocities roughly half those of the intermediate population.

Breathing motions of the different age populations
We therefore compute the breathing velocities for stars of different ages in three vertical slices: | | = [0, 400] pc, | | = [400, 800] pc, and | | = [800, 1200] pc. We calculate breath using Eqn. 4 separately for each age population. Fig. 6 shows breath in the ( , ) plane whereas Fig. 7 shows breath in the ( , ) plane. Strong breathing motions are present in all three age populations; breath has small amplitude near the mid-plane and increases with distance from it. This trend is present for all ages; at the largest heights the intermediate-age population has the strongest breath amplitude, but the youngest pop-ulation does not reach these heights so can only be compared in the intermediate layer. The intermediate population at the largest height shows that the compressive motions are closely associated to the peak density of the spiral in the outer region (excluding the central weak bar).
In summary, we find large-scale, non-zero vertical motions present both in the meridional plane as well as in the ( , ) plane at a fixed radial annulus. These non-zero bulk vertical motions appear preferentially at azimuthal locations which contain the strongest spiral structure. The resulting breathing velocity, breath , is small near the disc mid-plane, and increases with height, as found also by Debattista (2014). The breathing velocity in our model is largest for the intermediate age population at the largest heights, because the young population does not reach such heights; in the middle layer the young population has the largest breath , while the old population has the smallest.

A NOVEL TECHNIQUE TO DETECT THE BREATHING MOTIONS
In the previous section, we measured the distribution of breathing velocity in the ( , ) plane using Eq. 4. This requires calculating, at a particular spatial location, the mean vertical velocity, , of stars in both the upper and the lower layers (of width Δ ) of the disc, and then taking the difference of these two mean vertical velocities. For the simulated galaxy model we have the whole spatial coverage so that we can calculate quite robustly the spatial distribution of the resulting breathing velocities for stars at different vertical layers. However, for the Milky Way, this may not be the case due to the limitations in the spatial coverage and incompleteness of Gaia DR2.
Here, we develop a new technique for identifying breathing motions from the distribution of the mean vertical velocities ( ) as a function of height from the disc mid-plane ( = 0). First, we demonstrate this technique on our star-forming model, and then apply it to the Gaia DR2 data.
Above we showed that the breathing velocities increase away from the mid-plane. Therefore to leading order, the variation of with height can be fitted with a straight line. Breathing motions result in best-fit straight lines to the height variation of the mean vertical velocity that has a significantly non-zero slope. If the slope is positive, then it indicates an expanding (positive) breathing motion, while if the slope is negative then it indicates the presence of a compressing (negative) breathing motion.
To test this, we divide the stellar particles into two age populations: younger (1 − 6 Gyr), and older (6 − 12 Gyr) stellar particles. We exclude the very young stellar particles (age less than 1 Gyr) because of their limited vertical coverage. Fig. 8 shows a few examples of fitting a straight-line to the variation of with height. The straight-line is fitted via the package which uses the Levenberg-Marquardt algorithm. The non-zero slope indicates breathing motions with amplitude increasing with height from the mid-plane. The (absolute) value of the slope is larger for the younger stellar population than for the older stellar population. This implies that the breathing motions are stronger in the younger population than the older populations, in agreement with the findings of the previous section. We note that, at certain locations, there is a clear difference in the zero-point of the different age groups. When studied in detail, we found that these zero-point offsets arise because of weak bending waves, which have different amplitudes in the different age groups (Khachaturyants et al. submitted).
We compute the error on the slope ( ( )) and plot this in the   middle row of Fig. 9. For uniform comparison, we impose a simultaneous quality cut of / ( ) > 3 for both age populations, and choose only those bins where this quality-cut is met. The resulting distribution of the best-fit slopes is shown in the bottom panels of Fig. 9. Clear signatures of breathing motions (with amplitude increasing with height) are seen. The best-fit slopes for the younger population (1 − 6 Gyr) are, in general, larger than those for the older population (6 − 12 Gyr).
Lastly, we investigate how the slope varies as a function of the azimuthal angle for the young and the old populations relative to the density variation. To achieve this, we first choose a radial extent ranging from 6 kpc to 8 kpc (i. e., outside the CR of the spiral), where the model exhibits a strong spiral. It is evident from Fig. 1  (bottom panels) that the azimuthal locations of the density peaks vary as a function of radius in the chosen radial extent, because of the spiral's pitch angle. Therefore, to obtain a stronger signal of the slope, we rotate the stars in different 1 kpc-width radial bins, in such a way that the density peaks in our chosen radial extent coincide. Then, we recalculate the slope of the breathing velocity as a function of the rotated azimuthal angle ( ) for the young and the old stellar populations. This is shown in Fig. 10. We also calculate the total residual surface density (Σ( , )) as a function of the rotated azimuthal angle in this chosen radial extent. As evident from Fig. 10, the locations of the larger slopes of the breathing motions coincide with the azimuthal locations of the density minima while the azimuthal locations of the maximum density exhibit smaller values of the slope. In other words, the amplitudes of the breathing motions are larger in the inter-arm region. This trend is consistent with the earlier findings of Debattista (2014) that stars on the expanding side of the spiral always exhibit larger breathing motions when compared with those for the compressing side of the spiral. Debattista (2014) attributed this behaviour to the more abrupt density variation as stars leave the spirals compared with when they enter them. (The test particle and semi-analytic calculations of Faure et al. (2014) did not exhibit this behaviour because the spiral perturbation they used had a purely = 2 multiplicity.) Furthermore, the slope of the young stellar population is larger than the slope for the older stellar population, thereby demonstrating that indeed the breathing motion is stronger in the young stellar population than in the old stellar population.

COMPARISON WITH VERTICAL MOTIONS IN THE MILKY WAY
We now search for these signatures of spiral-driven breathing motions in the Gaia DR2 data as a function of stellar ages, and compare with the results obtained from the star-forming model.

Sample selection from the Sanders & Das dataset
We use the publicly available sample of stellar ages of Sanders & Das (2018) (hereafter, SD18), which provides a catalogue of distance, mass, and age for approximately 3 million stars in the Gaia DR2. This sample allows us to study the vertical motions for different age populations. Gaia DR2 suffers from distance biases from the parallax measurements (for detailed discussion, see e.g. Each stellar entry in the SD18 dataset is characterised by three flags that, depending on the nature of the analysis, can aid in filtering out problematic data. The entry characterises all issues related to the quality of the input data (spectroscopy, astrometry, photometry) from the above listed surveys and of the output data from the SD18 Bayesian frameworks (mass, age). If there are no issues in the input or output data, then = 0. The flag specifies if stellar entries are duplicates due to a crossmatch within a single survey or between multiple surveys. Within a single survey, the duplicate with the smallest vertical velocity uncertainty is kept, while duplicates between surveys are kept based on the entry and a survey hierarchy set in SD18 (APOGEE, GALAH, GES, RAVE-ON, RAVE, LAMOST and SEGUE). For example, if APOGEE and RAVE contain a duplicate entry with proper input and output data ( = 0) then, regardless of the relative errors, the APOGEE entry is kept and labelled as = 0. Lastly, the flag in SD18 encompasses the two prior flags and is 1 when = 0, = 0 and the star has a valid cross-match in Gaia DR2. For this work, we choose a sub-sample of stars from the SD18 dataset by applying a number of quality cuts. In particular, we choose only those stars for which distance ( ) is less than 10 kpc, 1/ < 10 kpc, = 1, and the ratio of parallax to error in parallax, / > 3. Since the parallax/distance bias is more significant for very young stars (for details see Sanders & Das 2018), we discard all the stars with ages less than 1 Gyr. Our selected subsample has a total of 3,147,069 stars.

Comparison of SD18 dataset with a bias-controlled dataset
Before we proceed to the calculation of breathing motions for the Milky Way using our selected sub-sample from the SD18 dataset, first we compare our sub-sample with a bias-controlled dataset for the Milky Way. We choose the Schönrich et al. (2019) dataset (hereafter, SME19) who corrected for the bias in the Gaia parallax, and derived the distances of all stars in the RVS (Radial Velocity Spectrograph) sample of the Gaia DR2 via a Bayesian framework (for details, see Schönrich et al. 2019).
Since the SME19 dataset does not include stellar ages, we do not include any age cuts on the SD18 dataset for this comparison. In both datasets we select stars with distance < 10 kpc and with the ratio of parallax to parallax error, / > 3. Furthermore, to reduce the contamination from both very young stars and halo stars (at larger heights), we further constrain the sample based on the radial action ( ) of the individual stars. The SD18 dataset includes the radial action of individual stars; for the SME19 dataset we calculate these ourselves. To calculate , we use the Stäckel fudge method (Binney 2012;Sanders & Binney 2016) from the AGAMA software library (Vasiliev 2019) in the potential of McMillan (2017). The Stäckel fudge function in AGAMA takes the Galactocentric 6D coordinates of the SME19 dataset and produces the full action-angle coordinates, including the values. We then compute the cumulative distribution function (CDF) of the radial action, ( ), of the SME19 dataset. The very young stars are, in general, on almost circular orbits and thus should have small values of . On the other hand, the halo stars at larger heights should have large values of . Therefore, we choose stars from the SME19 dataset, while rejecting the top and bottom 10 percent of stars in ( ). We denote this rejection process as Δ ( ) = 0.1. To make a uniform comparison, we apply the same radial action cut (Δ ( ) = 0.1) on the SD18 dataset using the same values of upper and lower as in the SME19 dataset. Fig. 11 shows the corresponding distributions of stars in the ( , ) plane and in the ( , ) plane, as well the distribution of mean vertical velocity, , in the ( , ) plane for both the SD18 and SME19 datasets. We split our selected samples of stars from both the datasets, into positive and negative Galactocentric azimuthal angles, and recalculate the distribution of , which is shown in Fig. 12. The distribution of the mean vertical velocity ( ) in the ( , ) plane, calculated for both the SD18 and SME19 datasets appear to match quite well. Further, we checked that a more conservative cut on , e.g. Δ ( ) = 0.15, does not alter the main findings presented here and in the subsequent sections.
We calculate the distributions of the bending velocity ( bend ) in kpc (using Eq. 5) for both the SD18 and SME19 datasets. The resulting bending velocities calculated from the two datasets are broadly in agreement; these results are not shown here but are evident in Fig. 11. We also compare the breathing motions, breath , between the two datasets. Because the breathing motions are significantly smaller than the bending motions, we start looking for the signature of the breathing motions in the Solar Neighbourhood in the following way. First, we calculate the distribution of the best-fit slope in the ( , ) plane around the Solar Neighbourhood, by fitting a straight line to the variation of with , in a similar way as done for the model. The corresponding distribution of the best-fit slope, and the associated error ( (slope)) in the ( , ) plane, for both the SD18 and the SME19 datasets are shown in Fig. 13 (see top and middle panels). The bottom panels of Fig. 13 show the bins for which a quality-cut of / ( ) > 3 is met simultaneously for these datasets. We find that the location where we get the most consistent signature of the breath is ∈ [7.6, 8.1] kpc and ∈ [0.9, 1.4] kpc. In this bin, the best-fit slope values, obtained from the two datasets agree with each other (within their error-bars). Moreover, we find that this bin shows the largest value for the best-fit slope, i. e., the signature of the largest breath . The presence of this 'consistent' signature of breath is not seen for any other bins shown in the bottom panels of Fig. 13. Hence, we will only consider this spatial bin for the subsequent study of the age-dependence of the breath in the Solar Neighbourhood. The corresponding slope of versus at the same spatial location using the two datasets are shown in Fig. 14. We note that there is an offset between the values calculated from both the datasets, especially near the disc mid-plane; this could possibly be due to the different zero-points used in the datasets. The values of the slopes obtained from the two datasets match quite well within their error-bars. In both cases, the / ( ) is found   Figure 13. Comparison of the breathing motions in the SD18 and SME19 datasets : distributions of the best-fit slope (top panels) and the associated errors (middle panels) in the ( , ) plane are shown for the SD18 (left panels) and the SME19 (right panels) datasets. Bottom panels show the same distribution of the best-fit slope, but only for those bins where a simultaneous quality cut of / ( ) > 3, for both the datasets, is met separately.
The magenta rectangles in the bottom panels denote the chosen spatial bin where we study the age variation of the breathing motions (see text for details). The Sun's position is also indicated in the bottom panels.
to be greater than 3 so that in both datasets the detection of breath is significant.

Vertical motions for different ages in the Milky Way
We now study the non-zero mean vertical motion, and the associated breathing motions, as a function of the stellar ages in the same volume as where we found the largest breathing motions. We can do this only using the SD18 dataset. Fig. 15 presents how the selected sample of stars from the SD18 dataset is divided into three equal-sized populations by cutting the cumulative distribution function, (age), of the stellar ages. The resulting sub-populations are: young (0 ≤ (age) < 1/3), intermediate (1/3 ≤ (age) < 2/3), and old (2/3 ≤ (age) < 1) stellar populations, or, in terms of the stellar ages, [1, 4.85] Gyr (young), side of the mid-plane are seen to move in opposite directions, i. e., with a non-zero breath at larger heights, while at moderate heights all populations exhibit bending motions (i. e., stars on either side of the mid-plane are moving in the same direction).
To measure the variation of the breathing velocity as a function of stellar age, we consider the same spatial location as before, i.e. ∈ [7.6, 8.1] kpc and ∈ [0.9, 1.4] kpc, where a strong (and consistent) breathing motion is present without an age cut. The vertical distributions and slopes of the three populations in this spatial bin are shown in Fig 18. There is an offset in the zero point among the different age groups, similar to what is seen in the model. We find a similar explanation as in the model, i. e., a difference in bending motion amplitudes in different age groups is giving rise to this offset. The presence of a non-zero slope in all three populations indicates that all of them take part in the breathing motion. The values of the best-fit slope vary with stellar age: the (absolute) value of the best-fit slope is largest for the intermediate stars (2.27 ± 0.25 km s −1 kpc −1 ) and smallest for the old stars (1.02 ± 0.19 km s −1 kpc −1 ) whereas the value of the best-fit slope for the young star falls in between (1.61 ± 0.06 km s −1 kpc −1 ). We check whether the stars at the larger distance from the disc mid-plane bias the fits by ignoring the two points at the largest height (| | = 1.5 kpc). The resulting bestfit slopes (−1.6 ± 0.1 km s −1 kpc −1 for young, −2.26 ± 0.6 km s −1 kpc −1 for intermediate, and −1.5 ± 0.04 km s −1 kpc −1 for old stellar populations) match, within the error-bars, with the slopes obtained including the largest heights, indicating that the large | | points are not biasing the fits.
To test the robustness of the best-fit slope to systematic biases in the Gaia distances, we lower the distances of all the stars in our selected sub-sample by 10 percent to simulate a parallax bias. After repeating the analysis of fitting a straight line to the variation of with height, the resulting change in the best-fit slope is negligible (less than 5 per cent) for all three age populations. This emphasises that the breathing motions are not sensitive to systematic errors of this type and accentuating the robustness of the detection of breathing motions in the Solar Neighbourhood.
So far, we have estimated the statistical errors on and also checked for the systematic errors on the best-fit slope. However, the uncertainties associated with the parallax and the proper motions, the quantities used to derive vertical motion of stars, also result in uncertainties in the calculation of the and the corresponding best-fit slope. To estimate this error on the best-fit slope, we perform a Monte Carlo analysis. First, we calculate the slope and the associated error of the best-fit straight line using the 6-D position-velocity of stars in our full selected sub-sample. Since there is a robust signal of breathing motions only in our selected spatial location in the Solar Neighbourhood, we restrict our Monte Carlo analysis for that spatial location only. Then, using the errors in parallax and the proper motion, we add random Gaussian errors in the parallax and the proper motion of each star in our sub-sample, and recalculate the 6-D position-velocity of stars in our selected sub-sample using G (Bovy 2015). Once we recalculate the 6-D position-velocity of stars, we repeat the exercise of fitting a straight-line to the variation of with height, , for all three stellar populations with different ages. This process is repeated 1000 times. The resulting distributions of best-fit slopes, calculated for that particular spatial bin, for all three stellar populations are shown in Fig. 19. The estimated errors in the best-fit slope measurements, obtained from the Monte Carlo analysis are ∼ 14 per cent for the young population, ∼ 16 per cent for the intermediate population, and ∼ 17 per cent for the old population.

SIGNATURES OF SPIRAL DRIVEN BREATHING MOTIONS
Using a star-forming simulation we have shown that a spiral density wave drives vertical breathing motions. The increasing amplitude of the breathing velocity with height above the mid-plane, previously reported by Debattista (2014), is consistent with the trend seen in the Milky Way (Gaia Collaboration et al. 2018a). This behaviour is an important clue to the source of breathing motions: it must be something local to affect stars near the mid-plane differently from those further up. The reason why spirals have this effect is easy to understand: on vertical scales small compared to the disc scale-length, the vertical gravitational field can be approximated, using Gauss's law, as a uniform field from an infinite sheet of mass corresponding to the mass enclosed within a cylinder of the same height and small radius. Both the simulation of Debattista (2014) and the one considered here show that the spiral density wave extends to large heights (see Fig. 3).
If the density of the disc is decomposed into an axisymmetric part and a spiral perturbation part: then the vertical velocities are in equilibrium with the axisymmetric part of the potential and 0 ( , ) gives rise to breath = 0. Since the azimuthally varying part, , is due to spirals, it extends to large height. The enclosed mass therefore increases with height as the vertical integral of ( , , ). Therefore the gravitational field increases with distance from the mid-plane, and breath increases with height, as we see in the simulations, and in the Milky Way (compare  2018b)). It is noteworthy that this vertical increase in breath is absent in the model of Faure et al. (2014), who used an analytic spiral potential with a very short scale height of 100 pc, much shorter than that of the unperturbed density distribution. As a result, the gravitational field reaches nearly constant value quite close to the mid-plane in their model. Consequently, breath shows a near-linear increase from the disc mid-plane, and after ∼ 0.7 0 , it starts to decline where 0 is the scale-height of the spiral perturbation. The vertical variation of breath in the Gaia data therefore indicates that the spiral structure in the Milky Way is vertically extended. The vertical increase of breath is different from the nearly constant bend as a function of height (compare both panels of Fig. C.6 in Gaia Collaboration et al. (2018b)). This indicates that the Milky Way's bending motions are probably not excited by spirals.
Moreover we examined the dependence of the breathing velocity on stellar age. The star-forming model exhibits a decreasing breath with increasing stellar age (see Fig. 10) indicating that increasing random motions decrease the vertical response of the stars. Fig. 18 shows that, using the SD18 dataset, the intermediate stellar population indeed yields a larger (absolute) value of the best-fit slope in the Solar Neighbourhood compared with the old stellar population, while the young stellar population has a slope that falls in between; however the difference between the young and intermediate age populations (1.8 ) is not statistically significant in this dataset. In other words, the breathing motion is stronger for the younger stellar populations (i.e., young and intermediate populations considered together) than that for the old stellar population in the Solar Neighbourhood.
For these reasons we therefore argue that the breathing motions seen in the Solar Neighbourhood are very likely to be driven by the spiral density wave structure. Widrow et al. (2014) proposed a scenario where a satellite galaxy, while plunging into the disc, can generate both bending and breathing motions, depending on the value of the satellite's velocity in the direction perpendicular to the disc. Such interactions also excite a strong spiral response within the disc which are probably partly responsible for the breathing motions seen in the simulations of Widrow et al. (2014). The vertical variation of the amplitude of the breathing modes excited by satellites is also unclear.

FUTURE PROSPECTS AND SUMMARY
Our requirement of stellar ages limits the number of stars available to us by almost a factor of two when compared with the Gaia (which contains 6,376,803 stars; Gaia Collaboration et al. 2018b). In addition, the azimuthal coverage in Gaia DR2 is not large (Δ ∼ 30 • ). This, in turn, prevents us from exploring the change from expanding to compressing breathing motions as a function of azimuthal angle, which would help constrain the extent of the spiral and its corotation radius (Faure et al. 2014;Debattista 2014). Gaia DR3 holds the potential to permit us to investigate the variation of breathing motions as a function of azimuth, in order to constrain the Milky Way's spiral structure better.
Our main findings are: • The spirals in the self-consistent star-forming model drive coherent vertical breathing motions, with amplitude which increases with the height from the disc mid-plane. We argue that this vertical variation is a consequence of the vertically extended spiral structure. Since Gaia DR2 also exhibits this increasing breathing amplitude with height, it implies that the Milky Way's spirals are also vertically extended.
• Prominent breathing motions are present in stellar populations of all ages in the model, with the strongest motions in the young stellar population, and weakest in the old one.
• Using Gaia DR2 with complementary age data from Sanders & Das (2018), we showed that the location near the Solar Neighbourhood with the smallest biases shows significant, strong breathing motions for all ages. In addition, the amplitude of the breathing motions varies with stellar age, with the oldest population displaying the weakest breathing amplitude (when compared with the intermediate and the young population, at 2.8 − 3.1 ). The difference between the young and intermediate age populations is smaller but less statistically significant (1.8 ). This behaviour is similar to the trend shown in the simulation model and is consistent with spirals being the driving mechanism of breathing modes.
The resemblance in the height and age variation of the amplitude of breathing velocity, in the star-forming model and in the location in the Milky Way showing prominent breathing motions from the Gaia DR2 dataset indicates that the observed breathing motions in the Milky Way are likely to be excited by spiral density waves. Due to the limited coverage and uncertainties in the Gaia parallax, we have restricted to the Solar Neighbourhood our search for breathing motions, but with the upcoming Gaia DR3, there will be a larger number of stars allowing us to improve the statistics in the Solar Neighbourhood.