The kinematic signature of the Galactic warp inGaiaDR1

A&A 601, A115 (2017) DOI: 10.1051/0004-6361/201629916 Astronomy & Astrophysics c ESO 2017 The kinematic signature of the Galactic warp in Gaia DR1 I. The Hipparcos subsample E. Poggio1, 2 , R. Drimmel2 , R. L. Smart2, 3 , A. Spagna2 , and M. G. Lattanzi2 1 2 3 Università di Torino, Dipartimento di Fisica, via P. Giuria 1, 10125 Torino, Italy e-mail: [email protected] Osservatorio Astrofisico di Torino, Istituto Nazionale di Astrofisica (INAF), Strada Osservatorio 20, 10025 Pino Torinese, Italy School of Physics, Astronomy and Mathematics, University of Hertfordshire, College Lane, Hatfield AL10 9AB, UK Received 17 October 2016 / Accepted 10 February 2017 ABSTRACT 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 long-lived 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 spectroscopically-identified 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 long-lived 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 long-lived 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 1. Introduction It has been known since the early HI 21-cm 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ópez-Corredoira 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 (Sanchez-Saavedra 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 line-of-sight. Being located within the disk of the Milky Way, systematic vertical motions will primarily manifest themselves to us in the direction perpendicular to our line-of-sight. 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 stellar-gaseous warp and the kinematics of their tracers, namely nearby red clump giants from Tycho2 and Cepheids with UCAC4 (Fourth US Naval Observatory Article published by EDP Sciences A115, page 1 of 14 A&A 601, A115 (2017) 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ópez-Corredoira et al. (2014) concluded that the data might be consistent with a long-lived warp, though they admit that smaller systematic errors in the proper motions are needed to confirm this tentative finding. Indeed, large-scale systematic errors in the ground-based proper motions compromise efforts to detect the Galactic warp. The first real hopes of overcoming such systematics came with global space-based 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 all-sky 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 long-lived warp. We choose the OB stars as they are intrinsically bright, thus can be seen to large distances, and are short-lived, 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 long-lived 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 pre-Gaia 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 A115, page 2 of 14 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 short-lived, 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 O-star 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 BJ magnitude <12. From the catalog presented in Sota et al. (2014), which is complete to BJ = 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 & Smith-Moore 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; Tycho-Gaia 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. 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 b (deg) E. Poggio et al.: The kinematic signature of the Galactic warp in Gaia DR1. I. 40 Table 1. Comparison of warp parameters for the models of Drimmel & Spergel (2001) and Yusifov (2004). 20 Rw (kpc) αw h0 (kpcαw −1 ) φw (◦ ) Drimmel & Spergel (2001), dust 7 2 0.073 0 Drimmel & Spergel (2001), stars 7 2 0.027 0 Yusifov (2004) 6.27 1.4 0.037 14.5 0 Notes. The radius Rw was scaled to account for different assumptions about the Sun – Galactic center distance in this work and in the considered papers. −20 −40 300 200 l (deg) 100 0 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). 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 spectro-photometric 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 l-b 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. proper motions (Sect. 3.5) of our two samples, including (or not) the expected effects of a stable (long-lived) non-precessing warp (Sect. 3.1). A comparison of the observed samples with the expectations from the different warp/no-warp 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. 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 flat-disk distribution with vertical exponential profile, then displace the z-coordinates by zw , where: zw (R, φ) = h(R) sin (φ + φw ). The warp phase angle φw determines the position of the line-ofnodes 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 h(R) = h0 (R − Rw )αw , 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 three-dimensional 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 (1) (2) where h0 and Rw 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 v̄z (R, φ). If we use the warped disk as described above and suppose that v̄R = 0 (i.e. that the disc is not radially expanding or collapsing), we obtain: v̄z (R, φ) = v̄φ h(R) cos (φ + φw ). R (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., A115, page 3 of 14 A&A 601, A115 (2017) −7 1 −6 µb (mas/yr) 0 −5 (B) −2 Mv −1 −4 (A) −3 −3 −2 360 270 180 l (deg) 90 0 −0.30 −0.25 −0.20 (B−V)0 −0.15 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◦ . 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. “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 zw = 0 will be the no-warp model, where Eq. (3) reduces to the trivial v̄z = 0. Figure 2 shows the prediction of a warp model, without errors, for the mean proper motions µb in the Galactic plane. For the no-warp 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 anti-center 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. giant fraction fg has been modeled as a function of the absolute magnitude as follows:   if M ≤ −7   1,  fg (M) =  −0.25M − 0.75, if − 7 ≤ M ≤ −4    0.25, if M ≥ −4, 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 A115, page 4 of 14 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. 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 birth-places. In our model we adopt R0 = 8.2 kpc as the Sun’s distance from the Galactic center, and a solar offset from the disk midplane of z0 = 25 pc, for galactocentric cylindrical coordinates (R, φ, z), as recommended by Bland-Hawthorn & Gerhard (2016). For the spiral arm geometry we adopt the model of Georgelin & Georgelin (1976), as implemented by Taylor & Cordes (1993), rescaled to R0 = 8.2 kpc, with the addition of a local arm described as a logarithmic spiral segment whose location is described by RLoc = RLoc,0 exp −(tan pφ), p being the arm’s pitch angle. The surface density profile across an arm is taken to be gaussian, namely: ρ ∝ exp (−da2 /w2a ), where da is the distance to the nearest arm in the R, φ plane, and wa = cw R is the arm half-width, with cw = 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′ |/hz ), where hz is the vertical scale height, z′ = z − zw − zLoc , zw (R, φ) HIP data (mv < 7.5) Simulation No Warp Simulation Warp 0 5 Stars / deg 10 15 E. Poggio et al.: The kinematic signature of the Galactic warp in Gaia DR1. I. Fig. 4. Modeled surface density of the OB stars. Sun’s position is indicated by the star. 1 https://www.r-project.org −20 −10 0 10 20 30 b (deg) 0.20 Fig. 5. Latitude distribution of the HIP2 OB stars. The red curve is a non-parametric 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 light-blue 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 dash-dot curves show the relative contributions of the major spiral arms (Sagittarius-Carina and Perseus) and the local arm, respectively, while the additional black solid curve is for a model without a warp. 0.10 0.05 Stars / deg 0.15 HIP data (mv < 7.5) Simulation 0.00 is the height of the warp as described in Sect. 3.1 and zLoc is a vertical offset applied only to the local arm. We generate the above spatial distribution in an iterative Monte-Carlo 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 (MV , (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 TGAS-like 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 hz = 70 pc and assuming zLoc = 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 RLoc,0 = 8.3 kpc, with a pitch angle of 6.5◦ and a half-width of 500 pc. (The curves in Figs. 5 and 6 are non-parametric fits to the distributions obtained through kernel density estimation with a Gaussian kernel, as implemented by the generic function density in R1 . 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). 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 Sagittarius-Carina arm, which we have not attempted to modify as we are primarily interested in the kinematics toward the Galactic anti-center. −30 360 315 270 225 180 135 90 45 0 l (deg) Fig. 6. Longitude distribution of the HIP2 OB stars. Meaning of the curves are as in the previous figure. 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. 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 co-rotation, as expected, such offsets are minimal. A115, page 5 of 14 A&A 601, A115 (2017) 5 0 0.2 0.4 5 6 Log ( N ) 100 0.6 Vmag Vmag 10 7 0.8 8 15 10 9 20 1 4 5 6 mv 7 10 25 8 0 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. 0.05 100.5 15 (B−V)obs (B−V) 1.0 20 1.5 25 Fig. 9. Fraction of HIP2 OB stars present in HIP-TGAS as a function of the observed color and the apparent magnitude. 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. 1.0 0.8 fHIP 0.6 0.4 0.2 3.5. Kinematics 0.0 5 6 7 8 9 10 11 12 VTmag Fig. 8. Fraction of Hipparcos completeness in function of apparent magnitude VT with respect to the Tycho-2 catalog. The dashed and the dotted line represent, respectively, the Hipparcos fraction for the stars above and below δ = −30◦ . 3.4. Completeness The completeness of the Hipparcos catalog decreases with apparent magnitude, as shown by the fraction of Tycho-2 stars in the Hipparcos catalog (Fig. 8). At VT ≤ 7.5 the HIP2 catalog is approximately 90% complete, reaching approximately 50% completeness at VT = 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 VT decreases in a similar way for the Johnson magnitude, that is fHIP (VT ) ≈ fHIP (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 A115, page 6 of 14 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 warp-induced 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: σ(1,2,3) = σ0(1,2,3) exp (R⊙ − R)/2hR , where hR = 2.3 kpc is the radial scale length and σ0(1,2,3) = (14.35, 9.33, 5.45) km s−1 are the three velocity dispersions in the solar neighborhood for the bluest stars (Dehnen & Binney 1998). A vertex deviation of lv = 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 Bland-Hawthorn & Gerhard (2016), we adopted Θ0 = 238 km s−1 for the circular rotation velocity at the Solar radius R0 . 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 E. Poggio et al.: The kinematic signature of the Galactic warp in Gaia DR1. I. -2 -4 -8 -6 μl (mas/yr) 0 2 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 non-zero 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). 360 270 180 l (deg) 90 0 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. are converted to Galactic coordinates, that is (µl ∗, µb ), where µl ∗ = µl cos b. Figure 10 shows the derived proper motions in Galactic longitude for both the data and simulations using the bivariate localconstant (i.e. Nadaraya-Watson) 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 dash-dotted 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 ∆VC = 6 km s−1 in the direction of Galactic rotation and ∆VR = 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 post-shock 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. 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 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 signal-to-noise 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. 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 A115, page 7 of 14 A&A 601, A115 (2017) 150 0 ¡0:1 50 log10(¾$) Frequency 100 ¡0:2 ¡0:3 ¡0:4 ¡0:5 0 ¡0:6 −3 −2 −1 ϖ (mas) 0 1 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). 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 fb ≈ 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 fbin = 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 A115, page 8 of 14 0 2 10 20 30 40 50 j¯j(deg) 60 70 80 90 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. 6000 j¯j < 40 40 < j¯j < 60 5000 j¯j > 60 4000 N −4 3000 2000 1000 0 ¡0:7 ¡0:6 ¡0:5 ¡0:4 ¡0:3 log10(¾$) ¡0:2 ¡0:1 0 Fig. 13. Distribution of the logarithm of the parallax uncertainties (mas) for the HIP-TGAS stars. Three subsets with different ecliptic latitude are shown. Table 2. Shape parameter k and scale parameter θ of the gamma distributions used to model the log10 (σ̟ ) distributions shown in Fig. 12. Ecliptic latitude |β| (deg) 0–40 40–60 60–90 k 1.5 1.2 1.1 θ 0.113 0.115 0.08 Offset –0.658 –0.658 –0.67 Notes. An additional offset is required in order to fit the distributions. 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. E. Poggio et al.: The kinematic signature of the Galactic warp in Gaia DR1. I. 0:5 0:4 0 0:3 ¡0:5 0:2 ¡1:0 1:0 0:1 0 ½ ( ¹®¤ ; ¹± ) 0:5 0:5 0 0:05 0:10 0:15 0:20 F ¢ ¾®(HIP2)=¢t (mas=yr) 0:25 Fig. 14. Comparison of the published error σµα∗ (TGAS) to the prediction based on Hipparcos uncertainties F σαH (m, β)/∆t 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). 0 ½ ( $ ; ¹®¤ ) ¾¹®(TGAS)( mas=yr) 1:0 ¡0:5 ¡1:0 1:0 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 σµα =   Cα F σαH (m, β)/∆t , where F is the correction factor applied to the Hipparcos astrometric uncertainties, as described in Sect. 3.6.1, σαH (m, β) 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 σµα / F σαH (m, β)/∆t for the stars of the TGAS(HIP2) sample. An analogous model is used for σµδ , with Cδ = 1.44. 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. 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 0 ½ ( $ ; ¹± ) 0:5 ¡0:5 ¡1:0 Fig. 15. Median correlations between parallaxes ̟ and proper motions µα , µδ in HIP-TGAS. 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 vz as a function of Galactocentic radius R that one would derive from the measured proper motions and spectro-photometric distances, comparing the data with what we expect from the no-warp and warp models. In the no-warp case, the true mean vertical velocities are zero, while a warp model predicts that they increase with R outside the radius Rw 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 no-warp 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 no-warp model, where the true vz (R) = 0 (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 A115, page 9 of 14 A&A 601, A115 (2017) 0 -10 -5 Vz - Vz(RSun) 5 HIP data (mv < 7.5) Simulation no warp Simulation with warp 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 R (kpc) Fig. 16. Bivariate Nadaraya-Watson 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 no-warp 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. the data was consistent with the no-warp 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 (Bailer-Jones 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. 4.2. Proper motion µb in function of Galactic longitude The aim of the present work is to determine whether a warp is favored over a no-warp model in either of our Hipparcos or Gaia A115, page 10 of 14 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 anti-center (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 no-warp 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 two-dimensional 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 R 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 Rthe 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 |li ) associated with each ith observed star according Q to each PDF. The likelihood associated N with the model is L = i=1 P(µb,i |li ), where N is the total number of stars in our dataset; for computational PN reasons, we used instead the log-likelihood ℓ = ln(L) = i=1 ln(P(µb,i |li ). 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 no-warp 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 log-likelihood 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), NB ≈ N(ln N)2 bootstrap resamples were generated. For each bootstrap resample, the log-likelihood was computed for the two models and the log-likelihood difference ∆ was calculated. Finally, after NB 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 no-warp model. 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 no-warp model, though the model of Yusifov (2004) cannot be excluded. E. Poggio et al.: The kinematic signature of the Galactic warp in Gaia DR1. I. Fig. 17. Distribution of the TGAS-HIP data in the l-µb plane (black dots), together with the probability density function P(µb |l) predicted by the no-warp model (left panel) and the warp model of Yusifov (2004; right panel). Table 3. Difference of the log-likelihoods of the warp and no-warp models ∆ ≡ ℓWARP − ℓNOWARP according to the warp parameters reported in Drimmel & Spergel (2001) (dust and stellar model) and Yusifov (2004). Warp model Drimmel & Spergel (2001), dust Drimmel & Spergel (2001), stars Yusifov (2004) ∆ –41.24 –5.47 –2.69 ̟ < 2 mas σ∆ P(∆ > 0) 9.72 0.00 3.93 0.09 6.99 0.35 (1 < ̟ < 2) mas ∆ σ∆ P(∆ > 0) –2.94 5.04 0.28 3.13 1.92 0.95 13.93 3.47 1.00 ∆ –60.35 –22.54 –10.48 ̟ < 1 mas σ∆ P(∆ > 0) 10.90 0.00 16.79 0.04 25.81 0.32 Notes. Log-likelihoods are calculated with the TGAS(HIP2) sample (̟ < 2 mas), containing 758 stars. We also show the results for the nearby ((1 < ̟ < 2) mas, 296 stars) and for the distant objects (̟ < 1 mas, 462 stars). The standard deviation σ∆ and the probability P(∆ > 0) are calculated using bootstrap resamples (see text). Table 4. Difference of the log-likelihoods of the warp (Yusifov 2004) and no-warp models ∆ ≡ ℓWARP − ℓNOWARP for the TGAS(HIP2) sample. Sample All mv < 7.5 ∆Q < ∆Q95% ∆Q < ∆Q90% No HIP2 binaries 95% Conf. Int. 90% Conf. Int. Nstars 758 310 749 690 672 730 692 ̟ < 2 mas ∆ σ∆ –2.69 6.99 0.68 4.44 –3.73 6.90 –7.35 6.56 –2.75 6.81 –0.86 5.94 2.79 4.91 P(∆ > 0) 0.35 0.57 0.30 0.13 0.34 0.44 0.72 N stars 296 129 293 257 267 289 273 (1 < ̟ < 2) mas ∆ σ∆ P(∆ > 0) 13.93 3.47 1.00 4.59 2.09 0.99 13.86 3.50 1.00 10.97 3.07 1.00 15.03 3.40 1.00 10.30 3.23 1.00 9.11 2.71 1.00 Nstars 462 181 456 433 405 444 420 ̟ < 1 mas ∆ σ∆ –10.48 25.81 –0.47 4.03 –11.54 25.89 –12.18 25.93 –11.85 26.27 –13.99 5.23 –10.85 4.13 P(∆ > 0) 0.32 0.45 0.31 0.30 0.31 0.00 0.00 Notes. Results are shown for the whole sample (All, mv < 8.5) and for the bright sample (mv < 7.5). We also report the results obtained removing the objects with high ∆Q, where ∆Q is the difference between the TGAS and Hipparcos proper motion (Lindegren et al. 2016). For the Hipparcos subset in TGAS, the 95th and 90th percentiles of the ∆Q distribution are ∆Q95% = 23 and ∆Q90% = 11. We also present the results excluding the stars labeled as binaries in van Leeuwen (2008). The 95% (90%) confidence interval is obtained considering the stars with proper motions µb between the 2.5th and the 97.5th (the 5th and the 95th) percentiles of the whole µb distribution (with 767 stars), without restricting to the range (−10 < µb < 5) mas/yr (see text). The standard deviation σ∆ and the probability P(∆ > 0) are calculated using bootstrap resamples (see text). However, on splitting our sample into distant and nearby subsamples, we find some of the warp models are clearly favored over the no-warp 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 A115, page 11 of 14 A&A 601, A115 (2017) Table 5. Difference of the log-likelihoods of the warp (Yusifov 2004) and no-warp models ∆ ≡ ℓWARP − ℓNOWARP for the HIP2 sample. HIP2 subset All All No HIP2 binaries mv < 7.5 mv < 7.5 No HIP2 binaries Nstars 989 838 498 404 ∆ –14.99 –11.86 –1.68 –2.29 σ∆ 6.50 5.74 4.47 4.21 P(∆ > 0) 0.01 0.02 0.35 0.29 Notes. Results are shown for the whole sample (All, mv < 8.5) and for the bright sample (mv < 7.5). We also present the results obtained removing the stars labeled as binaries in van Leeuwen (2008). The standard deviation σ∆ and the probability P(∆ > 0) are calculated using bootstrap resamples (see text). 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 short-lived/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 A115, page 12 of 14 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 large-scale 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 far-infrared 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 pre-Gaia 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 mVT = 8.5. Gaia DR1 TGAS astrometry is complete to about mVT = 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 multi-waveband 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 E. Poggio et al.: The kinematic signature of the Galactic warp in Gaia DR1. I. 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. Acknowledgements. 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. 2014-025R.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 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 Anderson, E., & Francis, C. 2012, Astron. Lett., 38, 331 Bahcall, J. N., & Soneira, R. M. 1980, ApJS, 44, 73 Bahcall, J. N., Casertano, S., & Ratnatunga, K. U. 1987, ApJ, 320, 515 Bailer-Jones, C. A. L. 2015, PASP, 127, 994 Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529 Bobylev, V. V. 2010, Astron. Lett., 36, 634 Bobylev, V. V. 2013, Astron. Lett., 39, 819 Bobylev, V. V. 2015, Astron. Lett., 41, 156 Brown, A. G. A., de Geus, E. J., & de Zeeuw, P. T. 1994, A&A, 289, 101 Burke, B. F. 1957, AJ, 62, 90 Chen, Y., Bressan, A., Girardi, L., et al. 2015, MNRAS, 452, 1068 Comeron, F., Torra, J., & Gomez, A. E. 1992, Ap&SS, 187, 187 de Zeeuw, P. T., Hoogerwerf, R., de Bruijne, J. H. J., Brown, A. G. A., & Blaauw, A. 1999, AJ, 117, 354 Dehnen, W., & Binney, J. J. 1998, MNRAS, 298, 387 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 Drimmel, R., & Spergel, D. N. 2001, ApJ, 556, 181 Drimmel, R., Smart, R. L., & Lattanzi, M. G. 2000, A&A, 354, 67 Drimmel, R., Cabrera-Lavers, A., & López-Corredoira, M. 2003, A&A, 409, 205 ESA 1997, The Hipparcos and Tycho catalogues, Astrometric and photometric star catalogues derived from the ESA Hipparcos Space Astrometry Mission, ESA SP, 1200 Feigelson, E. D., & Babu, G. J. 2012, Modern Statistical Methods for Astronomy (Cambridge, UK: Cambridge University Press) Flower, P. J. 1977, A&A, 54, 31 Freudenreich, H. T., Berriman, G. B., Dwek, E., et al. 1994, ApJ, 429, L69 Georgelin, Y. M., & Georgelin, Y. P. 1976, A&A, 49, 57 Gaia Collaboration (Brown, A. G. A., et al.) 2016, A&A, 595, A2 Gómez, F. A., Minchev, I., O’Shea, B. W., et al. 2013, MNRAS, 429, 159 Guijarro, A., Peletier, R. F., Battaner, E., et al. 2010, A&A, 519, A53 Houk, N. 1993, VizieR Online Data Catalog: III/51 Houk, N. 1994, VizieR Online Data Catalog: III/380 Houk, N., & Cowley, A. P. 1994, VizieR Online Data Catalog: III/31 Houk, N., & Smith-Moore, M. 1994, VizieR Online Data Catalog: III/133 Houk, N., & Swift, C. 2000, VizieR Online Data Catalog: III/214 Humphreys, R. M., & McElroy, D. B. 1984, ApJ, 284, 565 Kalberla, P. M. W., Dedes, L., Kerp, J., & Haud, U. 2007, A&A, 469, 511 Kerr, F. J. 1957, AJ, 62, 93 Laporte, C. F. P., Gómez, F. A., Besla, G., Johnston, K. V., & Garavito-Camargo, N. 2016, MNRAS, submitted [arXiv:1608.04743] Levine, E. S., Blitz, L., & Heiles, C. 2006, ApJ, 643, 881 Lindegren, L., & Kovalevsky, J. 1995, A&A, 304, 189 Lindegren, L., Lammers, U., Bastian, U., et al. 2016, A&A, 595, A4 López-Corredoira, M., Cabrera-Lavers, A., Garzón, F., & Hammersley, P. L. 2002, A&A, 394, 883 López-Corredoira, M., Abedi, H., Garzón, F., & Figueras, F. 2014, A&A, 572, A101 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 Marshall, D. J., Robin, A. C., Reylé, C., Schultheis, M., & Picaud, S. 2006, A&A, 453, 635 Martins, F., & Plez, B. 2006, A&A, 457, 637 Martins, F., Schaerer, D., & Hillier, D. 2005, A&A, 436, 1049 Michalik, D., Lindegren, L., Hobbs, D., & Lammers, U. 2014, A&A, 571, A85 Michalik, D., Lindegren, L., & Hobbs, D. 2015, A&A, 574, A115 Miyamoto, M., Yoshizawa, M., & Suzuki, S. 1988, A&A, 194, 107 Momany, Y., Zaggia, S., Gilmore, G., et al. 2006, A&A, 451, 515 Monari, G., Famaey, B., & Siebert, A. 2015, MNRAS, 452, 747 Monari, G., Famaey, B., Siebert, A., et al. 2016, MNRAS, 461, 3835 Oort, J. H., Kerr, F. J., & Westerhout, G. 1958, MNRAS, 118, 379 Peacock, J. A. 1983, MNRAS, 202, 615 Reed, B. C. 2001, PASP, 113, 537 Reid, M. J. 2008, in A Giant Step: from Milli- to Micro-arcsecond Astrometry, eds. W. J. Jin, I. Platais, & M. A. C. Perryman, IAU Symp., 248, 141 Reshetnikov, V., & Combes, F. 1998, A&A, 337, 9 Reylé, C., Marshall, D. J., Robin, A. C., & Schultheis, M. 2009, A&A, 495, 819 Robin, A. C., Reylé, C., & Marshall, D. J. 2008, Astron. Nachr., 329, 1012 Roeser, S., Demleitner, M., & Schilbach, E. 2010, AJ, 139, 2440 Sanchez-Saavedra, M. L., Battaner, E., & Florido, E. 1990, MNRAS, 246, 458 Scalo, J. M. 1986, Fund. Cosmic Phys., 11, 1 Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829 Sellwood, J. A. 2013, Dynamics of Disks and Warps, eds. T. D. Oswalt, & G. Gilmore, 923 Skiff, B. A. 2014, VizieR Online Data Catalog, 1 Smart, R. L., Drimmel, R., Lattanzi, M. G., & Binney, J. J. 1998, Nature, 392, 471 Sota, A., Maíz Apellániz, J., Walborn, N. R., et al. 2011, ApJS, 193, 24 Sota, A., Maíz Apellániz, J., Morrell, N. I., et al. 2014, ApJS, 211, 10 Taylor, J. H., & Cordes, J. M. 1993, ApJ, 411, 674 van Leeuwen, F. 2008, VizieR Online Data Catalog: III/311 Westerhout, G. 1957, Bull. Astron. Inst. Netherlands, 13, 201 Widrow, L. M., Gardner, S., Yanny, B., Dodelson, S., & Chen, H.-Y. 2012, ApJ, 750, L41 Widrow, L. M., Barber, J., Chequers, M. H., & Cheng, E. 2014, MNRAS, 440, 1971 Williams, M. E. K., Steinmetz, M., Binney, J., et al. 2013, MNRAS, 436, 101 Xu, Y., Newberg, H. J., Carlin, J. L., et al. 2015, ApJ, 801, 105 Yusifov, I. 2004, in The Magnetized Interstellar Medium, eds. B. Uyaniker, W. Reich, & R. Wielebinski, 165 Zacharias, N., Finch, C. T., Girard, T. M., et al. 2013, AJ, 145, 44 A115, page 13 of 14 A&A 601, A115 (2017) Table A.4. Median formal uncertainties for σµα∗ . 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. Table A.1. Median formal uncertainties for σα . mV 3–4 4–5 5–6 6–7 7–8 0–10 0.22 0.24 0.34 0.49 0.66 10–20 0.19 0.25 0.32 0.46 0.64 Ecliptic latitude (|β|, (deg)) 20–30 30–40 40–50 0.40 0.23 0.12 0.25 0.20 0.15 0.31 0.27 0.21 0.44 0.37 0.29 0.60 0.51 0.40 50–60 0.11 0.16 0.19 0.29 0.39 60–90 0.20 0.15 0.20 0.29 0.40 Table A.2. Median formal uncertainties for σδ . mV 3–4 4–5 5–6 6–7 7–8 0–10 0.16 0.16 0.23 0.32 0.44 10–20 0.12 0.18 0.22 0.32 0.44 Ecliptic latitude (|β|, (deg)) 20–30 30–40 40–50 0.31 0.19 0.12 0.18 0.17 0.17 0.23 0.24 0.23 0.33 0.32 0.32 0.44 0.44 0.45 50–60 0.13 0.17 0.21 0.31 0.42 60–90 0.26 0.16 0.21 0.29 0.40 50–60 0.13 0.19 0.26 0.38 0.53 60–90 0.14 0.16 0.22 0.31 0.45 Table A.3. Median formal uncertainties for σ̟ . mV 3–4 4–5 5–6 6–7 7–8 0–10 0.20 0.25 0.36 0.51 0.71 10–20 0.19 0.26 0.34 0.51 0.70 A115, page 14 of 14 Ecliptic latitude (|β|, (deg)) 20–30 30–40 40–50 0.19 0.21 0.16 0.25 0.23 0.24 0.35 0.34 0.33 0.49 0.47 0.47 0.69 0.66 0.66 mV 3–4 4–5 5–6 6–7 7–8 0–10 0.24 0.28 0.41 0.58 0.80 10–20 0.22 0.29 0.37 0.55 0.77 Ecliptic latitude (|β|, (deg)) 20–30 30–40 40–50 0.18 0.19 0.13 0.26 0.21 0.17 0.35 0.30 0.24 0.50 0.43 0.33 0.70 0.60 0.47 50–60 0.12 0.17 0.22 0.32 0.44 60–90 0.15 0.16 0.23 0.32 0.45 50–60 0.13 0.17 0.24 0.34 0.48 60–90 0.15 0.16 0.23 0.33 0.46 50–60 15 51 183 514 1321 60–90 18 61 196 577 1447 Table A.5. Median formal uncertainties for σµδ . mV 3–4 4–5 5–6 6–7 7–8 0–10 0.16 0.19 0.27 0.37 0.53 10–20 0.16 0.20 0.27 0.38 0.53 Ecliptic latitude (|β|, (deg)) 20–30 30–40 40–50 0.16 0.18 0.13 0.19 0.17 0.18 0.26 0.26 0.26 0.38 0.36 0.36 0.52 0.52 0.52 Table A.6. Number of stars in each bin. mV 3–4 4–5 5–6 6–7 7–8 0–10 20 69 209 565 1274 10–20 20 65 193 550 1287 Ecliptic latitude (|β|, (deg)) 20–30 30–40 40–50 20 12 17 63 57 61 209 169 167 581 538 522 1416 1372 1355