© The Authors 2024

1. Introduction

The Milky Way’s disc, similar to those of most other spiral galaxies, is warped. The outer disc is bent upwards in roughly the direction that the Sun is currently moving, and downwards on the opposite side. This was first discovered in 21 cm line measurements of the distribution of H I gas in the Galaxy (Kerr & Hindman 1957; Burke 1957; Westerhout 1957) before being found in the stellar population (Djorgovski & Sosin 1989; Freudenreich 1996) including the Cepheid population, where accurate distance estimation is possible to large distances (e.g. Chen et al. 2019; Skowron et al. 2019). Similarly, warps have been found to be common features in both the gas (e.g. Bosma 1981) and stars (e.g. Sánchez-Saavedra et al. 1990) of spiral galaxies.

The warp influences the movement of stars and leaves a signature in Milky Way stellar kinematics. Systematic vertical velocity trends, dependent on Galactic positions or angular momentum, have long been attributed to the warp (e.g. Dehnen 1998; Drimmel et al. 2000; Poggio et al. 2020; Cheng et al. 2020). For a warp that was not precessing, we would simply expect a net upward motion of stars moving towards the current higher end of the Milky Way’s warp (and downwards for stars moving towards the lower end) in a predictable way. However, a warp that is precessing relatively quickly around the galaxy does, at a certain radius, begin to overtake the stars orbiting around the centre of the Galaxy, and therefore the sign of the stars’ vertical velocity reverses. This is a trend very clearly seen for the Milky Way’s stars by Cheng et al. (2020).

Most kinematic models of the warp assume that the velocity of stars towards or away from the centre of the Galaxy is irrelevant. This assumption underlies “tilted ring” models of galactic warps, which are commonly used for analysis of both external galaxies (e.g. Rogstad et al. 1974; Briggs 1990) and the Milky Way (e.g. Dehnen et al. 2023). It is also an assumption that is often put into models that apply Jeans’ first equation to study the Milky Way’s warp (e.g. Drimmel et al. 2000; Poggio et al. 2020). However, a net radial velocity¹ for a set of stars does correspond to a net movement within the warped disc, which one would expect to correspond to a net motion vertically. As shown by Cheng et al. (2020), this naturally enters into Jeans’ first equation and can be treated quantitatively in the corresponding models.

With Gaia data release 2 (DR2: Gaia Collaboration 2016, 2018a) providing new insight into these stellar kinematics, the warp has been found to be precessing relatively rapidly in the prograde direction (the direction of the disc’s rotation) at rates of 10−14 km s⁻¹ kpc⁻¹ (Poggio et al. 2018; Cheng et al. 2020). We note that Chrobáková & López-Corredoira (2021) dispute this finding, though their estimate of $4_{-4}^{+6}$ km s⁻¹ kpc⁻¹ is consistent with either a stationary or rapidly precessing warp.

For a sample of Cepheids from Gaia data release 3 (DR3: Gaia Collaboration 2023a), Dehnen et al. (2023) found a similarly high prograde precession rate of 9 − 15 km s⁻¹ kpc⁻¹. This is a much smaller number of stars, but with much more accurate and precise distance estimates, and can be considered an independent data sample and analysis. It also specifically probes the warp in the population of young stars, as Cepheids have ages of a few hundred million years or less. This precession is unexpected in the context of the standard picture of a warp made up of tilted rings of stars, in which the precession is expected to be retrograde and much slower.

In most warp models, the line of nodes, ϕ_w, along which the warped disc intersects z = 0, is a parameter used to describe the orientation of the warp. A simple approximation that is suitable for relatively local samples, made for example by Poggio et al. (2018) and Cheng et al. (2020), is that ϕ_w is a constant. However, Briggs (1990) found that the outer parts of warps in external galaxies are twisted, and that ϕ_w forms a leading spiral in the outer part of the warp (this is known as “Briggs’ rule”). For more extended samples, the twist of the warp can be expected to be significant, and indeed both Chen et al. (2019) and Dehnen et al. (2023) found that their Cepheid sample’s warp followed Briggs’ rule beyond 12 kpc from the Galactic centre, but the latter study also found that the warp precession rate decreased outwards such that it could unwind the twisted warp in ∼100 Myr.

The cause of the Milky Way’s warp remains uncertain. Theories include that it is the torque that arises from a misalignment between the Galactic disc and the dark matter halo (e.g. Dubinski & Chakrabarty 2009), accretion of intergalactic matter onto the disc (e.g. Ostriker & Binney 1989), internal amplification of a small misalignment between the inner and outer disc (Sellwood & Debattista 2021), or the interaction with the Large Magellanic Cloud (e.g. Kerr 1957; Burke 1957) or with the Sagittarius dwarf galaxy (e.g. Bailin 2003). A better understanding of the warp’s behaviour and kinematics is needed for us to compare the predictions of these theories to observations.

With the increased data set provided by Gaia DR3, we are in a position to investigate the warp of the Milky Way over a large area, and to investigate the azimuthal variation of the warp’s kinematic signature. We focus purely on the kinematics of the sample because the physical distribution of the sample seen in the Gaia data is enormously affected by the dust extinction of the disc, which is not well understood. We therefore leave the study of the physical distribution of the warp to others, and will compare the results from our kinematic study to the physical warp found by these works. While doing this we are mindful that any limitations of the kinematic model will be translated into an alteration of the shape of the warp in the model. To keep a clear separation between results from kinematic models, and more direct measurements of the current physical shape of the warp, we use the subscript _phys when describing the latter.

Our ability to trace the kinematics of the warp across a large area of the disc allows us to investigate its properties in great detail. It allows us to, for the first time, investigate the twist of the warp’s line of nodes in the outer galaxy. We explore the precession rate and line of nodes of the warp, and how they vary with radius. We also examine the effect of including the radial velocity of stars in our kinematic models, following the pioneering study by Cheng et al. (2020). These authors put all the stars in a given radial range in a single bin, which was characterised by its centre. The influence of the radial velocity on the model vertical velocity is zero along the line of nodes, and the centres of these bins tend to lay along the Sun’s azimuth in the disc, which Cheng et al. (2020) had fixed as the line of nodes. Looking at the stars binned in azimuth gives us a much better lever arm to study the radial velocity’s effect on the vertical velocity in this model, and we show that this effect is significant.

In this paper we describe how we used Gaia DR3 data to study the kinematic signature of the Milky Way’s warp, and to measure its precession rate and line of nodes under different assumptions. Section 2 describes the data sample taken from Gaia DR3, and Sect. 3 describes the kinematic models we used to fit these data. In Sect. 4 we show how the various models fitted the measured average vertical velocities of stars across the Milky Way’s disc, before looking at the properties of these models in Sect. 5. Finally, we discuss these results in context and conclude the paper in Sects. 6 and 7, respectively.

2. Data

To study the kinematic signature of the Milky Way warp we need to know stars’ three-dimensional positions and velocities, which Gaia DR3 provides us to good accuracy for roughly 33 million stars. Since deriving the heliocentric distances from Gaia’s parallaxes alone leads to biased and increasingly uncertain estimates at large distances, we used the Bailer-Jones et al. (2021) photogeometric distances which combine Gaia measured parallaxes with the colour and magnitudes of the star, and a Bayesian prior. The ADQL query used to obtain the data is given in Appendix A. To prevent the very worst distance errors from affecting our results we excluded stars for which the relative parallax error is large (ϖ/σ_ϖ > 1), and for which the Gaia astrometric quality indicator RUWE > 1.4 (Lindegren et al. 2021).

We analysed the data in galactocentric cylindrical coordinates (R, ϕ, z), where the position and velocity of the Sun with respect to the Galaxy’s centre were taken to be the astropy default values². These are that the Sun is at a cylindrical galactocentric radius R_⊙ = 8.122 kpc, (GRAVITY Collaboration 2018) and z_⊙ = 20.8 pc above the Galactic plane, with velocity v_⊙ = (v_R, ⊙, v_ϕ, ⊙, v_z, ⊙) = (−12.9, −245.6, 7.78) km s⁻¹ (Drimmel & Poggio 2018; Reid & Brunthaler 2004). The Sun’s position in galactocentric azimuth is defined to be exactly ϕ_⊙ = 180°. Since we used a standard right-handed coordinate system, the Milky Way disc’s typical sense of rotation is negative, that is, disc stars are rotating towards lower angles in azimuth. Where we cite studies that have used different conventions for the Sun’s position in ϕ, and for the direction of rotation, the values we quote in this paper are converted into our coordinate system.

We removed contaminants from the sample that come from the Large and the Small Magellanic Clouds (LMC and SMC), and 47 Tucanae. All stars within 12° and 6° respectively of the centres of the LMC and SMC respective on-sky positions³ were removed. For 47 Tucanae we removed stars in the range of R ∈ (6.6, 6.9) kpc and ϕ ∈ (200° ,205° ) with z < −2500 kpc.

The stars in our sample were divided up into bins in galactocentric radii and azimuth. In radii, the bins were 500 pc intervals beginning at R = 4 kpc and extending outwards to the edge of our sample, while the azimuth intervals were 10° wide and extend from ϕ = 120° −240°. We disregarded any bin that contains fewer than 50 stars so that we could determine the average positions and velocities and their uncertainties from a sufficiently large sample. This leaves us with a sample of 26 891 917 stars in 327 bins. A map of the bins included in our sample, coloured by the number of stars in that bin, can be seen in Fig. 1.

Fig. 1. Distribution (left) and average vertical velocity (right) of our sample over the Galactic plane. This is binned, with each bin covering 0.5 kpc in R and 10° in ϕ, both in galactocentric coordinates, with stars in the disc rotating in the negative ϕ direction. We note that moving radially outwards from the position of the Sun ${\overline{v}}_z$ rises to around 10 km s−1 before falling to below zero, but that this behaviour disappears for ϕ ≳ 220°.

For each component of velocity, we used the median as our average value, which we preferred to the mean as it is more robust to outliers and contamination. We derived uncertainties on the average of each velocity component ${\overline{v}}_i$ using non-parametric bootstrapping. For each bin we resampled the velocities 1000 times, and determined the median value in each case. Our statistical uncertainty ${\sigma }_{{\overline{v}}_i,\mathrm{stat}}$ is then the standard deviation of the distribution of resampled medians.

As tests of the robustness of our results, we have looked at data samples with a more conservative cut of 1000 stars per bin, and tests in which we consider stars within limited ranges of R. In this latter case we looked at both 2 kpc annuli in 1 kpc steps with the outer R ranging from 13 to 16 kpc and at samples of stars with R less than a threshold value which ranges from 13 to 18 kpc. We discuss the outputs of these tests when describing our results.

The unprecedented number of stars in Gaia DR3 for which we have full six dimensional phase space measurements allowed us to probe the warp’s kinematic signature as a function of azimuth, and Fig. 1 also shows the average vertical velocity of the sample as a function of position in the Galaxy. It is immediately obvious that the stars in the range 10 ≲ R ≲ 15 kpc are moving upwards across most of the area covered, while those outside ∼15 kpc are moving downwards. We can also see that beyond ϕ ≈ 220° the stars in the range 10 ≲ R ≲ 15 kpc are not moving upwards. As we will discuss through the rest of this paper, these two clear features of the data would appear to indicate that the warp is precessing quickly, and that the line of nodes of the warp, as defined by kinematics, lies closer to ϕ ≈ 135° than the value measured from the positions of stars of around 180°.

3. The warp models

3.1. Warp shape

The mean vertical position of the plane of the disc in our model, z₀, depends on both R and ϕ, and we separated these dependencies out, adopting

$$\begin{array}{c}\hfill z_0(R,\varphi )=h(R)s(\varphi ,R),\end{array}$$(1)

where h(R) sets the greatest displacement due to the warp at a given R, while s(ϕ, R) is a function that takes values in the range −1 ≤ s ≤ 1.

For h(R), numerical experiment led us to use a piecewise-linear model, with the free parameters being the peak heights at 1 kpc steps from R = 7 to 16 kpc. Inside 6 kpc we assume zero warp, while beyond 16 kpc we extrapolate h(R) linearly from the last segment (or set it as constant if the final segment had a negative gradient – our results do not depend on this approximation). We denote these h₇, h₈, …, h₁₆. For our robustness tests where there is a cut in R we use the same h(R) model, but limit it to the radii where the model affects the smaller dataset.

It has been common in studies of the Milky Way’s warp to assume that the warp’s dependence on azimuth, s(ϕ, R) in Eq. (1), is a simple sinusoidal function with no dependence on R (see Drimmel et al. 2000; Cheng et al. 2020, and others). We followed this precedent and adopt the azimuthal dependence

$$\begin{array}{c}\hfill s(\varphi )=sin(\varphi -{\varphi }_w^{\prime }-{\omega }_{\mathrm{p}}t),\end{array}$$(2)

where ${\varphi }_w^{\prime }$ is the line of nodes, and ω_p is the warp’s precession rate. We present results using this model in this study. However, with the Gaia DR3 data we can go further than this, and we therefore investigated the radial change (or “twist”) of the line of nodes.

We therefore extended our model to allow both ${\varphi }_w^{\prime }$ and ω_p to vary linearly with R. Specifically, we took

$$\begin{array}{c}\hfill s(\varphi ,R)=sin(\varphi -{\varphi }_w^{\prime }(R)-{\omega }_{\mathrm{p}}(R)t),\end{array}$$(3)

where

$$\begin{array}{cc}& {\varphi }_w^{\prime }(R)={\varphi }_{w,0}^{\prime }+{\varphi }_{w,1}(R-R_{\mathrm{s}})\hfill \\ \hfill & {\omega }_{\mathrm{p}}(R)={\omega }_{\mathrm{p},0}+{\omega }_{\mathrm{p},1}(R-R_{\mathrm{s}}),\hfill \end{array}$$(4)

where R_s is a fixed at 12 kpc. This is purely a numerical convenience that ensures that ${\varphi }_{w,0}^{\prime }$ and ω_p, 0 took values that are easy to interpret (i.e. the values of the respective parameters at R = 12 kpc). We note that, given our formulation of the problem, typical values of ${\varphi }_w^{\prime }(R)$ correspond to the line of nodes on the far side of the Galactic disc. However, by symmetry the other line of nodes in our model is always at ${\varphi }_w={\varphi }_w^{\prime }+{180}^{°}$, so it is this value that we quote, along with ${\varphi }_{w,0}={\varphi }_{w,0}^{\prime }+{180}^{°}$ for models where the line of nodes varied with R.

3.2. Modelling the kinematic signature

Qualitatively, we can understand that stars moving from the lower part of the warp to the higher part will, today, be observed with a net upward motion. For a warp that varies sinusoidally with ϕ, this velocity will be a maximum at the line of nodes, and zero at the peak or trough of the warp. A first complication to this is that, for a precessing warp, the warp may be catching up with the stars. This can cause stars that appear to be heading towards the higher part of the warp to have a net downward motion. A subtler effect is that stars moving inwards or outwards should also follow the warp shape, and this could also affect the expected velocities of the stars.

We can create a quantitive model for the expected velocities in idealised models using Jeans’ first equation (Jeans 1915), which is an integral of the collisionless Boltzmann equation to give, effectively, the continuity equation. In cylindrical polar coordinates (R, z, ϕ) we can write it as

$$\begin{array}{c}\hfill 0=\frac{\partial n}{\partial t}+\frac{\partial (n{\overline{v}}_{\mathrm{R}})}{\partial R}+\frac{1}{R}\frac{\partial (n{\overline{v}}_{\varphi })}{\partial \varphi }+\frac{\partial (n{\overline{v}}_z)}{\partial z},\end{array}$$(5)

where n is the number density of stars in the disc. We assumed that n follows a double exponential profile for stars in the disc

$$\begin{array}{c}\hfill n(R,z\prime )=f(R)g(z\prime )=n_{0}exp(-\frac{|z^{\prime }|}{z_{\mathrm{h}}}-\frac{R}{R_{\mathrm{h}}}),\end{array}$$(6)

where the scale height z_h and scale length R_h are parameters that ultimately disappear at the end of the model derivation. The term z′ is the distance along the z-axis to the warped midplane of the disc at z₀ (Eq. (1)), and is given by

$$\begin{array}{c}\hfill z^{\prime }=z-z_0=z-h(R)s(\varphi ,R).\end{array}$$(7)

The standard approach to warp modelling in the Milky Way has been to assume a rigidly processing warp, and that the effect of v_R is negligible, yielding (Drimmel et al. 2000; Poggio et al. 2020)

$$\begin{array}{c}\hfill {\overline{v}}_{z,\mathrm{mod}}=(\frac{{\overline{v}}_{\varphi }}{R}-{\omega }_{\mathrm{p}})h(R)\frac{\partial s}{\partial \varphi }·\end{array}$$(8)

However, as we noted above, Cheng et al. (2020) showed how to include the effect of v_R motion in this equation, and we can follow their derivation, with the same assumptions that $\partial {\overline{v}}_{\varphi }/\partial \varphi =0$ and $\partial {\overline{v}}_z/\partial z=0$, to arrive at an expression for ${\overline{v}}_z$

$$\begin{array}{cc}\hfill {\overline{v}}_{z,\mathrm{mod}}=& (\frac{{\overline{v}}_{\varphi }}{R}-{\omega }_{\mathrm{p}})h(R)\frac{\partial s}{\partial \varphi }+{\overline{v}}_{\mathrm{R}}\frac{\partial h(R)}{\partial R}s(\varphi ,R)\hfill \\ \hfill & +{\overline{v}}_{\mathrm{R}}h(R)\frac{\partial s(\varphi ,R)}{\partial R}·\hfill \end{array}$$(9)

This differs from the equation derived by Cheng et al. (2020, their Eq. (9)), by the addition of the final term which includes $\frac{\partial s(\varphi ,R)}{\partial R}$. This new term was required when we introduced the twist to the line of nodes and precession rate as given in Eqs. (3) and (4).

Again, we note that the parameters of the model, such as h(R) and ϕ_w, are properties of a kinematic model which, if the model is correct, correspond to the physical height of the warp and its line of nodes. The model is inevitably incomplete, and there are clear degeneracies between h and ω_p implicit in Eq. (8). We draw a distinction between these model parameters and these properties measured in other studies by referring to the physical warp height and line of nodes as h_phys and ϕ_w, phys, respectively.

3.3. Models summary

The modelling approach described above leaves us with four different models that we fitted to the data. The names we give them through this study are:

1. Classical – Assume that we should neglect the v_R terms, and fix ω_p and ϕ_w as a single values (Eqs. (2) and (8)).

2. Classical twisted – Assume that we should neglect the v_R terms, but allow ω_p and ϕ_w to vary linearly with radius (Eqs. (3) and (8)).

3. v_R term – Include the v_R terms, and fix ω_p and ϕ_w as a single values (Eqs. (2) and (9)).

4. v_R term twisted – Include the v_R terms, and allow ω_p and ϕ_w to vary linearly with radius (Eqs. (3) and (9)).

For all of these models we treat ${\overline{v}}_z$ as the quantity of interest, for which we have measured values and their statistical uncertainties. The statistical uncertainties on ${\overline{v}}_{\varphi }$ and ${\overline{v}}_{\mathrm{R}}$ enter as, effectively, uncertainties in the model ${\overline{v}}_{z,\mathrm{mod}}$ velocity,

$$\begin{array}{c}\hfill {\sigma }_{{\overline{v}}_{z,\mathrm{mod}}}={\sigma }_{{\overline{v}}_{\varphi }}\frac{h(R)}{R}\frac{\partial s}{\partial \varphi }\end{array}$$(10)

and

$$\begin{array}{c}\hfill {\sigma }_{{\overline{v}}_{z,\mathrm{mod}}}={\sigma }_{{\overline{v}}_{\varphi }}\frac{h(R)}{R}\frac{\partial s}{\partial \varphi }+{\sigma }_{{\overline{v}}_{\mathrm{R}}}(\frac{\partial h(R)}{\partial R}s(\varphi ,R)+h(R)\frac{\partial s(\varphi ,R)}{\partial R}),\end{array}$$(11)

for the Classical and v_R term models, respectively. We added these in quadrature to the measurement uncertainty ${\overline{v}}_z$, along with a systematic uncertainty ${\sigma }_{{\overline{v}}_{z,\mathrm{sys}}}=0.7$ km s⁻¹ to give us

$$\begin{array}{c}\hfill {\sigma }_{{\overline{v}}_{z,\mathrm{tot}}}=\sqrt{{\sigma }_{{\overline{v}}_{z,\mathrm{stat}}}^2+{\sigma }_{{\overline{v}}_{z,\mathrm{mod}}}^2+{\sigma }_{{\overline{v}}_{z,\mathrm{sys}}}^2}.\end{array}$$(12)

The systematic uncertainty allows for the imperfections in both the model and the data, and prevents the fit being dominated by the few bins near the Sun that contain over a million stars (see Fig. 1) and therefore have tiny statistical uncertainties. The value of 0.7 km s⁻¹ was chosen such that simple models are, on average, around 1σ from the measurement in each bin. Changing this value by a factor of 2 in either direction did not change our results substantially.

To find the best fitting parameters of the model, and their associated uncertainties, we used a Markov chain Monte Carlo approach using the Goodman & Weare (2010) sampling scheme as applied in the Python package emcee (Foreman-Mackey et al. 2013). The log-probability that we sample was then simply

$$\begin{array}{c}\hfill lnp=-\frac{1}{2}{\displaystyle \sum \frac{{({\overline{v}}_{z,\mathrm{model}}-{\overline{v}}_{z,\mathrm{data}})}^2}{{\sigma }_{{\overline{v}}_{z,\mathrm{tot}}}^2}}+\mathrm{const},\end{array}$$(13)

which is equivalent to taking a flat prior on all of our parameters. We ran 50 walkers for 20 000 steps, discarding the first 20% as burn-in.

4. Fit to the data

The models we used in this study have ten parameters that describe the greatest height of the warp as a function of R (h₇, …, h₁₆) and either two or four which describe the position of the line of nodes and the pattern speed of the warp, with only the two twisted models having four parameters (Eq. (4)).

Appendix B gives the best fitting values of the parameters and their uncertainties in numerical form, while we present them graphically in this main body of the paper. In Figs. 2 and 3 we show the fits to the data. In these figures, as in all figures showing model parameters, we represent the uncertainty by randomly drawing 100 samples of the MCMC chain and plotting them as semi-transparent lines. The uncertainties quoted in Appendix B are taken from the 15.9 and 84.1 percentiles of the MCMC chain after burn-in. Experiments in which we limit our sample to bins with 1000 or more stars give results that are very similar to those shown here, except with larger uncertainties, and these results are shown as an additional table in the appendix.

Fig. 2. Fits to the ${\overline{v}}_z$ data for Classical models. Each panel shows data for 10° sectors of the Galaxy, and plots the average vz in 0.5 kpc bins, with associated uncertainties (statistical and systematic). The model fits with a single value of ϕw and ωp are plotted in black, while those of the twisted model are in light blue. Overall both models provide reasonable fits, without major differences.

Fig. 3. Fits to the ${\overline{v}}_z$ data for the vR term models. As in Fig. 2, each panel shows data for 10° sectors of the Galaxy, and plots the average vz in 0.5 kpc bins, with associated uncertainties (statistical and systematic). Note that the uncertainties are larger here because we include uncertainties propagated from the vR uncertainties. The models with a single value of ϕw and ωp are plotted in orange, while those of the twisted models are in green. Both models provided reasonable fits, but note in particular that model ${\overline{v}}_z$ values are all around zero at larger R.

All four kinematic models capture at least some of the features of the data, with the characteristic rise and fall of ${\overline{v}}_z$ as seen by Cheng et al. (2020) found from ϕ = 120° to ∼210°. For the Classical models, the data for ϕ ≳ 210 values were not perfectly fitted by the model. The model flattens towards zero in the outer radii, while the data is almost all below zero. This is somewhat improved for the twisted model, but not for the bins with the largest values of ϕ. For the v_R term models the model ${\overline{v}}_z$ flatten out to zero in the outer radii in a way that the data does not. The twisted versions of these models provide the best fit of all our models at ϕ > 210, but this comes at the cost of the worst fit of all the models at low ϕ.

In Fig. 4 we show a quantitative comparison of the quality of fit for the four models. We computed a quantity analogous to the chi-squared statistic per degree of freedom, ${\chi }_{\nu }^2$, for each model

$$\begin{array}{c}\hfill {\chi }_{\nu }^2=\frac{1}{N_{\mathrm{bin}}-N_{\mathrm{param}}}{\displaystyle \sum \frac{{({\overline{v}}_{z,\mathrm{model}}-{\overline{v}}_{z,\mathrm{data}})}^2}{{\sigma }_{{\overline{v}}_{z,\mathrm{tot}}}^2}},\end{array}$$(14)

Fig. 4. Quantitative comparison of the four models. We compute a quantity analogous to the chi-squared statistic per degree of freedom, ${\chi }_{\nu }^2$, for each model, and plot the distribution of these values for the final 50 000 steps of the MCMC chain. We note that the Classical models perform significantly better than the vR term.

where N_bin = 327 is the number of bins for our data and the number of parameters N_param = 12 or 14 for the different models. We emphasise that setting the systematic uncertainty that we apply to be 0.7 km s⁻¹ was guided by the premise that this would lead to ${\chi }_{\nu }^2$ ≈ 1 for simple models, so the absolute value here is far less important than the relative values for the different models. The Classical models perform significantly better than the v_R term models. We also see that the twisted models perform better than the single value models, but that the improvement is not as great as the difference between the v_R term and Classical models. These relative results are not changed if we change the systematic uncertainty by a factor of 2 in either direction.

In summary, while all four model provided fits to the kinematic data that seem reasonable to the eye, a quantitative comparison showed that the Classical models perform significantly better than the v_R term ones.

5. Properties of the models

We now look at the properties of the models, bearing in mind that the Classical models provide the better fit to the data.

5.1. Warp amplitude

In Fig. 5 we show the maximum amplitudes of the warp in our models, and how it compares with other measurements of the warp amplitude found in the literature. We divide this into results that were found purely from kinematics (labelled h) and from photometric data and star counts (labelled h_phys). The two Classical models reach ∼3 kpc at R ≈ 14 kpc, and then remain around the same height (with significant uncertainty) further out. The two v_R term models both stay at much smaller heights, reaching only ≲0.3 kpc by R ≈ 13 kpc and then noticeably decreasing further out, down to almost zero in the case of the untwisted model.

Fig. 5. Maximum warp amplitude as a function of radius for our four models (top panel), and for a selection of recent estimates for the older population – rather than in Cepheids or other young populations – in the literature (lower panel). Two of these studies looked at the lopsidedness of the warp and in these cases we give both the warp upwards and downwards, indicated by the two lines for Amôres et al. (2017, with the upward warp being the greater) and the uncertainty bar for the Red Giant Branch (RGB) results of Romero-Gómez et al. (2019, with the downward warp being the greater value).

The literature values we show come from studies which looked at the bulk of the population or specified older populations (López-Corredoira et al. 2002; Amôres et al. 2017; Cheng et al. 2020; Chrobáková et al. 2020; Romero-Gómez et al. 2019; Uppal et al. 2024), as opposed to those which looked at the gas or young populations, such as Cepheid variables, that tend to be less warped than all of these except the Chrobáková et al. (2020) estimate (h typically around 0.5 − 0.8 kpc at R = 14 kpc for Cepheids, e.g. Skowron et al. 2019; Dehnen et al. 2023; Cabrera-Gadea et al. 2024). We pick out specifically the results using red-giant branch (RGB) stars from Romero-Gómez et al. (2019). We chose results for 6 Gyr old stars from Amôres et al. (2017), corresponding to their second-oldest age bin. Their oldest age bin (8.3 Gyr) has significantly larger uncertainties and even steeper warps. Of the results shown, only ours and that from the Cheng et al. (2020) study comes purely from kinematic data, while the other studies used photometric data and star counts to trace the warp shape.

The Classical models have a warp amplitude that lies somewhere in the middle of these models, a stronger warp at ∼14 kpc than in all but the Amôres et al. (2017) case, but flattening out to be fairly similar to both the López-Corredoira et al. (2002) and Cheng et al. (2020). This flattening is not something one would expect for a typical warp, and does come with significant uncertainties. The v_R term models have a much smaller amplitude than all but the Chrobáková et al. (2020) warp, and do not look like a warp in the sense that they do not tend towards greater heights at greater radii.

Our Classical models are close to the results from Cheng et al. (2020), which in one sense is unsurprising, as we both fitted the kinematic data alone. However, one might have expected the results from our v_R term models to be more similar to theirs, since they also used these models in their analysis. The reason for this difference is probably that we analysed the data in several azimuthal bins at each radius, while Cheng et al. (2020) combined data at all azimuths into a single radial bin, thus reducing the effect of the v_R term on their model.

Our robustness tests in which we cut the data at different maximum values of R produced very similar height profiles to those found here for all models, and we show these in the appendix. Our robustness tests in which we look at 2 kpc wide rings in R had slightly low values of h at R = 16 in the Classical models, with large uncertainties. The results from the v_R term models looking at 2 kpc wide rings had greater h than those we see here, but this reflects a weakness of this piecewise-linear model applied to this limited sample, rather than a worrying feature of the data analysis. We discuss this, and show results from all of these robustness tests in Appendix B.

In summary, the Classical models, in addition to providing a better fit to the data, also provided a warp amplitude that is more in line with the literature. The v_R term models, as well as being a worse fit for the data, do not resemble warps as we would expect them. We draw the conclusion that the Classical models are much more likely to be a helpful representation of the warp in the Milky Way.

5.2. Line of nodes

The positions of the line of nodes for each of our four models are shown in Fig. 6. We plot it from R = 6 kpc outwards, as our model warp has zero amplitude inside this radius. In both the cases where the line of nodes is held at a fixed value, it lies around ϕ = 140°. This reflects the clear signature in the data (Fig. 1) that towards ϕ = 230° the vertical velocity is close to zero. The Classical twisted model has a line of nodes that forms a trailing spiral, while the v_R term twisted model’s line of nodes is a leading spiral that is much more tightly wound. In both cases these cross 140° at R ≈ 12 kpc, which is the radius where the velocity signal of the warp appears to be strongest.

Fig. 6. Position of the line of nodes for our Classical (left) and vR term models. The Classical twisted model has a line of nodes that forms a trailing spiral, while the vR term twisted model’s line of nodes is a leading spiral that is much more wound up.

This line of nodes is significantly different from that found in the literature using the current positions of stars. Most early 21st-century studies of the Milky Way’s warp place ϕ_w, phys very close to 180° (e.g. Drimmel et al. 2000; López-Corredoira et al. 2002; Levine et al. 2006) but some more recent work has placed it closer to 165° in both the Cepheid population (e.g. Skowron et al. 2019) and the older population (e.g. Momany et al. 2006; Amôres et al. 2017).

Both Chen et al. (2019) and Dehnen et al. (2023) find a line of nodes for samples of Cepheids that is close to 180° for R = 12 kpc, but then decreases almost linearly in ϕ (a leading spiral) to 140° around R = 16 kpc. This is measured using the positions of the Cepheids by both studies, while Dehnen et al. (2023) find the same result again when they also use the velocities to measure the warp’s properties from the angular momenta of the stars around the Galactic centre⁴. Inside 12 kpc, the warp measured using only the positions of the Cepheids is instead a trailing spiral in both studies. The warp measured by Dehnen et al. (2023) using angular momentum is generally too small inside 12 kpc to draw any conclusions about its line of nodes there.

Perhaps most relevantly for understanding our results, using Gaia DR2 Romero-Gómez et al. (2019) find that, for red giant stars, their estimates of the position of the line of nodes using the positions (ϕ_w, phys ∼ 180° −200°) and kinematics (ϕ_w ∼ 160° −170°) of the stars are rather different from one another. They attributed this to the warp being lopsided, with the line of nodes being offset from the maximum vertical velocity.

The robustness tests we conducted suggest that the trailing spiral of the line of nodes seen in the Classical twisted model is not a secure result from our data. In Fig. 7 we show two examples of results from these tests which illustrate this (more complete figures can be found in Appendix B). When we divided the data into 2 kpc annuli and fitted a Classical model to each, we found that the different lines of nodes appeared to be forming a leading spiral. When we varied the outermost radius of the data we fitted, we found that the twist of the line of nodes swings between trailing for the R < 13 and 14 kpc samples, to leading for the R < 15 and 16 kpc samples, then back to trailing for the R < 17 and 18 kpc samples. We therefore cannot make a strong statement about the twist of the line of nodes.

Fig. 7. Position of the line of nodes for two of our robustness tests. The left panel shows the results of our Classical model when we cut the data into 2 kpc wide rings from 11 − 13 kpc to 14 − 16 kpc. The right panel shows our Classical twisted model when we cut the data at different maximum radii from 13 to 18 kpc. To avoid clutter in the figures we only show the line of nodes for the best fitting model for each sample. It is worth noting that in the right panel the models for data with R < 15 or 16 kpc are so similar that they sit on almost perfectly top of each other. The different trends seen for the different samples make it hard to draw concrete conclusions about the twist of the line of nodes.

We can however say with some confidence that the kinematically estimated line of nodes based on this population is clearly leading that found in the positions of similar samples (or from Cepheids) and that this result is robust against selection of different subsets of our data.

5.3. Precession rate

Figure 8 shows the precession rate of the line of nodes for each of our four models. The Classical models had a precession rate of ω_p ≈ −13.5 km s⁻¹ kpc⁻¹ (prograde with the disc rotation). This would be equivalent of a full rotation around the galaxy in around half a Gyr. The Classical twisted model had negligible change in precession rate with radius.

Fig. 8. Precession of the line of nodes for our four models. Note that the disc is spinning in the negative sense, so the two Classical models are both precessing in the prograde sense, with very similar precession rates. The Classical twisted model has negligible change in precession rate with radius. The two vR term models have very different precession rates with the vR term twisted model having a dramatic change in precession rate with radius.

This relatively fast prograde precession of the warp is consistent with the results from Poggio et al. (2020) and Cheng et al. (2020). It is consistent with results for Cepheids from Dehnen et al. (2023) at R ≈ 12 kpc, but they found that the precession speed decreased to ∼6 km s⁻¹ kpc⁻¹ for R ≳ 13 kpc. Our Classical twisted model allows for this possibility, but does not find it.

The v_R term model had a retrograde precession of around 8 km s⁻¹ kpc⁻¹, while the v_R term twisted model had a dramatic change in precession rate with radius of ∼3 km s⁻¹ kpc⁻², changing sign from prograde at lower radii to retrograde for R ≳ 12 kpc. This seems rather unlikely, and we consider it be a consequence of the fitting algorithm hunting down a mathematically acceptable solution for a model where a physically meaningful one was not available.

Our robustness tests reinforce these results, with the tests in which we set a maximum R value all being consistent with these values. The tests with 2 kpc wide rings in R are the same for the Classical models, but the v_R models become more consistent with the Classical models. As is the case for the h(R) in these models, this reflects a weakness of this piecewise-linear model applied to this limited sample.

5.4. Results summary

The Classical models provided a better fit to the data, with an amplitude and shape consistent with typical warps, and a precession speed consistent with other measures in the literature. Their line of nodes did not lie where one would expect given results in the literature based on the positions of stars. The v_R term models provide a worse fit to the data, with an amplitude and shape that are not consistent with typical warps and, in the twisted case, a precession speed that varies dramatically with radius.

6. Discussion

6.1. The v_R term

Cheng et al. (2020) introduced the idea that one could include the velocity of stars radially in the disc (v_R) in the kind of Jeans modelling of the warp that we have used here. While they had a smaller dataset, and combined stars at all azimuths (ϕ) into a single bin at each radius, we have exploited the larger dataset from Gaia DR3 to look at the vertical velocity of the warp as a function of ϕ as well as R.

It is important to recognise that the magnitude of the v_R terms’ effect on the model vertical velocity (Eq. (9)) is 90° out of phase with the effect of the other term (the only term in the Classical models, Eq. (8)). The Classical models have greatest ${\overline{v}}_{z,\mathrm{mod}}$ at the line of nodes because their value depends on the derivative of warp height with ϕ, which is greatest there. The v_R terms, in contrast, have their maximum effect at the ϕ value where the warp height is maximum, because their value depends on the derivative of warp height with R. This means that the influence of the v_R terms is zero at the line of nodes.

Our results with the Classical models are very similar to those of Cheng et al. (2020), while those with the v_R term models are very different. The most likely explanation for this is that when the Cheng et al. (2020) study combined all azimuths into a single bin at each radius, they treated the stars as all being at the median azimuth of the bin. This will be close to 180° which they assumed was the location of the line of nodes. This means that the v_R terms will have had a small effect in their fits. In contrast, because we have fit the data in bins of ϕ at each radius, the v_R terms have a much greater effect on our fits across the disc.

We should be a little cautious, because as Gaia Collaboration (2023b) point out, v_R estimates are prone to systematic errors across the disc. We can expect that our use of the Bailer-Jones et al. (2021) Bayesian distances has reduced these errors, and since our robustness tests with much more local samples come to similar results, we are confident that including these v_R terms makes the fit to the data worse, even when we have allowed the properties to the warp to vary with radius in a way that, especially for the precession of the warp, seems physically implausible. It seems therefore that including these terms makes our model of the warp significantly worse.

The radial motions of stars are not the only things omitted from our Classical models that might be important for analysing large scale trends due to the warp. We ignore any possible variation of warp amplitude with time, and our variation with ϕ is a simple sinusoid. Cabrera-Gadea et al. (2024) found that for a sample of Cepheids the lopsidedness of the warp (parameterised as an m = 2 mode) was an important factor, while the change in amplitude was not. We defer investigation of these for the bulk population to future work.

As Box (1976) noted: all models are wrong, but some are useful. Our warp models cannot capture the full complexity of the warp, but they can be useful in understanding the warp’s properties. The assumption underlying our models is that the warp is fixed, except for its precession. We therefore require that a star moving radially “knows about” the height of the warp at the radii it is moving towards. The radial motions only exist because the warp is not fixed, and the outer parts of the disc have a complicated structure including spiral arms. Introducing these radial motions into our models is therefore a way of trying to capture the complexity of the warp. However, it seems that this is not a successful way of doing so.

6.2. The inconsistent line of nodes

The line of nodes in our Classical models lies near ϕ = 140° at R = 12 kpc. As discussed in Sect. 5 this is not the same as the line of nodes found by numerous studies of the warp in the Milky Way, whether of position of the bulk population or young populations such as Cepheids, which tend to be closer to 165° −180°.

From this we conclude that the observed structure and kinematics of the bulk population in the outer disc are inconsistent with the picture of a fixed, precessing, warped disc. This stands in contrast to the younger Cepheid population for which Dehnen et al. (2023) found warps defined by the positions or angular momenta of the stars were much more consistent – though not perfectly consistent, as they discuss. We note also that Cabrera-Gadea et al. (2024) found that the line of maximum z velocity for the Cepheid population, which for our model would be the line of nodes, trailed the line of nodes found from the positions of stars by ∼20°, the opposite of the trend we see in the older population.

Our approach within this paper has been to look only at the kinematics of the sample. This is a pragmatic choice, since it avoids dealing with the dust extinction within the Milky Way, but it has great scientific merit to look at the two probes (density distribution and kinematics) separately, and see if they are consistent. Doing so allows us to check our assumptions, and we can see that they do not hold up.

We are not the first to note that the line of nodes determined from kinematics is offset from the Galactic anticentre, with Poggio et al. (2018) noting that the maximum in ${\overline{v}}_z$ was at ϕ < 180° and Romero-Gómez et al. (2019) placing it in the range 160° −170°. Our study places the line of nodes even further from the anti-centre than these studies. This is appears to be because we fitted a model to the kinematics alone, and this increases the significance that ${\overline{v}}_z\approx 0$ at ϕ ≈ 220° even at R ∼ 12 kpc where we see the strongest positive ${\overline{v}}_z$ elsewhere in the disc. This is 90° away from the line of nodes in our warp model.

6.3. The disturbed Milky Way

The Milky Way’s disc is not a particularly calm and settled system. The spiral arms create large disturbances in the radial velocities of stars within the disc (e.g. Gaia Collaboration 2018b). Even before Gaia DR2, disturbance of the vertical velocity of stars across the disc had been and found and associated with “corrugations” or “bending” and “breathing” modes of the disc (Widrow et al. 2014; Schönrich & Dehnen 2018). One of the first discoveries with Gaia DR2 was the “phase-spiral” (Antoja et al. 2018) which is a disturbance seen in the density in the z − v_z plane of stars in the Solar neighbourhood. Further studies have shown that this is how far from the Solar neighbourhood this feature can be found (Xu et al. 2020), and that it grows particularly strong when looking specifically at stars with orbital guiding centre radii around 10 kpc (angular momentum ∼2300 km s⁻¹ kpc, Alinder et al. 2023).

Investigation of the outer disc by Gaia Collaboration (2021) and McMillan et al. (2022) showed that there is a bifurcation of the velocity distribution at radii near 12 kpc, which varies as we look around the disc. This is near the peak of the warp signal in our data.

None of this behaviour is captured by a simple warp model of the kind we have used here. Perhaps, therefore, we should not be surprised by the fact that our kinematic model of the warp is not consistent with the positions of stars in the outer disc. None-the-less, we have quantified this inconsistency and are left with quite a stringent constraint for any model of the creation of the Milky Way’s warp. A range of such models have been proposed in the recent literature, most focussing on the LMC and its effect on the dark-matter halo (e.g. Han et al. 2023), the Sagittarius dwarf (e.g. Poggio et al. 2021) or a combination of both (e.g. Laporte et al. 2018). To match what we observe in the Milky Way, these models need to have warps that appear to be precessing fast in the prograde sense, and the line of nodes in the kinematics of the bulk population must lead the line of nodes in the positions of the bulk population by a substantial amount.

In this sense, we also see a difference between the warp seen in the Cepheid population, which is young, and the bulk population, which is older. The results of Dehnen et al. (2023) and Cabrera-Gadea et al. (2024) suggest that the line of nodes in simple kinematic models of the Cepheid population are roughly consistent with that in position, or that they trail it, while our results have a line of nodes in the kinematic model that leads one found in position. It may be that the Cepheids, a dynamically cold population, had a different response to perturbation than the dynamically hotter bulk population. Perhaps a more likely explanation is that the Cepheids more closely represent the disturbance of the gas disc from which they formed, while the bulk population is dominated by stars that were already in the disc when it was disturbed, and therefore we see the response of the collisionless component of the disc to the disturbance.

The difference between the responses to a Sagittarius-like perturber of the gas and stellar discs in a Milky-Way-like galaxy was investigated by Tepper-García et al. (2022). They refer to the variation in ${\overline{v}}_z$ for stars and gas across the disc as a corrugation, and note that they are initially in phase, but move apart after 500 − 700 Myr. They argue that the difference between the two disturbances may help to determine the age of the perturbation. In this context the difference between the warp seen in the Cepheids and the bulk population may be a sign that the warp is due to a relatively old perturbation of the Milky Way’s disc, but we will defer more quantitive statements about this to future work.

7. Conclusions

In this study we have used Gaia DR3 data to study the vertical motion of stars affected by the Milky Way’s warp. The exquisite quality of the data allowed us to fit a kinematic model to observations spanning a significant fraction of the Milky Way’s disc. We have used this to investigate the position and behaviour of the warp, under the assumption that it is a fixed, precessing, warped disc.

The kinematic signature of the warp in the bulk population is consistent with a warp that is precessing at around 13 km s⁻¹ kpc⁻¹ in a prograde direction with a line of nodes around ϕ = 140°, which is 40° prograde of the Sun’s position. We cannot come to any robust conclusions about how this line of nodes varies with radius, because different subsamples of our data lead us to differing results.

Including the expected effects of radial motions on the continuity-equation-based model of the warp makes the fit to the data worse. We conclude that this is not a helpful addition to models of the warp’s kinematics.

The kinematic signature of the warp, under the assumption that it is fixed and precessing, does not line up with the position of the warp seen in star counts. We are forced to conclude that this model is not sufficient to describe the system, and that this provides a challenge to models of the formation of the warp in the Milky Way. The difference between this behaviour and that of the Cepheid population may provide a clue as to the nature and date of the last major perturbation of the Milky Way’s disc.

Acknowledgments

The authors are grateful to the anonymous referee for suggestions which strengthened this work, to Walter Dehnen for discussing his work on the warp with us, and providing us with a pre-publication copy of his paper. P.M. thanks Eloisa Poggio and Joss Bland-Hawthorn for helpful advice. The authors gratefully acknowledge support from project grants from the Swedish Research Council (Vetenskaprådet, Reg: 2017-03721; 2021-04153). Some of the computations in this project were completed on computing equipment bought with a grant from The Royal Physiographic Society in Lund. 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 research has made use of NASA’s Astrophysics Data System. This paper made use of the following software packages for Python: Numpy (Harris et al. 2020), AstroPy (Astropy Collaboration 2013, 2018, 2022), emcee (Foreman-Mackey et al. 2013), SciPy (Virtanen et al. 2020), Pandas (McKinney 2010), and Matplotlib (Hunter 2007).

References

Appendix A: ADQL query

Our sample is selected from the https://gea.esac.esa.int/archive/ table gaiadr3.gaia_source_lite, cross-matched with its table external.gaiaedr3_distance which contains the Bailer-Jones et al. (2021) distance estimates. We use the following ADQL query: We make the final selections to remove stars associated with 47 Tucanae, on parallax uncertainty and RUWE within our data processing.

SELECT source_id, ra, dec,
parallax, parallax_error, pmra, pmdec,
phot_g_mean_mag, bp_rp, ruwe,
r_med_photogeo, r_lo_photogeo, r_hi_photogeo,
radial_velocity, radial_velocity_error
FROM external.gaiaedr3_distance
JOIN gaiadr3.gaia_source_lite
USING (source_id)
WHERE 0=CONTAINS(POINT("ICRS",ra,dec),
    CIRCLE("ICRS",81.28,-69.78,12))
AND 0=CONTAINS(POINT("ICRS",ra,dec),
    CIRCLE("ICRS",13.1583,-72.8003,6))
AND radial_velocity IS NOT NULL
AND parallax IS NOT NULL

Appendix B: Best fitting parameter values and robustness tests

We present the best fitting parameters from all of our models in Table B.1. These correspond to the values that we illustrate graphically in Figs. 5, 6, and 8.

The results of our robustness test with a minimum of 1000 stars per spatial bin are in Tab. B.2, and it is clear that the results are consistent with the full sample, but with much larger uncertainties in the outer bins and in ω_p for the v_R term models. An increase in uncertainty is, of course, to be expected since there is less data to constrain the model.

Figures B.1 & B.2 show the results of our tests with samples selected by Galactocentric radius, either divided into 2 kpc wide annuli (Fig. B.1) or with an upper limit on R (Fig. B.2). In each case we show h(R), ϕ_w(R) and ω_w(R) for the best fitting model with indicative error bars that are symmetrical around the best fit.

Fig. B.1. Maximum warp heights (h – left group of panels), line of nodes (ϕw – centre group of panels), and warp pattern speed (ωp, right group of panels) as a function of radius for our four models, where the data have been binned into overlapping 2 kpc rings in R, as indicated. Note that the line of nodes appears to be twisting in a prograde direction in all cases.

Fig. B.2. Maximum warp heights (h – left group of panels), line of nodes (ϕw – centre group of panels), and warp pattern speed (ωp, right group of panels) as a function of radius for our four models, where the data are truncated at different radii, as indicated.

The maximum warp heights as a function of radius, h(R), are consistent with the values found for our main sample for the two Classical models, with a slight tendency in some cases for lower warp heights in the outer parts of the warp. These still appear to be consistent within the uncertainties. The two v_R term models are very much consistent for all the samples truncated at different R, as seen in Fig. B.2, but this is not the case for the samples divided into 2 kpc wide annuli, as seen in Fig. B.1 (lower left panels). In this case the warp heights are significantly higher than for v_R term models fitted to the other samples – much closer to the values found for the Classical models. It is noticeable that these fits are inconsistent with one another where they overlap, and that the radial gradient of each is negative. Because we have to start this model at a non-zero height (to allow for normal behaviour in a warp at these radii), the fitting algorithm is able to find solutions with a significant height but never with a positive gradient in h(R). This allows it to find these solutions without ‘paying the price’ of a positive h(R) gradient elsewhere in the disc that would presumably harm the fit with the data. Note that the v_R term models – unlike the Classical models – explicitly include the gradient of the warp height with radius in the model (eq. 9). These models are completely unrealistic physically, and should be thought of as a numerical difficulty introduced by using the piecewise linear model for h(R) for these samples.

The lines of nodes for the two Classical models applied to the different data samples tell an interesting, hard to disentangle, story which we discuss in Sect. 5.2. The samples divided into 2 kpc wide annuli (Fig. B.1) appear to show a leading spiral which joins up neatly and consistently for the twisted model. This is the exact opposite of that found with our full sample or when we set a stricter limit of 1000 stars per bin. We can gain a bit more insight by looking at results for the Classical twisted model when we truncate the sample at different radii (upper row, fourth panel of Fig. B.2, also shown as a polar plot in Fig. 7). Depending on the cut we take, we can either end up with a leading or trailing spiral of the line of nodes. This does not follow a simple trend, being leading if we truncate at R = 13 or 14 kpc, trailing if we truncate at 15 or 16 kpc, and leading again with truncations at larger radii.

These results lead us to conclude that our analysis is not capable of providing a robust answer to the question of the twist of the warp’s line of nodes.

The results for the pattern speed of the warp for the different data samples are, for the Classical models, consistent with our main samples’ values, sometimes with large uncertainties. The pattern speed found for 2 kpc annuli with the v_R terms are, as with the values of h(R), closer to those for the full sample with the Classical models. This is again due to the numerical difficulty introduced by using the piecewise linear model for h(R) for these samples – there is a strong correlation between h(R) and ω_p in the models (eqns. 8 & 9) so a problem with the description of h(R) transfers naturally to a problem with ω_p.

All Tables

All Figures

Fig. 1. Distribution (left) and average vertical velocity (right) of our sample over the Galactic plane. This is binned, with each bin covering 0.5 kpc in R and 10° in ϕ, both in galactocentric coordinates, with stars in the disc rotating in the negative ϕ direction. We note that moving radially outwards from the position of the Sun ${\overline{v}}_z$ rises to around 10 km s−1 before falling to below zero, but that this behaviour disappears for ϕ ≳ 220°.

In the text

Fig. 2. Fits to the ${\overline{v}}_z$ data for Classical models. Each panel shows data for 10° sectors of the Galaxy, and plots the average vz in 0.5 kpc bins, with associated uncertainties (statistical and systematic). The model fits with a single value of ϕw and ωp are plotted in black, while those of the twisted model are in light blue. Overall both models provide reasonable fits, without major differences.

In the text

Fig. 3. Fits to the ${\overline{v}}_z$ data for the vR term models. As in Fig. 2, each panel shows data for 10° sectors of the Galaxy, and plots the average vz in 0.5 kpc bins, with associated uncertainties (statistical and systematic). Note that the uncertainties are larger here because we include uncertainties propagated from the vR uncertainties. The models with a single value of ϕw and ωp are plotted in orange, while those of the twisted models are in green. Both models provided reasonable fits, but note in particular that model ${\overline{v}}_z$ values are all around zero at larger R.

In the text

	Fig. 4. Quantitative comparison of the four models. We compute a quantity analogous to the chi-squared statistic per degree of freedom, ${\chi }_{\nu }^2$, for each model, and plot the distribution of these values for the final 50 000 steps of the MCMC chain. We note that the Classical models perform significantly better than the vR term.
In the text

Fig. 5. Maximum warp amplitude as a function of radius for our four models (top panel), and for a selection of recent estimates for the older population – rather than in Cepheids or other young populations – in the literature (lower panel). Two of these studies looked at the lopsidedness of the warp and in these cases we give both the warp upwards and downwards, indicated by the two lines for Amôres et al. (2017, with the upward warp being the greater) and the uncertainty bar for the Red Giant Branch (RGB) results of Romero-Gómez et al. (2019, with the downward warp being the greater value).

In the text

	Fig. 6. Position of the line of nodes for our Classical (left) and vR term models. The Classical twisted model has a line of nodes that forms a trailing spiral, while the vR term twisted model’s line of nodes is a leading spiral that is much more wound up.
In the text

Fig. 7. Position of the line of nodes for two of our robustness tests. The left panel shows the results of our Classical model when we cut the data into 2 kpc wide rings from 11 − 13 kpc to 14 − 16 kpc. The right panel shows our Classical twisted model when we cut the data at different maximum radii from 13 to 18 kpc. To avoid clutter in the figures we only show the line of nodes for the best fitting model for each sample. It is worth noting that in the right panel the models for data with R < 15 or 16 kpc are so similar that they sit on almost perfectly top of each other. The different trends seen for the different samples make it hard to draw concrete conclusions about the twist of the line of nodes.

In the text

Fig. 8. Precession of the line of nodes for our four models. Note that the disc is spinning in the negative sense, so the two Classical models are both precessing in the prograde sense, with very similar precession rates. The Classical twisted model has negligible change in precession rate with radius. The two vR term models have very different precession rates with the vR term twisted model having a dramatic change in precession rate with radius.

In the text

	Fig. B.1. Maximum warp heights (h – left group of panels), line of nodes (ϕw – centre group of panels), and warp pattern speed (ωp, right group of panels) as a function of radius for our four models, where the data have been binned into overlapping 2 kpc rings in R, as indicated. Note that the line of nodes appears to be twisting in a prograde direction in all cases.
In the text

	Fig. B.2. Maximum warp heights (h – left group of panels), line of nodes (ϕw – centre group of panels), and warp pattern speed (ωp, right group of panels) as a function of radius for our four models, where the data are truncated at different radii, as indicated.
In the text