Structure, kinematics and time evolution of the Galactic Warp from Classical Cepheids

The warp is a well-known undulation of the Milky Way disc. Its structure has been widely studied, but only since Gaia DR2 has it been possible to reveal its kinematic signature beyond the solar neighbourhood. In this work we present an analysis of the warp traced by Classical Cepheids by means of a Fourier decomposition of their height ( 𝑍 ) and, for the first time, of their vertical velocity ( 𝑉 𝑧 ). We find a clear but complex signal that in both variables reveals an asymmetrical warp. In 𝑍 we find the warp to be almost symmetric in amplitude at the disc’s outskirts, with the two extremes never being diametrically opposed at any radius and the line of nodes presenting a twist in the direction of stellar rotation for 𝑅 > 11 kpc. For 𝑉 𝑧 , in addition to the usual 𝑚 = 1 mode, an 𝑚 = 2 mode is needed to represent the kinematic signal of the warp, reflecting its azimuthal asymmetry. The line of maximum vertical velocity is similarly twisted as the line of nodes and trails behind by ≈ 25 ◦ . We develop a new formalism to derive the pattern speed and change in amplitude with time (cid:164) 𝐴 of each Fourier mode at each radius, via a joint analysis of the Fourier decomposition in 𝑍 and 𝑉 𝑧 . By applying it to the Cepheids we find, for the 𝑚 = 1 mode, a constant pattern speed in the direction of stellar rotation of 9 . 2 ± 3 . 1 km/s/kpc, a negligible (cid:164) 𝐴 up to 𝑅 ≈ 14 kpc and a slight increase at larger radii, in agreement with previous works.


INTRODUCTION
The warp is an undulation in a galactic disc that makes its mean vertical height deviate from the mid plane in the outskirts of the galaxy.Between 40 − 50% of edge-on disc galaxies are found to be warped (Sanchez-Saavedra et al. 1990;Reshetnikov & Combes 1998), which implies that warps should be long lived phenomena or the formation mechanism a very recurrent one in the history of galactic discs.The Milky Way is not an exception, having a warp whose structure has been widely studied with different tracers like HI (Levine et al. 2006), dust (Marshall et al. 2006) as well as with different stellar populations (López-Corredoira et al. 2002;Romero-Gómez et al. 2019;Skowron et al. 2019a;Chen et al. 2019;Cheng et al. 2020;Chrobáková et al. 2020;Li et al. 2023).Although the Galactic warp has been known for a long time (Burke 1957), its origin is still puzzling.In order to elucidate the history and formation of the Milky Way's warp, it is important to characterise its main properties, as its structure and kinematics.
Classical Cepheids have proven exceptionally useful in tracing the structure and kinematics of the warp offering several key advantages ★ E-mail: mauro.cabrera@pedeciba.edu.uy to study the Galactic disc (Bobylev 2013;Skowron et al. 2019b,a;Chen et al. 2019).Being very young stars (with ages up to a few hundred million years, e.g.Catelan & Smith 2015) it is expected that they have recently inherited the warped structure of the HI gas where they have formed, while still having cold kinematics (vertical velocity dispersion < 5km s −1 , Chen et al. 2019) making it easier to observe the warp signal as secular dynamics has not had time to 'heat' or disturb it, as it would have for older populations (Binney & Tremaine 2008, Sec. 8.1).Also, belonging to such a young population means there is no contamination from any other Galactic component, e.g. the thick disc or halo, which means they exclusively trace the Galactic thin disc.In addition, Classical Cepheids are well-known standard candles (Leavitt 1908;Leavitt & Pickering 1912), offering extremely precise distance measurements (∼ 3% errors); they can be reliably identified based on their variability, making contamination from other stars negligible (e.g.Jayasinghe et al. 2019;Rimoldini et al. 2019Rimoldini et al. , 2023)); and being luminous (500 < / ⊙ < 20, 000, Catelan & Smith 2015), makes them observable throughout a large extent of the disc even with the optical surveys used to identify them at present (Udalski et al. 2018;Ripepi et al. 2023).Their only disadvantage is that they are relatively scarce, with fewer than 2500 Classical Cepheids in the deepest and most complete catalogues of the Galactic disc to date provided by OGLE-IV (Udalski et al. 2018) and the Third Gaia Data Release (DR3, Gaia Collaboration et al. 2020;Ripepi et al. 2023).For these reasons, Cepheids have been used to study the 3D structure of the warp in more detail than any other stellar population (Skowron et al. 2019b), providing evidence for a twisted line of nodes (Chen et al. 2019;Dehnen et al. 2023) and asymmetry in height between both extremes (Skowron et al. 2019a) similar to the HI warp (Levine et al. 2006).They have also been used to study the kinematics of the warp revealing its characteristic bulk vertical motion in the outskirts of the disc (Skowron et al. 2019a).
To describe the structure of the warp, several studies have shown Fourier series to be of great use due to their versatility to summarise any warp signal if enough modes are considered (e.g.Levine et al. 2006;Chen et al. 2019;Skowron et al. 2019a).These studies have focused on describing the structure of the warp, i.e. its mean height as a function of radius and azimuth.Using a catalogue of Classical Cepheids identified mainly with OGLE-IV and combined with Gaia DR2 astrometry, Skowron et al. (2019b) used a Fourier decomposition with up to 2 modes ( ≤ 2) and a fixed line of nodes (LON) to present the first map of the Galactic warp in the young population covering over half the disc.Chen et al. (2019), using a compilation of optical plus Wide-field Infrared Survey Explorer (WISE, Chen et al. 2018) Cepheid catalogues, studied the azimuthal dependence of the LON with radius finding it does not coincide with the Sun-Galactic Centre direction and that it presents a leading pattern, following Briggs's rules for HI warps in spiral galaxies (Briggs 1990).For the kinematics, Fourier series have been used to characterise the changes in mean vertical velocity (  ) in simulations (Chequers et al. 2018;Poggio et al. 2021), but insofar there have been no Fourier decomposition studies of the warp's kinematic signal with Cepheids (or any tracer) which can reflect and quantify its plausible azimuthal asymmetries and changes with radius.Previous studies with other -olderstellar populations have assumed the kinematic signal to be well represented by an  = 1 mode (Poggio et al. 2020;Cheng et al. 2020;Wang et al. 2020;Chrobáková & López-Corredoira 2021;Dehnen et al. 2023), as expected from a tilted rings model ( = 1 mode), but Romero-Gómez et al. (2019) argue this model is insufficient to explain the more complex kinematic signature they observed with Red Clump stars.
In this work we use a Fourier Decomposition method to study the structure and kinematics of the Galactic warp using Classical Cepheids as tracers.We use the Cepheid catalogue from Skowron et al. (2019b) combined with kinematic data from Gaia DR3 (Gaia Collaboration et al. 2020) to explore the dependence of the amplitudes and the azimuths of the modes as free parameters as a function of radius, which allows us to infer the position of the LON and Line of Maximum Vertical Velocity (LMV  ) for a general warp model that accounts for the lopsidedness of the warp.The new method we present here (Sec.5.2), based on a joint analysis of the Fourier series for  and   , allows us to infer the time evolution of the Fourier components of the warp: i.e. their pattern speed and instantaneous change in amplitude.The inference of the evolutionary terms of the Galactic warp has been tackled recently using different tracers, but mostly under the tilted rings model which assumes a symmetric warp.Poggio et al. (2020) and Cheng et al. (2020) have focused on inferring the pattern speed, while Wang et al. (2020) derived the change in amplitude.These works use general samples of stars with good quality Gaia DR2 and DR3 astrometry and available radial velocities; by not been focused on a specific tracer or having any age constraints, these parameters are representative of the general population of the disc as weighed by its star formation history, i.e. a stellar population of intermediate age (several Gyrs old).The recent work by Dehnen et al. (2023) derives both evolutionary terms for different radii for a sample of Cepheids via a tilted rings model, finding differential rotation and change in inclination of the rings.Our work uses the same stellar population to derive the time evolution parameters of the warp with a completely independent method.
The structure of the present paper is as follows.In Section 2 we present the Fourier decomposition method used to describe the Galactic warp's height and vertical velocity (Sec.2.2) and the method showing how these are combined to derive each mode's pattern speed and amplitude change (Section 2.3).In Section 2.4 we present the inference model used to estimate the warp model's parameters, including main conclusions from the inference validation performed using a mock catalogue.In Sec. 3 we describe the catalogue of Classical Cepheids used in this work.In Section 4 we apply the methods to this sample and summarise our results for the structure and kinematics of the Cepheid's warp (Sec.4.1) and those for the time evolution (Section 4.2).In Section 5 we discuss our results and compare with the previous literature.Our conclusions are summarised in Section 6.

Reference Frame
We begin by describing the coordinate system and reference frame we use throughout this paper.The origin of the reference frame is at the galactic center (GC), fixed with respect to an external inertial frame.Positions can be given in Cartesian, or cylindrical coordinates.The X-axis points from the GC away from the Sun, the Y-axis is parallel to the rotation velocity of the disc at the Sun position and the Z-axis is perpendicular to the Galactic plane forming a right-handed triad.In cylindrical coordinates we use the Galactocentric azimuthal angle  measured from the X-axis toward the Y-axis (i.e.opposite to Galactic rotation).In this coordinate system the Sun is at  ⊙ = 8.277 kpc (GRAVITY Collaboration et al. 2022),  ⊙ = 180 • and  ⊙ = 0.027 kpc (Chen et al. 2001).For velocities we use a Cartesian system whose origin is at rest with the GC and their axes parallel to the directions in which the X-Y-Z axes increase.This is an inertial system and thus does not rotate with the Galaxy, the Sun being along the negative X-axis only at present.This facilitates the kinematical and dynamical descriptions.We assume the Sun has Galactocentric cartesian velocity (  ,   ,   ) = (11.10,232.24, 7.25) km/s (Schönrich et al. 2010;Bovy 2015).

Fourier Decomposition of the Structure and Kinematics
We implement the Fourier decomposition method following Levine et al. (2006) and Chequers et al. (2018).The disc is divided into concentric Galactocentric rings, in each ring the mean behaviour as a function of the azimuth for  and   described by a Fourier sums up to  modes as (1) The amplitudes (   ,   ) and phases (  ,    ) are free parameters, obtained as a function of .In what follows we describe the method only for  (), being analogous for   ().The Fourier representation is flexible enough, even with a small number of modes, to describe many known warp shapes: for example a U-shaped warp will be mainly described by an  = 0 mode with increasing amplitude as a function of radius; an integral or S-shaped warp, will be mainly described by an  = 1 mode with increasing amplitude as a function of radius, and asymmetries will be mainly described by a combination of  = 1 and  ≥ 2 modes.
As we will see in Sec 2.4, it is also convenient to rewrite Eq. 1 in linear form with free parameters   ,   as where the transformation between   ,   and   ,   is given by (4)

Deriving Time Evolution
In this section we present a new formalism to derive the evolutionary terms from the warp, its pattern speed and the change in amplitude of each mode at each radii, disentangled from the motion of the stars.
From now on we denote the star's vertical height and vertical velocity as  and   (lowercase) and the Fourier fits to the warp as  and   (uppercase).
We begin by taking a ring at a radius  and considering a star that has no radial motion and constant angular velocity Ω, that it simply rotates around the Galactic Centre but following the warp of a razor thin disc.These assumptions are reasonable for dynamically cold populations such as the Cepheid stars we use in our analysis.Given that the stars follow the warp's shape, their height () at time  is given by the functional expression of the warp  (, ), which we can express as a Fourier series evaluated at the star's azimuth () as follows (5) We allow the amplitude and phase of each mode to evolve in time because we are interested in determining their instantaneous derivatives   and   .If we take the total derivative of () with respect to time we obtain the vertical velocity   of the star -not the warp-given by As expected, Eq. 6 involves terms regarding the time evolution of the warp (   and   ), and a term regarding the motion of the star due to its own angular velocity ( () = Ω).Now we want to link Eq. 6, which describes the velocity of just one star at azimuth (), to the  () and   () fits from the previous section, which describe the mean motion of all stars in the ring at a given time  0 (today).
In a razor thin disc the height of the disc at an arbitrary azimuth and the position of a star at the same azimuth must exactly coincide.Thus, it follows that the vertical height ( 0 ) and vertical velocity   ( 0 ) of a star at  0 and azimuth ( 0 ) =  0 must coincide with the Fourier fits ( ( 0 ),   ( 0 )) we obtained at that same  0 time.Taking  0 as today,   ( 0 ) =  ,0 and   =  ,0 , the amplitudes and phases obtained from the Fourier fits from Eqs. 1 and 2.
Similarly, the vertical velocity   ( 0 ) of the star must also coincide with the mean vertical velocity obtained from our Fourier fit   ( 0 ), evaluated at the star's azimuth.Setting   ( 0 ) =   ( 0 ) in the left hand side of Eq. 6 and expressing   ( 0 ) as the Fourier fit for   in its linear form (as shown in Eq. 3 for ) we obtain where    and    are the linear amplitudes resulting from the Fourier fits in velocity, calculated from the   and    obtained via Eq. 4. The terms sin( 0 −   ) and cos( 0 −   ) in the right hand side of Eq. 7 can be rewritten as Regrouping the terms as a function of  0 and using the orthogonality of the Fourier modes, we obtain that the amplitudes    ,    from the   fit are related to the amplitudes   and   from the  and the warp evolutionary terms as Solving this linear system of equations for   ( 0 ) and [ ( 0 ) −   ( 0 )] and writing back    ,    in terms of the amplitude and phase (  ,    ), the evolutionary terms of the warp are given by and Assuming that the -th mode has angular velocity   , then setting   =    +  ,0 in Eq. 12, we get the pattern speed for each mode as Therefore, having connected the Fourier fits in  and   at a given radius, Eqs. 13 and 14 describe how each pattern speed and amplitude change in time, allowing a reconstruct the time evolution of the warp as a function of radius.
We leave for a future work the publication of a more general framework that consider an azimuthal dependence not only of the vertical motion of the stars, but also their radial and azimuthal velocity, which would presumably result in a better inference of the time evolution of individual Fourier modes of the warp.1

The Inference
We have used Bayesian Inference to infer the   ,   (or   ,   ) that best describe the mean behaviour of the stars in a given ring when applying the methods from Secs.2.2 and 2.3 to a particular sample.Bayes' theorem (e.g. Sivia 2006) relates the Posterior distribution to the Likelihood (L) and the Prior () probability densities functions (PDF) as where  refers to the data, X to the model parameters and  stands for any other available information.In our case, X is a vector containing the linear amplitudes: For the model parameters we assume uniform priors, with sufficiently and arbitrarily large limits.Assuming that the observations are independent, the likelihood is expressed as the product of the individual likelihood of each single data point   , for which we assume a Gaussian distribution where we take  2  to be the square sum of the uncertainty in the measurement   and the intrinsic dispersion    of the variable, at that ring.This    is introduced to take into account the natural dispersion around the mean value that arises from the dynamics of the Galactic disc; in   , it measures the velocity dispersion and, in , it measures how thick the disc is at that ring.The intrinsic dispersion is not a free parameter in the fit.We estimate it as the mean dispersion in the variable of interest in equally spaced azimuthal bins, weighted by the number of stars in each bin because low number statistics dominate over observational errors.
In our case, because the model is linear in all parameters and we have assumed a uniform prior, the maximum a posteriori (MAP) X 0 coincides with the maximum of L and the posterior is exactly a Gaussian distribution with mean X 0 and covariance matrix  (see e.g.Sec. 1 in Hogg et al. 2010, for a detailed discussion).The posterior PDF can, therefore, be expressed as where X 0 is given by The covariance matrix  is the inverse of A: the matrix that contains in its entries the "projection" of each mode into the other ones (see Eq. A2 in Appendix A) weighted by the dispersion in the data, and the vector p has the "projection" of the data in each mode (see Eq. A4).
Because we use a Fourier series to represent a variable, one would expect the modes to be mutually independent and therefore not correlated.This is not usually the case.When we have discrete measurements, the modes are not mutually orthogonal unless the measurements are equally spaced in azimuth and have the same   .In this special case A is diagonal and, in consequence, the covariance matrix is too.This particular distribution allows the modes to be mutually independent.Naturally, we will never get this configuration from the data itself, but this method shows analytically that data that are more or less uniformly distributed in azimuth are preferred for a Fourier analysis of the whole disc: studies with a sparse and irregular azimuthal coverage will get modes that are not "fundamental", in the sense that they are not describing the modes of the warp itself.The effect gets worse with high-frequency modes  ≥ 2. This should be kept in mind when interpreting results for individual modes, nevertheless, it will not affect our conclusions on the description of the warp as a whole (the sum of the modes) in the regions well sampled by the data.
Finally, the disc is divided in rings such that we get a "continuum" view of how the modes and the warp change with the radius.To do so we take each ring to contain a fixed number  of stars out of the total  tot stars, the first ring starting with the star at the smallest radius.The second ring will start at the radius of the second star and have a width such that it also contains  stars, and so on for subsequent rings.This scheme implies that the rings will have a varying width, depending on the sample's radial distribution.We take the radius associated with each ring as the mean radius of the stars in it.This procedure allows us to have a continuous view, with all rings having the same number of stars  and, therefore, constant stochastic noise.It must be kept in mind, however, that only one out of every  consecutive rings will be independent.Changing the number of stars in each ring changes the smoothened parameters inferred as a function of the radius (the bigger N, the smoother it gets).Also, the change in N moves the mean radius of each ring, the tendency is that a bigger N makes the rings to move inwards (smaller radius), as expected for a density profile that decreases with radius.

Validation with simulations
Here we present our main conclusions about the performance of the methods described in the previous section, assessed by applying them to mock catalogues constructed from test particle simulations.As discussed in detail in Appendix.2, we used a test particle simulation of a warped disc from Romero-Gómez et al. (2019) to create a mock catalogue affected by the Gaia DR3 selection function (SF) and observational errors.A fiducial model, unaffected by the SF or by errors, is used as a baseline for comparison of the results of the Fourier Decomposition.The interested reader may find full details and discussion of these results in Appendix.2.
Our main results on how the SF affects the recovery of the warp as a whole, in different regions of the disc, are summarised here as follows: • For  the best sampled region, the quadrants I and III ( < 0 kpc) is recovered well (differences between the real and the recovered warp are smaller than    ) and the general tendency for all radii is recovered for both series summing up to  = 1 and  = 2.For  > 0 kpc (quadrants I and IV) the SF causes the warp to be exaggerated.This bias is reduced for outer radii as the main mode of the warp ( = 1) becomes greater than the intrinsic dispersion (see Fig. B2).
• For   the recovery is better than , although for the inner disc ( ⪅ 9 kpc) the recovery is poor for  > 0. The recovery in the sampled area is better than in  for both  = 1 and  = 2 (differences are smaller than     in most of the disc area, see Fig. B2).
• Main conclusion: The recovery of the full model (the Fourier sum) in both variables is robust in the well sampled regions for  > 10 kpc, i.e. second and third quadrants.In this region all the warp features are well recovered within the uncertainties given by the Posterior realisations.
We also tested how the individual modes are recovered, our main conclusion being: • The  = 0 mode is well recovered throughout the disc.
• The  = 1 mode tends to be overestimated in amplitude at the inner disc.In the outer regions ( > 10 kpc) its amplitude and phase are well recovered in both variables (i.e  and   ).
• The  = 2 mode captures the asymmetries and is well recovered where there is data, i.e. quadrants II and III and  > 10 kpc, but tends to be underestimated in amplitude and the general trend of the phase is poorly recovered in the whole disc.
• Main conclusion: The uncertainty on the recovery of the individual modes stems from the correlations between the modes which, in turn, appear as a consequence of the imperfect azimuthal coverage.The degeneracies introduced by the correlations mean that different combinations of amplitudes and phases for the individual modes can give the same sum model, in a finite azimuth range.A full azimuthal coverage would break this degeneracy and make the inference on the individual modes unique.The mode less affected by this degeneracies is the  = 1 mode due to its large amplitude.
The intrinsic dispersion in  is well recovered in the outer disc where the warp amplitude is larger than the dispersion.In the inner disc    tends to be off by 10%.For   , we find     is underestimated by 3% without dependency on the radius.
Given these results, we decide to include up to the  = 2 mode in the Fourier fits for this work because it offers the least biased recovery for the region of the disc where the warp is most prominent (i.e.outer radii).Reliable results for the inner region of the disc are limited to || ≳ 90 • , the region least affected by the SF with best coverage, where biases in the recovery are lowest.
We also tested the inference of the time evolution parameters   and   .We concluded that the recovery of the   for  = 2 is unreliable due to the biases and noise.For the  = 1 mode we conclude: •  1 is well recovered within its uncertainties particularly for the outer disc • The recovered  1 tends to be overestimated due to a slight overestimation in  1 , but the mean difference is ≈ 4 km s −1 kpc −1 .In the outer disc ( > 14 kpc) we recover the values of the fiducial model within the uncertainty.

THE CEPHEIDS SAMPLE
We use the catalogue of Milky Way Cepheids from Skowron et al. (2019a).The catalogue contains 2385 Classical (Type I) Cepheids identified mainly with the OGLE survey (for more details see Skowron et al. 2019b,a) with photometric distances computed based on mid-IR photometry from the Wide-field Infrared Survey Explorer (WISE) and the Spitzer Space Telescope, resulting in distance uncertainties of 3% on average.We cross-matched the Cepheid catalogue (at 1 arcsec tolerance) with Gaia DR3 to retrieve proper motions for these stars.Out of the 2381 Cepheids with Gaia proper motions, only 860 stars have radial velocities in DR3.In order to curate a homogeneous catalogue with full velocity information allowing us to compute   , we infer the missing line of sight velocity for all stars in the catalogue by assuming the Cepheid rotation curve from Ablimit et al. (2020) which has a slope of −1.33 km s −1 kpc −1 and takes the value 232 km s −1 at the solar radius2 .
We clean this sample by keeping stars with RUWE < 1.4,   ≤ 0.1 kpc and    ≤ 13 km s −1 .These upper bounds in  and   uncertainties guarantee a significant amount of stars whose uncertainties are at most of the order of    .To avoid clear outliers due to probable contaminants and the Magellanic Clouds we restrict the analysis to stars with || ≤ 2 kpc, |  | ≤ 30 km s −1 and 3 kpc <  < 18 kpc.These are very broad cuts that only remove very few (≈ 3%) clear outliers (5 sigmas) most of them due to the cut in   (only one star is removed for the cut in ).These constraints reduce the sample to a total of  tot = 1997 stars.

RESULTS
Here we present the results obtained by applying the methods described in Sec. 2 to the final sample with  = 2 and  = 200 stars in each ring.To calculate    in both variables we use 8 azimuth bins.The resulting amplitudes and phases as a function of radius for the best fitting (Maximum a Posteriori, MAP) models for  and   are provided in Table 1 and 100 posterior realisations are provided in Table 2. Figures C3 and C4 in the Appendix C show the amplitudes and phase (respectively) of each mode in  and in   as a function of the radius.

Structure and kinematics of the warp
In the following sections we analyse different features of the warp structure and kinematics.We analyse the full Fourier series obtained.Since the validation with simulations indicated results for the individual modes are prone to be biased due to correlations between the modes, we discuss and summarise this in Appendix C for the interested reader.

General structure of the warp
We show in the upper panels of Fig. 1 the results of three fits in  for different Galactocentric radii.Each panel shows, for rings of increasing radius, the Cepheids present in the ring, the best Fourier fit (black curve) and 500 random realisations from the Posterior PDF, the gray curves are fits to 200 bootstrapping realisation.The plots clearly show a growth in amplitude typical of an S-shaped warp, reaching a maximum of ≈ 1.1 kpc in the outskirts of the disc.The effect of the SF is evident, the azimuth range sampled increases with radius.Other features like the change of the warp as a function of  become clear in the second panel ( = 11.0 kpc), where a plateau is noticeable around  = 180 • .The third panel ( = 15 kpc) shows how from  ≈ 60 • to  ≈ 240 • the change in the warp between the extremes resembles a straight line more than a sinusoidal curve corresponding to a pure  = 1 mode would.This feature is correctly reproduced by the model thanks to the  = 2 mode; a simple tilted rings model ( = 1) cannot reflect it.The fits to the bootstrapping realisations shows that for  < 10 kpc the fits are affected by statistical noise as shown in the first panel of Fig. 1 at  = 8 kpc in the first and fourth quadrant of the galactic plane, this became more clear in the residuals plot in the bottom panels.For the outer radii the fits are less sensitive to statistical noise as we see in the second and third   panel where the posterior realisations coincide with the bootstrap realisations.For this reason we focus our analysis on the second a third quadrants.

𝑅 𝐴
Figure 2 shows   as a function of  for the same three rings shown in Fig. 1.The first panel ( = 8.0 kpc) of this figure, as well as in the previous one, shows how the few observed data points in regions most affected by the SF (e.g. ∼ 300 • ) strongly drive the fit in those regions.As discussed in Sec.C, this makes the inference unreliable for the inner disc at  < 10 kpc, except around the azimuth of the solar neighbourhood.Therefore, in what follows we will restrict our analysis to  ≥ 10 kpc.As radius increases (second and third panels) the amplitude of the warp in velocity grows but only mildly, as it is at most of the order of the intrinsic dispersion     ≈ 8 km s −1 even at the outer disc.This is in contrast with , where the amplitude of the warp exceeds the intrinsic dispersion by a factor of ≈ 3 in the outer disc.This low amplitude in comparison with     makes it harder to detect the kinematic signature of the warp, but at the outer disc it is clear there is a complex and asymmetrical behaviour, as seen in the third panel in Fig. 2. The bootstrap realisations for   show the same conclusions as in , but due to the low amplitude of the kinematic signal in the first and fourth quadrant the realisations show a greater dispersion than the posterior, illustrating that due to low number statistics noise is larger.For this reason we will focus the analysis of the kinematic signal to the second and third quadrants.

Asymmetries in height
First, we explore the asymmetries of the warp in height above and below the plane.The left panel of Figure 3 shows the difference between the maximum and minimum height reached by the warp above and below the plane in the North and South Galactic hemispheres respectively.Positive values in this plot, at any given radius, imply that the Northern extreme of the warp deviates more from the Galactic plane than the South.Up to  ≈ 12 kpc the northern extreme is larger than the southern, even within the uncertainties, showing an asymmetrical warp.This asymmetry decreases towards the outer disc, with the warp being almost symmetrical to within the uncertainties (≈ 100 pc) at  ≳ 13.5kpc.We should keep in mind that because of the SF, the extremes of the warp tend to be overestimated in the internal regions.However, a more accurate and reliable mea-surement of this asymmetry is expected at the outskirts of the disc from our validation tests (Sec.2.4.1).
Since we set the phases of each mode free, we can also track the azimuth of each extreme of the warp to explore the azimuthal asymmetry as a function of .The right panel of Fig. 3 explores the azimuthal asymmetry of the extremes of the warp as a function of radius by showing the smallest angular difference in the azimuths of the warp extremes in .In a simple tilted rings model of an S-shaped warp these extremes are always separated 180 • , even if the line of nodes is twisted.The plot clearly shows the extremes of the Cepheid warp are never diametrically opposed.The difference in azimuth starts at its lowest value of ≈ 120 • at  ≈ 10 − 11.5 kpc and increases up to ≈ 145 • at  ≈ 12.5 kpc after which it remains approximately constant.This is a robust feature that cannot be reproduced by an  = 1 warp, reinforcing the need for an  = 2 mode to describe the full warp.

Line of nodes and Line of Maximum 𝑉 𝑧
The overall behaviour of the best fitting (MAP) warp model for the Cepheids is shown in Fig. 4 in a face-on view of the disc with a colour scale indicating the mean height above/below the mid-plane.The line of nodes (from now on LON) and line of maximum vertical velocity (LMV  ) are indicated with the black and green lines respectively.A leading twist (i.e. in the direction of Galactic rotation) in both the LON and LMV  is evident, as well as an offset between the two.Fig. 5 shows the LON and LMV  azimuths (for  < 0) as a function of radius.The figure shows that the azimuth of the LON is well represented by the straight line (in the plane , ) with the parameters presented in Eq. 20, obtained from a fit to data in independent rings with  > 11.The LMV  also follows a linear tendency, well described by an almost constant azimuthal difference of 25.4 • with respect to the LON.

The time evolution of the warp traced by Cepheid
Here we present results for the pattern speed (from Eq. 14) and the change in amplitude with time (from Eq. 13) for the  = 1 mode obtained for the Cepheids.We ignore the  = 0 mode, since its pattern speed is ill-defined and its amplitude change is  0 sin(−  0 ) (this is shown in the right panel of Fig. C3).
Although the Fourier series for  and   have been fit with  = 2, we focus this analysis in the  = 1 mode, because it is the dominant mode of the warp and the recovery of the evolutionary terms for  = 2 are biased and noisy due to selection function effects (as shown in Sec.B1.2).From Sec. 2.4.1 we recall that, for our simulation,  1 and  1 are well-defined for  ≳ 12 kpc where  1 is non-zero and also well-defined (as shown in Sec.C).Therefore, we will restrict this part of the analysis to  ≳ 12 kpc.

Amplitudes
In Fig. 6 we present results for  1 as a function of .In the range  < 14.5 kpc the change in amplitude is negligible, for  > 15 kpc it shows a growing tendency3 reaching a maximum in the external disc of ≈ 5 km s −1 ≈ 5kpc Gyr −1 , this tendency is also present in the results by Dehnen et al. (2023).Based on our validation summarised in Sec.2.4.1 we expect these results to be unbiased over this radial range.

Pattern speed
Assuming the angular velocity Ω from the rotation curve by Ablimit et al. (2020), we obtained the pattern speed for the  = 1 mode from Eq. 14.Because in our reference frame the stars rotate in the direction in which  decreases, the angular velocities in the direction of stellar rotation are negative.To avoid confusion we present the angular velocities with their sign changed.Figure 7 presents, as a function of Galactocentric radius, the pattern speed of the  = 1 mode  1 (black curve), the angular velocity of the rotation curve Ω (red curve), the upper and lower limit given by the measurements by Poggio et al. (2020) (dotted blue lines) and results from Dehnen et al. (2023) (solid cyan line and dots).
We find that  1 decreases for  ≳ 11 kpc and shows a small oscillation for 13 < /kpc < 16.This overall behaviour, both the decrease and the oscillation, are observed by Dehnen et al. ( 2023) but at a slightly different radius.This difference may arise from their use of guiding radius and also because we use a mean radius to represent each ring, which tends to drive the results from the outer to the inner radii.We would need smaller uncertainties to ensure this oscillation is a physical phenomenon in the disc and not an artefact from our fits.However, the fact that it is also observed by Dehnen et al. (2023), with a sample that includes radial velocities, increases our confidence in the result.Our mean value observed for  > 12 kpc is in agreement with the results from previous works on measuring the pattern speed by Poggio et al. (2020) and Cheng et al. (2020), who assumed rigid body rotation for the warp.

Comparison with different warp observations
In Fig. 8 we compare different warp models in the literature to our results for  = 90 • (northern region) and  = 270 • (southern region) for  > 10 kpc.The various works cited here have different azimuthal and/or radial coverage, use different tracers, and have used different methods to fit for the warp.Table 3 summarises this information for the works presented in the figure.We have selected these works in order to compare against other dynamically young tracers like the gas, dust and OB stars.We also include results from a few warp models for dynamically older populations for which the time evolution of the warp has been inferred.Since the warp followed by the older population may differ from that of the young, in Sec.5.2 we will discuss the effect due to the assumed structure on the inference of the time evolution of the warp.
We begin by comparing our results against those from Skowron et al. (2019a), obtained for the same Cepheid sample as used here.Within the uncertainties the two coincide at almost all radii.The Skowron et al. (2019a) model behaves like an average smooth model around our results.The mean difference between both models for the northern region (for  ⪆ 10 kpc) is 0.054 kpc, and for the southern region is 0.043 kpc.This level of agreement is expected because we are using a subset of their sample, the differences being in how we model the warp.Skowron et al. (2019a) model the warp as a Fourier sum with  = 2 (as we do) but assume a constant phase for each mode as a function of radius (    = 0) and a second-degree polynomial for each amplitude (  () =   ( −   ) 2 where   is a constant) as a function of .Under these assumptions the resulting model has the form In consequence, there is a single Fourier sum that expresses the mean azimuthal behaviour of the warp at all  >   which is scaled We also compare our results with those from other warp models obtained for dynamically cold tracers like HI (Levine et al. 2006), Dust (Marshall et al. 2006), OB stars (Li et al. 2023), Cepheid (Chen et al. 2019) and pulsars (Yusifov 2004).Because Cepheids are a young population (< 500 Myr, e.g.Catelan & Smith 2015), they are expected to still retain the warp shape inherited from the gas and its star-forming regions, so the agreement among young tracers is expected.We also show the results from Amôres et al. ( 2017) selected for a young population with an age of 400 Myr compatible with that of Cepheids.In the northern region, within uncertainties, we found excellent agreement with all previous results for young tracers, and a clear disagreement with results from Amôres et al. (2017) inferred from star counts modelling using the Besançon Galactic model.The warp model from pulsars departs the most from ours, with a mean difference of 0.14 kpc (less than the intrinsic dispersion of Cepheids, see Fig. C3).For the HI model we found differences for  < 12 kpc, which may be due to the underestimation by the amplitude fitted to its own results by Levine et al. (2006) between 10 <  < 12. Compared to our results in the southern region, these works tend to underestimate the amplitude of the warp for  ⪆ 13 kpc.The warp traced by pulsars underestimates the height the most, compared to ours, with a maximum difference of 0.42 kpc.These differences may arise due to the symmetry imposed in the models for this radial range.The models from Levine et al. (2006), Chen et al. (2019) and Li et al. (2023) are strictly symmetric in this radial range, in consequence, the asymmetry given by the  = 2 mode between both regions can't be represented.The difference with the model from Marshall et al. (2006) may be due to its radial coverage which does not extend beyond  ∼ 13.
Although the degree of agreement in the southern region is not as good as in the north, its clear that all young tracers follow a similar warp (Chen et al. 2019;Skowron et al. 2019a;Li et al. 2023).The clear exception to this agreement is the result from Amôres et al. (2017).Although the disagreement with the Amôres et al. results in the south is not as strong as in the north, they still found a warp amplitude that is systematically lower than ours as well as all other works for Cepheids and similarly young tracers like dust and HI.
We now focus our attention on the intermediate population: Red Clump stars (López-Corredoira et al. 2002;Wang et al. 2020), A type stars (Ardèvol et al. 2023), K type stars (Cheng et al. 2020) and the full Gaia DR2 population (Chrobáková et al. 2020).Results from López-Corredoira et al. (2002) in the radial range  ⪅ 13 kpc spanned by its observations (thick part of the line) shows agreement with our results and, as with the young populations, the agreement is better for the northern region.However, extrapolating this warp model (thin part of the line) for the outer region of the disc would yield increasing differences that would grow up to the order of a few kpc.Also, the models by Cheng et al. (2020) and Wang et al. (2020) for a 1-3Gyr population are in agreement within uncertainties for the northern region in  ⪅ 12.In the southern region both models are in agreement with our results for  ⪅ 11.5 kpc, after this radius the differences increase up to several kpc in the outer regions.The warp model presented by Chrobáková et al. (2020) is in clear disagreement in both the northern/southern regions with all other warp models using similarly old tracers (like Cheng et al. 2020) and with ours and all other results for young tracers.As we will discuss in Sec.5.2, these differences in amplitudes between the models will become important in the determination of the pattern speed of the warp.Results from Ardèvol et al. (2023) for the kinematics of Atype stars population have shown a clear signal of the warp in the anticentre direction ( = 180), the increasing vertical velocity as a function of the radius from  ≈ 12 kpc, reaching ≈ 6 − 7 km s −1 at  = 14 kpc, similar to our results. 4he issue of the warp's dependence on age of the stellar tracer remains an open question.Older stellar populations like RGB stars, Red Clump stars and other tracers older than Cepheids may trace a similar warp considering the uncertainty in the parameters of each model and their validity range.Also, Cantat-Gaudin et al. (2020) reported that stellar clusters typically older than 1 Gyr trace the southern region of the warp similarly to the Cepheids.Thus, it is unclear whether there are significant discrepancies between the warps traced by older and younger populations.
Among previous results available at present, either the predictions of models with age dependency deviate significantly from the warp observed for bona fide young tracers like the Cepheids, as seen in the case of Amôres et al. (2017), or there is not enough discrepancy in the differences (considering uncertainties) to determine an age dependency for the warp, as in the case of results from Wang et al. (2020) shown in Fig 8 .The uncertainties in the parameters obtained by Wang et al. (2020) for all ages are large enough to allow for the agreement of all models from 1 to 12 Gyrs with our result with Cepheids.In particular the model for 9 Gyrs (not shown), an age completely incompatible with that of Classical Cepheids, is the one in best agreement with our results.Taking into account the restrictions present in some of the models regarding the asymmetry and radial dependence of the warp, each model's validity range in distance and azimuth, and the current precision of the observed warp using different tracers, it remains unclear whether or not there is an age dependency in the warp.

Asymmetries and deviations from the tilted rings model
Our results, as well as several previous ones, showed that a tilted rings model ( = () sin( − ()) or  = 1) does not explain many of the features observed at different radii in position and in kinematics.For , the presence of a plateau at 10 kpc ⪅  ⪅ 11 kpc and  ≈ 180 • shown in the second panel of Fig. 1, where the warp in  is already present, cannot be explained without an  = 2 mode.At that distance the tendency of the disc to warp towards the southern hemisphere is clear at  ≈ 240 • , still far enough in azimuth from the strong obscuration towards the bulge (|| < 90) to be an effect of the selection function.The bootstrapping realisations shown in Fig. 1 and Fig. 2 that the plateau is well recovered, for this reason we consider unlikely to be an artefact of statistical noise.The northern extreme lies in the first quadrant and so its inference is more affected by selection function effects due mainly to obscuration, hence, it is less well constrained than the southern extreme.Nevertheless, the extremes of the warp (in ) are found to be ≈ 120 • apart, while in a tilted rings model this difference must be 180 • by construction.The observed shape resembles the "S-Lopsided" warp model presented by Romero-Gómez et al. (2019).A better azimuthal coverage in the first quadrant and behind the bulge (currently unavailable due to extinction) would provide better constraints for this model.Our result is robust, however, since better coverage can only make the difference between the warp extremes even smaller if the northern extreme lies closer to the bulge.
In the kinematics, a static warp (i.e.  = 0 and   = 0) with a plateau would create a distinctive shape in   .If we consider a star rotating with angular velocity Ω following the shape presented in the second panel of Fig. 1, then, because the star rotates in the direction in which  decreases, after passing the minimum in  ≈ 300 • the star increases its vertical velocity until it reaches the plateau ( ≈ 180 • ) where   ≈ 0, then, on its way to the maximum  close to  ≈ 60 • the star gains   until a certain point after which its   decreases to zero when it reaches the maximum .This creates two maxima in   , one before the plateau and another one after it.A toy representation of a plateau would be  () =  1 sin() +  1 2 sin(2); for a star rotating with angular velocity Ω in a static plateau, this will give   () = Ω 1 (cos() +cos(2)), which shows the geometry described before.This shape is observed in the second panel in Fig. 2. We also take the ratio between the amplitudes of the modes  = 1 and  = 2 in  and   and found consistency with what is expected from the toy model (  1  2 ≈ 2 and  1  2 ≈ 1 around  ≈ 10.5).This peculiar signal was also observed in proper motions by Romero-Gómez et al. (2019) in the RGB population, who interpret it was a signal of the lopsidedness of the warp.As we see here it is actually a characteristic signal of the S-Lopsided model due to its plateau.An indirect evidence of the plateau is also illustrated in Fig. 5 by the large dispersion of the LON for  ≲ 11.5 kpc where the LON is ill-defined.For  > 11.5 kpc the plateau disappears, and the dispersion in Fig. 5 is sharply reduced as the disc is significantly inclined and the LON becomes welldefined.For  > 11.5 kpc other features that differ from a tilted rings model are still present, like the azimuthal asymmetry between the two extremes.The angular difference between them grows but never reaches 180 • , meaning that an  = 2 mode is needed to describe the galactic warp.In consequence the tilted rings (i.e. = 1) model is unable to accurately describe the observed azimuthal location of the warp extremes at any radius.
In Sec.4.1.2we presented our results of the asymmetry between the north and south extremes in the Cepheid's warp.Asymmetry between the height of the warp extremes, or lopsidedness, has also been reported by Chen et al. (2019) and Skowron et al. (2019a) for the Cepheids sample, by Levine et al. (2006) for the HI component and also by Romero-Gómez et al. (2019) for the OB and RGB populations.All these works seem to agree in the existence of an asymmetrical distribution, with the HI as the best exponent of this feature.In our results, the northern extreme is larger by ≈ 0.25 kpc at 11.5 < /kpc < 13 which declines to a mean difference ≈ 0.1 kpc for  > 13.5 kpc as shown in Fig. 3.For comparison the figure also shows the North/South asymmetry for the warp model obtained by Skowron et al. (2019a).This difference behaves like a mean trend of our result as a consequence of the assumption of constant phases for the modes and the polynomial radial dependence of the amplitudes.The observed asymmetry in the outer disc is similar to that found for the OB population at  ≈ 14 kpc by Romero-Gómez et al. ( 2019), but note that Romero-Gómez et al. ( 2019) report an amplitude for the warp traced by OB of 0.3 kpc, much lower than the 0.8 kpc we observe for the Cepheids warp.For the RGB stars, which are older than the Cepheids, Romero-Gómez et al. ( 2019) report a larger asymmetry (red line) but with the opposite sign.This would mean the RGB present a warp with similar amplitude to the Cepheids but larger at the southern extreme.As we show in the following discussion, this may be due to an azimuth dependency of the asymmetry measured.
We have also found an azimuthal dependency in the asymmetry.Fig. 8 shows how the warp at  = 90 • (north) and  = 270 • (south) presents a southern region with a larger departure from the midplane in the outskirts of the disc ( ≈ 15.5 kpc), but in Fig. 3 comparing the north and south extremes we get a northern warp height that is larger at all radii (see also Fig. 4), with a decline in asymmetry towards a minimum almost constant value at the outer disc ( > 13 kpc).
Here, the twisted LON can create misleading interpretations in the measurement of the asymmetry, depending on how this measurement is made.Because the LON is leading and closely centred around  ≈ 180 • (see Fig. 4), a sample which covers the region 90 • <  < 270 • will tend to cover mostly regions below the galactic plane (the mean  in azimuth between 90 • <  < 270 • is below the plane for  > 13 kpc).Hence, comparing  at symmetric azimuths  = 180 • ± Δ, rather than symmetric with respect to the LON, will tend to show a warp with larger amplitude below the plane.In Fig. 9 we show the result of taking differences between the absolute values in height above/below the plane at lines of sight symmetric with respect to the Sun-Anticentre line, i.e. | (180 • − Δ)| and | (180 • + Δ)| for different Δ.This clearly shows how measurements of the asymmetry at azimuths symmetric with respect to the Sun-Anticentre line will yield a different result than when the extremes of the warp are compared.This is a consequence of both the twisted LON and the extremes of the warp never being diametrically opposed (Fig. 3).If the effect introduced by the twist in the LON found with Cepheids is also present in other populations, then different values for the asymmetry may not be enough to ensure different warps for different populations if the azimuthal dependency of the LON is not taken into account.
The results for the RGB sample obtained by Romero-Gómez et al. ( 2019), who reported a warp larger in the south than in the north, are also shown in Fig. 9.The north/south extremes found by Romero-Gómez et al. ( 2019) are roughly symmetric with respect to the Sun-Anticentre line, so comparison in Fig. 9 is appropriate, and this shows their results are consistent with ours for various Δ at their measured distance of  ≈ 14 kpc.Results for the OB population that were shown in Fig. 3 to be in agreement with ours regarding the asymmetry are not shown here because they do not correspond to measurements made symmetric with respect to the Sun-Anticentre line.

Line of nodes and line of maximum 𝑉 𝑧
In Sec.4.1.3we presented results for the LON.As was previously reported by Chen et al. (2019) and more recently also by Dehnen et al. ( 2023), the LON in the warp traced by Cepheids is twisted in the direction of the stellar rotation, meaning a leading LON, as shown in Figs. 4 and 5.This leading LON is in accordance with Briggs's Third Rule for warps (Briggs 1990), which states that warp's LONs are straight for  <   0 and twist for  >   0 , where   0 is the Holmberg radius.Although these rules are derived for the warps traced by HI, they are expected to also apply for warps in the young population.Chen et al. (2019) estimate the Holmberg radius at   0 = 11.4 kpc, its LON and its   0 are plotted in Fig. 5 (cyan curve and dashed vertical line, respectively).We found better agreement between the   0 and our twist's starting radius, than that of Chen et al. (2019), which starts further out in the disc, as shown in Fig. 5.For  ⪆ 12.5 kpc, the LON obtained by Chen et al. (2019) is in quite good agreement with our results.We believe that the difference for  < 12.5 kpc between both works is because, fitting only with one mode, the  = 1 mode in Chen et al. (2019) has to represent the whole warp despite its asymmetries.For this reason, and its relatively low amplitude in  < 12.5 kpc, the  = 1 mode in Chen et al. (2019) behaves as a mean between our LON (the full fit) and the phase of our  = 1 mode (blue dots).The LON twist is also suggested by Romero-Gómez et al. (2019) to be present in the warp traced by RGB stars, but with an opposite direction, i.e. a trailing LON.However, Romero-Gómez et al. (2019) warn that this result may be driven by selection effects due to extinction.
In Sec.4.1.3we show how the LON and the LMV  have a similar twist but they do not overlap, having an almost constant phase offset of 25. • 4 between them.Both lines lie in the region of the disc best populated by our data (as seen in Fig. 4) and best recovered in our tests with simulations from Sec. B1.1.Therefore, we consider both lines to be robust and not affected by biases.The difference in phase between the two lines could be due to the presence of  ≥ 2 modes in the overall warp.Romero-Gómez et al. (2019) also found an offset between the LON and the maximum vertical proper motion for the RGBs and attributed it to the lopsidedness of the warp.According to Romero-Gómez et al. (2019, see their Fig. 8 and Sec. 5.1) the LMV  for the RGBs may lie at  ≈ 160 • − 170 • (they observe  , rather than   ), leading their LON (at  ≈ 180 • − 200 • ) and also ours, but with a twist opposite to our results with Cepheids.Again, this result for the RGBs may be subject to selection effects which may have affected the inference of the LON.
In a warp dominated by the  = 1 warp in both  and   , a change in amplitude with time could be responsible for the out-of-phase LON and LMV  .The phase offset  between the LON and LMV  is given by  1 −   1 + /2, so Eq. 13 translates into directly relating the phase offset with the amplitude change.There are several caveats, however.First, as we have shown, for Galactic Cepheids the  = 1 mode dominates the warp in  but not in   , in which the  = 2 mode has a comparable amplitude at all radii.
Second, an LMV  trailing the LON implies  < 0, and Eq.22 would require  1 < 0 in contradiction with our results and those from Dehnen et al. ( 2023) shown in Fig. 6, which show that  1 ≈ 0 up to  ∼ 14 kpc and  1 > 0 at larger radii.Therefore the evolution of the  = 1 mode alone cannot explain the observed phase offset between the LON and LMV  .We tested whether shifting the disc mid-plane can move the LON to coincide with the LMV  .To do so we would need to shift the stars by −240 pc, the mean vertical height of the stars along the LMV  .This would be too big a shift compared to the typical uncertainties of the position of the Sun above the Galactic plane (of the order of tenths of pc, see e.g.Chen et al. 2001), making this explanation unlikely.
In conclusion, we find the most plausible explanation for the phase offset between the LON and LMV  to be the presence of  > 1 modes which deviate the LMV  from the LON, meaning that lopsidedness would indeed be the main driver of this offset as suggested also by Romero-Gómez et al. (2019) for the RGB sample.Samples with larger azimuthal coverage and also with measured line of sight velocities may help to confirm modes with higher frequencies in  and   , and settle the reason behind this out-of-phase LON and LMV  .

The twist and velocity arcs
In Sec.4.1.3we showed the LON and LMV  are twisted and found these are well represented by straight lines in the plane ,  (Eq.20 for the LON and plus an offset for LMV  ).These parameters are also explored by Dehnen et al. (2023) who present two LONs as a function of guiding radius for the warp traced by Cepheids, obtained from two different methods (mean orbital plane and mean position plane) and find a rate of change in the LON of −14.7±0.7 deg/kpc (mean orbital plane) and −10.6 ± 0.8 deg/kpc (mean orbital position).Our result −12.7 ± 0.3 deg/kpc lies between the two values.
As the disc rotates and the LON wraps up and gets more and more twisted, the disc could appear to have ripples if the LON wraps around more than once around the disc.The rate of change of the LON with  can be associated with the inverse of a radial wavelength of these ripples.If we take a simple tilted rings model ( () ∝ sin(−())) and look at how it changes radially for a constant azimuth, e.g. = 0, then the warp will cross the plane at (  ) =  .Therefore, if the phase is described by  =  + ,  it can be easily associated with a wavelength by  LON = 2  5 .This wavelength is the radial distance between two warp peaks at a constant , if the LON and the amplitudes do not change its behaviour.For our LON we obtain  LON ≈ 28.4 kpc.In the left panel of Fig. 10 we show how this twisted LON creates long arcs in  for different azimuths.
The right panel of Figure 10 shows   as a function of  for different azimuthal cuts, in which the velocity is seen to create arcs whose peak changes in radius for different azimuths.These arcs are explained by the twisted LMV  together with the growth in amplitude in the kinematic signal as a function of .These arcs are a direct consequence of the twist in the LMV  because the kinematic signal does not decline, and also because the peaks of the arcs move outward as phi decreases, as expected for the leading LMV  .Of course the change in amplitude and the asymmetries play a role in the position of the maximum, but the main driver of this arcs is the twist in the LMV  .
These arcs in   have been observed by previous works using Gaia DR2 and DR3 with other stellar tracers (Cheng et al. 2020;Gaia Collaboration et al. 2021).Cheng et al. (2020) pointed out that these arcs in   are a consequence of the pattern speed in a tilted ring model, and indeed an arc can be created with just a constant pattern speed and a growing kinematic signal without a twisted LON.This is because the growing amplitude gets modulated by the factor (Ω − ) so   grows and then starts to decline as the co-rotation radius is achieved where   is null (if  = 0), creating an arc.But this explanation can't take into account the change of the arc shape as a function of the azimuth (as shown in the right panel of Fig. 10), which can only be due to the twisted LMV  , which is a consequence of the twisted LON and the combination of the different evolutionary terms of the warp modes, not only the pattern speed of the  = 1 mode.
Here we have shown these arcs in   are a direct consequence of the twisted LMV  .It is worth noticing that the LMV  does not track the line of maximum in   of the arcs presented in Fig. 10 due to the change in amplitude as a function of .
Because the LMV  has the same twist as the LON,  LON represents also the radial distance between two   peaks at a constant .Using the same  LON for the LMV  we may expect from a extrapolation of this oscillatory behaviour the minimum in   at the anti centre direction to be around  ≈ 28 kpc and the point of null   to be around  ≈ 21 kpc.Wang et al. (2023) used Gaia DR3 to map the disc population out to  ≈ 23 kpc.In their Fig. 3 the disc's vertical velocity goes from positive to negative values at  ≳ 20 kpc.These results seem to support our prediction, assuming  LON holds for the entire Gaia DR3 sample used by Wang et al. (2023).Future extended maps of   may prove helpful to explore whether this analysis also holds for  > 20 kpc and for other stellar populations.
Finally, Poggio et al. (2021) have shown with an N-body simulation of the Milky Way affected by the Sagittarius dwarf galaxy, that the  = 1 mode has prograde rotation if the Milky Way disc and Sagittarius are close to an interaction.After approximately a few Myr the prograde motion coherent with the  = 2 mode disappears.Perhaps the coherent rotation between the  = 1 and  = 2 modes close to an interaction found by Poggio et al. (2021) is the reason why the azimuth of the LON and the LMV  are well approximated by a monotonic (linear) dependence as a function of radius; without a coherent movement in the outskirts of the disc the LON may behave more erratically than is observed.

Time evolution
In Sec.4.2 we provided a new formalism to derive the time change of each mode's amplitude (   ) and pattern speed (  ) at each ring, from the Fourier decomposition of its  and   .This new formalism is free from assumptions on how the amplitude and phase of each mode depends on the radius.By applying it to the Cepheids we derived the pattern speed and change in the amplitude of the  = 1 mode for  > 12 kpc.The dominant mode in  for the warp is the  = 1 mode as expected for an S-type warp, so its evolution may drive most of the time evolution of the warp.In Fig. 7 we show that, within uncertainties,  1 shows a mean rotation of 9.2 ± 3.1 km/s/kpc6 and its error corresponds to the standard deviation of posterior realisations from the independent rings, in agreement with previous reports from Poggio et al. (2020) and Cheng et al. (2020), presented in Table. 4. Some oscillations are present, similar to the results obtained by Dehnen et al. (2023), the main difference being that for  < 12 kpc they found differential rotation slightly larger than we do, perhaps as a consequence of our use of the galactocentric radius as opposed to their use of the guiding centre.We note, however, our results from simulations (Appendix B) suggest  1 may be overestimated in this radial range.
Chrobáková & López-Corredoira ( 2021) present arguments about how the over-estimation in the amplitude of the warp leads to a overestimation in its pattern speed7 , therefore getting a lower amplitude of the warp will translate into a slower precession.This is well reflected by our Eq.14.However, the very low amplitude of the warp presented by Equation 14 also makes it clear why our results for  1 are similar to those of Poggio et al. (2020) and Cheng et al. (2020), despite different assumptions in the three warp models, like the amplitude or the fixed phase  1 .This equation shows that  1 depends on the difference between the phase of the mode in  and in   , therefore, it doesn't matter where they are located or if they are twisted, as long as the phase difference is the same.Because in the Milky Way the assumption that  1 −   1 ≈ −  2 seems to hold at least up to  ≈ 14 kpc, which is the same as assuming  1 ≈ 0, independently of which  1 the model adopts or if it is twisted or fixed  1 will not be influenced by this assumption as long the model adopts  1 ≈ 0. Also, the assumed amplitude affects  1 as was previously mentioned, but  1 gets saturated by over-estimations in  1 , because as  1 → ∞ then  1 → Ω (because the kinematic signature  1 makes  1 < Ω).This may be the reason why Cheng et al. (2020), with its larger amplitude, gets a larger  1 than ours, and also why Poggio et al. (2020) with the same kinematic signature gets a larger  1 , as it uses larger amplitudes.In this analysis of  1 we have left  1 constant because the kinematic amplitude of the warp seems similar for different tracers as shown by Gaia Collaboration et al. (2021).These could be the reasons why Poggio et al. (2020) and Cheng et al. (2020) get similar results to ours, even when they do not consider a twisted  1 and when their amplitudes are larger than ours.We should add that the difference in amplitude is not the only parameter that plays a role in this analysis, the rotation curve and the kinematic signal are not the same between the works cited and they can change the pattern speed measurements, so we expect that these differences to also play a role.
Previous works on the time evolution of the warp neglect the contribution by the change in amplitude to the warp's kinematics (Poggio et al. 2020;Cheng et al. 2020).Poggio et al. (2020) argue that the effect of  1 may be a second order effect in the kinematics.Our results shows empirically that the change in amplitude can be neglected at least up to  ≈ 15 kpc.Wang et al. (2020), finds the change in amplitude derived from the young population (≈ 1 Gy) to be null, which within uncertainties is consistent with our mean measurement up to the radial limit to which Wang et al. (2020) restricted its sample, i.e.  = 14 kpc.For  > 14 kpc we found  1 > 0, reaching a maximum  1 ≈ 5 km s −1 , this tendency is also observed in the change of the inclination in the tilted rings model by Dehnen et al. ( 2023) with similar values.
The prograde rotation of the  = 1 mode found with Cepheids is expected in the context of a disc embedded in a prolate halo as shown by Ideta et al. (2000) and Jeon et al. (2009).However, if this were the case, the prograde motion should be much slower (0.1 km/s/kpc to 1.5 km/s/kpc) than our result.
Although the  = 1 mode rotates almost rigidly, this does not guarantee a rigid rotation of the LON, because the  = 2 mode also plays a role in the LON evolution, and in   its amplitude is comparable to that of the  = 1 mode.Due to the poor recovery expected for the  = 2 mode (Sec.B1.2), a derivation of  2 and  2 with our data would be biased, so we cannot ensure the evolution of the LON or of the whole warp to be one with rigid rotation.The  = 1 mode also presents a growing amplitude for  > 15 kpc, as is also reported by Dehnen et al. (2023).For  < 14.5 kpc the changes in amplitude are insignificant within the uncertainties, therefore, we present a warp which, at first order, shows a stable behaviour for  < 14.5 kpc but still evolving in the outskirts of the disc.
In our derivation of   and   we have ignored the radial velocity and azimuthal changes in Ω. Considering a radial motion of 10 km/s even when radial velocities may seem to be slower (Cheng et al. 2020), we found that  1 may change by about 1 km/s and  1 by 2 km/s/kpc, which are in the order of the uncertainties.Also, that the radial bulk motion reported by Cheng et al. (2020) is inwards for  ⪆ 14 kpc, will translate into a decrease in the measurement of  1 unless it is considered.Therefore, the growth in amplitude for  ⪆ 14 kpc cannot be reduced by taking the radial motion into consideration (in fact, it should increase).These changes are smaller than the uncertainties in the results presented in this work, therefore we do not take them into account in our analysis.These features could be added to the analysis by considering a field of radial velocity and Ω described by Fourier sums at different radii.The extension of our formalism to account for the radial component will be presented in a future work.

CONCLUSIONS
In this work we have used the Skowron et al. (2019b) catalogue of Classical Cepheids to study the structure and kinematics of the Milky Way warp by means of Fourier Decomposition methods.These are the first results presented in the literature for the Fourier Decomposition of the warp in   .Our main results regarding the structure and kinematics of the warp are the following: • The warp is clearly lopsided, both in  and   .In , the amplitudes of the  = 1 and  = 2 modes are comparable up to  ∼ 13 kpc.At larger radii the  = 1 mode dominates, as found previously by Chen et al. (2019); Skowron et al. (2019b).In   , the amplitudes of the  = 1 and  = 2 modes are comparable at all radii.The  = 0 mode does not play a major role in the overall warp shape, we detect a bowl-like shape in the radial range 11.5 < /kpc < 13 with a maximum amplitude of ≈ 200 pc.In   the  = 0 mode is almost null for  > 10 kpc.
• The warp presents a plateau at 10 < /kpc < 11.The observed shape resembles that of the S-lopsided model from Romero-Gómez et al. (2019).The double peak observed in   at this radius is a kinematic signal associated with this plateau.It has also been observed in the proper motions of Red Clump stars by Romero-Gómez et al. (2019).
• The warp is clearly asymmetric up to  ∼ 13, with a Northern warp larger than the Southern warp.In the outer disc ( ≳ 13.5 kpc) the warp becomes symmetric to within uncertainties.
• The extremes of the Cepheid warp in  are never diametrically opposed.The difference in azimuth between the warp extremes is ∼ 120 • at  ∼ 10−11.5 kpc and increases up to 140 • at  ≈ 12.5 kpc, remaining constant at larger radii.
• The LON begins to twist at around  ≈ 11, which is close to the Holmberg radius for the Milky Way (11.4 kpc, Chen et al. 2019), in agreement with Briggs' rules (Briggs 1990).The LON's azimuth follows a linear relationship with radius, presented in Equation 20.We found a twist of −12.7 ± 0.3 • kpc .
• The LMV  does not coincide with the LON, but trails behind it with a constant offset of 25. • 4. We rule out that this offset is due to the change in amplitude with time of the  = 1 mode, and explain this offset as a consequence of the lopsidedness also present in the kinematics.
• The arcs in   as function of  observed in other stellar populations (Cheng et al. 2020;Gaia Collaboration et al. 2021) are also present in the Cepheids sample.We show these are a consequence of the twisted LMV  (see Figure 10).
We have also introduced a new formalism (Section 4.2), based on the joint analysis of the Fourier series in  and   , from which the pattern speed and instantaneous change in amplitude for each individual Fourier mode can be derived.By applying this formalism to the Fourier Decomposition obtained for the Cepheids in  and   , we derive the pattern speed and amplitude change of the  = 1 mode as a function of radius.Our main results are as follows: • The  = 1 mode shows a prograde differential rotation for 11 < (kpc) < 13 with  1 going from ∼ −20km s −1 kpc −1 at  ∼ 10 − 11 kpc to −9.18 km s −1 kpc −1 at  ∼ 13.Our results from simulations, however, suggest  1 may be overestimated in 11 < (kpc) < 13 this radial range.
Thanks to the precise measurements from Gaia DR3 and distances from Skowron et al. (2019a) to its sample of Cepheids, we can explore the complex signal of the warp in both its structure and kinematics.Future Cepheid samples with increased coverage in the first and fourth quadrants will contribute to better restrict the parameters of the warp.A better understanding of the warp kinematics is necessary to make more robust comparisons with simulations and with analytical models of its dynamics, which can lead to better constraints on the possible history of the warp and its role in the evolution of the Milky Way disc's dynamics.Furthermore, the complexity revealed may not be unique to the Galactic warp, understanding it will help also understand warps in external galaxies.
affected by the SF.All differences between the GT and SF models are much smaller than the corresponding velocity dispersion, which has a mean of 19.2 km/s throughout the disc.For both  and   the reduced chi square  2   shows that the GT model fits for  and   are good ( 2  ≈ 1 ∀ ).

B1.1 Azimuthal and radial biases
Figure B1 hinted the existence of regions in which the reconstruction of the warp given by the SF model lacks accuracy.We argue this is due to the correlation between modes introduced by not having a uniformly distributed sample in azimuth and by stochastic clumps in regions with fewer stars in the sample due to the SF.To illustrate this, Figure B2 shows, in each panel, a residual plot between the SF and the GT model in  and   (respectively top and bottom) normalised by the intrinsic dispersion given by the GT model (with fixed  = 2), for SF model fits with up to 1 (left) and 2 modes (right).Grey dots show the SF sample, the black star shows the Sun's position, the inner ring is  = 10 kpc where the warp begins and the outer ring is the radius in which the amplitude of  = 1 mode is bigger than the intrinsic dispersion in the variable.The red and blue colours correspond to over/under estimations by the SF model respectively.
In , for both  = 1 and  = 2 the differences are greater at  > 0 kpc in the inner region before the warp begins.The discrepancy is larger for  = 2 because when fitting with a higher number of modes, in the areas most affected by the SF the higher frequency modes tend to drive the fit towards the few data points available introducing spurious oscillations where there is less data.For both  = 1 and  = 2 the recovery is best for outer radii, where the warp amplitude is larger than the intrinsic dispersion.For  = 1, the differences start growing with radius due to the simulated warp's asymmetry which is not well represented by the Fourier series with  ≤ 1 modes, generating an  = 2 pattern in the differences.By contrast, the asymmetry is better captured by the series for  = 2 for which the discrepancies in the outer region are smaller.However, a hint of the  = 2 pattern in the differences still remains; this is due to a lower amplitude of the  = 2 mode recovered by the SF model (this is illustrated in left panel of Fig. B3).
For   we analyse the bottom panels in Fig. B2, the left plot for  = 1 and the right for  = 2.The differences between the SF and GT models are always smaller for   than the intrinsic dispersion in the whole disc, both for  = 1 and  = 2, in contrast with the recovery in  where differences exceed the intrinsic dispersion in the inner region.The best and worst recovery for   are found in the same regions as for  because the azimuthal distribution is the same for both samples; with the best recovery at negative , and the worst in the internal disc at positive .Finally, the differences between the SF and GT models are much lower in   than in .As also discussed in Romero-Gómez et al. (2019), this is expected because the SF creates exclusion zones in  due to high extinction near the Galactic plane, but doesn't in   because the correlation between  and   is weak for a given star.
Given these results, we decide to use  = 2 for Fourier fits for this work because it offers the least biased recovery for the region of the disc where the warp is most prominent (i.e.outer radii).Reliable results for the inner region of the disc are limited to 90 ≳  ≳ 270 kpc, the region least affected by the SF with best coverage, where biases in the recovery are lowest.

B1.2 Recovery of individual modes
So far we have analysed the recovery of the shape and kinematics of the warp as a whole, given by the sum of the  individual modes in the Fourier series.Now we will analyse how well each mode is recovered.
Each mode  is characterised by its amplitude   and its phase   in , and in   with   and    .In Fig. B3 we compare the amplitudes for  (left) and   (right) as a function of  recovered for the SF model (dark solid lines) against the values given by the GT model (pale solid lines) for each mode.The intrinsic dispersion as a function of radius is also plotted in each panel.
The left panel of Fig. B3 shows how for inner radii ( < 10 kpc) the disc is flat before the onset of the warp, as shown by the near zero amplitudes for all modes in the GT model.Particularly for  = 1, 2, the SF model finds non-zero amplitudes of the order of the intrinsic dispersion.Amplitudes are overestimated in the inner disc because the modes make the full Fourier series flat in the region less affected by the SF ( ≈ 180 • ), but it also tries to fit stochastic clumps far from the midplane at  ≈ 0 • where the SF has removed stars preferentially in the disc plane.At the outer parts of the disc, the  = 1 mode is overestimated by the SF model but the bias is reduced at the external part of the warp ( > 13kpc) where the  1 amplitudes are larger.The  = 2 mode is overestimated due to correlations with other modes when the whole fit of the series is driven by stochastic clumps at  ≲ 10 kpc, as for the  = 1 mode.The  = 0 mode is well recovered over the whole disc.
Some features observed for the amplitudes in  are present also in   .For example the amplitudes are not 0 km s −1 for  <  1 due to the sparse azimuthal coverage caused by the SF.For   , the amplitude of  = 1 is underestimated but the general trend is well recovered by the SF model for  >  1 = 10 kpc as in .The amplitudes of  = 0, 2 have differences between the SF and the GT model, also as in , which is expected because both amplitudes in the GT model are smaller than    , which makes them harder for the SF model to recover.
Similarly to Figure B3, in Figure B4 we compare how the phase of each mode in  (left) and   (right) is recovered by the SF model as a function of radius.We do not plot the phase for  = 0 because it can only take two possible values (−90 • and 90 • ).
For the inner disc at radii  ≲  1 = 10 kpc before the onset of the warp, it is normal that the phase is badly recovered for all modes because the (true) amplitudes are near zero at these radii and the phase becomes meaningless.For the  = 1 mode the phase for both  and   are very well recovered, with no significant bias, for  ≳ 12 kpc where  1 >    .For  = 2, the general trends are recovered for  ≳ 12 kpc, e.g. the twist in  and   showing the change of phase as a function of radius.However, we must be cautious in any particular analysis of  = 2 as an individual mode due to the lack of recovery by the SF model with this mode, its phase recovers some of its tendency but without accuracy.

B1.3 Intrinsic dispersion
Finally, we analyse the bias introduced by the SF to the intrinsic dispersion that our method calculates.To do so, we compute the fractional difference between the    obtained with the GT and SF samples.These differences for  (blue curve) and   (green curve) are plotted in Fig. B5 as a function of radius.The black vertical dotted line at 10 kpc indicates the beginning of the warp, the blue one when the mode  = 1 for  starts to be greater than    , the green one is the same as the blue but for   .First, for  the recovered    is increasingly overestimated at inner radii until the warp becomes greater than the disc's thickness; for larger radii, the recovered    decreases and is off just by 10% of the GT value.Both effects are due to the combination of the increased warp amplitude and the SF.The SF makes the stars near the plane very unlikely to be observed due to high extinction, while the stars away from the plane are less affected by it; since these stars are further away from the disc plane (because of the amplitude of the warp) this tends to inflate    for .This effect is expected to be smaller for a dynamically colder stellar population like Cepheids.For   , on the other hand, we find a mean underestimation of 3%, much smaller than for .We find the appearance of the warp signal in   has no effect in the ability to recover    .Overall the recovery of the intrinsic dispersion affects the inference on the amplitudes and phases in terms only of the dispersion of the posterior PDF, it does not introduce any systematic biases in the parameters themselves.

B1.4 Assumptions on the rotation curve
We tested how the assumed rotation curve may affect the inference in the simulations and with the real data.We did not find any systematic bias in the amplitudes, phases and intrinsic dispersion inferred from the mock catalogues when we used the   derived from the rotation curve, even when using different rotation curves.
In the case of the Cepheids, we tested whether changing the rotation curve offset by ±10 km s −1 could change our main results.We found that different offsets change the amplitudes of the   arcs by ∼ 1km s −1 but do not change the general trend of the kinematic signal of the warp.The changes in  1 and  1 due to changes in the rotation curve are insignificant in comparison with the uncertainties.

B2 Time Evolution
In this Section, we validate the inference of   and   by applying the formalism developed in Sec.2.3 to the simulated sample affected by the SF and comparing it to results for the GT model (as in Sec.B1.2).By doing this we're assuming that the formalism developed holds and will yield correct results for the GT model.Since the test particle simulation we are using has a fixed warp, we expect from this test to recover a constant amplitude and null pattern speed in the region at equilibrium with the potential (i.e. <  2 = 14 kpc).In the outer parts the warp would be expected to evolve with time as the stars relax in the potential.Because the warp model used in the test particle simulation is not constructed by definition as a Fourier series it is not straightforward to use this data to test the recovery of specific values of the time evolution parameters.More involved tests in this direction could be done in a future work to validate the method.
In what follows we analyse the difference between the parameters from both models for  > 10 kpc where the warp is present.We apply this formalism only to the  = 1 mode due to the bias and noisy recovery in the  = 2 mode parameters discussed in Sec.B1.2.
Figure B6 shows the results for the GT and SF models for  1 as a function of radius.For the GT model we get  1 = 0 for  > 12 kpc (black curve in Fig. B6).The variations observed in  1 at 10 < /kpc < 12 are expected in this region were the amplitude of the mode is still very low and it's pattern speed ill-defined.As the amplitude of  = 1 mode increases the pattern speed recovered for the SF model converges to results for the GT at the outermost radii.The mean overestimation in  1 for  > 12 kpc is of the order of 4 km s −1 /kpc, which is within the uncertainties given by the posterior realisations (gray dots).
Figure B7 shows the result for  1 for the GT (black curve) and the SF models (blue curve).The difference between the two for  < 12 kpc is due to the poor recovery in   1 as shown in Fig. B4.The mean difference for  < 12 kpc between the recovery with the SF and GT models is less than 2 km s −1 , which is within the uncertainties given by the posterior realisations (gray dots).For 12 < /kpc < 17, the general tendency for  1 is recovered within the uncertainty with not appreciable bias.The relatively large uncertainty in the recovery on  1 stems from small differences in  1 −   1 , which near /2 translate in large differences in  1 due to it dependence on the difference via a cosine function (Eq.13).The opposite happens for   because it depends on the difference via a sine function.

APPENDIX C: GOODNESS OF FIT AND RESULTS FOR INDIVIDUAL MODES
This appendix presents the results for the goodness of fits and summarises the results of Sec.4.1 for the individual modes.

C1 Goodness of fit
We have tested with the reduced Chi-square how meany mode where needed to do the fits in  and   for the Cepheids sample.Fig. C1 and Fig. C2 shows as a function of the radius the results for  and   .Clearly the result favoured the fits for  = 2 for both variables, showing the need of the  = 2 mode to reflect the asymmetries present in the warp.We have also computed the Bayesian information criteria (BIC, Ivezić et al. (2014)) for different radii and we found that  = 2 is always clearly the best model for  at all radii > 10 kpc.This is of special importance since  is more sensitive to biases due to the selection function problems.For   the fits with  = 2 is also the best case for the outer disc where the amplitude of the warp is significant.We have therefore chose the  = 2 model for both variables.

C2.1 Fits in 𝑧
In Fig. C3 (left panel), we present the results of the amplitudes for each mode and the intrinsic dispersion in  as a function of radius.Clearly the  = 1 mode (red) dominates the fit (it has a maximum of ≈ 1.1 kpc), as expected from an almost S-type like the Milky Way warp.The main mode that takes into account the asymmetries is  = 2 (violet), its amplitude begins to grow at  ≈ 10 kpc but never exceeds 250 pc.For  = 0 (yellow) we have a maximum of ≈ 200 pc.This mode can give asymmetry between both extremes of the warp, but its main purpose is to set the mean height in each ring, so it has the ability to represent radial ripples with no azimuthal dependence.For comparison, we plot the amplitudes for each mode from Skowron et al. (2019a) (dotted curves) obtained with exactly the same Cepheid sample but under the assumption of a monotonic dependency of   with  2 .For the  = 1 mode at  > 10 kpc both amplitudes are practically the same; for the other modes the amplitudes obtained by Skowron et al. (2019a) are similar to the mean behaviour of our results.
The wavy pattern in the amplitudes for  < 10 kpc should not be fully taken as real corrugations in the modes.In Sec.B1.2 we concluded that stochastic clumps in the  −  plane due to the SF generate correlations between the modes.This wavy pattern in  1 is removed if we take  = 1, so the wavy pattern is mainly due to correlations between  = 1 and  = 2.
In Fig. C4 we present the phases of the  = 1, 2 modes for  (top right for  = 1 and top left for  = 2) as a function of .First, lets consider  1 , our results and those of Skowron et al. (2019a) (dotted line) coincide in their general trends for the external region of the disc (> 10 kpc).For  = 1, a twist in the direction of the galactic rotation is well defined, beginning at  ≈ 13 kpc.For the internal region both phases are difficult to determine due to the low amplitude of the warp and because the azimuthal coverage is affected by the SF.For  2 there is more uncertainty than for  1 because  = 1 is better defined and dominates the warp.Within its uncertainty  2 agrees with the phase obtained by Skowron et al. (2019a) (red dotted line).For  > 10 kpc the phases, like the amplitudes, are better behaved than in the internal disc as we expected from Sec. B1.2.
Finally, given that we calculate the intrinsic dispersion for  in each ring, we can see how the disc traced by Cepheids becomes thicker at larger radius, as its shown with the black curve in the left panel of Fig. C3.This shows how the flare in this young population starts at around  ≈ 8 kpc with a height ≈ 100 pc to end up at a height ≈ 390 pc at  ≈ 15 kpc.Previous measurements on how thick the disc traced by Cepheid is (Skowron et al. 2019b;Chen et al. 2019) agree with our results for the scale and trend found from    .

C2.2 Fits in 𝑣 𝑧
The right panel of Fig. C3 presents amplitudes for the fits in vertical velocity as a function of galactocentric radius.For   the amplitudes show a smooth oscillating pattern.The important difference between the   and   is that in  the  = 1 mode dominates the warp at all radii; in   the kinematic signal of the warp is dominated by both  = 1 and  = 2, a result unexpected for a tilted rings model.The  = 1 mode in   starts to appear at  ≈ 12 kpc and at its maximum reaches an amplitude similar to the value of    ≈ 7.2 km s −1 .For  = 2 in   there is an oscillation, as in for  = 1 too.The amplitude of none of the kinematic modes never exceeds the intrinsic dispersion, by contrast to the warp in , in which they do.However, the amplitude of the oscillations in  = 1, 2 is larger than the uncertainty in each mode, making the result more significant.
For the phases,   1 rises for  > 11 kpc and   2 is nearly constant, declining for  > 14 kpc.Since the amplitudes in   for  = 1 and  = 2 are comparable, the connection between these behaviours and the twisting of the LMV  is not as straightforward as in  where  = 1 clearly dominates and the LON twist is evident in the decline of  1 for outer radii.
Finally, the intrinsic dispersion for   (black curve, Fig. C3 right panel) is found to be almost constant with radius at    ≈ 7.2 km s −1 .This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure 1.The upper panels are the vertical height  as a function of galactocentric azimuth  for three different rings with radius 8.0 kpc (left), 11.0 kpc (middle) and 15 kpc (right).Gray dots represent the Cepheids in each ring, the black solid line represents the Fourier curve MAP fit to the gray dots, and the oranges curves are 500 random realisations of the Posterior PDF.The gray curves are fits to 200 bootstrapping realisations.The bottom panels show the residuals between the best fit and each bootstrap realisation.

Figure 2 .
Figure 2. The upper panels are the vertical velocity   as a function of galactocentric azimuth  for three different rings with radius 8.0 kpc (left), 11.0 kpc (middle) and 15 kpc (right).Gray dots represent the Cepheids in each ring, the black solid line represents the Fourier curve MAP fit to the gray dots, and the oranges curves are 500 random realisations of the Posterior PDF.The gray curves are fits to 200 bootstrapping realisation.The bottom panels show the residuals between the best fit and each bootstrap realisation.

Figure 3 .
Figure 3. Left: Difference between the north and south extreme of the warp as a function of galactocentric radii from our results (black curve), the same is calculated for the warp model by Skowron et al. (2019a) (doted orange curve).Right: Least angular difference between the north and south extremes as a function of galactocentric radii from our results (black curve), the same is calculated for the warp model of Skowron et al. (2019a) (doted orange curve).The grey dots are 500 random realisation at each ring taken from the Posterior.

Figure 4 .
Figure 4. Face-on view of the best fitting (MAP) warp model for the Cepheids (gray dots).The colour scale represents the mean Z of the disc (blue above the plane and red below it).The line of nodes (LON, i.e.  = 0) is indicated with the black curve.The line of maximum vertical velocity (LMV  ) is indicated by the dark green curve.The different coloured lines correspond to lines of constant galactocentric azimuth.

Figure 5 .
Figure 5. Azimuth as a function of galactocentric radius for: the LON from this work (black curve), Chen et al. (2019) (cyan dashed curve) and Dehnen et al. (2023) (red curve); the LMV  (thick green curve) and  1 for  = 1 mode from our  fits (red dots).The gray and olive green dots are 500 realisations of the LON and the LMV  taken from the posterior at each ring, respectively.The vertical dashed line indicates the Holmberg radius for the Milky Way from Chen et al. (2019).The crosses indicate a sample of independent (disjoint) rings.

Figure 6 .
Figure 6.The change in amplitude  1 for the  = 1 mode as a function of the galactocentric radii from our fits (black curve) and the  1 by Dehnen et al. (2023) (cyan curve).Red dots indicate measures for independent rings.The gray dots around each  1 are 500 realisation taken from the posterior at each ring.

Figure 7 .
Figure 7. Minus the angular frequency as a function of the galactocentric radius for the pattern speed for the  = 1 mode  1 (black curve) from our fits, the angular velocity given by the rotation curve Ω (red curve) (Ablimit et al. 2020), the upper and lower limit found by Poggio et al. (2020), and results from Dehnen et al. (2023).The dots around  1 are 500 posterior realisations at each ring.

Figure 8 .
Figure 8. Vertical height of the warp as a function of galactocentric radius for slices at  = 90 • and  = 270 • , for this work and warp models in the literature summarised in Table 3.The shaded region represents the uncertainty in the warp model from Wang et al. (2020) (see text for more details).

Figure 10 .
Figure 10. (left panel) and   (right panel) as a function of galactocentric radius for nine evenly spaced azimuthal cuts from  = 120 • to 240 • .The range of azimuth is selected to be in the region of the disc more populated by data and less affected by the SF.In the right panel the black curve represent the mean from the coloured ones Chrobáková et al. (2020) seems unrealistic when compared to the rest of results from the literature, even compared to those with similar tracers asCheng et al. (2020).Chrobáková & López-Corredoira (2021) also present a warp model for the younger population in its sample, with very similar results as obtained for the total sample.This particular disagreement in amplitude with the models of the young populations may indicate that the model byChrobáková et al. (2020) may be significantly underestimating precession rate of the warp, as a consequence of the underestimated amplitude.

Figure B1 .
Figure B1.Each panel in the left column shows  as a function of  centred at four different Galactocentric azimuths (from top to bottom  = 0.0 • , 90.0 • , 180.0 • and 270.0 • ) with Δ = 5 • .The same is shown in the right column for   .The different solid curves show: the analytical prediction for the mean position of test particles in the warped potential (red curve), the ground truth model (black curve) and the selection function model (blue curve) for the variable .The gray dots represent the stars used for the GT model, and the blue ones represent the star used for the SF model.The mean absolute difference between the SF and GT models for each  is reported in the title.

Figure B2 .
Figure B2.Each plot shows the residuals between the fiducial model (  ) and the model recovered with the mock catalogue (  ), normalised by the intrinsic dispersion of the variable obtained in the fiducial model.The left and right panels show the residuals for  and   , respectively, for  = 1 (top) and  = 2 (bottom).In each of the plots the inner black ring at  = 10 kpc indicates where the warps starts, the outer ring is where the  = 1 mode begins to be greater in amplitude than the intrinsic dispersion.This happens at  = 12 kpc for  and  = 14 kpc for   .The pale gray dots are the stars in the mock catalogue used in the fit, and the black star indicates the position of the Sun.

Figure B3 .
Figure B3.Each panel show the amplitude for the  = 0 (yellow),  = 1 (blue),  = 2 (violet) modes and   (black) for the SF model (dark solid lines) and the GT model (pale solid lines) obtained for  (left panel) and   (right panel).Each amplitude and   is plotted as a function of radii.

Figure B4 .
Figure B4.Each panel shows the phase for the modes  = 1 (blue) and  = 2 (violet) for the SF model (solid dark line) and the GT model (solid pale line) obtained for  (left panel) and   (right panel), as a function of galactocentric radius .

Figure B5 .
Figure B5.The ratio between the difference in intrinsic dispersion by the GT model and the one obtained with the SF model, over the intrinsic dispersion of the GT model as a function of the galactocentric radii.The blue and green continuous curves correspond to the intrinsic dispersion of  and   , respectively.The black, blue, and green vertical dashed lines indicate the radius at which the warp begins, when the amplitude of the mode  = 1 is bigger than the intrinsic dispersion for  and for   , respectively.

Figure B6 .
Figure B6.Angular frequency as a function of the galactocentric radii of the pattern speed  1 for the  = 1 mode, calculated from Eq. 14 for the GT model (black curve) and for the SF model (red curve), the angular velocity of the stars Ω (green curve) is derived analytically from the Allen & Santillan (1991) potential.The red dots around each  1 are 500 realisation taken from the posterior at each ring for the SF model.

Figure B7 .
Figure B7.The change in amplitude  1 for the  = 1 mode as a function of the galactocentric radius calculated from Eq. 13 for the GT model (black curve) and for the SF model (red curve).The red dots around each  1 are 500 realisation taken from the posterior at each ring for the SF model.

Figure C1 .
Figure C1.Reduced Chi-square for the fits in  done with  = 1 (green curve) and  = 2 (blue curve) as a function of galactocentric radius.

Figure C2 .
Figure C2.Reduced Chi-square for the fits in   done with  = 1 (green curve) and  = 2 (blue curve) as a function of galactocentric radius.

Figure C3 .
Figure C3.Amplitudes of each mode for  (left panel) and   (right panel) as a function of galactocentric radius.The black curve shows the intrinsic dispersion for each radius in the respective variable.The dotted line for the amplitudes in  shows the results from Skowron et al. (2019a).The colour dots around each mode are 500 realisation taken from the posterior at each ring.

Figure C4 .
Figure C4.Each panel shows the phases of each mode as a function of galactocentric radius.The first two top panels are the results for  and the bottom two for   (left  = 1 and right  = 2).The doted line for the phases in  are the constant phases obtained by Skowron et al. (2019a).The colour dots around each mode are 500 realisation taken from the posterior at each ring.

Table 1 .
Amplitudes and phases as a function of radius for the best fitting (MAP) models for  and   .The first few rows of the table are shown to provide guidance regarding its form and content.The full version of the table is available in the electronic version.

Table 2 .
Amplitudes and phases as a function of radius for 100 posterior realisations for the  and   models.The first few rows of the table are shown to provide guidance regarding its form and content.The full version of the table is available in the electronic version.

Table 3 .
Skowron et al. (2019a)ture.The asterisk ( * ) indicates the model shown in Fig.8.  ) 2 .In our model, without these assumptions, we can represent how the azimuthal geometry of the warp changes with the radius, giving rise to the differences between both models.TheSkowron et al. (2019a)model has the ability to reproduce the mean asymmetries observed in the warp, but not the LON twist or azimuthal changes in the different modes, which affect where the maxima are located.