© ESO 2019

1. Introduction

Warped discs still represent a theoretical challenge. From surveys of edge-on galaxies, it is clear that ∼50−70% of spiral disc galaxies present stellar warped discs (e.g. Sanchez-Saavedra et al. 1990), suggesting that they are long-lived or repeatedly generated. Our Galaxy also presents a warped disc, first detected with 21 cm observations of the HI gas (Burke 1957; Westerhout 1957; Oort et al. 1958; Levine et al. 2006, among others), later in dust and stars using the Two Micron All Sky Survey infrared data (Freudenreich et al. 1994; Drimmel & Spergel 2001; López-Corredoira et al. 2002; Reylé et al. 2009; Amôres et al. 2017, among others), and more recently using Cepheids (Skowron et al. 2018; Chen et al. 2019). Most of these studies dedicate their efforts to determining the morphology of the warped disc, that is, the Galactocentric radius at which the disc starts bending, the phase angle of the line-of-nodes, its maximum amplitude, and its possible dependence upon the tracer.

Very few studies have focused on analysing the warp kinematically, in the sense of finding the effect of a warped disc on the kinematics of the stars. The first kinematic analyses were conducted using HIPPARCOS proper motions. Dehnen (1998) selected a set of kinematically unbiased main sequence and giants stars and found evidence of the stellar warp when plotting the vertical velocity as a function of the tangential velocity. Drimmel et al. (2000) plotted the vertical velocity of the stars as a function of the Galactocentric radius for a sample of OB stars, concluding that the kinematics of the stars towards the outer disc was inconsistent with the expectations from a long-lived non-precessing warp. Some years later, studies by Seabroke & Gilmore (2007) as well as Bobylev (2010, 2013) concluded the (then) available proper motion surveys do not allow complete studies of the Galactic warp. Using PPMXL proper motions and selecting disc Red Clump stars using 2MASS photometry, López-Corredoira et al. (2014) found kinematic evidence of the stellar warp, by plotting the vertical velocity as a function of the Galactocentric azimuth, concluding that their results cannot be reproduced by a population in statistical equilibrium. The quality of the Gaia Data Release 1 and the Tycho-Gaia Astrometric Solution (TGAS, Michalik et al. 2015; Brown et al. 2016) proper motions allowed Schönrich & Dehnen (2018) to review the study started by Dehnen (1998). Following a similar strategy, these authors find evidence of a kinematic signature of the Galactic stellar warp in the TGAS data in a cone towards the centre and anti-centre directions and estimate the onset of the warp at a guiding radius inside the Solar circle, R_g ≲ 7 kpc, in agreement with the previous work by Drimmel & Spergel (2001). Carrillo et al. (2018) and Schönrich & Dehnen (2018) also point out that the complexity of the data cannot be explained by a simple warp pattern, but also wave-like patterns, bending and breathing modes are reflected in the kinematic structure. The advent of the Gaia Data Release 2 (Gaia DR2), and the availability of proper motions for 1.3 billion sources (Gaia Collaboration 2018a), has allowed the study to be expanded to larger volumes and fainter sources. Gaia Collaboration (2018b) has already shown the complexity of the vertical motion and Poggio et al. (2018) review the data selection of OB-type and Giants stars performed in Gaia Collaboration (2018b) and show the kinematic evidence of the Galactic warp in vertical velocity using 2D maps.

The main challenge in our work is going one step forward. We want to use the spatial and kinematic data, together with models, to constrain the morphology of the Galactic warp. We assume this challenge by defining three Galactic warp models and using two characterisation methods developed by the authors specifically for this purpose. The complexity of the data leads us to consider not only simple symmetric tilted ring models, but also more complex asymmetric models, such as lopsided tilted rings or lopsided S-shaped warp models. We also use two methods to infer the structural parameters of the warp given the observables (positions and proper motions). The first is the LonKin method (Abedi et al. 2015), which provides information on the level of asymmetry of the warp and shape when applied to different warp models. The second is nGC3, a method from the Great Circle Cell Counts (GC3) family (Johnston et al. 1996; Mateu et al. 2011), which searches for over-densities of stars along great circle cells. The latter has already been used in simulations by Abedi et al. (2014) to assess the Gaia capabilities to derive the tilt angle, or amplitude of the warp, and the twist angle of the line-of-nodes. Here, we apply both methods to two different populations from Gaia DR2, namely OB-stars and red giant branch (RGB) stars (similarly but not equally selected as in Poggio et al. 2018). We define a strategy that allows us to select stars in both samples up to magnitude G = 20 mag in order to have enough statistics up to R ≤ 16 kpc and to disentangle the kinematic signature of the warp. Two tracer populations are necessary to study the dependency of morphology and kinematics on the age of the tracer, as recently suggested by Amôres et al. (2017). We study the characteristics of the Galactic warp in terms of spatial and kinematic data of the two tracer populations, with intrinsically different ages, in order to provide a first step towards understanding the origin of the Galactic warp. Additionally, the application of the methods to samples of test particles evolved with the different models specified above allows us to give an interpretation of the results in terms of shape, onset radius, and tilt angle.

The paper is organised as follows. In Sect. 2 we give a brief description of the characterisation methods (Sect. 2.1) used, we discuss the need to use an adequate distance estimator, and we give the final choice (Sect. 2.2); we also describe the Gaia mock catalogues used to compare with the data (Sect. 2.3). In Sect. 3, we describe the selection of the working samples from Gaia DR2 data and give an initial characterisation. In Sect. 4, we study the spatial density distribution of both samples (Sect. 4.1) and give the first spatial characteristics of the Galactic warp (Sect. 4.2). We then continue by studying the kinematic signature of the warp in Sect. 5. We first focus on the kinematics in 2D projection maps (Sect. 5.1) and in Sects. 5.2 and 5.3 we apply the LonKin and nGC3 methods to the two populations extracted from the Gaia DR2 data, respectively. In Sect. 6 we combine the spatial and kinematic distributions and we compare them with the models and data from the literature. Finally, in Sect. 7 we present our conclusions. This work is complemented with four Appendices: in Appendix A we give further details on the sample selection; in Appendix B we provide a full description of the warp models used; in Appendix C we detail the position and velocity transformations applied to a flat disc in order to warp it according to each of the three models; and in Appendix D we discuss the information the two methods provide when we apply them to synthetic data consisting of three different sets of particles, simulated with increasing complexity and reality.

2. Methods and data treatment

In order to analyse the warp signature in the disc kinematics, we use methods previously developed by the authors (Abedi et al. 2014), namely the LonKin method, and the GC3 method and its variations. Both methods are designed to provide structural and kinematic information of the warp from the use of positions, distances, and proper motions. Firstly, in Sect. 2.1, the methods are described and specify the requirements needed to apply them. Secondly, in Sect. 2.2, we use mock catalogues to study and define the appropriate choice for the distance estimator used throughout this paper.

2.1. Methods for warp detection and characterisation

The LonKin method looks for the signature of a possible warped disc in a plot of proper motion in Galactic latitude, μ_b, as a function of the heliocentric Galactic longitude l. The method as shown in this paper was developed in Abedi et al. (2015), and its main advantage is that it directly uses the Gaia observables. The method is a variation of those used in Drimmel et al. (2000) and López-Corredoira et al. (2014), which plot the vertical velocity, W, as a function of the Galactocentric radius and of the Galactocentric azimuth, respectively. Therefore, the main advantage of LonKin with respect to these is that it works in the space of the observables, using proper motions instead of the vertical velocity, W. In the vertical axis we plot the proper motion in latitude corrected by the reflection of the solar motion, that is, with respect to the Local Standard of Rest, μ_b, LSR. In the horizontal axis, we plot the heliocentric Galactic longitude, l, in segments of 20°. For each segment, we compute the median of the μ_b, LSR values and the lower (15.85-percentile) and upper (84.15-percentile) 1σ uncertainties on the estimation of the median as defined in Appendix A of Gaia Collaboration (2018b), to which we add in quadrature the uncertainty given by the solar velocity. Since the Galactic warp is a Galactocentric feature, we need to split the sample into cylindrical Galactocentric radial bins, for which we need to assume a Galactic solar radius and propagate the errors to the Galactocentric frame to obtain the cylindrical Galactocentric radius R. If the Galactic disc is flat, the median μ_b, LSR should be constant and equal to zero, but if it is warped, a particular variation will be introduced as a function of l. If the disc is symmetrically warped and the line-of-nodes is aligned with the Sun–Galactic centre line, the LonKin method predicts a maximum in μ_b, LSR in the anti-centre direction. If on the other hand the disc is asymmetric, that is, it is Lopsided, and the line-of-nodes, still on the Sun–Galactic centre line, the maximum in μ_b, LSR is no longer in the anti-centre direction, but is shifted towards longitudes coinciding with the maximum warp amplitude (see detailed discussion in Appendix D).

The GC3 methods (Johnston et al. 1996; Mateu et al. 2011) comprise different ways of searching for overdensities in great circle cells in the sky. In the most general version (mGC3) including full kinematic information – introduced in Mateu et al. (2011) – the method sweeps over the sky counting how many stars have position and velocities lying in a great circle within a given tolerance, each great circle being defined uniquely by its normal vector or pole. The all-sky sweep over all possible great circle cells results in a pole count map (PCM), a plot of the number of stars associated to each possible pole, and thus great circle cell, in the celestial sphere. To apply this method for the particular case of characterising the warp, PCMs are made for different Galactocentric radial bins and the star counts are made in great circles defined in a Galactocentric reference frame. If the Galactic disc is flat, the peak in stellar density should be located in the north Galactic pole of the PCM. If the disc is not flat, the peak of over-density should move in the PCM providing information on the tilt angle of the warp, as well as the azimuth (twist) of the line-of-nodes as a function of radius (see Abedi et al. 2014, for detailed examples). If the warp is lopsided, the signature in the PCM is not a single peak, but has a shape that depends on the warp model (see detailed discussion in Appendix D). In this work in particular, we use the nGC3 method of the family, which uses 3D position information and proper motions, without requiring line-of-sight velocities. As we have shown in Abedi et al. (2014), the use of the full velocity information in mGC3 severely limits the sample spatial coverage and does not produce a significant improvement in the results.

2.2. Selection of the distance estimator

The Lonkin and nGC3 methods both start off by binning the sample in Galactocentric distance. This is the most error-prone step in both methods (see Abedi et al. 2014), and therefore in order to reach as far as possible into the disc with the smallest possible distance bias it is crucial to have a well-behaved distance estimator. Thus, our analysis here is focused in estimating the Galactocentric distance bias for stars binned in consecutive Galactocentric rings, for different distance indicators under DR2 parallax error prescriptions (Lindegren et al. 2018). This will also allow us to estimate the maximum Galactocentric radius we can reach without introducing a significant bias.

We use the warped disc test particle simulation for RC stars up to G = 20 mag from Abedi et al. (2014) and explore as distance indicators: the distance computed as the reciprocal of the parallax (excluding stars with negative parallaxes); the Bayesian indicator with an exponentially decreasing space density (EDSD) prior proposed by Bailer-Jones (2015) and Astraatmadja & Bailer-Jones (2016) for values of the scale length parameter L = 1.35, 2.0 and 2.5 kpc; and the EDSD prior with a variable scale length L(l, b) dependent upon the line of sight (Bailer-Jones et al. 2018, hereafter BJ18). We use the posterior mode as a point estimator for the distance from here on, as in BJ18.

Figure 1 shows two summary statistics (median and mode) for the bias in the inferred Galactocentric radius r_gal, that is, the deviation of the median or mode of the true r_gal distributions for stars in a given inferred r_gal bin. In this plot we show the results obtained using stars with observed parallax error smaller than 50%, that is, f_obs = |Δϖ/ϖ| <  0.5. We checked that this f_obs threshold allows us to exclude stars with arbitrarily large fractional parallax errors (both positive and negative), while keeping the Galactocentric distance bias reasonably low. More details are given in Appendix A.1. We are aware that a change in the selected f_obs threshold would affect the distribution of true distances of the stars selected in the sample and therefore the optimal prior scale length.

Fig. 1. Median (top) and mode (bottom) bias of the true rgal distances of stars selected in bins of inferred rgal, for Gaia DR2 errors. Results for the reciprocal of the parallax are shown with the (yellow) solid+circle line. The prior scale lengths used with the Bayesian estimators are indicated with different colours as shown in the inset keys.

The figure shows the behaviour of the distributions is fairly stable up to r_gal ∼ 12−13 kpc, beyond which there is an abrupt change. For all estimators, the (median and mode) bias first increases up to that distance and then systematically plummets to large negative values. For any given estimator, it also shows that the mode of the true distance distribution in a given bin is less biased and more stable than the median, at all distances.

For the median r_gal, the best performance up to the threshold radius (∼13 kpc for RC stars) is obtained with the reciprocal of the parallax and the Bayesian EDSD prior with L = 2.5 kpc. At larger distances, the L = 2 kpc estimator performs slightly better than L = 2.5 kpc and much better than the reciprocal of the parallax. In general, the BJ18 L(l, b) model results most resemble those for the shortest scale length L = 1.35 kpc, for both the mode and median. This is consistent with the fact that values of L shorter than 1.35 kpc are the most common in the BJ18 model (see their Fig. 1).

For the mode r_gal, the best performance is consistently obtained with the L = 2 kpc estimator, which is effectively unbiased up to ∼14.5 kpc and showing the smallest bias (∼0.5 kpc) at 15.5 kpc. From r_gal >  13. kpc, the short scale estimators L = 1.35 kpc and L (BJ18) show a larger bias than for L = 2 kpc (∼0.3 kpc), half the corresponding bias expected for the median (∼0.6 kpc). Therefore, for the 0.5 kpc bins used in our analysis, the results for these two short scale length priors should not be overly different to those of the L = 2 kpc estimator when using the mode.

The threshold radius of ∼12−13 kpc that marks the change in behaviour for the bias depends upon the selected tracer, as this choice will set a particular dependence of the parallax errors with distance, via the intrinsic magnitude, the colour, and the distance modulus. Compared to the RC stars used in the tests presented here, we expect the threshold radius to be slightly larger for OB stars, around 16 kpc, since these stars are intrinsically brighter on average, based on similar tests conducted in Abedi et al. (2014) using the reciprocal of the parallax as an estimator. There will also be an additional contribution to the r_gal bias due to contaminant stars, which will have a different error distribution. However, this contribution is small: less than 10% of the stars in either of our samples (see discussion in Appendix A.2) and largely due to dwarf stars. They affect mostly the radial bins at r_gal <  10 kpc and contribute at most ∼0.15 kpc, in addition to the bias shown in Fig. 1, at any given radius.

As discussed in Luri et al. (2018) the best estimator depends upon the particular choice of the sample and its specific parallax error distribution. Our results show that there is no single estimator, of those considered here, that simultaneously outperforms all the others at all distance ranges and with all statistics. Overall, we find that for Gaia DR2 errors and f_obs <  0.5, the mode as the estimator with a prior scale length L = 2 kpc shows the best performance.

2.3. The use of Gaia mock catalogues

In this work we make use of Gaia mock catalogues in order to test the capabilities of the methods described above and the possible effects of Gaia selection function and astrometric errors in the characterisation of the Galactic warp. Details about the warp models, the coordinate transformations, and the generation of mock catalogues are given in Appendices B–D.1, respectively.

We generate three Gaia mock catalogues, one to test the null hypothesis, i.e. a catalogue with no warp, and a further two with imposed warp models, namely the Sine Lopsided and the S Lopsided model. We first generate a set of disc Red Clump initial conditions in a flat disc as in Abedi et al. (2014) and Romero-Gómez et al. (2015). This set of particles forms our null hypothesis. We then generate a test particle simulation by integrating the initial conditions using the same integration strategy as in Abedi et al. (2014) using the Sine Lopsided, a warp model that allows the test particle simulation to reach statistical equilibrium. On the contrary, we cannot ensure statistical equilibrium for the S Lopsided models, that is, we cannot guarantee that positions and velocities of a test particle simulation using the S Lopsided model will inherit those established by the model. We therefore show the results of the random realization of particles, that is, we apply the velocity transformation described in Appendix C to the initial flat conditions. The free parameters defining the warp disc are taken from Amôres et al. (2017) for a mean age population of 2.5 Gyr.

Finally, we generate a Gaia DR2 mock catalogue up to magnitude G = 20 with the extinction model from Drimmel et al. (2003) and using the prescription of the astrometric, photometric, and spectroscopic formal errors for Data Release 2¹. We simulate errors for a mission time of 22 months and for the bright stars (G <  13) we include a multiplicative factor of 3.6 to the error in parallax to match the distribution of uncertainties as a function of the G magnitude observed in Gaia DR2 data. The final mock catalogues have more than 2 million Red Clump star particles.

3. Selection of the two working samples in Gaia DR2

From the Gaia DR2 catalogue, we are interested in selecting two samples characterised by being intrinsically bright, in order to reach the outermost parts of the disc, and with different age ranges, in order to assess whether the structural and kinematic properties of the warp depend on the age of the population. In the HR diagram, we select a young population formed by mainly upper main sequence stars or OB-type stars (hereafter, OB sample) and an older population composed of all the stars in the RGB (hereafter, RGB sample). From Gaia DR2² we select stars up to magnitude G = 20, with an available parallax measurement and with f_obs <  0.5 (as discussed in the previous section), that is, an absolute value of the relative error in parallax of less than 50% (i.e., we keep stars with negative parallaxes). This first cut reduces the Gaia DR2 sample to 383 510 799 sources.

For these stars, we compute the distance using a Bayesian estimator with an exponentially decreasing space density prior with scale-length L = 2 kpc based on our analysis of the distance bias for different estimators (see Sect. 2.2). Once we have derived the distances, and in order to reduce the computational cost of the process, we make a second cut which consists in removing cool main sequence stars from the sample. To do that, we select stars with , following the extinction line, where is the absolute G magnitude of the star uncorrected for extinction; is given by ; and (G_BP − G_RP) is the observed colour. More details are given in Appendix A.2. This second cut reduces the sample to 86 814 618 sources. Since we want to characterise the Galactic warp, we only keep stars with Galactocentric spherical distance r_gal >  7 kpc and Galactocentric cylindrical distance R_gal <  16 kpc³. This second cut reduces the sample to 45 349 864 stars. For these stars, we compute the V-band absorption using the Drimmel extinction model (Drimmel et al. 2003) with re-scaling factors. Using A_V and (G_BP − G_RP), we compute A_G and the reddening E(G_BP − G_RP) according to a polynomial fit (J. M. Carrasco, priv. comm.). We refer to this sample as the “All-Stars” sample. Finally, the two tracers are selected as in Gaia Collaboration (2018b):

1. Young population, OB sample:

(1)

1. RGB sample:

(2)

where and (G_BP − G_RP)₀ = (G_BP − G_RP) − E(G_BP − G_RP). The OB sample consists of 1 860 651 stars and the RGB sample has 18 008 025 stars. This RGB sample is a mixture of stars with ages in the range 3−7 Gyr, while the OB sample is typically formed by young stars with ages ≤1 Gyr (e.g. Robin et al. 2012). Using the Gaia Universe Model Snapshot (Robin et al. 2012) with Gaia DR2 simulated errors, we estimated the contamination in the OB and RGB samples to be 9% and 8%, respectively (see Appendix A.2).

Figure 2 shows the resulting distributions in G magnitude and heliocentric distance, d, for both working samples. We emphasise here that both samples reach up to magnitude 20 and they extend well beyond the solar neighbourhood, having 1.1 × 10⁶ OB type stars and 13.6 × 10⁶ RGB stars up to d = 5 kpc.

Fig. 2. Distributions in magnitude G (top panel) and heliocentric distance (bottom panel) of the two working samples.

4. The spatial distribution of the OB and RGB samples

In this section, we analyse the spatial distribution of both working samples in order to detect the signature of the warp. We apply a coordinate transformation to the Galactocentric cartesian frame and we study the structural characteristics of both samples. Throughout this work, we use a right-handed Galactocentric reference frame. Its X-axis points along the Sun-Galactic centre line, away from the Sun. The Y-axis is orthogonal to the X-axis and on the Galactic plane in the direction of Galactic rotation. The Z-axis points towards the North Galactic pole. We also use cylindrical coordinates (R, θ, Z) with R being the Galatocentric distance and the azimuthal angle θ measured on the Galactic plane starting from the positive X-axis and towards the positive Y-axis. To perform these coordinate transformations from the Gaia DR2 data, we adopt a distance of the Sun to the Galactic centre of R_⊙ = 8.34 kpc (Reid et al. 2014) and the height of the Sun with respect to the Galactic plane of Z_⊙ = 0.027 kpc (Chen et al. 2001). We also adopt the circular velocity at the solar radius from Reid et al. (2014) of V_c(R_⊙) = 240 km s⁻¹⁴. We adopt the peculiar velocity of the Sun with respect to the Local Standard of Rest of km s⁻¹ (Schönrich et al. 2010).

4.1. Peculiar features in the stellar density distribution

In Fig. 3 we plot the stellar density map of the two samples in the (X,Y) cartesian Galactocentric projections. In general, both populations show the radial features characteristic of the high extinction lines (under-densities), which point to the presence of dust clouds in the direction of the line-of-sight, and of the error in distance of a given object (shown as elongated over-densities, also known as the finger-of-god effect). The effect of dust clouds is visible in both populations (OB and RGB samples), the most prominent clouds being those we find in the directions: 75° < l <  85°, 130° < l <  150°, 180° < l <  190°, and 250° < l <  270°, well in agreement with the Drimmel extinction map (Drimmel et al. 2003).

Fig. 3. Two dimensional stellar density map of the two working samples, namely the OB sample (left column) and the RGB sample (right column). The sizes of the bins are 0.32 kpc × 0.64 kpc for the OB sample and 0.16 kpc × 0.32 kpc for the RGB sample. We mark the Sun position with a star and the black dotted lines show circles at different Galactocentric radii, R = 8, 10, 12, 14, and 16 kpc. The colour bar is in log-scale and is the same for both panels. The Galaxy rotates clockwise.

Focusing on the stellar density map of the OB sample (see left panel of Fig. 3), the over-densities we observe clearly trace some known star forming regions. Zooming in at heliocentric distances closer than 1 kpc, we can see the Orion, Vela, and the Cygnus star forming regions, well described in Zari et al. (2018). At heliocentric distances in the range 1 <  d <  3 kpc, we see the Cassiopeia region (120° < l <  130°), the Outer Spur (230° < l <  250°), and the Carina star forming region (280° < l <  300°). In the latter is where we see the largest concentration of OB stars. This over-density is well traced by the early-type stars, including HII regions (e.g. Molina-Lera et al. 2016) or the classical Cepheids (Skowron et al. 2018), and it is often associated to the Carina–Sagittarius spiral arm of gas and dust (e.g. Martos et al. 2004). We observe that this Carina region over-density contributes to the Galactocentric rings 8 <  R <  10 kpc in the study of the Galactic warp.

In the right panel of Fig. 3, we plot the stellar density map corresponding to the RGB sample. Although it appears more homogeneous than the map for the OB stars, it presents a clear asymmetry in stellar density with respect to the centre–anticentre line. We show this feature in detail in Fig. 4, where we separate the RGB stars that are located above (left panel) and below (right panel) the Galactic plane. Although some of the under-densities visible in both panels have the characteristic radial trend of high-extinction lines, this significant over-density region at Galactocentric azimuths 135° < θ <  180° is clearly visible in both hemispheres, and is not expected in an axisymmetric Galactic potential. This over-density will require further detailed analysis in the future. Combining the information in Figs. 3 and 4, we conclude that the stellar density maps of the young OB and old RGB samples clearly show opposite over-densities with respect to the centre–anticentre line at Galactocentric radial bins in the range 8 <  R <  12 kpc.

Fig. 4. Two dimensional projections of the RGB sample. Left panel: Z >  0. Right panel: Z <  0. The size of the bins are 0.16 kpc × 0.32 kpc. The colour mapping scales logarithmically with number counts N. We mark the position of the Sun with a star and the grey solid lines show the Galactocentric azimuths: 135°, 180°, and 225°. The black dotted lines show circles at different Galactocentric radii, R = 8, 10, 12, 14, and 16 kpc.

4.2. Spatial signature of the Galactic warp

We now focus on the stellar signature of the warp detected at larger heliocentric distances by plotting the median Z for each of the samples (see Fig. 5). Both samples show a clear positive median in Z at Galactic azimuths 90° < θ <  135° and a negative median at 225° < θ <  270°, well in agreement with the expected behaviour of the Galactic warp (e.g. Drimmel et al. 2000; Levine et al. 2006). For intermediate azimuths, 135° < θ <  225°, where the signature of the line-of-nodes lies, the data show a high degree of complexity, with significant differences between the two populations. By comparing the warped mock catalogues with a Sine Lopsided warp model and a straight line-of-nodes in the Galactic centre–anti-centre direction (see middle panel of Fig. 6), we can associate a stripe with null median Z towards the anti-centre direction as the possible location of the line-of-nodes of the Galactic warp. We see that the extinction can modify the shape of the line-of-nodes, and therefore any solid conclusion about the position of the line-of-nodes requires a detailed correction for extinction. The RGB sample with more than 18 million stars shows a null median Z stripe inclined with respect to the Sun–anticentre direction towards θ ∼ 180−200°. The fact that the stripe is curved and with Galactic azimuth growing with Galactocentric radius could lead to the conclusion that the line-of-nodes is twisted. However, as mentioned above, the position could also be affected by the extinction. We emphasise that these results do not agree with previous studies. As an example, Momany et al. (2006) classically placed the line-of-nodes at θ <  180°. Regarding the young population, our OB sample presents a clumpy distribution of median Z, without a clear trace of the line-of-nodes. Some recent studies using classical Cepheids, place the line-of-nodes of the warp also at θ <  180° (Skowron et al. 2018; Chen et al. 2019). The Cepheids population is very young, so we would expect the OB sample to be comparable to the Skowron et al. (2018) and Chen et al. (2019) results. However, as mentioned above, the clumpy structure of the median Z we find for this young population does not allow us to support or rule out their conclusion.

Fig. 5. Two dimensional projections of the median Z for the two working samples, namely the OB sample (left column) and the RGB sample (right column). The sizes of the bins are 0.32 kpc × 0.64 kpc for the OB sample and 0.16 kpc × 0.32 kpc for the RGB sample. We only plot the bins with at least 15 sources. We mark the position of the Sun with a star. The grey solid lines show the Galactocentric azimuths: 135°, 180°, and 225°, and the black dotted lines show circles at different Galactocentric radii: R = 8, 10, 12, 14 and 16 kpc. The colour scale is different in each sample.

Fig. 6. Two dimensional projections of the median Z for the three mock catalogues, namely the flat-disc null hypothesis (left), the Sine Lopsided warped disc (middle), and the S Lopsided warped disc (right). The sizes of the bins are 0.16 kpc × 0.32 kpc. We mark the position of the Sun with a star. The grey solid lines show the Galactocentric azimuths: 135°, 180°, and 225°, and the black dotted lines show circles at different Galactocentric radii: R = 8, 10, 12, 14 and 16 kpc.

For comparison, in Fig. 6 we show the 2D projection of the median Z for the three mock catalogues, namely with a flat-disc null hypothesis, with a Sine Lopsided warp model, and with the S Lopsided warp model, from left to right, respectively. In the left panel, the disc imposed is flat and therefore the median Z is around zero, as seen in the large extent of the disc. There are small regions in which the median deviates slightly from zero because of the effect of the extinction and large astrometric errors. In the case of a mock catalogue including the effect of the Galactic warp (middle and right panels of Fig. 6), we see the expected median of Z for such a warped model, including some line-of-sight features corresponding to lines with high astrometric errors, namely l ∼ 50° and l ∼ 300−310°. These lines are also seen in Gaia DR2 data (see Fig. 5), but the effects of the extinction or Gaia DR2 errors do not mask the spatial signature of the Galactic warp. In the middle panel, for the SineLop model, we note the grey stripe that should correspond to the line-of-nodes bends in the anti-centre direction towards Galactic azimuths θ >  180°, as we also see in the Gaia DR2 data and mentioned above. On the other hand, in the SLop model (right panel), the grey stripe is broader as expected by construction of the SLop model (see Fig. B.1), thus it is not limited to a small region along the line-of-nodes.

Finally, in Fig. 7 we plot the density contours in the (Y,Z) projection. In this projection, the X-axis is perpendicular to the plot and we expect to see the up and down side of the warp towards the Y-positive and Y-negative axes, respectively. We show the median of the three samples, namely All-Stars (black), OB (blue), and RGB (red) samples. As expected, the median of the density in the (Y,Z) projection for the three samples is not flat, but it describes a warped profile. More importantly, the Gaia DR2 data definitively confirm that the amplitude of the warp is increasing with the age of the population.

Fig. 7. Contour density plots in the (Y,Z) projection of the All-Stars sample in grey scale. Solid lines show the median for the three samples: in black, blue, and red, the All-Stars, OB, and RGB samples, respectively. We note that the Y-axis is inverted. The line shows the median of the particles in bins of 2 kpc and the error bars show the lower and upper 1σ uncertainty in the position of the median.

The altitude of the warp at R = 14 kpc for the OB and RGB samples is given in the top rows of Table 1, where we also provide values given in other works found in the literature, for comparison. Our results show that the altitude of the Galactic warp at R = 14 kpc is of the order of 0.2 kpc and 1.0 kpc for the OB and RGB samples, respectively. These values are in agreement with those obtained by Drimmel & Spergel (2001) and López-Corredoira et al. (2002) using a Red Clump population, those by Reylé et al. (2009) using 2MASS star counts also dominated by the red giant population and, finally with those of Marshall et al. (2006) traced by dust. We also include the results derived using pulsars (Yusifov 2004) as it is generally accepted that progenitors of pulsars are young OB stars. The estimated altitude of the warp using the OB sample is smaller than that determined by the pulsars. Even though our results can be biased by the effects of the extinction, our results suggest that Amôres et al. (2017) values for all age populations are overestimated. Finally, we also want to mention that the RGB sample shows an inflection point at around R ∼ 13−14 kpc at the down side, similar to the one reported by Levine et al. (2006) when tracing the warp in the gas. Further work is required to analyse why such different tracers present a similar feature in this region.

Our RGB sample shows only a small degree of asymmetry in the sense that |Z(up)| ≲ |Z(down)|, in agreement with Marshall et al. (2006) and opposite to the values reported by Amôres et al. (2017), who, in the age range 2−7 Gyr, always found that |Z(up)| > |Z(down)|.

To conclude, the stellar warp is highly dependent on the age of the tracer population. Whereas no solid conclusion can be established from the OB sample, the huge amount of data of the Gaia DR2 RGB sample suggests a slightly lopsided warp in the sense that |Z(up)| ≲ |Z(down)|, with a line-of-nodes of the Galactic warp twisted towards Galactocentric azimuths θ ∼ 180−200° at R >  12 kpc.

5. The kinematic signature of the Galactic warp

In this section, we analyse the kinematic distribution of both working samples using three strategies. We start by projecting the median proper motion distribution into the Galactic plane (Sect. 5.1) and then we apply the two methods described in Sect. 2.1, namely the LonKin (Sect. 5.2) and the nGC3 PCM methods (Sect. 5.3), which extract and quantify the overall signal of the warp using both position and kinematic information. We show simplified theoretical expectations for three different warp models⁵ as a guide on how to interpret these results.

5.1. The kinematics in the (X,Y) plane

In Fig. 8, we plot the median of the proper motion in latitude corrected for the peculiar motion of the Sun projected into the (X,Y) plane and we analyse any substructures seen in this projection.

Fig. 8. Two dimensional maps of the median proper motion in the latitude direction corrected for the peculiar motion of the Sun, μb, LSR of the two working samples, namely the OB sample (left) and the RGB sample (right). The black dotted lines show circles at different Galactocentric radii, R = 8, 10, 12, 14 and 16 kpc. The sizes of the bins are 0.32 kpc × 0.64 kpc for the OB sample and 0.16 kpc × 0.32 kpc for the RGB sample. We only plot the bins with at least 15 sources. The grey solid lines show azimuths 135° and 225°, while the black horizontal line shows the centre–anti-centre direction. The position of the Sun is shown with a small black star.

The median μ_b, LSR distribution shows some radial features in the kinematic map, such as those at l ∼ 55°, 135°, 200°, and 220°. Such features are also present in Poggio et al. (2018), with similar but less extended samples. By using the mock catalogues, we study whether we can attribute these radial features to effects of the extinction. In Fig. 9 we show the median μ_b, LSR distribution for the mock catalogues with flat disc (null hypothesis, left panel) and with a Sine Lopsided and S Lopsided warped disc (middle and right panels). As expected, the proper motion in latitude has a median around zero in the case of a flat disc. According to a simple symmetric warp model (see Appendix B.1), as stars are moving from the down-side at θ ∼ 270° to the up-side at θ ∼ 90°, the expected distribution of the proper motion in latitude μ_b, LSR has maximum and positive values towards the anti-centre, decreasing from there towards both sides (Abedi et al. 2015). In the case of a Sine Lopsided disc, the maximum of the proper motion is shifted towards the side of larger amplitude (the up-side in the mock catalogues), while in the case of a S Lopsided disc the proper motion in latitude has a modulation, being almost zero along the line-of-nodes and maximum towards the up and down lines-of-sight. We can conclude that the effects of extinction or Gaia DR2 errors do not affect the predicted pattern.

Fig. 9. Two dimensional projections of the median μb, LSR for the three mock catalogues, namely the flat-disc null hypothesis (left column) and the Sine Lopsided warp disc (middle column) and the S Lopsided warp disc (right column). The size of the bins is 0.16 kpc × 0.32 kpc. We mark the Sun position with a star. The grey solid lines show the Galactocentric azimuths: 135°, 180°, and 225°, and the black dotted lines show circles at different Galactocentric radii, R = 8, 10, 12, 14, and 16 kpc.

The pattern we observe in the 2D projections in the OB and RGB samples (see Fig. 8) shows that, towards the anti-centre, the median μ_b, LSR does not have a uniform positive trend, but shows modulated changes. Whereas the median μ_b, LSR of the OB sample has a transition from negative to positive values, the reverse is observed in the RGB sample. This is a clear first signature of the age dependency of the vertical motion. Its possible relation to a bending mode (Gaia Collaboration 2018b) is discussed in Sect. 6.

These maps also allow us to detect any possible asymmetry with respect to the Sun–anti-centre line, which might be related to the lopsidedness of the Galactic warp. For the RGB sample (right panel of Fig. 8) we detect a significant stripe with maximum positive μ_b, LSR towards Galactic azimuths θ ∼ 160−170°. Furthermore, this shift of about 20−30° from the Sun–anti-centre line is also visible throughout the map, causing the null median values to be asymmetric (further discussed in Sect. 5.2). For the OB sample (left panel of Fig. 8), the kinematic structure is more clumpy, but again there is a clear shift with respect to the Sun–anti-centre line, with lower negative median μ_b, LSR towards the down side of the warp (l ∼ 270−310°).

Finally, we report the existence of a feature at l ∼ 100−120° at the outermost radii R >  12 kpc in the RGB sample. We refer to it as the “blob”. The lack of OB stars in this area prevents us from deciphering whether or not there exists a counterpart of this feature in the young population. We suggest that the blob is not related to the warp and we discuss it further in Sect. 6.

5.2. The warp as seen by the LonKin method

Here we show the LonKin result for the OB and RGB, together with a toy theoretical prediction. Figure 10 shows the results for the OB sample (blue lines) and the RGB sample (red lines) in Galactocentric radial bins of increasing distance from top to bottom panels. Only the longitude bins for which there are at least 300 stars are shown. Figure D.2 shows the model predictions for a random realisation of three warp models, described in detail in Appendix D.

Fig. 10. LonKin method applied to the OB sample (blue lines) and RGB sample (red lines). We divide into three panels according to Galactocentric rings, from the inner to the outer disc, from top to bottom. In each panel, different distance bins are indicated with solid, dotted and dashed lines consistently from nearest to farthest. The shaded vertical region in the bottom panel corresponds to the “blob” (see text for details). The line shows the median of the particles in bins of 20° and the error bars show the lower and upper 1σ uncertainty in the position of the median, including the uncertainty in the solar velocity in quadrature.

As already hinted at in the previous section, Gaia DR2 data reveal clear differences between the two populations. As for the OB stars, the inner annuli show an almost flat trend with a median proper motion of about −0.2 mas yr⁻¹, while for the RGB, the median is slightly positive. In the middle annuli, the proper motion increases, as expected in a simple warp model, showing a relatively planar trend for the OB stars in the outer disc, from l = 90° to l = 270°, and a sharp decrease towards the inner disc. Conversely, the RGB stars clearly show that the maximum in the proper motion is not achieved in the anti-centre direction, while it shows two local maxima at about l = 140° and l = 220°, the one in the second quadrant being slightly higher. As we move towards the outer annuli (bottom panel), the OB star distribution becomes clearly asymmetric (see top left panel of Fig. 3), and the LonKin only provides information in the third and fourth quadrant, showing a maximum in the proper motion around l = 220°. The proper motion for the RGB sample shows a clearly different trend from that of the OB stars. The trend that we could initially only just make out in the middle annuli is clear here, with the maxima moving slightly towards the inner disc. Also, the median proper motion in the anti-centre direction is almost zero for the radial bin 14 <  R <  15 kpc, and becomes negative in the outermost ring. We can also distinguish a significant local minimum at l ∼ 120° corresponding to the “blob” already detected in the 2D maps shown above (see Fig. 8). This minimum in proper motion in latitude could be related to a negative median Z at the same region (see right panel of Fig. 5) due to high extinction and the fact that the distribution of RGB stars is not homogeneous above and below the disc. It might also be caused by some type of substructure present in the southern hemisphere.

For comparison, in Fig. 11 we show the trend in proper motion predicted by the three mock catalogues, namely the flat disc (dotted lines), the SineLop (dashed lines), and the SLop (solid lines) warped disc models. As expected, the proper motion in latitude for a flat disc is approximately null, while for a SineLop warped disc it increases towards the outer parts and is not symmetric with respect the line-of nodes (l = 180°). For the SLop warped disc it shows a modulation around zero towards the line-of-nodes, and has two relative maxima towards the largest amplitudes. As mentioned in Sect. 2.3, the SineLop and SLop models are fixed using the free parameters given by Amôres et al. (2017) with ⟨age⟩ = 2.5 Gyr, even though we have found these to be over-estimated (Sect. 4). The amplitude of the curves in Fig. 10 supports our findings from Sect. 4 that the warp has a low amplitude. The fact that the proper motion for the RGB data has a double-peak tendency and that it becomes negative towards the anti-centre at large Galactocentric distances (14 <  R <  15 kpc) indicates that an S Lopsided model would be the correct choice as a starting model (see further discussion in Sect. 6).

Fig. 11. LonKin method applied to the mock catalogues without warp (dotted lines), with the SineLop warp (dashed lines) and with the SLop warp (solid lines) models, including extinction, Gaia selection function, and Gaia DR2 errors. The line shows the median of the particles in bins of 20° and the error bars show the lower and upper 1σ uncertainty in the position of the median and they include in quadrature the uncertainty in the solar velocity.

The LonKin method also allows us to estimate the initial radius of the warp. We find that the initial radius for the OB sample is 12 <  R <  13 kpc, while for the RGB sample it is 10 <  R <  11 kpc. We again confirm an age dependence of the warp characteristics, as already shown in Sect. 4.

In Fig. 12 we check whether the results of applying the LonKin method shown in Fig. 10 depend on the choice of the distance estimator. We show the results of applying the LonKin method to the radial annuli 13 <  R <  14 with L = 2 kpc (solid line) and with L = 1.35 kpc (dotted line) for the OB (blue line) and the RGB (red lines) samples. We note that, qualitatively, both sets of curves are consistent within their respective error bars, which overlap in most longitude bins. Therefore, our results are robust against the choice of the distance estimator.

Fig. 12. Comparison of results for the LonKin methods for the radial bin R = 13−14 kpc obtained for two different prior scale-lengths, L = 2 kpc (solid lines) and L = 1.35 kpc (dotted lines) in red for the RGB sample and in blue for the OB sample.

5.3. nGC3 signature of the warp

Here we show the nGC3 result for the OB and RGB, together with a toy theoretical prediction. The nGC3 PCMs for the OB and RGB samples are shown in the left and right panels of Fig. 13, respectively, for r_gal bins of 0.5 kpc in width. Qualitatively, the two sets of PCMs for both tracers have a similar behaviour: there is a peak at the north Galactic pole, that is, at the centre of the PCM, that remains fixed at all r_gal and another peak in the lower-left quadrant of the PCM, whose co-latitude increases with r_gal. This latter peak is the clear signature of a warped disc whose amplitude increases with radius. Comparing the PCMs for the two samples, we see that the amplitude of the warp for the RGB stars is larger than for the OB stars at a given radius, in the sense that the centroid of the peak that traces the warp (cross) reaches co-latitudes close to 4° for the RGBs compared to for the OBs, at r_gal ∼ 14 kpc.

Fig. 13. nGC3 PCMs for the OB (left panels) and RGB samples (right panels). The corresponding rgal range is indicated at the top right of each panel. The cross indicates the centroid (weighed by pole counts) of the main off-pole peak of the PCM. The maps were produced using a grid spacing and tolerance of . The colour scales linearly with star counts, indicating larger values with darker colours. In all panels, dotted circles and radial lines correspond to parallels at co-latitudes increasing by 2° and meridians at azimuths every 30°, respectively.

These tilt values estimated from the PCMs are compatible with those shown by the Y-Z density plot in Fig. 7 at the same distance. The behaviour of the tilt angle as a function of r_gal is summarised in Fig. 14, where the larger amplitude of the tilt angle for the RGB compared to the OB is evident, at all radii. The figure also compares these results with those obtained with a shorter prior scale length (L = 1.35 kpc); for the two tracers, both sets are in excellent agreement, showing our results are robust against the choice of the distance estimator.

Fig. 14. Comparison of results for the nGC3 method obtained for two different prior scale lengths, L = 2 kpc (crosses and filled circles) and L = 1.35 kpc (open squares and open diamonds) in red for the RGB sample and in blue for the OB sample. We show the inclination angles of the main nGC3 PCM peak (crosses in Fig. 13) as a function of rgal. The plateau in the tilt angle observed for the RGBs (rgal ≲ 9 kpc) is due to contamination from the central peak affecting the centroid when the two peaks are close.

Figure 15 shows the predictions given by the three simple warp models applied to the same mock catalogues as in Fig. 11, at the Galactocentric ring r_gal = 13−14 kpc. The simultaneous signals observed at the north pole and at increasing latitude seem incompatible with the tilted-ring-like SineLop model and closer to the prediction of the SLop model. Although the SineLop model could explain the peak whose co-latitude increases with distance, it cannot account for the signal observed in the north pole and remains fixed with distance, while additionally requiring the line of nodes to be rotated by as much as ∼30°. On the other hand, the SLop model predicts qualitatively similar signals simultaneously: the fixed signal at the pole, and the signal at increasing co-latitude that traces the increase of the mean amplitude of the warp with distance. Figure 15 shows this prediction for the SLop model in the right panel, with signal at azimuth 240°, similar to the observed peak in the PCM⁶. This suggests the SLop model can do a better job at explaining the data, although we do caution that the details of the observed signal do not exactly match the predictions of the model, possibly due to a combination of model mismatch and differences in the extinction and selection function of the real and mock catalogues.

Fig. 15. nGC3 PCMs for a Gaia DR2 mock catalogue (with errors, extinction, and selection function) of test particles at rgal between 13 and 14 kpc, for the flat disc (left) the Sine Lopsided warp (middle) and the S Lopsided warp (right). We note the asymmetry in the PCM signature of the SLop model, which is due to the stars lost in the I and IV quadrant when the effect of extinction and the selection function are included. The corresponding PCMs for the SineLop and SLop models without the effect of errors, extinction, and selection function is shown for comparison in Fig. D.5.

Also, the inclined shape of the peak that traces the warp (centroid indicated with a cross) seems very similar to the PCM signature shown for mock catalogue PCMs in the SLop model. We do stress that the theoretical predictions given by the simple warp models do not fully reproduce the features in the observed PCMs, though the S Lopsided model is the one that qualitatively better matches the data, as already suggested by the LonKin method.

We emphasise that for the SLop model the mean tilt angle derived from the PCM is not trivially related to the up and down tilt angles of the model due to the overlap of the two signatures (see discussion in Appendix D). Deriving values for these parameters would require detailed modelling. Looking at the change of slope of the tilt angle as a function of r_gal in Fig. 12 we estimate the starting radius of the warp to be around 12−12.5 kpc for the OB sample and around 10.5−11 kpc for the RGB sample, in agreement with the reported values using the LonKin method.

6. The richness and complexity of the vertical motion

From Sects. 4 and 5, we clearly confirm the presence of the Galactic warp in Gaia DR2 data. The projections of the positions to the plane perpendicular to the Sun–anti-centre line (Fig. 7) of the OB and RGB samples reveal that the amplitude of the warp is different according to the age of the population, the warp in the RGB sample being more prominent than in the OB sample.

The simple theoretical predictions we have shown are intended to shed some light and help us interpret the results in terms of different warp models. The 2D projection into the Galactic plane suggests that the line-of-nodes in the RGB sample is shifted towards Galactocentric azimuth θ >  180° and twisted in this direction from R >  12 kpc, though this shift could be partially induced by the extinction. The results of applying both the LonKin and nGC3 PCM methods, together with the kinematic 2D projections, suggest that we can indeed rule out that the Galactic warp is symmetric. Both the SineLop and SLop models predict that when the warp is lopsided, the maximum proper motion does not coincide with the line-of-nodes, as we have clearly observed with the LonKin method and in the kinematic 2D maps, with both tracers. Poggio et al. (2018), using similar tracers but with smaller radial and magnitude coverage, and Chen et al. (2019), using Cepheids, note that the maximum vertical velocity they observe in their maps is not along the Sun–anti-centre line and say this might be due to the Sun not being in the line-of-nodes. Here, we trace that kinematic signature further out in the disc and give an alternative interpretation based on our model predictions: that the offset between the line-of-nodes and the line of maximum vertical proper motion is due to the lopsidedness of the warp of the disc.

To be sure that the features we see in the kinematic maps and are not an artefact due to the contribution of the rotation curve to the proper motion in latitude, we make the following test. By assuming a constant rotation curve, with rotational velocity V_c(R) = 240 km s⁻¹, we assign at each star in the RGB sample the velocity due to the flat rotation curve at the observed distance of each star and we compute the median proper motion in latitude due only to this effect. The result of this test is shown in the left panel of Fig. 16. The offset between the line-of-nodes and the line of maximum vertical proper motion is not reproduced in this test, meaning that the rotation curve does not contribute to the outcomes we claim are due to the warp. We do not see any of the other effects discussed below. To highlight this, in the right panel we also show the result of subtracting this contribution from the observed map of Fig. 8. The resulting plot shows the same features as our original plot (Fig. 8, right), confirming that they cannot be accounted for by a combined effect of the rotation curve and selection function of our sample.

Fig. 16. Contribution of a flat rotation curve (with the selection function of our sample) to the proper motion in latitude. Left panel: given the spatial distribution of the RGB sample we compute the contribution to the proper motion in latitude with respect to the LSR of a flat rotation curve. Right panel: we subtract from the proper motion in latitude the contribution due to a flat rotation curve. Bin size is the same as in the RGB sample in Fig. 8. Black and grey lines and the black small star are the same as in Fig. 8.

Using both the LonKin and PCM nGC3 methods we can establish a dependence on age for the starting radius of the warp, with the starting radius of the OB sample being farther out than that of the RGB sample. More precisely, the starting radius for the young OB population would be around 12−12.5 kpc and that of the older RGB population around 10.5−11 kpc. This trend is in agreement with Amôres et al. (2017), though they predict slightly smaller onset radii of about 10 kpc and 9 kpc for the age bins 0−1 Gyr and 3−7 Gyr, respectively, and taking into account that the authors place the Sun at R_⊙ = 8 kpc. Using the OB stars of the TGAS-HIPPARCOS subsample, which extend up to 3 kpc from the Sun, Poggio et al. (2017) find that the proper motions of the nearby OB stars are consistent with the signature of a kinematic warp, while those of the more distant stars (parallax < 1 mas) are not. The authors also suggest that additional phenomena may cause systematic vertical motions that are masking the expected warp signal. Bobylev (2010) uses Tycho-2 nearby (0.3 <  d <  1 kpc) Red Clump stars to infer the parameters of the local warp from their kinematics. Already at this close distance to the Sun, the author detects the signal of the warp in the deformation tensor. Poggio et al. (2018), analysing similar young and old samples up to G <  16 mag covering a smaller area in the disc, find that both samples show similar vertical motions, without providing any indication of the starting radius. Gaia Collaboration (2018b) show different vertical velocity maps for the Giant sample from those of the OB sample, again without any mention of starting radius, but indeed signalling differences between both tracers. Schönrich & Dehnen (2018) and Huang et al. (2018), using Gaia-TGAS and LAMOST-TGAS data, respectively, find the signature of the warp already at the solar neighbourhood, without studying different age tracers. Similarly, Kawata et al. (2018) select a stripe along the centre–anti-centre direction using Gaia DR2 from R = 5−12 kpc and find a signal of the warp when plotting the vertical velocity as a function of radius.

Furthermore, the presence of two relative maxima with different amplitude and a negative dip in proper motions in the LonKin plots for the RGB sample favour the SLop model as a more suitable description of the warp. The fact that the LonKin method favours a SLop shape is corroborated by the nGC3 PCM method. The over-density in the PCMs is not point-like, as would be expected from a symmetric warp; instead, it is elongated and rotated with respect to the meridians, which together suggest lopsidedness. In addition, at all radii there is significant signal in the PCMs at the north Galactic pole, a feature characteristic of the SLop model alone, not observed in the SineLop model. The signature observed in the Gaia DR2 data PCMs most resembles what we see in the mock catalogue with a SLop model (see right panel of Fig. 15). As stated above, the simple tilted-ring models used here do not reproduce every feature revealed by the kinematic data, but they suggest that the SLop model could be an appropriate starting point to model the Galactic warp. The SLop model is designed such that the disc is flat not just along the line-of-nodes, like in the SS and SineLop models, but also in a region around it, out to roughly the warp starting radius (see Fig. B.1 for a schematic representation of the three warp models).

Figures 8, 10, and 13 clearly show the high complexity of the Galactic disc structure and it is clear that in addition to it being detected in density in Fig. 7, the kinematics are in agreement with the presence of a Galactic warp. The large amount of data and its high accuracy reveal that the kinematic maps (Fig. 8) are not dominated by a single phenomenon, as we might expect from a Galaxy that has not yet reached equilibrium (Poggio et al. 2017; Antoja et al. 2018). The 2D projection of the median proper motion in latitude corrected from the peculiar motion of the Sun, μ_b, LSR, reveals complex and non-uniform trends. The radial heliocentric features are clearly related to extinction effects. However, the data also highlights wave-like trends that resemble those of a bending mode, apart from that given by the Galactic warp. The bending mode, as defined recently by Widrow et al. (2014) or Gaia Collaboration (2018b), is the mean of the vertical component of the velocity measured at two parallel layers above and below the Galactic plane. In this sense, the projection of median proper motion in latitude can be viewed as a measure of the bending mode at (large) heliocentric distances where the line-of-sight velocity does not dominate the motion. Thus, from the analysis of the proper motion in latitude we can get a glimpse of the effect of the Galactic warp (via the LonKin and nGC3 PCM methods) and the bending modes (via the 2D projection) at large distances. The features with median negative proper motion alternated with median positive values in the 2D map clearly have a wave-like pattern reminiscent of a bending mode (as previously found by (Schönrich & Dehnen 2018).

Throughout this analysis we have detected two significant features in the structure and kinematic characterisation of the RGB sample: Firstly, an over-density of RGB stars above the Galactic plane or, possibly, an under-density of RGB stars below the Galactic plane towards Galactocentric azimuths 135° < θ <  180°. Secondly, in a similar but more constrained azimuthal range and in the outer disc, R >  12 kpc, the so-called blob, a minimum in proper motion in latitude, which is not predicted by the theoretical models and is not consistent with a wave-like pattern either. We wonder whether these two feature are related. This sudden drop in proper motion is located near the region where there is also a lack of stars in the northern hemisphere at the same longitude range (see right panel of Fig. 5). There are not enough stars in the OB sample in that region to draw this conclusion. A similar gap in this quadrant is also seen in the open cluster distribution in the disc reported by Cantat-Gaudin et al. (2018, priv. comm.), while Bisht et al. (2016) find a bend of the Galactic disc towards the southern latitude using young star clusters with a heliocentric distance larger than 2 kpc in the longitude range l = 130° −180°. The origin of the TriAnd overdensity (Majewski et al. 2004; Deason et al. 2014) has for a long time been thought to reside in this same longitude range and below the Galactic plane; recently, using the chemical composition, this overdensity has been related to the Galactic disc, extending well beyond 20 kpc from the Galactic centre. This is a much more distant feature than the blob, but our lack of data in the range from 15 to 20 kpc prevents us from concluding on whether or not there is a connection between the two. From the studies mentioned, it seems that this region shows different features. We therefore suggest that the blob we detect as an overdensity below the Galactic plane with negative median proper motion is a physical feature in the Galactic disc.

7. Conclusions

We used two different tracers from Gaia DR2 data, namely a young bright population, the OB sample (1.8 million sources), and the RGB sample (18 million sources), in order to study their kinematics and relate them to the Galactic warp structure. To achieve this goal, we used different and complementary analyses that provide structural (Sect. 4) and kinematic (Sect. 5) information on the warp. Our results using spatial data only can be summarised as follows.

1. The RGB sample presents a clear asymmetry in stellar density with respect to the centre–anti-centre line at Galactocentric azimuths 135° < θ <  180° clearly visible in both hemispheres, and not expected in an axisymmetric Galactic potential.

2. The RGB sample shows a null median-Z stripe, the line-of-nodes, rotated with respect to the Sun–anti-centre direction towards θ  ∼  180−200° and twisted from R >  10 kpc. This suggests that the line-of-nodes of the Galactic warp is twisted.

3. The Gaia DR2 data definitively confirm that the altitude of the warp increases with the age of the population. The altitude of the Galactic warp at R = 14 kpc is of the order of 0.2 kpc and 1.0 kpc for the OB and RGB samples, respectively.

4. The RGB sample reveals a slightly lopsided Galactic warp in the sense |Z(down)| − |Z(up)|∼250 pc. No solid conclusion can be established from the OB sample, for which lopsidedness seems to be much smaller.

Combining spatial and kinematic information, we make the following conclusions.

1. The median vertical proper motion, μ_b, LSR, values towards the anti-centre show a clear vertical modulation. Whereas the OB sample has a transition from negative to positive values, the reverse is observed in the RGB sample. This is a clear first signature of the age dependency of the bending vertical motion.

2. The offset between the line-of-nodes (θ ∼ 180−200°) and the line of maximum vertical proper motion, detected in a significant stripe towards Galactic azimuths θ ∼ 160−170°, can be naturally associated by our models to the lopsidedness of the warp of the disc.

3. The fact that the vertical proper motion for the RGB data has a double-peak tendency in the LonKin and that it becomes negative towards the anti-centre at large Galactocentric distances (14 <  R <  15 kpc) together with the signature found in the PCM, are again clear indications of lopsidedness and suggests that an S Lopsided model would be the correct choice as a starting model.

4. Both the LonKin and the nGC3 PCM methods also allow us to estimate the initial radius of the warp. We find that the initial radius for the OB sample is 12 <  R <  13 kpc, while for the RGB sample this is 10 <  R <  11 kpc, again with a clear age dependency.

5. The data also highlight wave-like trends that resemble those of bending and breathing modes.

The kinematic maps of both samples reveal a complex disc that is not in equilibrium, in agreement with recent studies (e.g. Antoja et al. 2018) that reflect multiple phenomena. New models with a higher degree of complexity and methods that are flexible enough to be able to capture the complexity of the disc kinematics are required.

Acknowledgments

The authors thank the anonymous referee for a thorough review that helped to improve the manuscript. The authors wish to thank E. Masana and J. M. Carrasco for their help in the computation of A_G for the full Gaia DR2 sample. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This work was supported by the MINECO (Spanish Ministry of Economy) – FEDER through grant ESP2014-55996-C2-1-R and MDM-2014-0369 of ICCUB (Unidad de Excelencia “María de Maeztu”), the European Community’s Seventh Framework Programme (FP7/2007-2013) under grant agreement GENIUS FP7 – 606740. We also acknowledge the team of engineers (GaiaUB-ICCUB) in charge to set up and maintain the Big Data platform (GDAF) at University of Barcelona. CM is grateful for the hospitality and support from ICCUB and IA-UNAM, where part of this research was carried out. We thank the PAPIIT program of DGAPA/UNAM for their support through grant IG100319. The authors acknowledge the use of TOPCAT (Taylor 2005) throughout the course of this investigation.

References

Appendix A: Details of the sample selection

A.1. Cut in f_obs

Here we show how the proposed cut in relative error in parallax, namely f_obs, is appropriate in order to obtain a distance estimator with a minimum bias in distance at large Galactocentric distances. In Fig. A.1 we show the same as in Fig. 1 but without performing any cut on f_obs. The bias in distance at large Galactocentric distances is larger when using any of the scale-lengths of the exponentially decreasing prior, either using the median or the mode.

Fig. A.1. Same as Fig. 1 but without any cut in fobs.

Indeed, without performing any cut in f_obs, the number of OB and RGB stars increases up to 3948339 and 32506963, respectively, which corresponds to an increase of 47% and 55%, respectively. However, the bias in distance obtained at large Galactocentric distances is larger. We therefore decided to keep the proposed cut at f_obs <  0.5.

A.2. The observed HR diagram

We describe here in detail all the steps we performed until we obtained the de-reddened Hertzsprung-Russel diagram (HR diagram) from which we defined the All-Stars sample. As mentioned in the main text, first we select from Gaia DR2 all stars up to magnitude G = 20 with an available parallax measurement and with f_obs <  0.5. This first cut reduces the Gaia DR2 sample to 383 510 799 sources. For these stars, we compute the distance using a Bayesian estimator (see Sect. 2.2) and , the absolute G magnitude of the star uncorrected for absorption. We now want to remove cool main sequence stars from the sample. To do so, we infer what would be the extinction line of Gaia DR2 data from the relation between the absorption in G, A_G, and the reddening, E(G_BP − G_RP) from Arenou et al. (2018), which is about 1.95 and will determine the slope of the extinction line. In Fig. A.2 we show the observed HR diagram, as a function of the observed colour (G_BP − G_RP), and in black dotted lines two extinction curves for 5 mag of absorption. Subsequently, a hypothetical star with , observed colour (G_BP − G_RP) = 4, and absorption A_G ≤ 5 would have an extinction curve like the black dotted one in Fig. A.2. Therefore, we perform a cut parallel to the extinction line , and with a zero point chosen so that a typical star with observed colour (G_BP − G_RP) = 0 would have to include the whole OB sample of the main sequence. This cut is shown as a red curve in Fig. A.2. This second cut reduces the sample to 86 814 618 sources. From this sample we remove the stars in the Large and Small Magellanic Clouds, the Andromeda Galaxy (M 31), and five globular clusters visible in the (l, b) projected density: NGC 6205, NGC 6341, NGC 7078, NGC 1851, NGC 7089 (Gaia Collaboration 2018c), which are distant objects in the halo that may have large uncertainties in parallax that can create artificial over-densities in the vertical space.

Fig. A.2. Observed HR diagram not corrected for absorption or extinction: as a function of the observed colour. The red line shows the cut in and the two black dotted and dot-dashed lines show two examples of the extinction line of 5 mag.

From the All-Stars sample, we select the OB and RGB samples using Eqs. (1) and (2), which is the same recipe as in Gaia Collaboration (2018b). In order to assess the possible contamination in each sample, we use the Gaia mock catalogue of GUMS (Robin et al. 2012) with the errors in parallax scaled to Gaia DR2 performances as in Luri et al. (2018). Since GUMS contains the spectral type, we are able to assess the possible contamination our strategy may include. Therefore, we perform the same cuts as described above and select the OB and RGB samples again as in Eqs. (1) and (2). We estimate that 91% of the stars in the OB sample are of the spectral type O, B, and A0, and therefore only a 9% of the OB sample are contaminants. The RGB sample includes 92% sub-giants, giants, and super-giants, and therefore only 8% of the stars in the RGB sample are contamination, entirely due to dwarf stars.

Appendix B: Warp models

We describe the three warp models we use in this study to compare with data: one symmetric as in Abedi et al. (2014), the “Simple Symmetric warp” model, and two lopsided models first used in this study: the “Sine Lopsided” and the “S Lopsided” models. In Fig. B.1 (first to third row), we show a schematic representation in 3D (left column) and in projection along the line-of-nodes (right column) of the three warp models and an example of the twisted S Lopsided model (fourth row), whose mathematical forms are described in detail in the following subsections.

Fig. B.1. Schematic plot of a warped disc according to the three warp models (see text): First row: simple Symmetric warp. Second row: Sine Lopsided warp. Third row: S Lopsided warp. Fourth row: S Lopsided warp with a twisted line-of-nodes. The red line shows the line-of-nodes and the orange sphere marks the position of the Sun. The blue line, perpendicular to the line-of-nodes, shows the maximum amplitude of the warp. The amplitude values for the schematic plot have been increased for the sake of clarity.

B.1. The Simple Symmetric warp

The Simple Symmetric warp (SS) is already described in detail in Abedi et al. (2014). It consists of splitting the Galactic disc into Galactocentric rings and tilting each one by an angle ψ as a function of the Galactocentric radius. The expression for the tilt angle is:

(B.1)

where the function f(R) denotes the dependence of the tilt angle with radius and we express it as a power law:

(B.2)

where R₁ and R₂ are the initial and final radius of the warp, and the exponent α fixes the shape of the power law.

We notice that once the warp parameters R₁, R₂, and α have been fixed, the tilt angle only depends on the Galactocentric radius R, and has a constant maximum amplitude at R ≥ R₂ of ψ₂. This defines a symmetric warp with respect to the line-of-nodes, i.e. the maximum amplitude of the warp is the same at corresponding points on both sides of the line-of-nodes. According to observational evidence, the up (down) side of the Galactic warp is towards the northern (southern) hemisphere (e.g. Oort et al. 1958; López-Corredoira et al. 2002; Levine et al. 2006).

B.2. The Sine Lopsided model

The Sine Lopsided model (SineLop) is a variation of the Simple Symmetric warp model. It introduces a dependence on the Galactocentric azimuth, θ, while keeping the same radial dependence as the Simple Symmetric warp model. Therefore, the maximum amplitude of the warp changes on either side of the line-of-nodes. The tilt angle is given by:

(B.3)

where A = (ψ_up + ψ_down)/2, B = (ψ_up − ψ_down)/2 and ψ_up and ψ_down are the maximum amplitudes of the warp at the up and down warps, respectively, at the azimuth perpendicular to the line-of-nodes. We note that using this expression, rings are no longer planar, that is a warped ring when seen edge-on is not a flat line, but is slightly bent at the line-of-nodes, as illustrated by the edge-on view shown in Fig. B.1 (second row, right). This curvature is caused by the sine dependence on the azimuth, which allows us to introduce an arbitrary up or down maximum amplitude asymmetry in a smooth (continuous and differentiable) manner across the line of nodes, which is not attainable by simply setting different tilt angles at either side of the line of nodes in the tilted ring model of the previous section. Hence, the expression of Eq. (B.3) provides a lopsided disc (see second row panels of Fig. B.1). In the warp of the Galactic disc, the amplitude in the up side is larger than that in the down side (Levine et al. 2006; Marshall et al. 2006).

B.3. The S Lopsided model

The S lopsided model (SLop) is designed to have an S-shape when seen edge-on, that is, it is flat at the line-of-nodes. It is a variation of the Sine Lopsided model with a different dependence on the Galactocentric azimuth. We first define an angle function given by:

(B.4)

where, as previously, ψ_up and ψ_down are two constant angles that are the maximum amplitudes of the warp at the up and down sides, respectively.

The tilt angle for this model is given by:

(B.5)

B.4. Twisted versions

In all warp models, the line-of-nodes can either be a straight line aligned with the Sun–Galactic centre line or can be twisted by an angle ϕ:

(B.6)

where it increases with radius from the x-axis (see last row of Fig. B.1.

Appendix C: Warp transformations

The purpose of this section is to develop the equations that provide the spatial and velocity coordinates for particles in a warped disc. We apply a warp transformation to the spatial coordinates of particles in a flat disc, and then find the corresponding transformation for the velocities. As we see below, velocities require a different transformation.

C.1. The spatial transformation

As mentioned in the previous section, to warp a flat disc we use a rotation by an angle ψ with respect to the positive x–axis (see Fig. C.1):

(C.1)

where (x′,y′,z′) are the warped coordinates and (x, y, z) are the Cartesian coordinates in the original flat disc. Therefore, the x–axis will be the line of nodes (our warp does not include twisting, which implies a non-straight line of nodes).

In the case of a symmetric warp, the tilt angle ψ is a function of Galactocentric radius only, thus generating a family of concentric, flat, tilted rings. In the general case of a lopsided warp, the tilt angle will also be a function of the azimuthal, Galactocentric angle θ (measured with respect to the positive x–axis), thus introducing the bending needed to produce the asymmetric warp. In the following section we work with the latter, more general case.

C.2. The velocity transformation

One may think that velocities are transformed using the same transformation used for positions, but this is not the case. Velocities are spatial displacements per time, and so they are transformed as:

(C.2)

where summation over repeated indices is implied.

In order to compute the Λ′ matrix, first we obtain the derivatives of the tilt angle in Cartesian coordinates. We note that the dependence of the tilt angle is given by:

(C.3)

The components of the Jacobian of the transformation between polar and Cartesian coordinates are:

(C.4)

The partial cartesian derivatives of the tilt angle are then:

(C.5)

(C.6)

Fig. C.1. Y-Z projection of a tilted circular orbit at R = 14 kpc using each of the three warp models: simple warp (red solid line), simple lopsided warp (blue short dashed line), and S lopsided warp (green long dashed line).

With these, we can now proceed to get the individual components of the Λ′ matrix using Eq. (C.1):

If we now define the following functions: P ≡ z cos(ψ) + y sin(ψ),  Q ≡ y cos(ψ) − z sin(ψ), the Λ′ matrix can be written in a very compact way as:

(C.7)

where the partial derivatives of ψ are given by Eqs. (C.5) and (C.6).

We note that in the case of a fixed tilt angle (inclined plane), or if we consider pure circular rotation along flat, tilted rings, the partial derivatives are null and the Λ′ matrix coincides with the spatial transformation matrix of Eq. (C.1). But in the general case, the ψ-derivatives make the velocity transformation different from the spatial one.

The transformed (but unscaled) velocity components are then given by:

(C.8)

(C.9)

(C.10)

Now, since the velocity transformation redirects the velocity vectors, but should not change their magnitudes, we need to scale the velocity components making sure that the v_x remains constant. In other words, we need to satisfy:

(C.11)

This implies:

(C.12)

which means that the scaling for the transformed velocity components is:

(C.13)

where the scaling factor is:

(C.14)

It only remains to specify (∂ψ/∂R) and (∂ψ/∂θ) for each of the three warp models used.

For the Simple Symmetric warp model, the tilt angle derivatives are simply:

(C.15)

where f(R) is as in Eq. (B.2).

For the Sine Lopsided model, the derivatives of ψ_SL are:

(C.16)

Finally, for the S Lopsided model, the derivatives are:

(C.17)

where g(θ) is as in Eq. (B.4).

In the transformation outlined here, the twist is not included for reasons of brevity and clarity. To include a twist, an additional rotation around the positive z-axis by an angle ϕ should be added to the rotation matrix (Eq. (C.1)) according to Eq. (B.6). The resulting combined rotation matrix should be used to compute the new version of the Λ′ matrix.

Appendix D: Prediction of the kinematic signature of the warp as seen by the Lonkin and nGC3 methods

In this section, we discuss the expected theoretical signatures of the three warp models described in Appendix B for the LonKin and the nGC3 PCM methods.

To predict and understand the kinematic signatures of the different warp models we use three different test scenarios of increasing complexity: a circular orbit, a random realization of warp models, and a test particle simulation.

We fix the free parameters of the warp models to the values derived from Amôres et al. (2017) for a population of Red Clump stars (hereafter, RC): R₁ = 10.1 kpc, R₂ = 14 kpc, α = 1.1 and , and . In the case of the twisted line-of-nodes, we fix ϕ_max = 20° and γ = 1.1. In the remainder of this section we assume these numerical values for the warp model parameters.

D.1. The three test scenarios

The first scenario is the simplest case of a planar circular orbit warped following each of the three warp models. We take a circular orbit at R_gal = R₂ = 14 kpc, with circular velocity of V_c(R₂) = 215 km s⁻¹, assuming the rotation curve of Allen & Santillan (1991). The equations we use to warp positions and velocities are detailed in Appendix C.

The second scenario is a random realisation obtained by taking an initial set of particles in a planar disc and applying position and velocity transformations to each of the particles. As the initial set, we take the set of disc Red Clump initial conditions as in Abedi et al. (2014) and Romero-Gómez et al. (2015). These are a subset of the RGB stars, well defined in the HR diagram, that are generated according to a disc population of RC stars with the corresponding velocity dispersion⁷ and with the density in statistical equilibrium with the Miyamoto-Nagai disc potential (see Appendix A of Romero-Gómez et al. 2015). We note that by construction, the realisation has non-zero radial velocity dispersion.

The third scenario considered is the test particle simulation, for which we take the same initial conditions for a disc RC population as we did for the realisation. The number of particles in the simulation is chosen to have a realistic surface number density of RC stars in the Solar neighbourhood (Czekaj et al. 2014; Romero-Gómez et al. 2015). The integration strategy is the same as in Abedi et al. (2014), where we integrated the initial conditions to obtain a test particle simulation in statistical equilibrium with the imposed potential. This strategy consists of integrating the particles in a flat disc potential for two periods (of an orbit at R_gal = 14 kpc), then integrating for another five periods while the disc potential is changed towards the chosen warp model, and in the third step, integrating two further periods to allow particles to settle into the final warped disc potential. We want to stress that the success of this strategy to obtain a set of particles in statistical equilibrium depends on the warp model. We have confirmed that it is successful for the Simple Symmetric warp and the Sine Lopsided models, but the test particle simulation with the S Lopsided model imposed does not reach statistical equilibrium, even if we increase the number of periods in which the potential is warped with time. In this case, we cannot ensure that the positions and velocities of the test particles using the S Lopsided model will inherit those established by the model.

We generate a Gaia DR2 mock catalogue up to magnitude G = 20 using the prescription of the astrometric, photometric, and spectroscopic formal errors for Data Release 2⁸. We simulate errors for a mission time of 22 months and, for the bright stars (G <  13), we include a multiplicative factor of 3.6 to the error in parallax to match the distribution of uncertainties as a function of the G magnitude observed in Gaia DR2 data.

D.2. Theoretical signature in the LonKin method

First, we consider the simplest case of a warped circular orbit. In Fig. D.1 we show how a warped circular orbit at R_gal = 14 kpc looks when we apply the LonKin method using the three warp models described in Appendix B, with and without twist of the line-of-nodes. As expected from Abedi et al. (2015), the proper motion in latitude, μ_b, LSR, has a maximum in the anti-centre direction when the orbit is tilted using the Simple Symmetric warp model. As soon as the symmetry is broken, that is, in the lopsided models, the maximum is no longer in the anti-centre and the proper motion in latitude has a different behaviour depending on the model used to warp the orbit. The Sine Lopsided model simply displaces the maximum towards the longitude of the maximum amplitude of the warp. When the orbit is warped using the S Lopsided model, the proper motion in latitude becomes zero at the anti-centre direction. We note that in this model, the disc is flat at the line-of-nodes. The main characteristic of the S Lopsided model is the presence of two relative maxima at l around 90° and 270°.

Fig. D.1. LonKin method applied to a warped circular orbit at Rgal = 14 kpc using each of the three warp models: Simple Symmetric warp, SS (black solid line), the Sine Lopsided warp, SL (red solid line), and the S Lopsided warp, SLop (blue solid line). The twisted version of the three models is shown by dotted, dashed, and dotted-dashed lines, respectively.

In Fig. D.2, we apply the LonKin method to a random realisation of particles inside the cylindrical Galactocentric radial bin 13−14 kpc. As expected, the curve corresponding to the Simple Symmetric warp is essentially flat in the anti-centre direction and symmetric with respect to the line-of-nodes, while the signature of the Sine Lopsided warp displaces the maximum towards the longitudes of maximum warp amplitude and the S Lopsided warp shows two clear maxima at l ∼ 90° and l ∼ 270°. We note that warping the ensemble of particles using this last model, the proper motion in latitude is no longer zero at the line-of-nodes. This fact is due to the velocity dispersion of the generated particles.

Fig. D.2. LonKin method applied to random realisations of the reference models, at the Galactocentric ring 13 <  R <  14 kpc. Simple Symmetric warp (SS) in black, Sine Lopsided (SineLop) in red, and S Lopsided (SLop) in blue. The line shows the median of the particles in bins of 20° and the error bars (shaded areas) show the lower and upper 1σ uncertainty in the position of the median.

Finally, in Fig. D.3 we show the result of applying the LonKin method to the test particle simulations of disc RC stars. We divide Fig. D.3 in three panels according to the Galactocentric radial bin of the particles. From top to bottom: 9 <  R <  10 and 10 <  R <  11 kpc; 11 <  R <  12 and 12 <  R <  13 kpc; and 13 <  R <  14, 14 <  R <  15 and 15 <  R <  16 kpc. The warp models used are the Sine Lopsided (left) and the S Lopsided (right). We stress here that only the test particle simulation using the Sine Lopsided model is in statistical equilibrium, while the result of the test particle simulation using the S Lopsided model will not follow the trend predicted for such model (compare with Figs. D.1 and D.2) because it is not in statistical equilibrium. As for the Sine Lopsided test particle simulation, we note that the median proper motion in latitude increases towards the outer disc and the outermost curve shows the expected shape for this model (compare with Figs. D.1 and D.2), namely the maximum of the curve is shifted towards the maximum amplitude of the warp. Thus, if the disc of the Milky Way is warped following a Sine Lopsided model, we should expect a signature in the LonKin method qualitatively similar to that in the left panel of Fig. D.3. As expected, the results regarding the test particle simulation not being in statistical equilibrium using the S Lopsided model have lost the predicted signature for this model in the previous scenarios. The proper motion in latitude does not show the two relative maxima. Since the simulation is not in equilibrium, it is expected that the median values of the proper motion in latitude per longitude bin resemble those for a symmetric warp: we simply recover the increase of the proper motion with respect to the Galactocentric radius (Eqs. (B.1) and (B.2)), but it clearly does not show a lopsided behaviour.

Fig. D.3. LonKin method applied to RC test particle simulations. Left panel: using the Sine Lopsided model. Right panel: using the S Lopsided model. From top to bottom: different radial Galactocentric bins, specified in the legend, expanding from the inner to the outer disc. Horizontal dashed line shows the zero-axis while the vertical dashed line shows the anti-centre direction at l = 180°. As in Fig. D.2, the error bars show the lower and upper 1σ uncertainty of the median.

D.3. Theoretical signature in the nGC3 method

The nGC3 PCM signature for the first test scenario, a circular orbit, is shown in Fig. D.4 for each of the three warp models. The Simple Symmetric warp model produces a single peak (black dot) at a co-latitude corresponding to the tilt of the warp. Lopsided models produce double-symmetric contours, since each side of the warp has a different amplitude ψ_up/down. For the Sine Lopsided model the PCM signature (red line) is a double or bow-tie contour with extremes at co-latitudes corresponding to the up and down warp amplitudes ψ_up/down. For the S Lopsided model, there is also a double contour (blue line) but, in this case, each peak extends from the north pole to the co-latitude corresponding to the maximum tilt angle of the up/down warp amplitudes. This is clearly illustrated by how the circles of the S Lopsided signature (blue) touch the extremes of the Sine Lopsided signature bow-tie (red). Notice also that the signature goes back to the PCM center in between the two loops. This is because the warp disappears next to the line-of-nodes for this model.

Fig. D.4. nGC3 PCM signatures for a circular orbit at rgal = 14 kpc, for each of the three warp models: Simple Symmetric warp (black), the Sine Lopsided warp (red), and the S Lopsided warp (blue). The dotted circles correspond to parallels at co-latitudes of 2°,4°,6°,8° and 10°, with the outermost one (solid black line) corresponding to 10°.

The PCM signature for the second test scenario, a random realisation of particles with r_gal in the range 13–14 kpc, is shown in Fig. D.5. The figure shows that the signatures for the different warp models have similar overall features as in Fig. D.4, with some differences for the S Lopsided model. The Simple Symmetric model produces a single well-defined peak; the SineLop model produces a bow-tie signature caused by two overlapping peaks, barely resolved in this plot; and the SLop model produces two overlapping peaks, in this case with a triangular shape extended in azimuth. For the first two models, the Simple Symmetric and SineLop, the main PCM peaks stay along the ϕ = 270° meridian, which corresponds to the pole’s azimuth of the (straight) line of nodes assumed in this example. On the contrary, for the SLop model, even though the line of nodes is straight, there is signal at azimuths approximately ±30° away from the ϕ = 270° azimuth of the line-of-nodes. This means that for the SLop model the azimuth of PCM regions where there is signal does not give the twist angle straight away; instead, the twist angle is given by the azimuth of the meridian along which the signature shows reflection symmetry.

Fig. D.5. nGC3 PCMs for a random realisation of particles in the rgal range 13–14 kpc, for each of the three warp models: Simple Symmetric warp (left), the Sine Lopsided warp (centre) and the S Lopsided warp (right). The colour scales linearly with star counts, indicating larger values with darker colours. Dotted circles and radial lines correspond respectively to parallels at co-latitudes increasing by 2° and meridians at azimuths every 30°.

The PCM signature for the last test scenario, the test particle simulations, is shown in Fig. D.6. In this figure we show the change in PCM signature going, from top to bottom, from the error-free test particles to the mock catalogue described in Appendix D.1 with Gaia DR2 simulated errors and selection function. For the SineLop model, the left column of the figure shows the bow-tie signature expected from the previous tests, which is only mildly distorted in the last row due to the effect of the observational errors. For the SLop model the signature differs from that of the ensemble of particles, even in the error-free set (top row). As discussed in Sect. 2.1, the test particle set of this model is in an impulsive regime as the SLop warped potential does not allow the particles to reach statistical equilibrium. Therefore, the PCMs shown here are not to be taken as characteristic of the SLop model itself, but as a guide to what the pole count signature might look like for a population of stars that is not in statistical equilibrium with the SLop model.

Fig. D.6. nGC3 PCMs for a Gaia DR2 mock catalogue of test particles at rgal between 13 and 14 kpc, for each of the two lopsided warp models: the Sine Lopsided warp (left) and the S Lopsided warp (right). Top: no errors and no Selection Function. Middle: no errors and Gaia DR2 selection function (G <  20). Bottom: Gaia DR2 errors and selection function (G <  20).

All Tables

All Figures

	Fig. 1. Median (top) and mode (bottom) bias of the true rgal distances of stars selected in bins of inferred rgal, for Gaia DR2 errors. Results for the reciprocal of the parallax are shown with the (yellow) solid+circle line. The prior scale lengths used with the Bayesian estimators are indicated with different colours as shown in the inset keys.
In the text

	Fig. 2. Distributions in magnitude G (top panel) and heliocentric distance (bottom panel) of the two working samples.
In the text

Fig. 3. Two dimensional stellar density map of the two working samples, namely the OB sample (left column) and the RGB sample (right column). The sizes of the bins are 0.32 kpc × 0.64 kpc for the OB sample and 0.16 kpc × 0.32 kpc for the RGB sample. We mark the Sun position with a star and the black dotted lines show circles at different Galactocentric radii, R = 8, 10, 12, 14, and 16 kpc. The colour bar is in log-scale and is the same for both panels. The Galaxy rotates clockwise.

In the text

Fig. 4. Two dimensional projections of the RGB sample. Left panel: Z >  0. Right panel: Z <  0. The size of the bins are 0.16 kpc × 0.32 kpc. The colour mapping scales logarithmically with number counts N. We mark the position of the Sun with a star and the grey solid lines show the Galactocentric azimuths: 135°, 180°, and 225°. The black dotted lines show circles at different Galactocentric radii, R = 8, 10, 12, 14, and 16 kpc.

In the text

Fig. 5. Two dimensional projections of the median Z for the two working samples, namely the OB sample (left column) and the RGB sample (right column). The sizes of the bins are 0.32 kpc × 0.64 kpc for the OB sample and 0.16 kpc × 0.32 kpc for the RGB sample. We only plot the bins with at least 15 sources. We mark the position of the Sun with a star. The grey solid lines show the Galactocentric azimuths: 135°, 180°, and 225°, and the black dotted lines show circles at different Galactocentric radii: R = 8, 10, 12, 14 and 16 kpc. The colour scale is different in each sample.

In the text

Fig. 6. Two dimensional projections of the median Z for the three mock catalogues, namely the flat-disc null hypothesis (left), the Sine Lopsided warped disc (middle), and the S Lopsided warped disc (right). The sizes of the bins are 0.16 kpc × 0.32 kpc. We mark the position of the Sun with a star. The grey solid lines show the Galactocentric azimuths: 135°, 180°, and 225°, and the black dotted lines show circles at different Galactocentric radii: R = 8, 10, 12, 14 and 16 kpc.

In the text

Fig. 7. Contour density plots in the (Y,Z) projection of the All-Stars sample in grey scale. Solid lines show the median for the three samples: in black, blue, and red, the All-Stars, OB, and RGB samples, respectively. We note that the Y-axis is inverted. The line shows the median of the particles in bins of 2 kpc and the error bars show the lower and upper 1σ uncertainty in the position of the median.

In the text

Fig. 8. Two dimensional maps of the median proper motion in the latitude direction corrected for the peculiar motion of the Sun, μb, LSR of the two working samples, namely the OB sample (left) and the RGB sample (right). The black dotted lines show circles at different Galactocentric radii, R = 8, 10, 12, 14 and 16 kpc. The sizes of the bins are 0.32 kpc × 0.64 kpc for the OB sample and 0.16 kpc × 0.32 kpc for the RGB sample. We only plot the bins with at least 15 sources. The grey solid lines show azimuths 135° and 225°, while the black horizontal line shows the centre–anti-centre direction. The position of the Sun is shown with a small black star.

In the text

Fig. 9. Two dimensional projections of the median μb, LSR for the three mock catalogues, namely the flat-disc null hypothesis (left column) and the Sine Lopsided warp disc (middle column) and the S Lopsided warp disc (right column). The size of the bins is 0.16 kpc × 0.32 kpc. We mark the Sun position with a star. The grey solid lines show the Galactocentric azimuths: 135°, 180°, and 225°, and the black dotted lines show circles at different Galactocentric radii, R = 8, 10, 12, 14, and 16 kpc.

In the text

Fig. 10. LonKin method applied to the OB sample (blue lines) and RGB sample (red lines). We divide into three panels according to Galactocentric rings, from the inner to the outer disc, from top to bottom. In each panel, different distance bins are indicated with solid, dotted and dashed lines consistently from nearest to farthest. The shaded vertical region in the bottom panel corresponds to the “blob” (see text for details). The line shows the median of the particles in bins of 20° and the error bars show the lower and upper 1σ uncertainty in the position of the median, including the uncertainty in the solar velocity in quadrature.

In the text

Fig. 11. LonKin method applied to the mock catalogues without warp (dotted lines), with the SineLop warp (dashed lines) and with the SLop warp (solid lines) models, including extinction, Gaia selection function, and Gaia DR2 errors. The line shows the median of the particles in bins of 20° and the error bars show the lower and upper 1σ uncertainty in the position of the median and they include in quadrature the uncertainty in the solar velocity.

In the text

	Fig. 12. Comparison of results for the LonKin methods for the radial bin R = 13−14 kpc obtained for two different prior scale-lengths, L = 2 kpc (solid lines) and L = 1.35 kpc (dotted lines) in red for the RGB sample and in blue for the OB sample.
In the text

Fig. 13. nGC3 PCMs for the OB (left panels) and RGB samples (right panels). The corresponding rgal range is indicated at the top right of each panel. The cross indicates the centroid (weighed by pole counts) of the main off-pole peak of the PCM. The maps were produced using a grid spacing and tolerance of . The colour scales linearly with star counts, indicating larger values with darker colours. In all panels, dotted circles and radial lines correspond to parallels at co-latitudes increasing by 2° and meridians at azimuths every 30°, respectively.

In the text

Fig. 14. Comparison of results for the nGC3 method obtained for two different prior scale lengths, L = 2 kpc (crosses and filled circles) and L = 1.35 kpc (open squares and open diamonds) in red for the RGB sample and in blue for the OB sample. We show the inclination angles of the main nGC3 PCM peak (crosses in Fig. 13) as a function of rgal. The plateau in the tilt angle observed for the RGBs (rgal ≲ 9 kpc) is due to contamination from the central peak affecting the centroid when the two peaks are close.

In the text

Fig. 15. nGC3 PCMs for a Gaia DR2 mock catalogue (with errors, extinction, and selection function) of test particles at rgal between 13 and 14 kpc, for the flat disc (left) the Sine Lopsided warp (middle) and the S Lopsided warp (right). We note the asymmetry in the PCM signature of the SLop model, which is due to the stars lost in the I and IV quadrant when the effect of extinction and the selection function are included. The corresponding PCMs for the SineLop and SLop models without the effect of errors, extinction, and selection function is shown for comparison in Fig. D.5.

In the text

Fig. 16. Contribution of a flat rotation curve (with the selection function of our sample) to the proper motion in latitude. Left panel: given the spatial distribution of the RGB sample we compute the contribution to the proper motion in latitude with respect to the LSR of a flat rotation curve. Right panel: we subtract from the proper motion in latitude the contribution due to a flat rotation curve. Bin size is the same as in the RGB sample in Fig. 8. Black and grey lines and the black small star are the same as in Fig. 8.

In the text

	Fig. A.1. Same as Fig. 1 but without any cut in fobs.
In the text

	Fig. A.2. Observed HR diagram not corrected for absorption or extinction: as a function of the observed colour. The red line shows the cut in and the two black dotted and dot-dashed lines show two examples of the extinction line of 5 mag.
In the text

Fig. B.1. Schematic plot of a warped disc according to the three warp models (see text): First row: simple Symmetric warp. Second row: Sine Lopsided warp. Third row: S Lopsided warp. Fourth row: S Lopsided warp with a twisted line-of-nodes. The red line shows the line-of-nodes and the orange sphere marks the position of the Sun. The blue line, perpendicular to the line-of-nodes, shows the maximum amplitude of the warp. The amplitude values for the schematic plot have been increased for the sake of clarity.

In the text

	Fig. C.1. Y-Z projection of a tilted circular orbit at R = 14 kpc using each of the three warp models: simple warp (red solid line), simple lopsided warp (blue short dashed line), and S lopsided warp (green long dashed line).
In the text

	Fig. D.1. LonKin method applied to a warped circular orbit at Rgal = 14 kpc using each of the three warp models: Simple Symmetric warp, SS (black solid line), the Sine Lopsided warp, SL (red solid line), and the S Lopsided warp, SLop (blue solid line). The twisted version of the three models is shown by dotted, dashed, and dotted-dashed lines, respectively.
In the text

	Fig. D.2. LonKin method applied to random realisations of the reference models, at the Galactocentric ring 13 <  R <  14 kpc. Simple Symmetric warp (SS) in black, Sine Lopsided (SineLop) in red, and S Lopsided (SLop) in blue. The line shows the median of the particles in bins of 20° and the error bars (shaded areas) show the lower and upper 1σ uncertainty in the position of the median.
In the text

Fig. D.3. LonKin method applied to RC test particle simulations. Left panel: using the Sine Lopsided model. Right panel: using the S Lopsided model. From top to bottom: different radial Galactocentric bins, specified in the legend, expanding from the inner to the outer disc. Horizontal dashed line shows the zero-axis while the vertical dashed line shows the anti-centre direction at l = 180°. As in Fig. D.2, the error bars show the lower and upper 1σ uncertainty of the median.

In the text

	Fig. D.4. nGC3 PCM signatures for a circular orbit at rgal = 14 kpc, for each of the three warp models: Simple Symmetric warp (black), the Sine Lopsided warp (red), and the S Lopsided warp (blue). The dotted circles correspond to parallels at co-latitudes of 2°,4°,6°,8° and 10°, with the outermost one (solid black line) corresponding to 10°.
In the text

Fig. D.5. nGC3 PCMs for a random realisation of particles in the rgal range 13–14 kpc, for each of the three warp models: Simple Symmetric warp (left), the Sine Lopsided warp (centre) and the S Lopsided warp (right). The colour scales linearly with star counts, indicating larger values with darker colours. Dotted circles and radial lines correspond respectively to parallels at co-latitudes increasing by 2° and meridians at azimuths every 30°.

In the text

	Fig. D.6. nGC3 PCMs for a Gaia DR2 mock catalogue of test particles at rgal between 13 and 14 kpc, for each of the two lopsided warp models: the Sine Lopsided warp (left) and the S Lopsided warp (right). Top: no errors and no Selection Function. Middle: no errors and Gaia DR2 selection function (G <  20). Bottom: Gaia DR2 errors and selection function (G <  20).
In the text