Issue 
A&A
Volume 601, May 2017



Article Number  A115  
Number of page(s)  14  
Section  Galactic structure, stellar clusters and populations  
DOI  https://doi.org/10.1051/00046361/201629916  
Published online  15 May 2017 
The kinematic signature of the Galactic warp in Gaia DR1
I. The Hipparcos subsample
^{1} Università di Torino, Dipartimento di Fisica, via P. Giuria 1, 10125 Torino, Italy
email: poggio@oato.inaf.it
^{2} Osservatorio Astrofisico di Torino, Istituto Nazionale di Astrofisica (INAF), Strada Osservatorio 20, 10025 Pino Torinese, Italy
^{3} School of Physics, Astronomy and Mathematics, University of Hertfordshire, College Lane, Hatfield AL10 9AB, UK
Received: 17 October 2016
Accepted: 10 February 2017
Context. The mechanism responsible for the warp of our Galaxy, as well as its dynamical nature, continues to remain unknown. With the advent of high precision astrometry, new horizons have been opened for detecting the kinematics associated with the warp and for constraining possible warp formation scenarios for the Milky Way.
Aims. The aim of this contribution is to establish whether the first Gaia data release (DR1) shows significant evidence of the kinematic signature expected from a longlived Galactic warp in the kinematics of distant OB stars. As the first paper in a series, we present our approach for analyzing the proper motions and apply it to the subsample of Hipparcos stars.
Methods. We select a sample of 989 distant spectroscopicallyidentified OB stars from the new reduction of Hipparcos, of which 758 are also in the first Gaia data release (DR1), covering distances from 0.5 to 3 kpc from the Sun. We develop a model of the spatial distribution and kinematics of the OB stars from which we produce the probability distribution functions of the proper motions, with and without the systematic motions expected from a longlived warp. A likelihood analysis is used to compare the expectations of the models with the observed proper motions from both Hipparcos and Gaia DR1.
Results. We find that the proper motions of the nearby OB stars are consistent with the signature of a kinematic warp, while those of the more distant stars (parallax <1 mas) are not.
Conclusions. The kinematics of our sample of young OB stars suggests that systematic vertical motions in the disk cannot be explained by a simple model of a stable longlived warp. The warp of the Milky Way may either be a transient feature, or additional phenomena are acting on the gaseous component of the Milky Way, causing systematic vertical motions that are masking the expected warp signal. A larger and deeper sample of stars with Gaia astrometry will be needed to constrain the dynamical nature of the Galactic warp.
Key words: Galaxy: kinematics and dynamics / Galaxy: disk / Galaxy: structure / proper motions
© ESO, 2017
1. Introduction
It has been known since the early HI 21cm radio surveys that the outer gaseous disk of the Milky Way is warped with respect to its flat inner disk (Burke 1957; Kerr 1957; Westerhout 1957; Oort et al. 1958), bending upward in the north (I and II Galactic quadrants) and downward in the south (III and IV Galactic quadrants). The Galactic warp has since been seen in the dust and stars (Freudenreich et al. 1994; Drimmel & Spergel 2001; LópezCorredoira et al. 2002; Momany et al. 2006; Marshall et al. 2006; Robin et al. 2008; Reylé et al. 2009). Our Galaxy is not peculiar with respect to other disk galaxies: more than 50 percent of spiral galaxies are warped (SanchezSaavedra et al. 1990; Reshetnikov & Combes 1998; Guijarro et al. 2010). The high occurrence of warps, even in isolated galaxies, implies that either these features are easily and continuously generated, or that they are stable over long periods of time. In any case, the nature and origin of galactic warps in general are still unclear (Sellwood 2013).
While many possible mechanisms for generating warps in disk galaxies have been proposed, which is actually at work for our own Galaxy remains a mystery. This is due to the fact that while the shape of the Galactic warp is known, its dynamical nature is not; vertical systematic motions associated with the warp are not evident in radio surveys that only reveal the velocity component along the lineofsight. Being located within the disk of the Milky Way, systematic vertical motions will primarily manifest themselves to us in the direction perpendicular to our lineofsight. More recent studies of the neutral HI component (Levine et al. 2006; Kalberla et al. 2007) confirm that the Galactic warp is already evident at a galactocentric radius of 10 kpc, while the warp in the dust and stellar components are observed to start inside or very close to the Solar circle (Drimmel & Spergel 2001; Derriere & Robin 2001; Robin et al. 2008). Thus, if the warp is stable, the associated vertical motions should be evident in the component of the stellar proper motions perpendicular to the Galactic disk.
A first attempt to detect a kinematic signature of the warp in the proper motions of stars was first made using OB stars (Miyamoto et al. 1988). More recently, a study of the kinematic warp was carried out by Bobylev (2010, 2013), claiming a connection between the stellargaseous warp and the kinematics of their tracers, namely nearby red clump giants from Tycho2 and Cepheids with UCAC4 (Fourth US Naval Observatory CCD Astrograph Catalog, Zacharias et al. 2013) proper motions. Using red clump stars from the PPMXL survey (Positions and Proper Motions “Extra Large” Catalog, Roeser et al. 2010), LópezCorredoira et al. (2014) concluded that the data might be consistent with a longlived warp, though they admit that smaller systematic errors in the proper motions are needed to confirm this tentative finding. Indeed, largescale systematic errors in the groundbased proper motions compromise efforts to detect the Galactic warp. The first real hopes of overcoming such systematics came with global spacebased astrometry. However, using Hipparcos (ESA 1997) data for OB stars, Smart et al. (1998) and Drimmel et al. (2000) found that the kinematics were consistent neither with a warp nor with a flat unwarped disk.
Before the recent arrival of the first Gaia data release (Gaia Collaboration 2016, DR1), the best allsky astrometric accuracy is found in the new reduction of the Hipparcos catalog (van Leeuwen 2008, HIP2), which improved the quality of astrometric data by more than a factor of two with respect to the original Hipparcos catalog. For the Hipparcos subsample the Gaia DR1 astrometry is improved further by more than an order of magnitude. The primary aim of this work is to assess whether either the HIP2 or the new Gaia astrometry for the OB stars in the Hipparcos shows any evidence of the systematics expected from a longlived warp. We choose the OB stars as they are intrinsically bright, thus can be seen to large distances, and are shortlived, so are expected to trace the motions of the gas from which they were born. We select stars with spectral types of B3 and earlier. For B3 stars, stellar evolutionary models (e.g. Chen et al. 2015) predict masses ranging from 7 to 9 M_{⊙}, corresponding to a MS time of about 10–50 Myr for solar metallicity. We find that the kinematics of this young population do not follow the expected signal from a longlived stable warp.
Our approach is to compare the observations with the expectations derived from a model of the distribution and kinematics of this young population of stars, taking into full account the known properties of the astrometric errors, thereby avoiding the biases that can be introduced by using intrinsically uncertain and biased distances to derive unobserved quantities. In Sect. 2 we describe the data for our selected sample of OB stars from HIP2 and from Gaia DR1. In Sect. 3 we present the model developed to create mock catalogs reproducing the observed distributions. In Sect. 4 we report the results of comparing the proper motion distributions of our two samples with the probability distribution of the proper motions derived from models with and without a warp. In the last sections we discuss the possible implications of our results and outline future steps.
2. The data
In this contribution we will analyze both the preGaia astrometry from the new Hipparcos reduction (van Leeuwen 2008), as well as from Gaia DR1 (Gaia Collaboration 2016). First, our approach here to analyzing the proper motions is significantly different from that used previously by Smart et al. (1998) and Drimmel et al. (2000) for the first Hipparcos release; any new results based on new Gaia data cannot be simply attributed to better data or better methods. Also a study of the Hipparcos error properties is necessary for understanding the astrometric error properties of the Hipparcos subsample in Gaia DR1 that we consider here, because of the intrinsic connection, by construction, between the astrometry of the new Hipparcos reduction and Gaia DR1 (Michalik et al. 2015). Finally, the two samples are complementary, as the Hipparcos sample can be considered more complete and is substantially larger than the Hipparcos subsample of OB stars in Gaia DR1 with superior astrometry, as explained below.
We select from the new Hipparcos catalog (hereafter HIP2) the young OB stars, because of their high intrinsic luminosity. Moreover, being shortlived, they are expected to trace the warped gaseous component. However, the spectral types in the HIP2 are simply those originally provided in the first Hipparcos release. In the hope that the many stars originally lacking luminosity class in the Hipparcos catalog would have by now received better and more complete spectral classifications, we surveyed the literature of spectral classifications available since the Hipparcos release. Most noteworthy for our purposes is the Galactic Ostar Spectroscopic Survey (GOSSS; Maíz Apellániz et al. 2011; Sota et al. 2011, 2014), an ongoing project whose aim is to derive accurate and selfconsistent spectral types of all Galactic stars ever classified as O type with B_{J} magnitude <12. From the catalog presented in Sota et al. (2014), which is complete to B_{J} = 8 but includes many dimmer stars, we imported the spectral classifications for the 212 stars that are present in the HIP2 catalog. Thirteen of these HIP2 sources were matched to multiple GOSSS sources, from which we took the spectral classification of the principle component. Also worth noting is the Michigan catalog of HD stars (Houk & Cowley 1994; Houk 1993, 1994; Houk & SmithMoore 1994; Houk & Swift 2000), which with its 5th and most recent release now covers the southern sky (δ < 5°), from which we found classifications for an additional 3585 OB stars. However, these two catalogs together do not cover the whole sky, especially for the B stars. We therefore had to resort to tertiary sources that are actually compilations of spectral classifications, namely the catalog of Stellar Spectral Classifications (4934 stars; Skiff 2014), and the Extended Hipparcos Compilation (3216 stars; Anderson & Francis 2012). In summary, we have spectral classifications for 11 947 OB stars in Hipparcos.
We select from the HIP2 only those stars with spectral type earlier than B3, with an apparent magnitude V_{J} ≤ 8.5, and with Galactic latitude b < 30°, resulting in 1848 OB stars. From this sample of HIP2 stars we define two subsamples: a HIP2 sample whose measured Hipparcos parallax is less than 2 mas, and a TGAS(HIP2; TychoGaia Astrometric Solution) sample consisting of those HIP2 stars that appear in the Gaia DR1 whose measured TGAS parallax is less than 2 mas. The cut in parallax, together with the cut in Galactic latitude, is done to remove local structures (such as the Gould Belt). Our HIP2 sample contains 1088 stars (including 18 stars without luminosity class), while our TGAS(HIP2) sample contains only 788 stars. This lack of HIP2 in TGAS stars is largely due to the completeness characteristics of DR1, discussed further in Sect. 3.4 below.
Fig. 1 Our final sample of Hipparcos OB stars on the sky, plotted in Galactic coordinates. The dashed line shows the orientation of the Gould belt according to Comeron et al. (1992). Colored points indicate the stars that are identified members of the OB associations Orion OB1 (red), Trumpler 10 (purple), Vela OB2 (blue), Collinder 121 (green) and Lacerta OB1 (cyan). 
Notwithstanding the parallax cut we found that there were HIP2 stars in our sample that are members of nearby OB associations known to be associated with the Gould Belt. We therefore removed from the HIP2 sample members of the Orion OB1 association (15 stars, as identified by Brown et al. 1994) and, as according to de Zeeuw et al. (1999), members of the associations Trumpler 10 (two stars), Vela OB2 (seven stars), and Lacerta OB1 (nine stars), all closer than 500 pc from the Sun. We also removed 33 stars from the Collinder 121 association as it is thought to also be associated with the Gould Belt. With these stars removed, we are left with a final HIP2 sample of 1022 stars. It is worth noting that the superior Gaia parallaxes already result in a cleaner sample of distant OB stars: only 21 members of the above OB associations needed to be removed from the TGAS(HIP2) sample after the parallax cut. Figure 1 shows the position of the stars in our two samples in Galactic coordinates.
Due to the above mentioned parallax cut, our sample mostly contains stars more distant than 500 pc. Though our analysis will only marginally depend on the distances, we use spectrophotometric parallaxes when a distance is needed for the HIP2 stars, due to the large relative errors on the trigonometric parallaxes. Absolute magnitudes and intrinsic colors are taken from Martins et al. (2005) and Martins & Plez (2006) for the O stars, and from Humphreys & McElroy (1984) and Flower (1977) for the B stars.
We note that the OB stars with V_{J} ≤ 7.5 (approximately 90% complete) beyond the Solar Circle (90° <l< 270°) show a tilt with respect to the Galactic plane that is consistent with a Galactic warp: a robust linear fit in the lb space (l normalized to 180°−l) yields a slope of 0.049 ± 0.007. However, given the possible effects of patchy extinction, it would be risky to make any detailed conclusions about the large scale geometry of the warp from this sample with heliocentric distances limited to a few kiloparsecs.
3. The model
Here we describe the model used to produce synthetic catalogs and probability distribution functions of the observed quantities to compare with our two samples of Hipparcos OBstars described above, taking into account the error properties of the Hipparcos astrometry (for the HIP2 sample) and of the Hipparcos subsample in Gaia DR1 (for the TGAS(HIP2) sample), and applying the same selection criteria used to arrive at our two samples. The distribution on the sky and the magnitude distribution of the HIP2 sample to V = 7.5 (556 stars), assumed to be complete, are first reproduced in Sects. 3.2 and 3.3, using models of the color–magnitude and spatial distribution of the stars, and a threedimensional extinction model. The completeness of the Hipparcos, and of TGAS with respect to Hipparcos, is described in Sect. 3.4 and is used to model our two samples down to V = 8.5. Then a simple kinematic model for the OB stars is used to reproduce the observed distribution of proper motions (Sect. 3.5) of our two samples, including (or not) the expected effects of a stable (longlived) nonprecessing warp (Sect. 3.1). A comparison of the observed samples with the expectations from the different warp/nowarp models is presented in Sect. 4.
The model that we present here is purely empirical. Many parameters are taken from the literature, while a limited number have been manually tuned when it was clear that better agreement with the observations could be reached. We therefore make no claim that our set of parameters are an optimal set, nor can we quote meaningful uncertainties. The reader should thus interpret our choice of parameters as an initial “first guess” for a true parameter adjustment, which we leave for the future when a larger dataset from Gaia is considered. In any case, after some exploration, we believe that our model captures the most relevant features of the OB stellar distribution and kinematics at scales between 0.5−3 kpc.
Comparison of warp parameters for the models of Drimmel & Spergel (2001) and Yusifov (2004).
3.1. Warp
In this Section we describe our model for the warp in the stellar spatial distribution and the associated kinematics. In Sect. 3.3 below we construct a flatdisk distribution with vertical exponential profile, then displace the zcoordinates by z_{w}, where: (1)The warp phase angle φ_{w} determines the position of the lineofnodes of the warp with respect to the Galactic meridian (φ = 0). The increase of the warp amplitude with Galactocentric radius, h(R), is described by the height function (2)where h_{0} and R_{w} are the warp amplitude and the radius at which the Galactic warp starts, respectively. The exponent α_{w} determines the warp amplitude increase. Table 1 reports three different sets of warp parameters taken from the literature and used later in our analysis in Sect. 4.2. Assuming that the Galaxy can be modeled as a collisionless system, the 0th moment of the collisionless Boltzmann equation in cylindrical coordinates gives us the mean vertical velocity . If we use the warped disk as described above and suppose that (i.e. that the disc is not radially expanding or collapsing), we obtain: (3)Equations (1)–(3)assume a perfectly static warp. It is of course possible to construct a more general model by introducing time dependencies in Eqs. (1)and (2), which will result in additional terms in Eq. (3), including precession or even an oscillating (i.e., “flapping”) amplitude. For our purpose here, to predict the expected systematic vertical velocities associated with a warp, such time dependencies are not considered.
This above model we refer to as the warp model, with three possible sets of parameters reported in Table 1. Our alternative model with z_{w} = 0 will be the nowarp model, where Eq. (3)reduces to the trivial .
Figure 2 shows the prediction of a warp model, without errors, for the mean proper motions μ_{b} in the Galactic plane. For the nowarp model (here not shown), we expect to have negative μ_{b} values symmetrically around the Sun as the reflex of the vertical component of Solar motion, progressively approaching 0 with increasing heliocentric distance. For a warp model, a variation of μ_{b} with respect to Galactic longitude is introduced, with a peak toward the anticenter direction (l = 180°). Figure 2 also shows that a variation of the warp phase angle has a rather minor effect on the kinematic signature, nearly indistinguishable from a change in the warp amplitude.
Fig. 2 According to the warp model, the true μ_{b} in the Galactic plane as a function of Galactic longitude at heliocentric distances of 0.5 kpc (A) and 1.5 kpc (B). For each set of curves, the thick line represents the case with warp phase φ_{w} = 0° and the two thin curves show φ_{w} = ± 20°. 
3.2. Luminosity function
There are different initial luminosity functions (ILF) in the literature for the upper main sequence (Bahcall & Soneira 1980; Humphreys & McElroy 1984; Scalo 1986; Bahcall et al. 1987; Reed 2001). Given the uncertainties in the ILF for intrinsically bright stars (absolute magnitude M< −3), we assume N(M) ∝ 10^{αM}, and use the value α = 0.72 that we find reproduces well the apparent magnitude distribution (Fig. 7) with the spatial distribution described in Sect. 3.3. We use a main sequence color–magnitude relation consistent with the adopted photometric calibrations (see Sect. 2). Absolute magnitudes M are randomly generated consistent with this ILF then, for a given absolute magnitude, stars are given an intrinsic color generated uniformly inside a specified width about the main sequence (Fig. 3), which linearly increases as stars get fainter. According to an assumed giant fraction (see below), a fraction of the stars are randomly labeled as giant. The absolute magnitude of these stars is incremented by −0.5 mag, and their color is generated uniformly between the initial main sequence color and the reddest value predicted by our calibrations, (B−V)_{0} = −0.12. The giant fraction f_{g} has been modeled as a function of the absolute magnitude as follows:in order to roughly reproduce the fraction of giants in the observed catalog. We caution that this procedure is not intended to mimic stellar evolution. Instead, we simply aim to mimic the intrinsic color–magnitude distribution (i.e., Hess diagram) of our sample.
Fig. 3 Color–magnitude diagram used to generate synthetic intrinsic colors. The dark and light orange regions show, respectively, the main sequence and the giant regions. The density of the two regions (here not shown) depends on the ILF and on the giant fraction. 
3.3. Spatial distribution and extinction
Since we wish to model the distribution of OB stars on scales larger than several hundred parsecs, we use a mathematical description of this distribution that smooths over the inherent clumpy nature of star formation, which is evident if we consider the distribution of young stars within 500 pc of the Sun. On these larger scales it is nevertheless evident that the OB stars are far from being distributed as a smooth exponential disk, but rather trace out the spiral arms of the Galaxy, being still too young to have wandered far from their birthplaces. In our model we adopt R_{0} = 8.2 kpc as the Sun’s distance from the Galactic center, and a solar offset from the disk midplane of z_{0} = 25 pc, for galactocentric cylindrical coordinates (R,φ,z), as recommended by BlandHawthorn & Gerhard (2016). For the spiral arm geometry we adopt the model of Georgelin & Georgelin (1976), as implemented by Taylor & Cordes (1993), rescaled to R_{0} = 8.2 kpc, with the addition of a local arm described as a logarithmic spiral segment whose location is described by R_{Loc} = R_{Loc,0}exp−(tanpφ), p being the arm’s pitch angle. The surface density profile across an arm is taken to be gaussian, namely: , where d_{a} is the distance to the nearest arm in the R,φ plane, and w_{a} = c_{w}R is the arm halfwidth, with c_{w} = 0.06 (Drimmel & Spergel 2001). An overview of the modeled surface density distribution is shown in Fig. 4. The stars are also given an exponential vertical scale height ρ ∝ exp(−  z′  /h_{z}), where h_{z} is the vertical scale height, z′ = z−z_{w}−z_{Loc}, z_{w}(R,φ) is the height of the warp as described in Sect. 3.1 and z_{Loc} is a vertical offset applied only to the local arm.
Fig. 4 Modeled surface density of the OB stars. Sun’s position is indicated by the star. 
We generate the above spatial distribution in an iterative MonteCarlo fashion. Ten thousand positions in (x,y) coordinates are first generated with a uniform surface density to a limiting heliocentric distance of 11 kpc, and with an exponential vertical profile in  z′ . The relative surface density Σ(x,y) is evaluated at each position according to our model described above, and positions are retained if u < Σ(x,y)/max(Σ), where u is a uniform random deviate between 0 and 1. Each retained position is assigned to a (M_{V},(B−V)_{0}) pair, generated as described in the previous section. The extinction to each position is then calculated using the extinction map from Drimmel et al. (2003) and the apparent magnitude is derived. Based on the apparent magnitude, a fraction of the stars with V ≤ 8.5 are randomly retained consistent with the completeness model of Hipparcos, while for a TGASlike catalog an additional random selection of stars is similarly made as a function of the observed apparent magnitude and color, as described in the following section. This procedure is iterated until a simulated catalog of stars is generated matching the number in our observed sample.
Good agreement with the HIP2 distribution in Galactic latitude was found adopting a vertical scale height h_{z} = 70 pc and assuming z_{Loc} = 25 pc, as shown in Fig. 5. Figure 6 compares the modeled and observed distributions in Galactic longitude, which is dominated by the local arm. This observed distribution is reproduced by placing the local arm at a radius of R_{Loc,0} = 8.3 kpc, with a pitch angle of 6.5° and a halfwidth of 500 pc. (The curves in Figs. 5 and 6 are nonparametric fits to the distributions obtained through kernel density estimation with a Gaussian kernel, as implemented by the generic function density in R^{1}. The smoothing bandwith is fixed for all the curves in the same figure, with values of 2.5° and 15° for the latitude and the longitude distribution, respectively).
Fig. 5 Latitude distribution of the HIP2 OB stars. The red curve is a nonparametric fit to the selected HIP2 sample, the red dashed curve shows the additional contribution of the Orion OB1 association, while the red dotted curve shows the added contributions of the Trumpler 10, Vela OB2, Collinder 121, and Lacerta OB1 associations. The blue curve and lightblue shaded area shows the average and 2σ confidence band of the simulated longitude distribution, based on 30 simulated instances of the sample. The black dotted and dashdot curves show the relative contributions of the major spiral arms (SagittariusCarina and Perseus) and the local arm, respectively, while the additional black solid curve is for a model without a warp. 
Fig. 6 Longitude distribution of the HIP2 OB stars. Meaning of the curves are as in the previous figure. 
Figure 7 shows the resulting apparent magnitude distribution, as compared to the HIP2 sample. Comparing the observed and the simulated longitude distributions in Fig. 6, we note that our model fails to reproduce well the observed distribution in the longitude range l = 300−360 degrees. This probably reveals a deficit in the geometry adopted for the SagittariusCarina arm, which we have not attempted to modify as we are primarily interested in the kinematics toward the Galactic anticenter. It should also be noted that, for both the longitude and latitude distributions, the presence or absence of a warp (modeled as described in Sect. 3.1) has very little effect.
Fig. 7 Apparent magnitude distribution for the data (histogram) and the simulations (black dots). The error bars show 2σ uncertainty, calculated with 30 simulated samples. 
Our approach assumes that the Hess diagram is independent of position in the Galaxy. We recognize this as a deficit in our model, as the spiral arms are in fact star formation fronts, in general moving with respect to Galactic rotation. We thus expect offsets between younger and older populations, meaning that the Hess diagram will be position dependent. However, if the Sun is close to corotation, as expected, such offsets are minimal.
3.4. Completeness
Fig. 8 Fraction of Hipparcos completeness in function of apparent magnitude V_{T} with respect to the Tycho2 catalog. The dashed and the dotted line represent, respectively, the Hipparcos fraction for the stars above and below δ = −30°. 
The completeness of the Hipparcos catalog decreases with apparent magnitude, as shown by the fraction of Tycho2 stars in the Hipparcos catalog (Fig. 8). At V_{T} ≤ 7.5 the HIP2 catalog is approximately 90% complete, reaching approximately 50% completeness at V_{T} = 8.5. The fact that the Hipparcos catalog was based on an input catalog built from then extant groundbased surveys results in inhomogenous sky coverage. In particular we find a north/south dichotomy at Declination δ ≈ −30°. We assume that the completeness fraction of Hipparcos stars in function of V_{T} decreases in a similar way for the Johnson magnitude, that is f_{HIP}(V_{T}) ≈ f_{HIP}(V_{J}).
As already noted in Sect. 2, TGAS contains only a fraction of the stars in Hipparcos. We find that the completeness of TGAS with respect to Hipparcos is strongly dependent on the observed magnitude and color of the stars: 50% completeness is reached at V_{J} = 6.5 mag and B−V = 0, with the brightest and bluest stars missing from TGAS. This incompleteness is a result of the quality criteria used for constructing TGAS and of the difficulty of calibrating these stars due to their relative paucity. Figure 9 shows a map of the TGAS(HIP2) completeness as a function of apparent magnitude and color. The completeness reaches a maximum plateau of about 80%; however, this is not uniform across the sky. Due to the scanning strategy of the Gaia satellite, and the limited number of months of observations that have contributed to the DR1, some parts of the sky are better covered than others. This results in a patchy coverage, which we have not yet taken into account. However, this random sampling caused by the incomplete scanning of the sky by Gaia is completely independent of the stellar properties, so that our TGAS(HIP2) sample should trace the kinematics of the stars in an unbiased way.
Fig. 9 Fraction of HIP2 OB stars present in HIPTGAS as a function of the observed color and the apparent magnitude. 
3.5. Kinematics
Now that the spatial distribution has been satisfactorily modeled, we can address the observed distribution of proper motions. We point out that the kinematic model described in this section is independent from the inclusion (or not) of the warpinduced offset in latitude proper motions (Sect. 3.1). We adopt a simple model for the velocity dispersions along the three main axes of the velocity ellipsoid: , where h_{R} = 2.3 kpc is the radial scale length and are the three velocity dispersions in the solar neighborhood for the bluest stars (Dehnen & Binney 1998). A vertex deviation of l_{v} = 30° is implemented, as measured by Dehnen & Binney (1998) for the bluest stars, although we find that it has no significant impact on the proper motion trends.
As recommended by BlandHawthorn & Gerhard (2016), we adopted Θ_{0} = 238 km s^{1} for the circular rotation velocity at the Solar radius R_{0}. Given that current estimates of the local slope of the rotation curve vary from positive to negative values, and that our data are restricted to heliocentric distances of a few kpc, we assume a flat rotation curve. In any case, we have verified that assuming a modest slope of ± 5 km s^{1}/kpc does not significantly impact the expected trend in proper motions.
After local stellar velocities (U,V,W) of the stars are generated, proper motions are calculated assuming a Solar velocity of v_{⊙} = (U_{⊙},V_{⊙},W_{⊙}) = (11.1,12.24,7.25) km s^{1} (Schönrich et al. 2010). Observed proper motions in (α,δ) are derived by adding random errors as per an astrometric error model, described in the Sect. 3.6. Finally, the proper motions in equatorial coordinates are converted to Galactic coordinates, that is (μ_{l} ∗ ,μ_{b}), where μ_{l} ∗ = μ_{l}cosb.
Figure 10 shows the derived proper motions in Galactic longitude for both the data and simulations using the bivariate localconstant (i.e. NadarayaWatson) kernel regression implemented by the npregbw routine in the np R package with bandwidth h = 45°. The solid black line shows the trend obtained for the simulation with the above listed standard parameters. Our simple model of Galactic rotation fails to reproduce the observations, even if we assume Θ_{0} = 220 or 260 km s^{1} (upper and lower dashdotted black lines, respectively). We also tried modifying the (U_{⊙},V_{⊙},W_{⊙}) components of the solar motion (equivalent to adding a systematic motion to the Local Standard of Rest), but without satisfactory results. We finally obtained a satisfactory fit by assigning to the stars associated with the Local Arm an additional systematic velocity of ΔV_{C} = 6 km s^{1} in the direction of Galactic rotation and ΔV_{R} = 1 km s^{1} in the radial direction. Such a systematic velocity could be inherited from the gas from which they were born, which will deviate from pure rotation about the Galactic center thanks to postshock induced motions associated with the Local Arm feature. Similar, but different, systematics may be at play for the other major spiral arms, which we have not tried to model given the limited volume that is sampled by this Hipparcos derived dataset. In any case, the addition (or not) of these systematic motions parallel to the Galactic plane does not significantly influence the the proper motions in Galactic latitude, as discussed in Sect. 4.
Fig. 10 μ_{l} in function of Galactic longitude for the data (red curve) together with the 95% bootstrap confidence band (pink shaded area). The three black dotted curves show the trend obtained with simulations with circular velocity 260, 238 and 220 km s^{1}, respectively, from the lowest to the highest curve. Simulations with an additional velocity to the Local Arm (see text) produce the blue curve, for which the light blue shaded area shows the 2σ confidence band, calculated with 30 simulated catalogs. 
As is well known, the International Celestial Reference Frame (ICRF) is the practical materialization of the International Celestial Reference System (ICRS) and it is realized in the radio frequency bands, with axes intended to be fixed with respect to an extragalactic intertial reference frame. The optical realization of the ICRS is based on Hipparcos catalog and is called Hipparcos Celestial Reference Frame (HCRF). The reference frame of the new reduction of Hipparcos catalog was identical to the 1997 one, aligned with the ICRF within 0.6 mas in the orientation vector (all three components) and within 0.25 mas/yr in the spin vector ω (all three components) at the epoch 1991.25 (van Leeuwen 2008). It is evident that a nonzero residual spin of the HCRF with respect to the ICRF introduces a systematic error in the Hipparcos proper motions. Depending on the orientation and on the magnitude of the spin vector, the associated systematic proper motions can interfere or amplify a warp signature and must therefore be investigated and taken into account (Abedi et al. 2015; Bobylev 2010, 2015). In the following section, when modeling the HIP2 sample, we consider the effects of such a possible spin, adding the resulting systematic proper motions to the simulated catalogs following Eq. (18) of Lindegren & Kovalevsky (1995), and using the residual spin vector (ω_{x},ω_{y},ω_{z}) ≃ (−0.126, + 0.185, + 0.076) mas yr^{1} as measured by Gaia (Lindegren et al. 2016).
3.6. Error model
Our approach to confronting models with observations is to perform this comparison in the space of the observations. Fundamental to this approach is having a proper description (i.e. model) of the uncertainties in the data. For this purpose we construct an empirical model of the astrometric uncertainties in our two samples from the two catalogs themselves. Below we first describe the astrometric error model for the HIP2 sample, based on the errors in the HIP2 catalog, and then that of the TGAS(HIP2) sample based on the error properties of the Hipparcos subsample in Gaia DR1. We note that, while we are here principally interested in the proper motions, we must also model the uncertainties of the observed parallaxes ϖ since we have applied the selection criteria ϖ < 2 mas to arrive at our OB samples, and this same selection criteria must therefore be applied to any synthetic catalog to be compared to our sample.
3.6.1. Hipparcos error model
It is known that the Hipparcos astrometric uncertainties mainly depend on the apparent magnitude (i.e. the signaltonoise ratio of the individual observations) and on the ecliptic latitude as a result of the scanning law of the Hipparcos satellite, which determined the number of times a given star in a particular direction on the sky was observed. These dependencies are not quantified in van Leeuwen (2008), which only reports the formal astrometric uncertainties for each star. To find the mean error of a particular astrometric quantity as a function of apparent magnitude and ecliptic latitude we selected the stars with (B−V) < 0.5 from the HIP2 catalog, consistent with the color range of our selected sample of OB stars. We then binned this sample with respect to apparent magnitude and ecliptic latitude and found, for each bin, the median errors for right ascension α, declination δ, parallax ϖ, proper motions μ_{α ∗}, and μ_{δ}. The resulting tables are reported in the Appendix, which gives further details on their construction.
Fig. 11 Histogram of the observed parallax distribution. The dashed and the solid curves show, respectively, the synthetic distributions with F = 1 and 1.5 (see text for explanation). 
However, before using these formal HIP2 uncertainties to generate random errors for our simulated stars, we first evaluated whether the formal errors adequately describe the actual accuracy of the astrometric quantities. For this purpose the distribution of observed parallaxes is most useful, and in particular the tail of the negative parallaxes, which is a consequence of the uncertainties since the true parallax is greater than zero. In fact, using the mean formal uncertainties in the parallax to generate random errors, we are unable to reproduce the observed parallax distribution in our sample (see Fig. 11). Assuming that our model correctly describes the true underlying distance distribution, we found that the formal HIP2 uncertainties must be inflated by a factor of F = 1.5 to satisfactorily reproduce the observed distribution. This factor F is then also applied to the mean formal uncertainties of the other astrometric quantities. We note that this correction factor is larger than that implied from an analysis of the differences between the Hipparcos and Gaia DR1 parallaxes. (See Appendix B of Lindegren et al. 2016.)
To better fit the HIP2 proper motion distributions we also take into consideration stellar binarity. Indeed, approximately f_{b} ≈ 20% of stars of our sample have been labeled as binary, either resolved or unresolved, in the HIP2 catalog. For these stars, the uncertainties are greater than for single stars. Therefore, we inflated the proper motion errors for a random selection of 20% of our simulated stars by a factor of f_{bin} = 1.7 to arrive at a distribution in the errors comparable to the observed one.
Finally, we also performed similar statistics on the correlations in the HIP2 astrometric quantities published by van Leeuwen (2008), using the four elements of the covariance matrix relative to the proper motions (see Appendix B of Michalik et al. 2014). We found that the absolute median correlations are less than 0.1, and therefore we do not take them into account.
3.6.2. TGAS(HIP2) error model
A detailed description of the astrometric error properties of the TGAS subset in Gaia DR1 is described in Lindegren et al. (2016). However, on further investigation we found that the error properties of the subset of 93 635 Hipparcos stars in Gaia DR1 are significantly different with respect to the larger TGAS sample. In particular, we found that the zonal variations of the median uncertainties seen with respect to position on the sky are much less prominent for the Hipparcos stars in DR1, and are only weakly dependent on ecliptic latitude. The parallax errors with respect to ecliptic latitude are shown in Fig. 12. Meanwhile the TGAS(HIP2) parallax errors show no apparent correlation with respect to magnitude or color. Figure 13 shows the distribution of parallax uncertainties for three ecliptic latitude bins, which we model with a gamma distribution having the parameters reported in Table 2. However, Lindegren et al. (2016) have warned that there is an additional systematic error in the parallaxes at the level of 0.3 mas. In this work we only use the parallaxes to split our sample into two subsets, and we have verified that adding an additional random error of 0.3 mas to account for these possible systematic errors does not affect our results.
Fig. 12 Logarithm of the parallax errors (mas) in function of ecliptic latitude for the Hipparcos subsample in TGAS. The points show the medians, while the error bars show the 10th and the 90th percentiles of the distribution. 
Fig. 13 Distribution of the logarithm of the parallax uncertainties (mas) for the HIPTGAS stars. Three subsets with different ecliptic latitude are shown. 
The errors for the proper motions also show a weak dependence on ecliptic latitude, as well as additional dependence with respect to magnitude. Indeed, we find that the proper motion errors for the Hipparcos subset in Gaia DR1 are strongly correlated the Hipparcos positional errors, as one would expect, given that the Hipparcos positions are used to constrain the Gaia DR1 astrometric solutions (Michalik et al. 2015). We use this correlation to model the proper motion errors of the Hipparcos subsample in DR1. Figure 14 show the agreement which results when we take as our model , where F is the correction factor applied to the Hipparcos astrometric uncertainties, as described in Sect. 3.6.1, is the Hipparcos error in right ascension, interpolated from Table A.1 in the Appendix, and Δt is the difference between the Gaia (J2015) and Hipparcos (J1991.25) epoch. The adopted coefficient C_{α} = 1.42 is the median of for the stars of the TGAS(HIP2) sample. An analogous model is used for σ_{μδ}, with C_{δ} = 1.44.
Fig. 14 Comparison of the published error σ_{μα ∗}(TGAS) to the prediction based on Hipparcos uncertainties for each star of the Hipparcos subset in TGAS (see text). The dashed line represents the bisector. The solid line has null intercept and coefficient C_{α} = 1.42, which is used to calibrate our error model (see text). 
Finally, in contrast to the correlations in the HIP2 sample, we find that the correlations in DR1 between the astrometric quantities of the Hipparcos subsample vary strongly across the sky, but are significantly different from the complete TGAS sample, shown in Fig. 7 of Lindegren et al. (2016). Figure 15 shows the variation across the sky of the correlations between the parallaxes and the proper motions. As we can see, the correlations between proper motion components are relevant. To take this into account, we generate the synthetic proper motion errors from a bivariate Gaussian distribution with a covariance matrix which includes the σ_{μα} and σ_{μδ} predicted by the above described model and the correlations from the first map in Fig. 15.
Fig. 15 Median correlations between parallaxes ϖ and proper motions μ_{α}, μ_{δ} in HIPTGAS. 
4. Comparison between models and data
In this section we first attempt to derive the systematic vertical motions from the observed proper motions and compare these with those predicted from the model, to demonstrate the weaknesses of this approach. We then compare the observed proper motions of our two samples with the proper motion distributions derived from models with and without a warp.
4.1. Mean vertical velocity in function of Galactocentric radius
We first consider the mean vertical velocity as a function of Galactocentic radius R that one would derive from the measured proper motions and spectrophotometric distances, comparing the data with what we expect from the nowarp and warp models. In the nowarp case, the true mean vertical velocities are zero, while a warp model predicts that they increase with R outside the radius R_{w} at which the Galactic warp starts (see Eq. (1)). Figure 16 shows the mean vertical velocities, after removing the solar motion, for the data and the nowarp and warp simulated catalogs. The simulated catalogs include the modeled errors, as described in Sect. 3.6. Taking into account both distance and proper motion errors, the observed trend is biased toward negative velocities with increasing distance. This bias is particularly evident with the nowarp model, where the true (dashed line in Fig. 16), but similarly affects the warp model. One might be tempted to proceed to compare models to the data in this space of derived quantities, assuming the error models are correct, but this approach gives the most weight to the data at large distances, that is, those with the highest errors and the most bias. Indeed, from Fig. 16 one might quickly conclude that the data was consistent with the nowarp model, based however on trends that are dominated by a bias in the derived quantities.
A better approach is to compare the data to the models in the space of the observations, that is, the mean proper motions as a function of position on the sky, thereby avoiding the biases introduced by the highly uncertain distances. That is, it is better to pose the question: which model best reproduces the observations? Our approach will make minimum use of the spectrophotometric distances and parallaxes to avoid the strong biases introduced when using distances with relatively large uncertainties to arrive at other derived quantities, as in the example above. Indeed, whether based on spectrophotometric data or parallaxes, distance is itself a derived quantity that can suffer from strong biases (BailerJones 2015). Nevertheless, distance information is useful. We have already used the distance information contained in the parallaxes for selecting our two samples with the criteria of ϖ< 2 mas (see Sect. 2), while in Sect. 4.2, we will split our HIP and TGAS(HIP) samples into nearby and distant sources.
Fig. 16 Bivariate NadarayaWatson regression estimator of stellar vertical velocities as a function of Galactocentric radius, using a bandwidth of h = 0.5 kpc. The same regression bandwidth has been used for the data (red), the nowarp model (black) and warp model of Yusifov (2004) (blue). The 95% bootstrap confidence band is shown for the data. For each of the two models, the non parametric regressions are performed for 20 simulated catalogs, obtaining the curves as the mean values and the shaded areas as the 95% uncertainty. 
4.2. Proper motion μ_{b} in function of Galactic longitude
Fig. 17 Distribution of the TGASHIP data in the lμ_{b} plane (black dots), together with the probability density function P(μ_{b}  l) predicted by the nowarp model (left panel) and the warp model of Yusifov (2004; right panel). 
The aim of the present work is to determine whether a warp is favored over a nowarp model in either of our Hipparcos or Gaia DR1 samples. In Sect. 3.1 we showed that the Galactic warp, if stable, results in a distinct trend in the proper motions perpendicular to the Galactic plane μ_{b} with respect to Galactic longitude, with higher (i.e. more positive) proper motions toward the Galactic anticenter (see Fig. 2). Our approach is to compare the distribution of the observed proper motions with the expectations derived from a warp model and a nowarp model, taking into full account the known properties of the astrometric errors. To achieve this we adopt the approach of calculating the likelihood associated with a given model as the probability of the observed data set arising from the hypothetical model (as described in Peacock 1983).
Given an assumed model (i.e. parameter set), we generate 500 thousand stars and perform a twodimensional kernel density estimation in l−μ_{b} space to derive the conditional probability density function (PDF), P(μ_{b}  l), the probability of observing a star with proper motion μ_{b} at a given Galactic longitude, where P(μ_{b}  l) is constrained to satisfy ^{∫}P(μ_{b}  l) dμ_{b} = 1. The motivation for using P(μ_{b}  l) is that we want to assign the probability of observing a given value of μ_{b} independent of the longitude distribution of the stars, which is highly heterogeneous. That is, we wish to quantify which model best reproduces the observed trend of the proper motions with respect to longitude. In any case, we also performed the below analysis using P(l,μ_{b}), imposing the normalization ^{∫}P(l,μ_{b}) dldμ_{b} = 1, and obtain similar results. Figure 17 shows the conditional PDFs for the warp/nowarp models for the TGAS(HIP2) sample. The PDFs for the HIP2 sample (not shown) are very similar, as the proper motion distribution is dominated by the intrinsic velocity dispersion of our sample rather than by the proper motion errors.
Once the PDFs for the two different models are constructed, we found the probability P(μ_{b,i}  l_{i}) associated with each ith observed star according to each PDF. The likelihood associated with the model is , where N is the total number of stars in our dataset; for computational reasons, we used instead the loglikelihood . Also for practical reasons we applied a cut in μ_{b}, considering only the range (− 10 <μ_{b}< 5) mas/yr when calculating ℓ, reducing our HIP2 dataset to 989 stars, and our TGAS(HIP2) sample to 791 stars. Below we will confirm that this clipping of the data does not impact our results by considering alternative cuts on μ_{b}.
For a given sample, the difference between the loglikelihoods of a warp model and the nowarp model (i.e. the ratio of the likelihoods), Δ ≡ ℓ_{WARP}−ℓ_{NOWARP}, is found as a measure of which model is more likely. We performed a bootstrap analysis of the loglikelihood to quantify the significance level of the obtained Δ. Bootstrap catalogs were generated by randomly extracting stars N times from the observed set of N stars of the dataset (resampling with replacement). As suggested by Feigelson & Babu (2012), N_{B} ≈ N(lnN)^{2} bootstrap resamples were generated. For each bootstrap resample, the loglikelihood was computed for the two models and the loglikelihood difference Δ was calculated. Finally, after N_{B} resamples, the standard deviation σ_{Δ} of the distribution of Δ is determined, while the integral of the normalized Δ distribution for Δ > 0 gives P(Δ > 0), the probability of the warp model being favored over the nowarp model.
Difference of the loglikelihoods of the warp and nowarp models Δ ≡ ℓ_{WARP}−ℓ_{NOWARP} according to the warp parameters reported in Drimmel & Spergel (2001) (dust and stellar model) and Yusifov (2004).
Difference of the loglikelihoods of the warp (Yusifov 2004) and nowarp models Δ ≡ ℓ_{WARP}−ℓ_{NOWARP} for the TGAS(HIP2) sample.
Difference of the loglikelihoods of the warp (Yusifov 2004) and nowarp models Δ ≡ ℓ_{WARP}−ℓ_{NOWARP} for the HIP2 sample.
Table 3 collects the results for the TGAS(HIP2) dataset for the three different warp models whose parameters are given in Table 1, for the full dataset as well as for the two subsets of distant (ϖ < 1 mas) and nearby (2 >ϖ > 1 mas) stars. For the full dataset none of the warp models are favored over a nowarp model, though the model of Yusifov (2004) cannot be excluded. However, on splitting our sample into distant and nearby subsamples, we find some of the warp models are clearly favored over the nowarp model for the nearby stars.
In Tables 4 and 5 we show the results of our analysis, for both the HIP2 and TGAS(HIP2) samples, considering various subsamples with alternative data selection criteria to test for possible effects due to incompleteness and outliers. Separate PDFs were appropriately generated for the selections in magnitude and parallax. Here we show the results using the Yusifov (2004) model, as it is the most consistent with the data, as indicated by the the maximum likelihood measurements shown in Table 3. Again, the chosen warp model is not clearly favored nor disfavored until we split our sample into distant and nearby subsamples.
Various selection criteria were applied to investigate the role of possible outliers in biasing the outcome. We tried to remove the stars identified as binaries in the Hipparcos catalog (van Leeuwen 2008) and, for the TGAS subset, the objects with a high difference between the Gaia and Hipparcos proper motion (Lindegren et al. 2016). We also removed the high proper motion stars (i.e., the tails in the μ_{b} distributions), to exclude possible runaway stars or nearby objects with significant peculiar motions. As shown in Tables 4 and 5, the exclusion of these possible outliers does not change our findings, confirming that the warp model (Yusifov 2004) is preferred for the nearby objects, but rejected for the distant stars.
A further test was performed, removing the most obvious clumps in the l,b, and in the l,μ_{b} space (for example the one centered on l ≈ 80° and μ_{b} ≈ 2.5 mas/yr, see Fig. 17), to study the effect of the intrinsic clumpiness of the OB stars. We also removed the stars which are part of the known OB association Cen OB2 according to de Zeeuw et al. (1999). The obtained results (here not shown) are very similar to the ones in Tables 4 and 5.
5. Discussion
We have used models of the distribution and kinematics of OB stars to find the expected distribution of proper motions, including astrometric uncertainties, for two samples of spectroscopically identified OB stars from the new Hipparcos reduction and Gaia DR1. The resulting PDFs of the proper motions perpendicular to the Galactic plane produced by models with and without a warp are compared to the data via a likelihood analysis. We find that the observed proper motions of the nearby stars are more consistent with models containing a kinematic warp signature than a model without, while the more distant stars are not. Given that the warp signal in the proper motions is expected to remain evident at large distances (see Sect. 3.1), this result is difficult to reconcile with the hypothesis of a stable warp, and we are forced to consider alternative interpretations.
Keeping in mind that our sample of OB stars is tracing the gas, one possibility is that the warp in the gas starts well beyond the Solar Circle, or that the warp amplitude is so small that no signal is detectable. However, most studies to date suggest that the warp in the stars and in the dust starts inside or close to the Solar Circle, (Drimmel & Spergel 2001; Derriere & Robin 2001; Yusifov 2004; Momany et al. 2006; Robin et al. 2008), while the warp amplitude in the gas at R≈ 10 kpc (Levine et al. 2006) is consistent with warp models of sufficiently small amplitudes. Indeed, Momany et al. (2006) found an excellent agreement between the warp in the stars, gas and dust using the warp model of Yusifov (2004), the same model that we used in Sect. 4.2. Another possible scenario is that the warp of the Milky Way is a shortlived/transient feature, and that our model of a stable warp is not applicable. This hypothesis would be consistent with the finding that the warp structure may not be the same for all Milky Way components, as argued in Robin et al. (2008), but in contradiction to the findings of Momany et al. (2006) cited above. Finally, our expected kinematic signature from a stable warp could be overwhelmed, or masked, by other systematic motions. Evidence has been found for vertical oscillations (Widrow et al. 2012; Xu et al. 2015), suggesting the presence of vertical waves, as well as kinematic evidence of internal breathing modes (Williams et al. 2013) in the disk. Both have been attributed as being possibly caused by the passage of a satellite galaxy (Gómez et al. 2013; Widrow et al. 2014; Laporte et al. 2016), while breathing modes could also be caused by the bar and spiral arms (Monari et al. 2015, 2016). If such effects as these are present, then sampling over a larger volume of the Galactic disk will be necessary to disentangle the kinematic signature of the largescale warp from these other effects. Also, a comparison of the vertical motions of young stars (tracing the gas) and a dynamically old sample could also confirm whether the gas might possess additional motions due to other effects.
We have compared the proper motions of our two samples with the expectations from three warp models taken from the literature. Among these the model based on the farinfrared dust emission (Drimmel & Spergel 2001) can be excluded based on the kinematic data from Gaia DR1 that we present here. In addition, the Drimmel & Spergel (2001) dust warp model also predicts a vertical motion of the Local Standard of Rest of 4.6 km s^{1}, which would result in a vertical solar motion that is clearly inconsistent with the measured proper motion of Sag A^{∗} (Reid 2008). This calls into question the finding of Drimmel & Spergel (2001) that the warp in the dust and stars are significantly different. However, the proper motion data do not strongly favor the other two warp models, that of Yusifov (2004) and Drimmel & Spergel (2001) based on the stellar nearinfrared emission. As pointed out in Sect. 3.1, the local kinematic signature produced by a warp model with φ_{w} ≠ 0 is quite similar to that of a warp model with φ_{w} = 0 of smaller amplitude. In short, the parameters of even a simple symmetric warp cannot be constrained from local kinematics alone. In any case, we stress that the observed kinematics of the most distant OB stars are not consistent with any of the warp models.
6. Conclusion and future steps
Our search for a kinematic signature of the Galactic warp presented here is a preliminary study that adopts an exploratory approach, aimed at determining whether there is evidence in the Gaia DR1 and/or in the preGaia global astrometry of the new Hipparcos reduction. While unexpected, our finding that distant OB stars do not show evidence of the kinematic signature of the warp is in keeping with the previous results of Smart et al. (1998) and Drimmel et al. (2000), who analyzed the original Hipparcos proper motions using a simpler approach than employed here.
We point out that this work only considers a small fraction of the Gaia DR1 data with a full astrometric solution, being restricted to a subset of Hipparcos stars in DR1 brighter than m_{VT} = 8.5. Gaia DR1 TGAS astrometry is complete to about m_{VT} = 11, and potentially will permit us to sample a significantly larger volume of the disk of the Milky Way than presented here. In future work we will expand our sample to a fainter magnitude limit, using selection criteria based on multiwaveband photometry from other catalogs. We will also compare the kinematics of this young population to an older population representative of the relaxed stellar disk.
Understanding the dynamical nature of the Galactic warp will need studies of both its structural form as well as its associated kinematics. Gaia was constructed to reveal the dynamics of the Milky Way on a large scale, and we can only look forward to the future Gaia data releases that will eventually contain astrometry for over a billion stars. We expect that Gaia will allow us to fully characterize the dynamical properties of the warp, as suggested by Abedi et al. (2014, 2015), and allow us to arrive at a clearer understanding of the nature and origin of the warp. At the same time, Gaia may possibly reveal other phenomena causing systematic vertical velocities in the disk of the Milky Way.
Acknowledgments
This work has been partially funded by MIUR (Ministero dell’Istruzione, dell’Universitá e della Ricerca), through PRIN 2012 grant No. 1.05.01.97.02 “Chemical and dynamical evolution of our Galaxy and of the galaxies of the Local Group” and by ASI under contract No. 2014025R.1.2015 “Gaia Mission – The Italian Participation to DPAC”. E. Poggio acknowledges the financial support of the 2014 PhD fellowship programme of INAF. This work has made use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, http://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. The authors thank the referee for comments and suggestions that improved the overall quality of the paper. They also thank Scilla degl’Innocenti for the helpful discussions.
References
 Abedi, H., Figueras, F., Aguilar, L., et al. 2014, in EAS Pub. Ser., 67, 237 [Google Scholar]
 Abedi, H., Figueras, F., Aguilar, L., et al. 2015, in Highlights of Spanish Astrophysics VIII, eds. A. J. Cenarro, F. Figueras, C. HernándezMonteagudo, J. Trujillo Bueno, & L. Valdivielso, 423 [Google Scholar]
 Anderson, E., & Francis, C. 2012, Astron. Lett., 38, 331 [NASA ADS] [CrossRef] [Google Scholar]
 Bahcall, J. N., & Soneira, R. M. 1980, ApJS, 44, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Bahcall, J. N., Casertano, S., & Ratnatunga, K. U. 1987, ApJ, 320, 515 [NASA ADS] [CrossRef] [Google Scholar]
 BailerJones, C. A. L. 2015, PASP, 127, 994 [NASA ADS] [CrossRef] [Google Scholar]
 BlandHawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bobylev, V. V. 2010, Astron. Lett., 36, 634 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Bobylev, V. V. 2013, Astron. Lett., 39, 819 [NASA ADS] [CrossRef] [Google Scholar]
 Bobylev, V. V. 2015, Astron. Lett., 41, 156 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, A. G. A., de Geus, E. J., & de Zeeuw, P. T. 1994, A&A, 289, 101 [NASA ADS] [Google Scholar]
 Burke, B. F. 1957, AJ, 62, 90 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, Y., Bressan, A., Girardi, L., et al. 2015, MNRAS, 452, 1068 [NASA ADS] [CrossRef] [Google Scholar]
 Comeron, F., Torra, J., & Gomez, A. E. 1992, Ap&SS, 187, 187 [NASA ADS] [CrossRef] [Google Scholar]
 de Zeeuw, P. T., Hoogerwerf, R., de Bruijne, J. H. J., Brown, A. G. A., & Blaauw, A. 1999, AJ, 117, 354 [NASA ADS] [CrossRef] [Google Scholar]
 Dehnen, W., & Binney, J. J. 1998, MNRAS, 298, 387 [NASA ADS] [CrossRef] [Google Scholar]
 Derriere, S., & Robin, A. C. 2001, in The New Era of Wide Field Astronomy, eds. R. Clowes, A. Adamson, & G. Bromage, ASP Conf. Ser., 232, 229 [Google Scholar]
 Drimmel, R., & Spergel, D. N. 2001, ApJ, 556, 181 [NASA ADS] [CrossRef] [Google Scholar]
 Drimmel, R., Smart, R. L., & Lattanzi, M. G. 2000, A&A, 354, 67 [NASA ADS] [Google Scholar]
 Drimmel, R., CabreraLavers, A., & LópezCorredoira, M. 2003, A&A, 409, 205 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 ESA 1997, The Hipparcos and Tycho catalogues, Astrometric and photometric star catalogues derived from the ESA Hipparcos Space Astrometry Mission, ESA SP, 1200 [Google Scholar]
 Feigelson, E. D., & Babu, G. J. 2012, Modern Statistical Methods for Astronomy (Cambridge, UK: Cambridge University Press) [Google Scholar]
 Flower, P. J. 1977, A&A, 54, 31 [NASA ADS] [Google Scholar]
 Freudenreich, H. T., Berriman, G. B., Dwek, E., et al. 1994, ApJ, 429, L69 [NASA ADS] [CrossRef] [Google Scholar]
 Georgelin, Y. M., & Georgelin, Y. P. 1976, A&A, 49, 57 [NASA ADS] [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2016, A&A, 595, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gómez, F. A., Minchev, I., O’Shea, B. W., et al. 2013, MNRAS, 429, 159 [NASA ADS] [CrossRef] [Google Scholar]
 Guijarro, A., Peletier, R. F., Battaner, E., et al. 2010, A&A, 519, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Houk, N. 1993, VizieR Online Data Catalog: III/51 [Google Scholar]
 Houk, N. 1994, VizieR Online Data Catalog: III/380 [Google Scholar]
 Houk, N., & Cowley, A. P. 1994, VizieR Online Data Catalog: III/31 [Google Scholar]
 Houk, N., & SmithMoore, M. 1994, VizieR Online Data Catalog: III/133 [Google Scholar]
 Houk, N., & Swift, C. 2000, VizieR Online Data Catalog: III/214 [Google Scholar]
 Humphreys, R. M., & McElroy, D. B. 1984, ApJ, 284, 565 [NASA ADS] [CrossRef] [Google Scholar]
 Kalberla, P. M. W., Dedes, L., Kerp, J., & Haud, U. 2007, A&A, 469, 511 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kerr, F. J. 1957, AJ, 62, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Laporte, C. F. P., Gómez, F. A., Besla, G., Johnston, K. V., & GaravitoCamargo, N. 2016, MNRAS, submitted [arXiv:1608.04743] [Google Scholar]
 Levine, E. S., Blitz, L., & Heiles, C. 2006, ApJ, 643, 881 [NASA ADS] [CrossRef] [Google Scholar]
 Lindegren, L., & Kovalevsky, J. 1995, A&A, 304, 189 [NASA ADS] [Google Scholar]
 Lindegren, L., Lammers, U., Bastian, U., et al. 2016, A&A, 595, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 LópezCorredoira, M., CabreraLavers, A., Garzón, F., & Hammersley, P. L. 2002, A&A, 394, 883 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 LópezCorredoira, M., Abedi, H., Garzón, F., & Figueras, F. 2014, A&A, 572, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Maíz Apellániz, J., Sota, A., Walborn, N. R., et al. 2011, in Highlights of Spanish Astrophysics VI, eds. M. R. Zapatero Osorio, J. Gorgas, J. Maíz Apellániz, J. R. Pardo, & A. Gil de Paz, 467 [Google Scholar]
 Marshall, D. J., Robin, A. C., Reylé, C., Schultheis, M., & Picaud, S. 2006, A&A, 453, 635 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Martins, F., & Plez, B. 2006, A&A, 457, 637 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Martins, F., Schaerer, D., & Hillier, D. 2005, A&A, 436, 1049 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Michalik, D., Lindegren, L., Hobbs, D., & Lammers, U. 2014, A&A, 571, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Michalik, D., Lindegren, L., & Hobbs, D. 2015, A&A, 574, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Miyamoto, M., Yoshizawa, M., & Suzuki, S. 1988, A&A, 194, 107 [NASA ADS] [Google Scholar]
 Momany, Y., Zaggia, S., Gilmore, G., et al. 2006, A&A, 451, 515 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Monari, G., Famaey, B., & Siebert, A. 2015, MNRAS, 452, 747 [NASA ADS] [CrossRef] [Google Scholar]
 Monari, G., Famaey, B., Siebert, A., et al. 2016, MNRAS, 461, 3835 [NASA ADS] [CrossRef] [Google Scholar]
 Oort, J. H., Kerr, F. J., & Westerhout, G. 1958, MNRAS, 118, 379 [NASA ADS] [CrossRef] [Google Scholar]
 Peacock, J. A. 1983, MNRAS, 202, 615 [NASA ADS] [Google Scholar]
 Reed, B. C. 2001, PASP, 113, 537 [NASA ADS] [CrossRef] [Google Scholar]
 Reid, M. J. 2008, in A Giant Step: from Milli to Microarcsecond Astrometry, eds. W. J. Jin, I. Platais, & M. A. C. Perryman, IAU Symp., 248, 141 [Google Scholar]
 Reshetnikov, V., & Combes, F. 1998, A&A, 337, 9 [NASA ADS] [Google Scholar]
 Reylé, C., Marshall, D. J., Robin, A. C., & Schultheis, M. 2009, A&A, 495, 819 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Robin, A. C., Reylé, C., & Marshall, D. J. 2008, Astron. Nachr., 329, 1012 [NASA ADS] [CrossRef] [Google Scholar]
 Roeser, S., Demleitner, M., & Schilbach, E. 2010, AJ, 139, 2440 [NASA ADS] [CrossRef] [Google Scholar]
 SanchezSaavedra, M. L., Battaner, E., & Florido, E. 1990, MNRAS, 246, 458 [NASA ADS] [Google Scholar]
 Scalo, J. M. 1986, Fund. Cosmic Phys., 11, 1 [NASA ADS] [EDP Sciences] [Google Scholar]
 Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829 [NASA ADS] [CrossRef] [Google Scholar]
 Sellwood, J. A. 2013, Dynamics of Disks and Warps, eds. T. D. Oswalt, & G. Gilmore, 923 [Google Scholar]
 Skiff, B. A. 2014, VizieR Online Data Catalog, 1 [Google Scholar]
 Smart, R. L., Drimmel, R., Lattanzi, M. G., & Binney, J. J. 1998, Nature, 392, 471 [NASA ADS] [CrossRef] [Google Scholar]
 Sota, A., Maíz Apellániz, J., Walborn, N. R., et al. 2011, ApJS, 193, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Sota, A., Maíz Apellániz, J., Morrell, N. I., et al. 2014, ApJS, 211, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Taylor, J. H., & Cordes, J. M. 1993, ApJ, 411, 674 [NASA ADS] [CrossRef] [Google Scholar]
 van Leeuwen, F. 2008, VizieR Online Data Catalog: III/311 [Google Scholar]
 Westerhout, G. 1957, Bull. Astron. Inst. Netherlands, 13, 201 [Google Scholar]
 Widrow, L. M., Gardner, S., Yanny, B., Dodelson, S., & Chen, H.Y. 2012, ApJ, 750, L41 [NASA ADS] [CrossRef] [Google Scholar]
 Widrow, L. M., Barber, J., Chequers, M. H., & Cheng, E. 2014, MNRAS, 440, 1971 [NASA ADS] [CrossRef] [Google Scholar]
 Williams, M. E. K., Steinmetz, M., Binney, J., et al. 2013, MNRAS, 436, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Xu, Y., Newberg, H. J., Carlin, J. L., et al. 2015, ApJ, 801, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Yusifov, I. 2004, in The Magnetized Interstellar Medium, eds. B. Uyaniker, W. Reich, & R. Wielebinski, 165 [Google Scholar]
 Zacharias, N., Finch, C. T., Girard, T. M., et al. 2013, AJ, 145, 44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Hipparcos astrometric errors
Tables A.1–A.5 show the median formal errors of right ascension α, declination δ, parallax ϖ, proper motion components μ_{α} and μ_{δ} in function of apparent magnitude and ecliptic latitude for the HIP2 stars. They were obtained considering the entire HIP2 catalog (as given by van Leeuwen 2008) excluding the stars redder than (B−V) = 0.5. We also excluded stars for which there was a claim of binarity in van Leeuwen (2008), taking account for binary systems after the single stars errors are generated, as described in the text. To construct the tables, we binned the resulting sample of 15 197 HIP2 stars with respect to apparent magnitude and ecliptic latitude and found the median errors for each bin. Table A.6 shows the number of objects in each bin.
Median formal uncertainties for σ_{α}.
Median formal uncertainties for σ_{δ}.
Median formal uncertainties for σ_{ϖ}.
Median formal uncertainties for σ_{μα ∗}.
Median formal uncertainties for σ_{μδ}.
Number of stars in each bin.
All Tables
Comparison of warp parameters for the models of Drimmel & Spergel (2001) and Yusifov (2004).
Difference of the loglikelihoods of the warp and nowarp models Δ ≡ ℓ_{WARP}−ℓ_{NOWARP} according to the warp parameters reported in Drimmel & Spergel (2001) (dust and stellar model) and Yusifov (2004).
Difference of the loglikelihoods of the warp (Yusifov 2004) and nowarp models Δ ≡ ℓ_{WARP}−ℓ_{NOWARP} for the TGAS(HIP2) sample.
Difference of the loglikelihoods of the warp (Yusifov 2004) and nowarp models Δ ≡ ℓ_{WARP}−ℓ_{NOWARP} for the HIP2 sample.
All Figures
Fig. 1 Our final sample of Hipparcos OB stars on the sky, plotted in Galactic coordinates. The dashed line shows the orientation of the Gould belt according to Comeron et al. (1992). Colored points indicate the stars that are identified members of the OB associations Orion OB1 (red), Trumpler 10 (purple), Vela OB2 (blue), Collinder 121 (green) and Lacerta OB1 (cyan). 

In the text 
Fig. 2 According to the warp model, the true μ_{b} in the Galactic plane as a function of Galactic longitude at heliocentric distances of 0.5 kpc (A) and 1.5 kpc (B). For each set of curves, the thick line represents the case with warp phase φ_{w} = 0° and the two thin curves show φ_{w} = ± 20°. 

In the text 
Fig. 3 Color–magnitude diagram used to generate synthetic intrinsic colors. The dark and light orange regions show, respectively, the main sequence and the giant regions. The density of the two regions (here not shown) depends on the ILF and on the giant fraction. 

In the text 
Fig. 4 Modeled surface density of the OB stars. Sun’s position is indicated by the star. 

In the text 
Fig. 5 Latitude distribution of the HIP2 OB stars. The red curve is a nonparametric fit to the selected HIP2 sample, the red dashed curve shows the additional contribution of the Orion OB1 association, while the red dotted curve shows the added contributions of the Trumpler 10, Vela OB2, Collinder 121, and Lacerta OB1 associations. The blue curve and lightblue shaded area shows the average and 2σ confidence band of the simulated longitude distribution, based on 30 simulated instances of the sample. The black dotted and dashdot curves show the relative contributions of the major spiral arms (SagittariusCarina and Perseus) and the local arm, respectively, while the additional black solid curve is for a model without a warp. 

In the text 
Fig. 6 Longitude distribution of the HIP2 OB stars. Meaning of the curves are as in the previous figure. 

In the text 
Fig. 7 Apparent magnitude distribution for the data (histogram) and the simulations (black dots). The error bars show 2σ uncertainty, calculated with 30 simulated samples. 

In the text 
Fig. 8 Fraction of Hipparcos completeness in function of apparent magnitude V_{T} with respect to the Tycho2 catalog. The dashed and the dotted line represent, respectively, the Hipparcos fraction for the stars above and below δ = −30°. 

In the text 
Fig. 9 Fraction of HIP2 OB stars present in HIPTGAS as a function of the observed color and the apparent magnitude. 

In the text 
Fig. 10 μ_{l} in function of Galactic longitude for the data (red curve) together with the 95% bootstrap confidence band (pink shaded area). The three black dotted curves show the trend obtained with simulations with circular velocity 260, 238 and 220 km s^{1}, respectively, from the lowest to the highest curve. Simulations with an additional velocity to the Local Arm (see text) produce the blue curve, for which the light blue shaded area shows the 2σ confidence band, calculated with 30 simulated catalogs. 

In the text 
Fig. 11 Histogram of the observed parallax distribution. The dashed and the solid curves show, respectively, the synthetic distributions with F = 1 and 1.5 (see text for explanation). 

In the text 
Fig. 12 Logarithm of the parallax errors (mas) in function of ecliptic latitude for the Hipparcos subsample in TGAS. The points show the medians, while the error bars show the 10th and the 90th percentiles of the distribution. 

In the text 
Fig. 13 Distribution of the logarithm of the parallax uncertainties (mas) for the HIPTGAS stars. Three subsets with different ecliptic latitude are shown. 

In the text 
Fig. 14 Comparison of the published error σ_{μα ∗}(TGAS) to the prediction based on Hipparcos uncertainties for each star of the Hipparcos subset in TGAS (see text). The dashed line represents the bisector. The solid line has null intercept and coefficient C_{α} = 1.42, which is used to calibrate our error model (see text). 

In the text 
Fig. 15 Median correlations between parallaxes ϖ and proper motions μ_{α}, μ_{δ} in HIPTGAS. 

In the text 
Fig. 16 Bivariate NadarayaWatson regression estimator of stellar vertical velocities as a function of Galactocentric radius, using a bandwidth of h = 0.5 kpc. The same regression bandwidth has been used for the data (red), the nowarp model (black) and warp model of Yusifov (2004) (blue). The 95% bootstrap confidence band is shown for the data. For each of the two models, the non parametric regressions are performed for 20 simulated catalogs, obtaining the curves as the mean values and the shaded areas as the 95% uncertainty. 

In the text 
Fig. 17 Distribution of the TGASHIP data in the lμ_{b} plane (black dots), together with the probability density function P(μ_{b}  l) predicted by the nowarp model (left panel) and the warp model of Yusifov (2004; right panel). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.