An intuitive 3D map of the Galactic warp's precession traced by classical Cepheids

TL;DR

This paper uses classical Cepheids to create an accurate 3D model of the Milky Way's warp, revealing its orientation, spiral pattern, and relation to the gas disk, thus improving understanding of Galactic structure.

Contribution

It presents the first detailed 3D Galactic warp model based on Cepheids, showing the warp's orientation, spiral pattern, and consistency with gas disk morphology.

Findings

Warp's line of nodes is offset from the Galactic Center-Sun line.

The warp exhibits a leading spiral pattern consistent with Briggs' rule.

The stellar disk extends to at least 20 kpc, matching the gas disk morphology.

Abstract

The Milky Way's neutral hydrogen (HI) disk is warped and flared. However, a dearth of accurate HI-based distances has thus far prevented the development of an accurate Galactic disk model. Moreover, the extent to which our Galaxy's stellar and gas disk morphologies are mutually consistent is also unclear. Classical Cepheids, primary distance indicators with distance accuracies of 3-5%, offer a unique opportunity to develop an intuitive and accurate three-dimensional picture. Here, we establish a robust Galactic disk model based on 1339 classical Cepheids. We provide strong evidence that the warp's line of nodes is not oriented in the Galactic Center-Sun direction. Instead, it subtends a mean angle of 17.5 \pm 1 (formal) \pm 3 (systematic) and exhibits a leading spiral pattern. Our Galaxy thus follows Briggs' rule for spiral galaxies, which suggests that the origin of the warp is associated with torques forced by the massive inner disk. The stellar disk traced by Cepheids follows the gas disk in terms of their amplitudes; the stellar disk extends to at least 20 kpc. This morphology provides a crucial, updated map for studies of the kinematics and archaeology of the Galactic disk.

Tables

Table 1. Table 1: Parameters of the linear and power-law warp models applied. RMSE: Root mean square error.

$z_w=a{\left(R-R_w\right)}^b\sin \left(\varphi -{\varphi }_w\right)$
	$R_w$ (kpc)	${\varphi }_w$ (∘)	$a$	$b=1$	RMSE
All	$9.26\pm 0.19$	$17.4\pm 1.2$	$0.148\pm 0.006$		0.256
WISE	$9.31\pm 0.20$	$18.0\pm 1.6$	$0.146\pm 0.006$		0.215
Optical	$9.01\pm 0.40$	$16.1\pm 1.7$	$0.148\pm 0.011$		0.292
	$R_w$ (kpc)	${\varphi }_w$ (∘)	$a$	$b$	RMSE
All	$7.72\pm 0.27$	$17.5\pm 1.0$	$0.060\pm 0.012$	$1.33\pm 0.08$	0.210
WISE	$6.93\pm 0.72$	$18.6\pm 1.4$	$0.028\pm 0.015$	$1.61\pm 0.19$	0.188
Optical	$7.85\pm 0.27$	$15.8\pm 1.3$	$0.093\pm 0.019$	$1.14\pm 0.09$	0.225

Table 2. Table 2: Adopted infrared extinction laws and possible biases affecting the distance modulus.

	$\alpha$	$A_J/A_{K_s}$	$A_H/A_{K_s}$	$A_{W1}/A_{K_s}$	$A_{W2}/A_{K_s}$	$A_{\left[3.6\right]}/A_{K_s}$	$A_{\left[4.5\right]}/A_{K_s}$	$A_{\left[5.8\right]}/A_{K_s}$	$A_{\left[8.0\right]}/A_{K_s}$	Bias (mag)
ref. 9	2.05	3.005	1.717	0.506	0.340	0.478	0.341	0.234	0.321	0
	1.61	2.438	1.501	0.657	0.551	0.626	0.549	0.489	0.519	0.092
ref. 46	1.79	2.720	1.599	0.591	0.463	0.553	0.461	0.389	0.426	0.044
ref. 45	1.66	2.660	1.545			0.553	0.451	0.334	0.372	0.042
ref. 47	1.66	2.50	1.54			0.560	0.430	0.430	0.430	0.045

Full text

{affiliations}

Key Laboratory for Optical Astronomy, National Astronomical Observatories, Chinese Academy of Sciences, 20A Datun Road, Chaoyang District, Beijing 100012, China

Kavli Institute for Astronomy and Astrophysics, Peking University, Yi He Yuan Lu 5, Hai Dian District, Beijing 100871, China

Department of Astronomy, China West Normal University, Nanchong 637009, China

Department of Physics and Astronomy, Macquarie University, Balaclava Road, Sydney, NSW 2109, Australia

Research Centre for Astronomy, Astrophysics and Astrophotonics, Macquarie University, Balaclava Road, Sydney, NSW 2109, Australia

International Space Science Institute–Beijing, 1 Nanertiao, Zhongguancun, Hai Dian District, Beijing 100190, China

We have compiled samples of classical Cepheids from the Wide-field Infrared Survey Explorer (WISE) catalogue of periodic variables[8] (our ‘WISE Cepheid sample’) as well as from a number of optical surveys (collectively referred to as our ‘optical Cepheid sample’). We will discuss both samples separately, because the catalogues’ optical and infrared passbands are characterised by significantly different photometric and extinction sensitivities. Highly accurate Cepheid distances can be estimated using their well-established wavelength-dependent period–luminosity relations. To mitigate the influence of extinction in the Galactic plane and of photometric uncertainties at infrared wavelengths, we adopted the ‘infrared multi-passband optimal distance method’[9] to determine accurate Cepheid distances. Contaminants, including Type-II Cepheids, long-period eclipsing binaries and quasi-periodic variables were removed using Gaia Data Release 2 parallaxes[10]. Cepheids located in areas centered on the Magellanic Clouds were also excluded. Careful sample selection resulted in a tally of 2330 classical Cepheids for further analysis.

Distances were converted to 3D $XYz$ and spherical $R\phi z$ coordinates by adopting a reference frame centered on the Galactic Center and a solar Galactocentric distance $R_{0}=8.0$ kpc. Here, $\phi$ is the Galactocentric angle in the anticlockwise direction (aligned with the disk’s rotation axis) with respect to the solar position ($\phi_{\odot}=0^{\circ}$). Since Gaia parallaxes are reliable within $\sim$5 kpc, we only selected Cepheids within the volume $R<20,|z|<2$ kpc to avoid significant contamination by Type-II Cepheids. Our downselection included 1459 Cepheids with distances accurate to $<$5%, corresponding to a distance modulus standard deviation $<$0.108 mag. Cepheids located clearly away from the best-fitting warp model ($\Delta>1$ kpc) were also excluded (see below). Our final sample contained 1339 Cepheids.

Figure 1a shows the 3D distribution of our final Cepheid sample. It covers two-thirds of the disk. The bluish-violet ‘optical’ Cepheids are distributed in the solar neighbourhood and the Galactic anticenter direction. The red WISE Cepheids are predominantly located on either side of our optical sample. The Galactic warp is clearly visible, particularly its downward deviation from the Galactic plane. To model the warp, we first adopted the commonly used model, $z_{\rm w}=a(R-R_{\rm w})^{b}\sin(\phi-\phi_{\rm w})$, where $z_{\rm w}$, $R_{\rm w}$ and $\phi_{\rm w}$ are the warp height, its onset radius and the line of nodes (LON), respectively. Linear and power-law nonlinear least-squares fits were applied to all sample Cepheids. For $R<9.0$ kpc, a power-law description is most appropriate (see the grid in Figure 1), whereas the distribution becomes linear for $R>9.0$ kpc: see Figure 5 and Table 1. In addition, our WISE Cepheid sample is contaminated by fewer overtone Cepheids than the optical sample (see below). To limit the impact of possible selection effects, we fitted our model to the WISE and optical Cepheid samples both separately and simultaneously. Owing to their more extended distribution across the disk, parameters based on the WISE Cepheids are intrinsically more accurate. Nevertheless, all model results are internally consistent given the uncertainties. This implies that selection effects are thus minor or indeed negligible.

We carefully checked our results based on Markov-chain Monte Carlo (MCMC) simulations using the MATLAB MCMC toolbox[11]. The resultant excellent agreement, both as regards the parameter values and the well-defined Gaussian distributions (see Figure 7) implies that our results are indeed robust. The LON angle $\phi_{\rm w}$ is insensitive to the adopted type of model because of the weak correlation between $\phi_{\rm w}$ and $(a,b,R_{\rm w})$. We also considered a warp model expressed in spherical coordinates[12], but this only marginally affected the LON angle, because our Cepheids are located close to the Galactic plane.

To explore whether our derived stellar warp morphology agrees with that of the Galaxy’s Hi gas warp, at least regarding their $z$ heights, we projected the Cepheids’ 3D distribution onto the plane defined by the maximum height of the warp, assuming $\sin(\phi-\phi_{\rm w})=1$. Figure 2 compares our model with both the Hi warp and the distributions defined by a number of other tracers[1, 13, 14, 15]. The Cepheid warp is in excellent agreement with the Hi warp out to $R\simeq$ 15 kpc ($\Delta|z|<50$ pc). At larger radii, the Hi warp’s Fourier m = 0 and m = 2 modes start to dominate. This enhances its amplitude, particularly in the northern warp. Despite the detection of a similar $m=2$ amplitude for our Cepheid sample, $W_{2}=0.14\pm 0.03$ kpc kpc*-1*, that sample’s rather different LON angle $\phi_{\rm W2}=174\pm 10^{\circ}$—which deviates from the Hi LON by \sim$$45^{\circ}—precludes us from assessing this type of behaviour in the Cepheid warp (see Figure 7).

Although the warp traced by pulsars is generally comparable with our Cepheid warp, the roughly 30% uncertainties affecting pulsar distances render any details unreliable. The warps traced by Two Micron All Sky Survey[16] (2MASS) red clump stars agree with neither the Cepheid nor the Hi warps. This may be owing to sampling incompleteness at $R>13$ kpc: since red clump stars are some four magnitudes fainter than Cepheids, stellar crowding, the Sun’s location close to the warp and background glare towards the Galactic Center imply that the former are more likely detected close to the Galactic anticenter. This problem is compounded if the warp’s LON deviates from the Sun–Galactic Center direction (see below). Evolution of the warp’s morphology was initially suggested based on the shallow 2MASS data and a population synthesis model[17]. However, the recent warp kinematics map[18] based on Gaia proper motions, combined with a flat rotation curve, suggests that the old and young stellar populations exhibit similar warp features.

The warp’s LON, combined with theoretical predictions, can help us constrain the warp’s origin. However, the Milky Way’s LON has never been studied using tracers with distance accuracies better than 20%. The statistical and systematic uncertainties in the distances unavoidably caused reduced accuracies in both the LON’s mean value and its shape. The LON of the Milky Way’s Hi warp is closely aligned with the Galactic Center–Sun direction[19], $\phi_{\rm w}=0^{\circ}$. Similarly, observations of dust and stars with the Cosmic Background Explorer (COBE)[13] and of red clump stars with 2MASS[14] yielded $\phi_{\rm w}=0^{\circ}$ and $\phi_{\rm w}=-5\pm 5^{\circ}$ (statistical error), respectively. Although a 2MASS-based study of red clump stars and red giants[20] found $\phi_{\rm w}\sim 15^{\circ}$, its authors did not provide an assessment of the uncertainties nor of their selection effects. Our application of accurate distances implies that the Milky Way’s global LON deviates significantly from the Sun–Galactic Center direction, $\phi=17.5^{\circ}\pm 1.0^{\circ}$ (formal, including propagation errors) $\pm 3.0^{\circ}$. The latter, systematic error is introduced by the combination of the Sun’s height above the Galactic midplane, $z=25$ pc[21], and the difference between the current-best Galactocentric distance, $R_{\odot}=8.3$ kpc[22], and that adopted here. Note that we determined the mean LON angle and its uncertainty assuming that the LON does not vary with Galactocentric radius.

To ascertain whether the LON is stable as a function of galactocentric radius, we subdivided our Cepheid samples using two selection cuts: (i) objects within 1.0 kpc-wide bins and (ii) equal numbers of Cepheids (95) in each bin; 1.0 kpc is the optimal bin size according to the Freedman–Diaconis rule[23]. The power-law warp model was fitted to both samples, adopting fixed $(R_{\rm w},a,b)$. We also estimated the LONs and their uncertainties: see Figure 3 (top). Both trends are globally similar for $R=9$–16 kpc, which thus suggests that the general trend is not affected by problems associated with insufficient sample sizes. The blue dots are also plotted in polar coordinates to allow comparison with the disk’s rotation (see Figure 3, bottom). To validate the derived LON, we performed MCMC simulations of $\phi$ for different conditions. They included fixed and free $(a,b,R_{\rm w})$ parameters, 5% and 10% limiting distance accuracy cuts and resampling of the Cepheids in the northern and southern warps. We also estimated the propagated uncertainties associated with the objects’ distances. All resulting LONs are mutually consistent (see Figure 8). We adopted the largest values of the statistical and systematic errors as our final uncertainties.

We also attempted a kinematic analysis of the Cepheid sample, adopting the proper motions and radial velocities from Gaia Data Release 2. Since two-thirds of the Cepheids do not have radial velocity measurements, we evaluated their radial velocities assuming a flat rotation curve[24], $v_{\rm c}=240$ km s*-1*. The typical uncertainty was based on the scatter in the rotation velocities of the other one-third of our sample with such measurements, $\Delta v_{\rm c}=13$ km s*-1*. We converted the 3D velocities to $(v_{r},v_{\phi},v_{z})$ in Galactocentric coordinates, assuming[25] $(U,V,W)=(11.1,12.24,7.25)$ km s*-1*. The mean uncertainty $\langle\sigma_{v_{z}}\rangle$ is around 4.2 km s*-1*; only Cepheids with uncertainties $\sigma_{v_{z}}<5$ km s*-1* are plotted in the $XY$–$v_{z}$ diagram. The spatial LON agrees globally with the $v_{z}$ maxima distribution.

A clear increase of the LON is apparent at $12\leq R\leq 15$ kpc in Figure 3 (top). In theory[5], the retrograde precession rates of the outer disk caused by, respectively, the massive inner disk and external torques scale as $R^{-4}$ and $R$. This represents the first clue that the Galactic warp traces a leading spiral pattern, which validates the notion that the origin of the outer disk’s pattern is predominately induced by torques associated with the massive inner disk. However, near $R=15.5$ kpc the LON appears to twist, possibly because of either external forcing of the misaligned outer halo or satellite accretion. In addition, within $R=12$ kpc the decrease of the LON with radius is likely caused by a decrease of the rotational speed (see the diagonal ridge in the $v_{\phi}-R$ diagram[26]).

Although LON precession of the Galactic warp has not yet been detected, it has been reported for 12 other spiral galaxies[4]. Those latter galaxies approximately follow Briggs’ rule: the LON remains straight within $R_{25}$ and advances in the rotation direction from around the Holmberg radius, $R_{\rm{Ho}}$. The LON traced by Cepheids conforms with this rule. Quantitatively, for a radial thin-disk scalelength of $R_{\rm d}=2.6\pm 0.5$ kpc[21], the Milky Way’s $R_{25}=3.0R_{\rm d}$ and $R_{\rm{Ho}}=4.4R_{\rm d}$[27] are located at 7.8 and 11.4 kpc, respectively. These radii agree well with the onset radius of the warp and the leading spiral pattern. Although the uncertainty in the scalelength is significant, the agreement of the $R_{\rm{Ho}}/R_{25}$ ratio supports a similar warp pattern in the Milky Way as observed for Briggs’ spiral galaxies.

Finally, we estimated the $z$-height residuals, $\Delta|z|$: see Figure 4. The clear flare seen in the $\Delta|z|$ residuals confirms the high reliability of both the data and our warp model. To quantify the parameters of the flare, we estimated the scaleheight based on the top-10 percentile of Cepheids in 1 kpc bins. The Cepheid flare agrees well with the Hi flare in the region of overlap. In detail, the Cepheid flare is smoother in the inner disk, whereas the Hi flare is better defined in the outer disk because of the decreasing completeness of our Cepheid sample at those radii. Three of the five previously confirmed Cepheids[28] which are located behind the Galactic Center in the flare region are shown as magenta stars. If we assume that the flare morphology behind the Galactic Center is similar to that on the near side, these three Cepheids appear on the far end of the Cepheid (or gas) flare scaleheight at these radii.

{methods}

0.1 Cepheid sample selection and the optimal distance method

The ‘WISE Cepheids’ were detected based on the full five-year WISE all-sky survey. They are affected by incompleteness for long periods ($P>10$ days) and in crowded regions in the inner disk (because of the WISE observation model). The completeness of Cepheids in the WISE variables catalogue[8] is approximately 80% (with respect to the optical Cepheid catalogue[31]) in the magnitude range of interest. Significant incompleteness of the Cepheids in the WISE variables catalogue occurs at long periods ($P>10$ days) and for low amplitudes (${{\rm Amp}_{W1}}<0.2$ mag). WISE is ineffective in detecting overtone Cepheids, which are characterised by half the amplitudes of fundamental Cepheids. For example, based on 9649 classical Cepheids in the Magellanic Clouds[32], the mean $I$-band amplitudes of the fundamental-mode and first-overtone Cepheids are 0.47 mag and 0.20 mag for $P>2$ days (the period range where both types of Cepheids overlap). Statistically, the fraction of first-overtone Cepheids is only 3.4% (8 out of 237) among the Magellanic Cloud Cepheids rediscovered in the WISE variables catalogue. Because of the even smaller number of photometric detections in the Milky Way, no known Galactic overtone Cepheids in the optical Cepheid Catalogue[31] were rediscovered in the WISE catalogue.

Our optical Cepheid sample is based on detections in optical passbands. The sample was compiled based on both Cepheid catalogues[30, 31] and variable star catalogues[33, 34, 35, 36, 37]. WISE Cepheids were not double counted. Since not all of these catalogues clearly separate fundamental-mode and overtone Cepheids, the optical sample may be affected by distance problems caused by unrecognised overtone Cepheids. Nevertheless, the dominant contaminants are Type-II Cepheids, eclipsing binaries and rapidly rotating stars. By virtue of accurate Gaia parallaxes at distances within 5 kpc, most of these contaminants have been excluded (see below).

Distances to our selected Cepheids were determined using the infrared multi-passband optimal distance method[9]. We adopted the 2MASS $JHK_{\rm s}$, Spitzer Space Telescope[38] [3.6], [4.5], [5.8] and [8.0], and WISE[39] $W1$ and $W2$ filters. In each band, the distance modulus was estimated using ${\rm DM}_{\lambda}=\langle m_{\lambda}\rangle-M_{\lambda}-A_{\lambda}$. Here, $\langle m_{\lambda}\rangle$ is the mean apparent magnitude, $M_{\lambda}$ the absolute magnitude determined from the Galactic Cepheid period–luminosity relations[3, 40] and $A_{\lambda}$ the extinction given by the infrared extinction law[9] and $A_{K_{\rm s}}$. The $K_{\rm s}$-band extinction was adjusted to achieve a weighted average distance modulus with the smallest possible standard deviation. The weights were based on the total uncertainties in ${\rm DM}_{\lambda}$, which include the photometric error, deviations of single-epoch magnitudes from the mean magnitude, as well as period and extinction uncertainties. Weights were set to 0 in bands without detections. The extinction uncertainty $\sigma_{\lambda,{\rm ext}}$ increases as $A_{K_{\rm s}}$ increases, so the extinction uncertainty dominates the weights for objects affected by higher extinction. In fact, this method yields the optimal distance based on a balance of extinction and photometric errors. The near-infrared distance is usually determined as ${\rm DM}=\langle m_{K_{\rm s}}\rangle-M_{K_{\rm s}}-R_{K_{\rm s}}\times E(H-K_{\rm s})$ [41, 42]. If we adopt a weight of unity in $H,K_{\rm s}$ and 0 in any other band, both methods become identical.

For the full sample, the statistical error is the larger of the propagated error and the internal fitting error of the optimal distance. The propagated error includes photometric uncertainties, deviations from the mean magnitude, the intrinsic scatter in the period–luminosity relations and period uncertainties. For single-epoch 2MASS photometry, the mean uncertainties associated with conversion to mean magnitudes are 0.100, 0.082 and 0.076 mag in $JHK_{\rm s}$, respectively, if we adopt the average full amplitudes $0.345\pm 0.091$, $0.286\pm 0.087$ and $0.265\pm 0.087$ mag, respectively, based on 275 fundamental-mode Cepheids[43]. Some Cepheids observed with the Spitzer Space Telescope have observations obtained during two epochs; the adopted error is 0.05 mag. The amplitude relations in near-infrared bands are based on template light curves[43]. Mid-infrared amplitudes are assumed to be no larger than those in the $K_{\rm s}$ band. If these uncertainties are independent in each band, the final, propagated uncertainties are $\sigma_{1}=\sqrt{1/(\sum{1/{\sigma_{\lambda}}^{2}})}$. The intrinsic scatter in the period–luminosity relations for different wavelengths is notindependent; therefore, the scatter in the $W1$ band is adopted here (0.082 mag).

As regards the systematic error, the main contributors are the zero point of the period–luminosity relation and the choice of extinction law. The zero-point uncertainty of the infrared period–luminosity relation is around 0.033 mag[44, 9]. The uncertainty in the near-infrared extinction law is the main contributor to the distance error in the Galactic plane; it can be up to 15% for heavily obscured stars[45]. We adopted the infrared extinction law determined using Cepheids in the Galactic Center direction (see Table 2, first row). It is comparable to the disk’s extinction law based on red clump stars[46, 47, 48] if we adopt the same near-infrared extinction index $\alpha$. Since $\alpha$ could be variable, we estimate the mean bias in the distance modulus for different infrared extinction laws (see Table 2). Note that if $A_{W1}/A_{K_{\rm s}}$ and $A_{W2}/A_{K_{\rm s}}$ are not available[46, 48], the relative extinction values pertaining to the nearby [3.6] and [4.5] bands are adopted. Half of the difference in the distance modulus associated with adopting either $\alpha=1.61$ [49] or $\alpha=2.05$ was treated as the error in the extinction; the mean deviation was 0.046 mag. In all of these statistical and systematic uncertainties, uncertainties caused by the intrinsic scatter in the period–luminosity relations and the optimal distance fitting dominate.

0.2 Exclusion of Contaminants

The main contaminants, Type-II Cepheids and long-period contact binaries, were excluded based on parallax determinations from Gaia Data Release 2. Reliable parallaxes were selected by requiring $\varpi>0.2$ mas, $\sigma_{\varpi}/\varpi<0.2$ and $G<16$ mag, where $\varpi$ and $\sigma_{\varpi}$ are the Gaia parallaxes and their uncertainties, respectively, and G denotes Gaia G band magnitudes. False Cepheids were excluded based on the large differences between parallaxes derived from the period–luminosity relation distances and the actual Gaia parallaxes, $|\varpi_{\rm PL}-\varpi|>3\sigma_{\varpi_{t}}$. Here, $\sigma_{\varpi_{t}}$ is the square root sum of the parallax error and the photometric distance error. We did not correct for possible systematic offsets in the Gaia parallaxes, since any such offset is small compared with the other uncertainties we need to deal with. Nevertheless, we tested implementation of a correction of $-46$ $\mu$as[50]. The number of objects in our final sample only decreased by 29 (corresponding to Cepheids with parallax differences in the range 3–$5\sigma_{\varpi_{t}}$). This has a negligible influence on the resulting mean warp parameters and the LON: the difference associated with adopting the corrected or uncorrected parallaxes is less than 10% of the statistical uncertainty. Type-II Cepheids are typically 2–3 mag fainter than classical Cepheids for a given period, whereas long-period contact binaries are at least 4 mag fainter. If Type-II Cepheids at a true distance of 5 kpc were mistaken for classical Cepheids, distances of 12.5–20 kpc would be estimated, somewhat depending on the pulsation period. Long-period contact binaries at 5 kpc would be placed at distances in excess of 30 kpc if they were assumed to be classical Cepheids. Since Gaia parallaxes are reliable out to distances of order 5 kpc, they can be used as independent distance tracers to exclude contaminants. This thus ensures the integrity of our Cepheid sample within approximately 15 kpc.

Type-II Cepheids at distances of 5–8 kpc and $z$ heights $|z|<0.8$ kpc ($R=15$–20 kpc, $|z|<2.0$ kpc if treated as Type-I Cepheids) may not be umambiguously excluded based on the Gaia parallaxes. In young environments, the ratio of Type-II to Type-I Cepheids is small, however. Thanks to the warp feature, this ratio could be estimated. Since the warp is not obvious for $R<10$ kpc, we can assume a symmetrical distribution of Type-II Cepheids at positive and negative $z$ heights. If they are treated as Type-I Cepheids, half of the Type-II Cepheids would appear at $z=-z_{\rm w}$. In other words, a false warp will be produced by these remnant Type-II Cepheids. Based on this idea, the number of Cepheids located within 0.5 kpc in $z$ height of a false warp and 1.0 kpc away from the real warp are considered contaminations. The percentage of contaminants is $2n(-z_{\rm w})/n(z_{\rm w})=8$%, where $n(z_{\rm w})$ and $n(-z_{\rm w})$ are the numbers of Cepheids located in the real and the falsely negative warp in the raw sample, respectively. The 4% Type-II Cepheids in the falsely negative warp were excluded by the selection cut, whereas another (negligible) 4% (5 objects) may remain mixed in with our final sample.

0.3 Validation of the warp model

We considered both linear and power-law models to model the warp. Figure 5 shows that the power-law model is better than the linear model at radii up to $R=7$–9 kpc. For objects at $R>9$ kpc, the two models are comparable. This means that the linear model is not suitable for Cepheids at $R<9$ kpc. Therefore, the results for the linear model were determined using Cepheids at Galactocentric distances greater than 9 kpc. MCMC simulations were performed to verify the warp model and investigate correlations among the parameters. In Figure 6, the Gaussian distributions and the similar values validate the results of the nonlinear least-squares method.

We converted $R\phi z$ to $R\phi\theta$ to investigate the warp model in spherical coordinates and adopted $\psi=a\psi_{\rm w}(R-R_{\rm w})^{b}\sin(\phi-\phi_{\rm w})$[12], where $\psi$ is the tilt angle. The $\phi$ values thus determined, both the mean value and the corresponding values as a function of radius, based on this warp model are almost the same as those of our above results. This means that the spatial distortion caused by adopting spherical coordinates is small. This can be understood based on two arguments. First, the tilt angle of the Cepheid warp is small ($4^{\circ}$ at a distance of 17 kpc; see Figure 2), so the spatial distortion is negligible. Second, $\phi$ is almost independent of $R,z$ (see Figure 6), so that adoption of $R,z$ or $R,\theta$ has little influence on $\phi$.

To test whether or not the $m=2$ warp model is realistic, we rearranged the model to read $z_{\rm w}=W_{0}+W_{1}(R-R_{\rm w})\sin(\phi-\phi_{\rm w})+W_{2}(R-R_{\rm W2})\sin(2\phi-\phi_{\rm W2})$[1]. $W_{0},W_{1},W_{2}$ are the $z$ amplitudes of the $m=0,1,2$ modes, respectively, and $\phi_{\rm w},\phi_{\rm W2}$ are the LON angles for the $m=1,2$ modes. For convenience, $R_{\rm W2}$ was adopted as the sample’s minimum Galactocentric distance. $W_{2}$ is around $0.01\pm 0.01$ and $0.02\pm 0.02$ kpc kpc*-1* for Cepheids in the range $R<10,10<R<12$ kpc, which means that the $m=2$ warp model is not obvious at $R<12$ kpc. Analysis of 146 Cepheids at $R>15$ kpc shows an obvious $m=2$ amplitude, $W_{2}=0.14\pm 0.03$ kpc kpc*-1*, similar to that of the Hi model, $W_{2}({\rm H\sc I})=0.12$ kpc kpc*-1*, whereas $\phi_{\rm W2}$ is rather different. For the gas model, the line of maxima for the $m=2$ mode is roughly aligned with the lines of maxima of the $m=1$ mode, which is different from that of the Cepheids: the lines of maxima for the two modes deviate by approximately $45^{\circ}$. The best values from the MCMC simulation are shown in Figure 7.

0.4 Validation of the warp’s LON

To verify the apparent precession trend of the warp’s LON shown in Figure 3, different conditions which would affect the result are considered. The LON in each radial bin is first tested based on MCMC simulations with fixed and free $a,b,R_{\rm w}$ parameters. The 50, 16 and 84 percentiles in the probability distribution are adopted as the median value and the corresponding errors. A comparison with the results of our nonlinear least-squares fitting method is shown in Figure 8a, b. Since the Cepheid sample was selected by imposing a limiting distance accuracy of 5%, a test was done using a different sample containing Cepheids with distance accuracies better than 10%. The warp’s LONs in different radial bins were again analysed using MCMC simulations. The trend is shown in Figure 8c. The excellent agreement of LON trends in Figure 8a, b, c means that the LON is robust among different methods and samples.

The effect of the propagation of distance uncertainties was quantified by means of Monte Carlo simulations. In addition to the statistical uncertainties, we also simulated the deviations caused by inclusion of 10% overtone Cepheids and of 2.6% systematic distance uncertainties. Based on 2000 realisations, the mean LONs and their standard deviations are shown in Figure 8d. Another possible systematic effect may be caused by the unequal distributions of Cepheids in the northern and southern warps, or by their spatial clumpiness. Since the number of Cepheids in the northern warp is half that in the southern warp, we resampled the southern warp to consider equal numbers. To avoid Cepheid clumpiness, we did not include Cepheids at $|\phi|<10^{\circ}$ in our tests. The sample was randomly selected 1000 times; the mean LONs and their standard deviation are plotted in Figure 8e. Figure 8 shows that all LON trends agree well, which serves as strong validation of the precession trend. The final uncertainties in the LONs are based on the largest of the statistical uncertainties and the systematic deviations.

0.5 Kinematics of the Cepheid warp

As shown in the kinematic map based on upper main-sequence stars and giants[18], the maximum median value $v_{z}$ is around 7.5 km s*-1*. Indeed, stars around the warp LONs have higher absolute velocities, $v_{z}$. Investigation of the LON of the kinematic warp requires even higher accuracies for both distance and velocity measurements. Limited by the larger mean uncertainty $\langle\sigma_{v_{z}}\rangle=4.2$ km s*-1*, only Cepheids with uncertainties $\sigma_{v_{z}}<5$ km s*-1* appear in Figure 9. We therefore only use the kinematic map as an additional tool to characterise the spatial warp. The blue and red data points are Cepheids with obviously positive and negative $v_{z}$, respectively. They are indeed reliable given their high signal-to-noise ratios $>$ 3. The red dots representing the possible kinematic LONs evidently confirm a tilted LON and agree well with the spatial LONs (considering the prevailing uncertainties). This tilted LON traced by Cepheids is also consistent with that traced by upper main-sequence stars and giants[18]. In addition, more negative $v_{z}$ Cepheids are located around and beyond the maximum amplitude direction of the southern warp, which confirms that the orientation of the LON is on the left-hand side of the solar direction, $\phi>0$.