Gaia kinematics reveal a complex lopsided and twisted Galactic disc warp

TL;DR

This study uses Gaia DR2 data to analyze the complex structure and kinematics of the Milky Way's warp, revealing age-dependent warp features, lopsidedness, and wave-like patterns in the stellar disk.

Contribution

It provides the first detailed kinematic models of the Galactic warp that incorporate age dependency, lopsidedness, and wave-like features based on Gaia data.

Findings

Age-dependent warp height: 0.2 kpc for OB, 1 kpc for RGB.

Warp onset radius: 12-13 kpc for OB, 10-11 kpc for RGB.

Detection of lopsidedness and wave-like bending modes.

Abstract

There are few warp kinematic models of the Galaxy able to characterise structure and kinematics. These models are necessary to study the lopsidedness of the warp and the twisting of the line-of-nodes of the stellar warp, already seen in gas and dust. We use the \Gaia~Data Release 2 astrometric data up to $G=20$mag to characterise the structure of the Galactic warp, the vertical motions and the dependency on the age. We use two populations up to galactocentric distances of $16$kpc, a young (OB-type) and an old (Red Giant Branch, RGB). We use the nGC3 PCM and LonKin methods based on the Gaia observables, together with 2D projections of the positions and proper motions in the Galactic plane. We confirm the age dependency of the Galactic warp, both in positions and kinematics, being the height of the Galactic warp of about $0.2$kpc for the OB sample and of $1.$kpc for the RGB at a galactocentric distance of $14$kpc. Both methods find that the onset radius is $12\sim 13$kpc for the OB sample and $10\sim 11$kpc for the RGB. From the RGB sample, we find from galactocentric distances larger than $10$kpc the line-of-nodes twists away from the Sun-anticentre line towards galactic azimuths $\sim 180-200^{\circ}$ increasing with radius, though possibly influenced by extinction. The RGB sample reveals a slightly lopsided stellar warp with $\sim 250$pc between the up and down sides. The line of maximum of proper motions in latitude is systematically offset from the line-of-nodes estimated from the spatial data, which our models predict as a kinematic signature of lopsidedness. We also show a prominent wave-like pattern of a bending mode different in the OB and RGB, and substructures that might not be related to the Galactic warp nor to a bending mode. GDR2 triggers the need for complex kinematic models, flexible enough to combine both wave-like patterns and an S-shaped lopsided warp.[abridged]

Figures

Tables

Table 1. Table 1: Values of the warp up ( Z u ​ p subscript 𝑍 𝑢 𝑝 Z_{up} ) and down ( Z d ​ o ​ w ​ n subscript 𝑍 𝑑 𝑜 𝑤 𝑛 Z_{down} ) altitude with respect to the Galactic plane at galactocentric radius R = 14 𝑅 14 R=14  kpc expressed in kpc. If the authors report an average value or find no difference between the up and down side, we specify it with the symbol ∼ similar-to \sim in the Z d ​ o ​ w ​ n subscript 𝑍 𝑑 𝑜 𝑤 𝑛 Z_{down} column. We refer to Drimmel & Spergel ( 2001 ) as D01, López-Corredoira et al. ( 2002 ) as LC02, Yusifov ( 2004 ) as Y04 using ATNF Pulsar catalogue (Manchester et al. 2005 ) , Marshall et al. ( 2006 ) as M06, Reylé et al. ( 2009 ) as M09, Amôres et al. ( 2017 ) as A17.

Source	$Z_{up}$	$Z_{down}$	Wavelength
OB (this work)	0.23	-0.19	Gaia optical
RGB (this work)	0.97	-1.22	Gaia optical
D01	1.34	$\sim$	COBE/DIRBE NIR & FIR
LC02	1.23	$\sim$	2MASS NIR
Y04	0.62	-0.58	ATNF pulsars
M06 (dust)	0.74	-0.68	2MASS NIR
R09	0.50	$\sim$	2MASS NIR
A17 (0.5Gyr)	3.3	-0.4	2MASS NIR
A17 (2.5Gyr)	2.1	-1.2	2MASS NIR
A17 (4.0Gyr)	2.0	-1.7	2MASS NIR
A17 (6.0Gyr)	3.9	-2.7	2MASS NIR

Equations

$$M_G$$

$$(G_{BP}-G_{RP}{)}_0$$

$$M_G$$

$$(G_{BP}-G_{RP}{)}_0$$

$${\psi }_{SS}(R;R_1,R_2,\alpha ,{\psi }_2)={\psi }_{2}f(R;R_1,R_2,\alpha )$$

$$f(R;R_1,R_2,\alpha )=\{\begin{array}{ll}0,& R\le R_1\\ {\left(\frac{R-R_1}{R_2-R_1}\right)}^{\alpha },& R_1<R<R_2\\ 1,& R\ge R_2\end{array}$$

$${\psi }_{SL}(R,\theta ;R_1,R_2,\alpha ,{\psi }_{up},{\psi }_{down})=f(R;R_1,R_2,\alpha )[A+B\sin (\theta )]$$

$$g(\theta )=\{\begin{array}{ll}{\psi }_{up},& 0\le \theta <180^{\circ }\\ {\psi }_{down},& 180^{\circ }\le \theta <360^{\circ }\end{array}$$

$${\psi }_{SLop}(R,\theta ;R_1,R_2,\alpha ,{\psi }_{up},{\psi }_{down})=f(R;R_1,R_2,\alpha )g(\theta ){\sin }^2(\theta )$$

$$\varphi (R;R_1,R_2,\gamma ,{\varphi }_{max})=\{\begin{array}{ll}0,& R\le R_1\\ {\varphi }_{max}{\left(\frac{R-R_1}{R_2-R_1}\right)}^{\gamma },& R_1<R<R_2\\ {\varphi }_{max},& R\ge R_2\end{array}$$

$$\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ z^{\prime }\end{array}\right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos (\psi )& -\sin (\psi )\\ 0& \sin (\psi )& \cos (\psi )\end{array}\right)\left(\begin{array}{c}x\\ y\\ z\end{array}\right)=\left(\begin{array}{c}x\\ y\cos (\psi )-z\sin (\psi )\\ z\cos (\psi )+y\sin (\psi )\end{array}\right)$$

$$v_i^{\prime }=\frac{d{x}_i^{\prime }}{dt}=\frac{\partial x_i^{\prime }}{\partial x_j}\frac{d{x}_j}{dt}={{\mathbf{\Lambda }}^{\prime }}_{i,j}{v}_j$$

$$\psi =\psi (R,\theta ),\text{with}R=(x^2+y^2{)}^{1/2},\theta =\arctan y/x$$

$$\frac{\partial R}{\partial x}=\frac{x}{R},\frac{\partial R}{\partial y}=\frac{y}{R},\frac{\partial \theta }{\partial x}=-\frac{y}{R^2},\frac{\partial \theta }{\partial y}=\frac{x}{R}$$

$$\frac{\partial \psi }{\partial x}=\frac{\partial \psi }{\partial R}\frac{\partial R}{\partial x}+\frac{\partial \psi }{\partial \theta }\frac{\partial \theta }{\partial x}=\frac{x}{R}\frac{\partial \psi }{\partial R}-\frac{y}{R^2}\frac{\partial \psi }{\partial \theta }$$

$$\frac{\partial \psi }{\partial y}=\frac{\partial \psi }{\partial R}\frac{\partial R}{\partial y}+\frac{\partial \psi }{\partial \theta }\frac{\partial \theta }{\partial y}=\frac{y}{R}\frac{\partial \psi }{\partial R}+\frac{x}{R^2}\frac{\partial \psi }{\partial \theta }$$

$${{\mathbf{\Lambda }}^{\prime }}_{1,1}=(\partial x^{\prime }/\partial x)=1,$$

$${{\mathbf{\Lambda }}^{\prime }}_{1,2}=(\partial x^{\prime }/\partial y)=0,$$

$${{\mathbf{\Lambda }}^{\prime }}_{1,3}=(\partial x^{\prime }/\partial z)=0$$

$${{\mathbf{\Lambda }}^{\prime }}_{2,1}=\frac{\partial }{\partial x}[y\cos (\psi )-z\sin (\psi )]=-\frac{\partial \psi }{\partial x}[z\cos (\psi )+y\sin (\psi )]$$

$${{\mathbf{\Lambda }}^{\prime }}_{2,2}=\frac{\partial }{\partial y}[y\cos (\psi )-z\sin (\psi )]=\cos (\psi )-\frac{\partial \psi }{\partial y}[z\cos (\psi )+y\sin (\psi )]$$

$${{\mathbf{\Lambda }}^{\prime }}_{2,3}=\frac{\partial }{\partial z}[y\cos (\psi )-z\sin (\psi )]=-\sin (\psi )$$

$${{\mathbf{\Lambda }}^{\prime }}_{3,1}=\frac{\partial }{\partial x}[z\cos (\psi )+y\sin (\psi )]=\frac{\partial \psi }{\partial x}[y\cos (\psi )-z\sin (\psi )]$$

$${{\mathbf{\Lambda }}^{\prime }}_{3,2}=\frac{\partial }{\partial y}[z\cos (\psi )+y\sin (\psi )]=\sin (\psi )+\frac{\partial \psi }{\partial y}[y\cos (\psi )-z\sin (\psi )]$$

$${{\mathbf{\Lambda }}^{\prime }}_{3,3}=\frac{\partial }{\partial z}[z\cos (\psi )+y\sin (\psi )]=\cos (\psi )$$

$${\mathbf{\Lambda }}^{\prime }=\left(\begin{array}{ccc}1& 0& 0\\ -P(\partial \psi /\partial x)& \cos (\psi )-P(\partial \psi /\partial y)& -\sin (\psi )\\ Q(\partial \psi /\partial x)& \sin (\psi )+Q(\partial \psi /\partial y)& \cos (\psi )\end{array}\right)$$

$$v_x^{\prime }=v_x$$

$$v_y^{\prime }=-v_{x}P(\partial \psi /\partial x)+v_y[\cos (\psi )-P(\partial \psi /\partial y)]-v_z\sin (\psi )$$

$$v_z^{\prime }=v_{x}Q(\partial \psi /\partial x)+v_y[\sin (\psi )+Q(\partial \psi /\partial y)]+v_z\cos (\psi )$$

$$v_x^2+v_y^2+v_z^2=(v_x^t{)}^2+(v_y^t{)}^2+(v_z^t{)}^2\mathrm{and}{v}_x=v_x^t$$

$$v_y^2+v_z^2=(v_y^t{)}^2+(v_z^t{)}^2$$

Full text

*Aims. *We use the vertical information coming from the Gaia  Data Release 2 astrometric data up to $G=20$ mag to characterise the structure of the Galactic warp, the related vertical motions and the dependency of the Galactic warp on the age.

*Methods. *We analyse two populations up to galactocentric distances of $16$ kpc, a young bright sample mainly formed by OB stars and an older one of Red Giant Branch (RGB) stars. We use two methods (the Pole Count Maps of Great Circle bands and galactic longitude - proper motion in latitude lines) based on the Gaia   observables, together with two dimensional projections of the positions and proper motions in the Galactic plane.

*Results. *This work confirms the age dependency of the Galactic warp, both in positions and kinematics, being the height of the Galactic warp of the order of $0.2$ kpc for the OB sample and of $1.0$ kpc for the RGB at a galactocentric distance of $14$ kpc. Both methods find that the onset radius of the warp is $12\sim 13$ kpc for the OB sample and $10\sim 11$ kpc for the RGB. From the RGB sample, we find from galactocentric distances larger than $10$ kpc the line-of-nodes twists away from the Sun-anticentre line towards galactic azimuths $\sim 180-200\degr$ increasing with radius, though possibly influenced by extinction. Also, the RGB sample reveals a slightly lopsided stellar warp with $\sim 250$ pc difference between the up and down sides. The line of maximum of proper motions in latitude is systematically offset from the line-of-nodes estimated from the spatial data, which our warp models predict as a kinematic signature of lopsidedness. We also show a prominent wave-like pattern of a bending mode different in the OB and RGB samples. Both positions and kinematics also reveal substructures that might not be related to the large scale Galactic warp nor to the bending mode.

Conclusions. Gaia  Data Release 2 data reveals a high degree of complexity in terms of both positions and velocities that triggers the need for complex kinematic models, flexible enough to combine both wave-like patterns and an S-shaped lopsided warp.

1 Introduction

Warped discs still represent a theoretical challenge nowadays. From surveys of edge-on galaxies, it is clear that $\sim 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 determine 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}\lesssim 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 (hereafter, GaiaDR2), and the availability of proper motions for $1.3$ billion sources (Gaia Collaboration et al. 2018a), has allowed to expand the study to larger volumes and fainter sources. Gaia Collaboration et al. (2018c) 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 et al. (2018c) and show the kinematic evidence of the Galactic warp in vertical velocity using two-dimensional 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 to consider not only simple symmetric tilted rings models, but 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. In this work, we apply both methods to two different populations from GaiaDR2, namely OB-stars and Red Giant Branch 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\leq 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 give 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), and we 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 GaiaDR2 data and give an initial characterisation. In Sect. 4, we study the spatial density distribution of both samples (Sect. 4.1) and we 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 GaiaDR2 data, respectively. In Sect. 6 we combine the spatial and kinematic distributions and we compare it 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; 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 Great Circle Cell Counts (hereafter, 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. First, in Sect. 2.1, the methods are described and specify the requirements needed to apply them. Second, 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, $\mu_{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 uses directly 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 the Local Standard of Rest, $\mu_{b,LSR}$. In the horizontal axis, we plot the heliocentric galactic longitude, $l$, in segments of $20^{\circ}$. For each segment, we compute the median of the $\mu_{b,LSR}$ values and the lower ($15.85$-percentile) and upper ($84.15$-percentile) $1\sigma$ uncertainties on the estimation of the median as defined in the appendix A of Gaia Collaboration et al. (2018c), to which we add in quadrature the uncertainty given by the solar velocity. Since the Galactic warp is a galactocentric feature, we require 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 $\mu_{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 $\mu_{b,LSR}$ in the anti-centre direction. If, on the other hand, the disc is asymmetric, this is, it is Lopsided, and the line-of-nodes, still on the Sun-Galactic Centre line, the maximum in $\mu_{b,LSR}$ is no longer in the anti-centre direction, but shifted towards longitudes coinciding the maximum warp amplitude (see detailed discussion in appendix D).

The family of Great Circle Cell Counts (hereafter, GC3) methods (Johnston et al. 1996; Mateu et al. 2011) comprises 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 (hereafter, 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 radii (see Abedi et al. 2014, for detailed examples). If the warp is lopsided, the signature in the PCM is not a single peak, but it has a shape that depends on the warp model (see detailed discussion in appendix D). In this work in particular, we will 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), so, in order to reach as far as possible in 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_{\mathrm{gal}}$, i.e. for stars in a given inferred $r_{\mathrm{gal}}$ bin, how deviated is the median or mode of their true $r_{\mathrm{gal}}$ distributions. In this plot we show the results obtained using stars with observed parallax error smaller than 50%, i.e. $f_{obs}=|\Delta\varpi/\varpi|<0.5$. We have 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.

The figure shows the behaviour of the distributions is fairly stable up to $r_{\mathrm{gal}}\sim 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 plummets to large negative values systematically. For any given estimator, it also shows 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_{{\mathrm{gal}}}$, the best performance up to the threshold radius ($\sim 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 resemble the most 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_{{\mathrm{gal}}}$, the best performance is consistently obtained with the $L=2$ kpc estimator, which is effectively unbiased up to $\sim 14.5$ kpc and showing the smallest bias ($\sim 0.5$ kpc) at $15.5$ kpc. From $r_{{\mathrm{gal}}}>13.$ kpc, the short scale estimators $L=1.35$ kpc and $L$(BJ18) show a larger bias than for $L=2$ kpc ($\sim 0.3$ kpc), half the corresponding bias expected for the median ($\sim 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 differ too much from those of the $L=2$ kpc estimator when using the mode.

The threshold radius of $\sim 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_{\mathrm{gal}}$ bias due to contaminant stars, which will have a different error distribution. However, this contribution is small, $<10\%$ of the stars in either of our samples (see discussion in A.2) and largely due to dwarf stars. They affect mostly the radial bins at $r_{\mathrm{gal}}<10$ kpc and contribute at most $\sim 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 GaiaDR2 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, C and D.1, respectively.

We generate three Gaia mock catalogues, one to test the null hypothesis, i.e. a catalogue with no warp, the other 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); Romero-Gómez et al. (2015). This set of particles will form 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, i.e. 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 thus show the results of the random realization of particles, i.e. 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 GaiaDR2 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 2111A fortran code to generate the Gaia errors is provided in https://github.com/mromerog. 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 GaiaDR2 data. The final mock catalogues have more than 2 million Red Clump star particles.

3 Selection of the two working samples in GaiaDR2

We are interested in selecting from the GaiaDR2 catalogue 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. We select in the HR diagram of a young population formed by mainly upper main sequence stars or OB-type stars (hereafter, OB sample) and of an older population composed of all the stars in the red giant branch (hereafter, RGB sample). From GaiaDR2222We use the GaiaDR2 catalogue within the Big Data infrastructure known as the Gaia Data Analytics Framework (GDAF) cluster in the Universitat de Barcelona (Tapiador et al. 2017), firstly used for a scientific purpose in Mor et al. (2018). It runs on Apache Spark (https://spark.apache.org/), which is an engine coming from business science suited to deal with large surveys. 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), i.e. an absolute value of the relative error in parallax less than 50$\%$ (i.e., we keep stars with negative parallaxes). This first cut reduces the GaiaDR2 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 $M_{G}^{\prime}$ ¡ 1.95*($G_{BP}$-$G_{RP}$)+2., following the extinction line, where $M_{G}^{\prime}$ is the absolute G magnitude of the star uncorrected for extinction; $M_{G}^{\prime}$ is given by $M_{G}^{\prime}=G-5\log 10(d)+5$; 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_{\mathrm{gal}}>7$ kpc and galactocentric cylindrical distance $R_{\mathrm{gal}}<16$ kpc333The LonKin methods works in galactocentric cylindrical bins, while the PCM works in galactocentric spherical bins.. 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, private communication). We refer to this sample as the ”All-Stars” sample. Finally, the two tracers are selected as in Gaia Collaboration et al. (2018c):

Young population, OB sample:

[TABLE] 2. 2.

RGB sample:

[TABLE]

where $M_{G}=M_{G}^{\prime}-A_{G}$ and ($G_{BP}$-$G_{RP}$)0=($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 typically is formed by young stars with age $\leq 1$ Gyr (e.g. Robin et al. 2012) . Using the Gaia Universe Model Snapshot (Robin et al. 2012) with GaiaDR2 simulated errors, we have 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\times 10^{6}$ OB type stars and $13.6\times 10^{6}$ RGB stars up to $d=5$ kpc.

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,\theta,Z)$ with $R$ the Galatocentric distance and the azimuthal angle $\theta$ measured on the Galactic plane starting from the positive X-axis and towards the positive Y-axis. To perform these coordinate transformations from the GaiaDR2 data, we adopt a distance of the Sun to the Galactic centre of $R_{\odot}=8.34$ kpc (Reid et al. 2014) and the height of the Sun with respect to the Galactic plane of $Z_{\odot}=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_{\odot})=240$ km s*-1444We use these values for consistency and continuation of the Gaia Collaboration et al. (2018c) work. We have checked the more recent values for the Solar radius and circular velocity at the solar radius from Gravity Collaboration et al. (2018) do not change at all the results of the paper.. We adopt the peculiar velocity of the Sun with respect to the local standard of rest of $(U_{\odot},V_{\odot},W_{\odot})=(11.1\begin{subarray}{c}+0.69\\ -0.75\end{subarray},12.24\begin{subarray}{c}+0.47\\ -0.47\end{subarray},7.25\begin{subarray}{c}+0.37\\ -0.36\end{subarray})$ km s-1*(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), being the most prominent the ones we find in the directions: $75\degr<l<85\degr$, $130\degr<l<150\degr$, $180\degr<l<190\degr$ and $250\degr<l<270\degr$, well in agreement with the Drimmel extinction map (Drimmel et al. 2003).

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 one kpc, we could see the Orion, Vela and the Cygnus star forming regions, well described in Zari et al. (2018). At heliocentric distances between $1<d<3$ kpc, we see the Cassiopeia region ($120\degr<l<130\degr$), the Outer Spur ($230\degr<l<250\degr$) and the Carina star forming region ($280\degr<l<300\degr$). In the latter, it is where we have a larger 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\leavevmode\nobreak\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 apparently looks 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 in detail this feature 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\degr<\theta<180\degr$ is clearly visible in both hemispheres, and 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 between 8<R<12\leavevmode\nobreak\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\degr<\theta<135\degr$ and a negative median at $225\degr<\theta<270\degr$, 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\degr<\theta<225\degr$, where the signature of the line-of-nodes lies, the data shows a high degree of complexity, with significant differences between the two populations. By comparing with 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 could 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 so 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 $\theta\sim 180-200\degr$. The fact that the stripe is curved and with galactic azimuth growing with galactocentric radius could induce to think that the line-of-nodes is twisted, though as mentioned above, the position could also be affected by the extinction. We emphasise that these results do not agree with previous works. As an example, Momany et al. (2006) classically placed the line-of-nodes at $\theta<180\degr$. 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 $\theta<180\degr$ (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.

For comparison, in Fig. 6 we show the two dimensional projection of the median Z for the three mock catalogues, namely with a flat disc – null hypothesis, with a Sine Lopsided warp model and 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 warped model, including some line-of-sight features, which correspond to lines with high astrometric errors, namely $l\sim 50\degr$ and $l\sim 300-310\degr$. These lines are also seen in GaiaDR2 data (see Fig. 5), but the effect of the extinction or GaiaDR2 errors do not mask the spatial signature of the Galactic warp. In the middle panel, for the SineLop model, we note the gray stripe that should correspond to the line-of-nodes bends in the anti-centre direction towards galactic azimuths $\theta>180$°, as we also see in the GaiaDR2 data and mentioned above. On the other hand, in the SLop model (right panel), the gray stripe is broader as expected by construction of the SLop model (see Fig. 19 in appendix B), 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/down side of the warp towards the Y-positive/Y-negative axis, 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 important, the GaiaDR2 data definitively confirms that the amplitude of the warp is increasing with the age of the population.

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), López-Corredoira et al. (2002) using a Red Clump population, also 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 over-estimated. Finally, we also want to mention that the RGB sample shows an inflection point at around $R\sim 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)|\lesssim|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 GaiaDR2 RGB sample suggests a slightly lopsided in the sense $|Z(up)|\lesssim|Z(down)|$, with a line-of-nodes of the Galactic warp twisted towards galactocentric azimuths $\theta\sim 180-200\degr$ 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 by applying 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 models555The warp models used are described in detail in appendix B. 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 Sun’s peculiar motion projected into the (X,Y) plane and we analyse any substructures seen in this projection.

The median $\mu_{b,LSR}$ distribution shows some radial features in the kinematic map, such as those at l $\sim 55\degr,\,135\degr,\,200\degr,\,220\degr$. 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 $\mu_{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 $\theta\sim 270\degr$ to the up-side at $\theta\sim 90\degr$, the expected distribution of the proper motion in latitude $\mu_{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-sights. We can conclude that the effect of extinction or GaiaDR2 errors do not affect the predicted pattern.

The pattern we observe in the two dimensional projections in the OB and RGB samples (see Fig. 8) shows that, towards the anti-centre, the median $\mu_{b,LSR}$ has not a uniform positive trend, but shows modulated changes. Whereas the median $\mu_{b,LSR}$ of the OB sample have 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 et al. 2018c) 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 $\mu_{b,LSR}$ towards galactic azimuths $\theta\sim 160-170\degr$. Even more, this shift of about $20-30\degr$ from the Sun–anti-centre line is also visible all through the map, making 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 $\mu_{b,LSR}$ towards the down side of the warp (l $\sim 270-310\degr$).

Finally, we report the existence of a feature at l $\sim 100-120\degr$ 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 does not allow us to know whether 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 22 shows the model predictions for a random realisation of three warp models, described in detail in appendix D.

As already hinted in the previous section, GaiaDR2 data reveals 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*-1*, 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 the OB stars a quite planar trend in the outer disc, from l=$90\degr$ to l=$270\degr$, and a sharp decrease towards the inner disc, whereas 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\degr$ and l=$220\degr$, the one in the second quadrant being slightly higher. As we move towards the outer annuli (bottom panel), the OB stars 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\degr$. The proper motion for the RGB sample shows a clearly different trend from that of the OB stars. The trend that we could initially glimpse 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$\sim 120\degr$, this corresponds to the ”blob” already detected in the two dimensional 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 about null, while for a SineLop warped disc it increases towards the outer parts and it is not symmetric with respect the line-of nodes ($l=180$ °), and for the SLop warped disc it shows a modulation being around zero towards the line-of-nodes and it 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 a right choice as a starting model (see further discussion in Sect. 6).

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 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. 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.

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_{\mathrm{gal}}$ bins of $0.5$ kpc width. Qualitatively the two sets of PCMs for both tracers have a similar behaviour: there is a peak at the north Galactic pole, i.e. the centre of the PCM, that remains fixed at all $r_{\mathrm{gal}}$ and another peak in the lower-left quadrant of the PCM, whose co-latitude increases as $r_{\mathrm{gal}}$ does. 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\degr$ for the RGBs compared to $\sim 2\aas@@fstack{\circ}5$ for the OBs, at $r_{\mathrm{gal}}\sim 14$ kpc.

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_{\mathrm{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.

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_{\mathrm{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 it the line of nodes to be rotated by as much as $\sim 30\degr$. 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 warp’s mean amplitude with distance. Figure 15 shows this prediction for the SLop model in the right panel, with signal at azimuth $240\degr$, similar to the observed peak in the PCM666For more details, see also the theoretical PCM prediction (and its discussion) for the different warp models provided in Appendix D.3 and Figs. 24, 25 and 26.. 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 match exactly the model’s predictions, possibly due to a combination of model mismatch and differences in the extinction and selection function of the real and mock catalogues.

Also, the inclined shaped of the peak that traces the warp (centroid indicated with a cross) seems most alike 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 the LonKin method points to.

We emphasise that for the SLop model the mean tilt angle derived from the PCM is not trivially related to the up/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_{\mathrm{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 GaiaDR2 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, being the warp in the RGB sample more prominent than in the OB sample.

The simple theoretical predictions we have shown are intended to shed some light and help interpreting the results in terms of different warp models. The two dimensional projection into the Galactic plane suggests that the line-of-nodes in the RGB sample is shifted towards galactocentric azimuth $\theta>180\degr$ 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 two dimensional 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 disc’s warp.

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*-1*, 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 in the conclusions we claim are due to the warp. We do not see any of the other effects we discuss below. To highlight this, in the right panel we also show the result of subtracting this contribution to 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.

Using both the LonKin and PCM nGC3 methods we can establish a dependence on age for the starting radius of the warp, being that of the OB sample 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 for the older RGB around $10.5-11$ kpc. This trend is in agreement with Amôres et al. (2017), though they predict onset radii slightly smaller, 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_{\odot}=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 et al. (2018c) show different vertical velocity maps for the Giant sample from those of the OB sample, again without venturing any starting radius, but indeed signalling differences between both tracers. Schönrich & Dehnen (2018); 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 GaiaDR2 from $R=5-12$ kpc and finds 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 an 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, coupled, 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 GaiaDR2 data PCMs resembles the most what we see in the mock catalogue with an SLop model (see right panel of Fig. 15). As stated above, the simple tilted-ring-models used here do not reproduce every feature the kinematic data reveal, 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 the warp starting radius roughly (see Fig. 19 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 the 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 reached yet equilibrium (Poggio et al. 2017; Antoja et al. 2018). The two dimensional projection of the median proper motion in latitude corrected from the Sun’s peculiar motion, $\mu_{b,LSR}$, reveals complex and non-uniform trends. The radial heliocentric features are clearly related to extinction effects. But the data also highlights wave-like trends that remind those of a bending mode, apart from that given by the Galactic warp. The bending mode, as defined recently by Widrow et al. (2014); Gaia Collaboration et al. (2018c), 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 glimpse the effect of the Galactic warp (via the LonKin and nGC3 PCM methods) and the bending modes (via the two dimensional projection) at large distances. The features with median negative proper motion alternated with median positive values in the two dimensional map clearly have a wave-like pattern reminiscent of a bending mode (as previously found by Schönrich & Dehnen 2018).

Throughout all this analysis we have detected two significant features in the structure and kinematic characterisation of the RGB sample: first, 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\degr<\theta<180\degr$. Second, 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 it 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, private communication), while Bisht et al. (2016) find a bend of the Galactic disc towards the southern latitude using young star clusters with an heliocentric distance larger than $2$ kpc in the longitude range $l=130\degr-180\degr$. In this same longitude range and below the galactic plane, there has been long debated the origin of the TriAnd overdensity (Majewski et al. 2004; Deason et al. 2014) which recently using the chemical composition 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 whether there is a connection or not between the two. From all these works, it seems that this region shows different features. So we 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 have used two different tracers from GaiaDR2 data, namely a young bright population, the OB sample (1.8 Million sources), and the Red Giant Branch, 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 have used different and complementary analyses that provide structural (Sect. 4) and kinematic (Sect. 5) information of the warp. Our results using spatial only data are the following:

The RGB sample presents a clear asymmetry in stellar density with respect to the centre–anti-centre line at galactocentric azimuths $135\degr<\theta<180\degr$ clearly visible in both hemispheres, and not expected in an axisymmetric Galactic potential. 2. 2.

The RGB sample shows a null median-Z stripe, the line-of-nodes, rotated with respect to the Sun–anti-centre direction towards $\theta\sim 180-200\degr$ and twisted from $R>10$ kpc. This suggests that the line-of-nodes of the Galactic warp is twisted. 3. 3.

The GaiaDR2 data definitively confirms 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. 4.

The RGB sample reveals a slightly lopsided Galactic warp in the sense $|Z(down)|-|Z(up)|\sim 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 conclude:

The median vertical proper motion, $\mu_{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. 2.

The offset between the line-of-nodes ($\theta\sim 180-200\degr$) and the line of maximum vertical proper motion, detected in a significant stripe towards galactic azimuths $\theta\sim 160-170\degr$, can be naturally associated by our models to the lopsidedness of the disc’s warp. 3. 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 a right choice as a starting model. 4. 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 is $10<R<11$ kpc, again with a clear age dependency. 5. 5.

The data also highlights wave-like trends that remind those of bending and breathing modes.

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

Acknowledgements

The authors thank the anonymous referee for a thorough review that helped improve the manuscript. The authors wish to thank E. Masana and J.M. Carrasco their help in the computation of $A_{G}$ for the full GaiaDR2 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.

Appendix A Details on 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 have a distance estimator with a minimum bias in distance at large galactocentric distances. In Fig. 17 we show the same as in Fig. 1 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.

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, decide 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 have performed until we obtained the de-reddened Hertzsprung-Russel diagram (HR diagram) from which we have defined the All-Stars. As mentioned in the main text, first we select from GaiaDR2 all stars up to magnitude G=20, with an available parallax measurement and with $f_{obs}<0.5$. This first cut reduces the GaiaDR2 sample to $383,510,799$ sources. For these stars, we compute the distance using a Bayesian estimator (see Sect. 2.2) and $M_{G}^{\prime}=G-5\log 10(d)+5$, 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 GaiaDR2 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 it will determine the slope of the extinction line. In Fig. 18 we show the observed HR diagram, $M_{G}^{\prime}$ as a function of the observed colour ($G_{BP}$-$G_{RP}$), and in black dotted lines two extinction curves for $5$ mag of absorption. So a hypothetical star with $M_{G}^{\prime}=5$, observed colour ($G_{BP}$-$G_{RP}$)$=4$ and absorption $A_{G}\leq 5$ would have an extinction curve like the black dotted one. Therefore, we perform a cut parallel to the extinction line $M_{G}^{\prime}$ ¡ 1.95*($G_{BP}$-$G_{RP}$)+2., and with a zero point chosen so that a typical star with observed colour ($G_{BP}$-$G_{RP}$)=0. would have $M_{G}^{\prime}=2$ to include all the OB sample of the main sequence. This cut is shown as a red curve in Fig. 18. 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 (M31), and five globular clusters visible in the (l,b) projected density: NGC 6205, NGC 6341, NGC 7078, NGC 1851, NGC 7089 (Gaia Collaboration et al. 2018b), which are distant objects in the halo that may have large uncertainties in parallax that can create artificial over-densities in the vertical space.

From the All-Stars, we select the OB and RGB using the same recipe as in Gaia Collaboration et al. (2018c) (see Eqns. 1 and 2). 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 GaiaDR2 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 described above and we select the OB and RGB samples again as in Eqns. 1 and 2. We estimate that $91\%$ of the stars in the OB sample are of the spectral type O, B and A0, so only a $9\%$ of the OB sample are contaminants. The RGB sample includes $92\%$ of sub-giants, giants and super-giants, so 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”, and two lopsided models first used in this work: the ”Sine Lopsided” and the ”S Lopsided” models. In Fig. 19 (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.

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 $\psi$ as a function of the galactocentric radius. The expression for the tilt angle is:

[TABLE]

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

[TABLE]

where $R_{1}$ and $R_{2}$ are the initial and final radius of the warp, and the exponent $\alpha$ fixes the shape of the power law.

We notice that once the warp parameters $R_{1}$, $R_{2}$, and $\alpha$ have been fixed, the tilt angle only depends on the galactocentric radius $R$, and has a constant maximum amplitude at $R\geq R_{2}$ of $\psi_{2}$. 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 model. It introduces a dependence on the galactocentric azimuth, $\theta$, 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:

[TABLE]

where $A=(\psi_{up}+\psi_{down})/2$, $B=(\psi_{up}-\psi_{down})/2$ and $\psi_{up}$ and $\psi_{down}$ are the maximum amplitudes of the warp at the up and down warps, respectively, at the azimuth perpendicular to the line-of-nodes. Note that, using this expression, rings are not planar any more, i.e. a warped ring when seen edge-on is not a flat line, but it is slightly bent at the line-of-nodes, as illustrated by the edge-on view shown in Fig. 19 (second row, right). This curvature is caused by the sine dependence on the azimuth, that allows us to introduce an arbitrary up/down maximum amplitude asymmetry in a smooth (continuous and differentiable) manner across the line of nodes, which is not attainable by just 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. 5 provides a lopsided disc (see second row panels of Fig. 19). In the Galactic disc’s warp, the amplitude in the up side is larger than 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:

[TABLE]

where, as previously, $\psi_{up}$ and $\psi_{down}$ are two constant angles which are the maximum amplitudes of the warp at the up and down sides, respectively.

The tilt angle for this model is given by:

[TABLE]

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 it can be twisted an angle $\phi$:

[TABLE]

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

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 will see, 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 $\psi$ with respect to the positive $x$–axis (see Fig. 20):

[TABLE]

Where $(x^{\prime},y^{\prime},z^{\prime})$ are the warped coordinates and $(x,y,z)$ are the Cartesian coordinates in the original flat disc. So 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 $\psi$ 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 be also a function of the azimuthal, galactocentric angle $\theta$ (measured with respect to the positive $x$–axis), thus introducing the bending needed to produce the asymmetric warp. In the next subsection 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:

[TABLE]

where summation over repeated indices is implied.

In order to compute the ${\mathbf{\Lambda}}^{\prime}$ 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:

[TABLE]

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

[TABLE]

The partial cartesian derivatives of the tilt angle are then:

[TABLE]

[TABLE]

With these, we can now proceed to get the individual components of the ${\mathbf{\Lambda}}^{\prime}$ matrix using equation (9):

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

If we now define the following functions: $P\equiv z{\cos(\psi)}+y{\sin(\psi)},\quad Q\equiv y{\cos(\psi)}-z{\sin(\psi)}$, the ${\mathbf{\Lambda}}^{\prime}$ matrix can be written in a very compact way as:

[TABLE]

where the partial derivatives of $\psi$ are given by equations (13) and (14).

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 ${\mathbf{\Lambda}}^{\prime}$ matrix coincides with the spatial transformation matrix of equation (9). But in the general case, the $\psi$-derivatives make the velocity transformation different from the spatial one.

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

[TABLE]

[TABLE]

[TABLE]

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:

[TABLE]

This implies:

[TABLE]

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

[TABLE]

where the scaling factor is:

[TABLE]

It only remains to specify $(\partial\psi/\partial R)$ and $(\partial\psi/\partial\theta)$ for each of the three warp models used.

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

[TABLE]

where $f(R)$ is as in Eq. 4.

For the Sine Lopsided model, the derivatives of $\psi_{SL}$ are:

[TABLE]

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

[TABLE]

where $g(\theta)$ is as in Eq. 6.

In the transformation outlined here, the twist is not included. This for reasons of brevity and clarity. To include a twist, an additional rotation around the positive z-axis by an angle $\phi$ should be added to the rotation matrix (Eq. 9) according to equation (Eq. 8). The resulting combined rotation matrix should be used to compute the new version of the $\mathbf{\Lambda}’$ 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 Sect. 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_{1}=10.1$ kpc, $R_{2}=14$ kpc, $\alpha=1.1$ and $\psi_{2}=7\aas@@fstack{\circ}5$, $\psi_{up}=7\aas@@fstack{\circ}5$ and $\psi_{down}=4\aas@@fstack{\circ}25$. In the case of the twisted line-of-nodes, we fix $\phi_{max}=20\degr$ and $\gamma=1.1$. In the remainder of this section we will 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{\mathrm{gal}}=R_{2}=14$ kpc, with circular velocity of $V_{c}(R_{2})=215$ km s*-1*, 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 positions and velocities 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); Romero-Gómez et al. (2015). They are a subset of the red giant branch stars, well defined in the HR diagram. They are generated according to a disc population of RC stars with the corresponding velocity dispersion777As in Romero-Gómez et al. (2015), we fix the velocity dispersion of RC K-giant stars at the Sun position to be $30.3$ km s*-1*, $23.6$ km s*-1and $16.6$ km s-1*in the radial, tangential and vertical directions, respectively (Binney & Tremaine 2008, and references therein), and a constant scale-height of $300$ pc (Robin & Creze 1986). 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{\mathrm{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 it is successful for the Simple Symmetric 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 GaiaDR2 mock catalogue up to magnitude G = 20 using the prescription of the astrometric, photometric and spectroscopic formal errors for Data Release 2888A fortran code to generate the Gaia errors is provided in https://github.com/mromerog. 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 GaiaDR2 data.

D.2 Theoretical signature in the LonKin method

First, we consider the simplest case of a warped circular orbit. In Fig. 21 we show how a warped circular orbit at $R{\mathrm{gal}}=14$kpc looks when we apply the LonKin method using the three warp models described in Sect. B, with and without twist of the line-of-nodes. As expected from Abedi et al. (2015), the proper motion in latitude, $\mu_{b,LSR}$, has a maximum in the anti-centre direction when the orbit is tilted using the Simple Symmetric model. As soon as the symmetry is broken, i.e. in the lopsided models, the maximum is no longer in the anti-center 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 warp’s maximum amplitude. When the orbit is warped using the S Lopsided model, the proper motion in latitude becomes zero at the anti-centre direction. 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\degr$ and $270\degr$.

In Fig. 22, 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\sim 90\degr$ and $l\sim 270\degr$. 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.

Finally, in Fig. 23 we show the result of applying the LonKin method to the test particle simulations of disc RC stars. We divide Fig. 23 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. 21 and 22) 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. 21 and 22), 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. 23. 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 bins resemble those for a symmetric warp: we simply recover the increase of the proper motion with respect to the galactocentric radius (Eqs. 3 and 4), but it clearly does not show a lopsided behaviour.

D.3 Theoretical signature in the nGC3 method

The nGC3 PCM signature for the first test scenario, a circular orbit, is shown in Fig. 24 for each of the three warp models. The Simple Symmetric model produces a single peak (black dot) at a co-latitude corresponding to the warp’s tilt. Lopsided models produce double-symmetric contours, since each side of the warp has a different amplitude $\psi_{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/down warp amplitudes $\psi_{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.

The PCM signature for the second test scenario, a random realisation of particles with $r_{\mathrm{gal}}$ in the range 13–14 kpc, is shown in Fig. 25. The figure shows that the signatures for the different warp models have similar overall features as in Fig. 24, with some differences for the S Lopsided model. The SS 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 SS and SineLop, the main PCM peaks stay along the $\phi=270\degr$ 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 $\pm 30\degr$ away from the $\phi=270\degr$ azimuth of the line-of-nodes. This means that, for the SLop model, the azimuth of PCM regions were 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.

The PCM signature for the last test scenario, the test particle simulations, is shown in Fig. 26. 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 Sect. D.1 with GaiaDR2 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, this model’s test particle set 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.