KINEMATICAL AND CHEMICAL VERTICAL STRUCTURE OF THE GALACTIC THICK DISK. I. THICK DISK KINEMATICS^*,^†

The interstellar reddening E(B−V) was derived for each target from the Schlegel et al. (1998) maps, and the K_s magnitude and (J − K_s) color were de-reddened by means of the transformations of Cardelli et al. (1989). The distance to each star was then estimated through a color–absolute-magnitude relation calibrated on 47 Tucanae, a disk globular cluster (Zinn 1985) whose stellar population is very similar to the Galactic thick disk in both age and metallicity (Wyse & Gilmore 2005). The fit of the cluster red giant branch yields the relation (L. Miller 2008, private communication)

$$\begin{equation} K_{s}= -7.886\cdot (J-K_{s}) + 16.302, \end{equation} \tag{ 1 }$$

which were translated into absolute magnitude and de-reddened color assuming (m − M)_V = 13.50 ± 0.08 and E(B − V) = 0.024 ± 0.004 for the cluster distance modulus and reddening, respectively (Gratton et al. 2003).

Inspecting the 2MASS photometric data of 47 Tuc used to derive Equation (1), we found that the uncertainty on M_K is of the order of ∼0.2 mag. This is only a marginal contribution to the final error in distance, because the main source of uncertainty is the relatively wide range of age and metallicity covered by thick disk stars. Indeed, the 2MASS photometric errors have only negligible impact, because they do not exceed 0.03 mag for our fainter targets (K_s ∼ 10.6). We estimated the effect of the age and metallicity distribution on the derived absolute magnitudes by means of Yale-Yonsei isochrones (Yi et al. 2003). We assumed a scatter of 0.3 dex in metallicity that should include the bulk of thick disk stars (Carney et al. 1989) when excluding the scarcely populated low- and high-metallicity tails (Beers & Sommer-Larsen 1995; Bensby et al. 2007), as well as a scatter of 2 Gyr in age (Bensby et al. 2003; Feltzing et al. 2003; Reddy et al. 2006). They were considered uncorrelated, because the age–metallicity relation for the thick disk is very weak (Bensby et al. 2007). We finally estimated the error on distance to be ∼20%, quadratically summing all the relevant contributions.

Thin disk stars do not follow the age and metallicity distribution assumed to derive Equation (1), and their distances should be systematically biased. Indeed, younger, more metal-rich red giants are intrinsically fainter than our estimate. Comparing the absolute magnitudes calculated from Equation (1) with Yale-Yonsei isochrones between 2 and 8 Gyr and metallicity following the age–metallicity relation of Haywood (2001) and Rocha-Pinto et al. (2000), we found that the distance of thin disk stars would be overestimated by 10%–20%. This systematic error is small, but it has a relevant consequence on the contamination of the sample, which is artificially stretched to larger heights from the plane.

Our sample is contaminated by the thin disk, because the targets were selected through a color cut that efficiently eliminates low-metallicity stars, but which excludes only a small fraction of metal-rich objects, as can be deduced from Figure 1 of Girard et al. (2006). In Figure 2 we plot an estimate of the fraction of thin disk stars in the sample, assuming 300 and 900 pc as thin and thick disk scale height, respectively, and a thick disk local normalization of 12% (Jurić et al. 2008). We also calculated the contamination in the case of a distance overestimate of 15% for thin disk objects. The curves can be considered an upper limit, because the color cut, unaccounted for in the calculation, could have removed part of the contamination. In the same figure we indicate the results of a rough empirical estimate, obtained fitting a double Gaussian to the distribution of the vertical velocity component (W) in five bins of 50 stars. The dispersions of the two Gaussian components were kept fixed, and the only free parameter of the fit was the number of thin disk stars in the bin. The error bars show the results when varying σ_W in the range 14–16 and 35–40 km s⁻¹ for the thin and the thick disk, respectively. This test was performed only for |Z| ≈ 1–2 kpc, where both populations contribute by more than 20% of the objects, i.e., 10 stars in each bin. The empirical error bars agree well with the curve obtained when the distance bias is considered, although the observed contamination fall-off with |Z| seems steeper than the modeled one.

Zoom In Zoom Out Reset image size

In conclusion, our sample is contaminated by thin disk stars, whose distances are overestimated. This affects even the estimate of their kinematics, which is also biased. Therefore, we will not analyze the results for the thin disk, and we will limit our study to |Z| ⩾ 1.3 kpc, to avoid the strong contamination (⩾50%) of the nearest subsample.

Girard et al. (2006) estimated that the halo contamination in the sample should be about 8%; hence, we would expect 41 halo members among the 514 targets with |Z| ⩾ 1.3 kpc. The cut in W removed 10 probable halo contaminants and, according to the statistics of Chiba & Beers (2000), 18 of the excluded low-metallicity targets are expected to be halo stars. The cuts in U and V could also have removed some halo objects, but the effect of this cut is harder to quantify, because an unknown fraction of the outliers is probably objects with problematic proper motion or wrong distance. The residual halo contamination in our sample of 412 targets should therefore be lower than ∼3% (13 stars), and it can be safely neglected.

The radial velocity (RV) for all the targets was measured with a cross-correlation (CC) technique (Tonry & Davis 1979) as implemented in the IRAF fxcor task. The procedure was described in detail in Moni Bidin (2009), and we give only a brief summary here. The spectrum of each star was cross-correlated with three template stars observed in the same run, namely, HD 180540, HD 223559, and HD 213893 (this last replaced by HD 224834 for Coralie data), encompassing the temperature range of the targets. Their RV was fixed by a CC with a synthetic spectrum generated with the SPECTRUM⁶ code (Gray & Corbally 1994), because we found poor agreement between the available literature sources. The analysis of the 17 standard stars acquired during observations, and of the solar spectra collected each night, revealed that the RV zero point thus defined was biased by 0.3–0.7 km s⁻¹, depending on the instrument, and this offset was removed.

The three measurements were averaged, although they never differed by more than 0.2 km s⁻¹. The CC was restricted to the range 4600–6800 Å (5000–6800 Å for MIKE data), to avoid the poor-quality blue section and the telluric bands at longer wavelengths. All RVs were corrected to heliocentric velocities, and then the position of the telluric band at 6875 Å was used to correct the RVs for instrumental effects, mainly caused by a displacement on CCD between the lamp and science spectra and an off-center position of the star inside the slit (see, for example, the analysis of Moni Bidin et al. 2006). Corrections of up to 2 km s⁻¹ were applied, but with little scatter (∼0.5 km s⁻¹) within each observing night.

The final RV errors were estimated as the quadratic sum of the five relevant uncertainties introduced in the procedure: the CC and wavelength calibration error, the uncertainty on the zero point and its offset, and the error on the estimate of the correction for instrumental effects. The resulting errors were typically in the range 0.4–0.7 km s⁻¹. The final RVs of all the 824 stars will be made available at the CDS Web site.⁷

We found 211 previous RV measurements for the stars in our sample. Our results agree excellently with literature sources: the mean difference (in the sense ours−literature) is 0.4 ± 2.7 km s⁻¹, where the uncertainty is the rms of the differences. The comparison with the works that share with us more than five stars in common is given in Table 2. The mean RV difference is always of the order of 1 km s⁻¹, except for Majewski et al. (2004), whose RVs are higher than ours by a non-negligible amount (∼5 km s⁻¹). However, the number of stars in common is too small to conclude that this offset is significant.

Table 2. Mean RV Difference between this Work and the Quoted Reference for the Stars in Common

Reference	No. of Stars	$\overline{\Delta _\mathrm{RV}}$
		(km s−1)
Flynn & Freeman (1993)	145	0.8 ± 2.5
Zwitter et al. (2008)	9	−1.2 ± 1.5
Beers & Sommer-Larsen (1995)	9	0.4 ± 1.8
Majewski et al. (2004)	8	−4.7 ± 2.1
Duflot et al. (1995)	6	−1.1 ± 3.8
Jones (1972)	6	−0.7 ± 0.9
Evans (1970)	6	1.8 ± 4.2

Download table as:  ASCIITypeset image

The proper motion, RV, and distance of each target were combined to derive its (U, V, W) velocity components in the Galactic cylindrical reference frame, where U is positive toward the Galactic center, V is directed in the sense of Galactic rotation, and W points toward the north Galactic pole. The error on these velocities was derived propagating the uncertainty on the observed quantities. We assumed a proper-motion error of 3 mas yr⁻¹ for all the stars, as this value is more realistic than the uncertainties quoted in the SPM3 catalog (T. Girard 2009, private communication; see also Girard et al. 2006). The velocities were corrected subtracting the solar peculiar motion (U_☉, V_☉, W_☉) = (10.0, 5.1, 7.2) km s⁻¹ (Dehnen & Binney 1998). Schönrich et al. (2010) recently proposed the revised values (U_☉, V_☉, W_☉) = (11.0, 12.2, 7.3) km s⁻¹, but we preferred to adopt the older ones, for sake of continuity with previous works. In any case, the definition of the solar motion does not affect the velocity dispersions or the off-diagonal terms of the dispersion matrix (Equation (2)), while the effects on the mean values are discussed in Section 5.1.

The sample was then divided into several bins, in which we calculated the mean velocities, the dispersions (σ_U, σ_V, σ_W), and the non-diagonal terms of the dispersion matrix:

$$\begin{equation} \sigma ^2_{\alpha \beta }=\frac{1}{(N-1)}\Sigma _i (v_{\alpha,i}-\overline{v_{\alpha }})(v_{\beta,i}-\overline{v_{\beta }}), \end{equation} \tag{ 2 }$$

where the sum is extended to all the stars in the bin and α, β = U, V, W. The results are given in Table 3. The bins were defined by the requirement that their centers were equally spaced by 0.1 kpc from |Z| = 1.5 to 4.5 kpc, to uniformly sample the variation with |Z| of the kinematical quantities, while the width was imposed by the fixed number of stars per bin (see below). This implied a partial overlap of the bins, increasing with distance from the plane owing to the decreasing density of observed stars. However, Moni Bidin et al. (2010) have shown that the binning definition does not alter the results, and in fact our results are very similar to that work, despite the very different binning schemes. The bins with $\overline{\vert Z\vert }\ge 2.4$ kpc comprised 50 stars each, while at lower heights, where the number of observed stars is larger, the bin size was increased to 80 targets for $2.1\le \overline{\vert Z\vert }\le 2.4$ kpc and 100 targets for $\overline{\vert Z\vert }\le 2.1$ kpc. We thus ensured a good statistic of thick disk stars in bins contaminated by the thin disk.

Table 3. Measured Kinematical Quantities and Associated Uncertainties in Each Bin

|Z|	$\overline{U}$	$\overline{V}$	$\overline{W}$	σU	σV	σW	σ2UW	σ2UV	σ2VW
(kpc)	(km s−1)	(km s−1)	(km s−1)	(km s−1)	(km s−1)	(km s−1)	(km2 s−2)	(km2 s−2)	(km2 s−2)
1.5	−31 ± 4	−46 ± 4	17 ± 1	81 ± 5	57 ± 5	38 ± 1	570 ± 310	400 ± 470	10 ± 210
1.6	−20 ± 4	−55 ± 4	5 ± 1	77 ± 5	59 ± 5	37 ± 1	840 ± 290	−180 ± 460	−20 ± 220
1.7	−28 ± 4	−61 ± 5	3 ± 1	79 ± 5	64 ± 5	38 ± 1	930 ± 310	−130 ± 510	−40 ± 250
1.8	−24 ± 4	−69 ± 5	8 ± 1	79 ± 5	61 ± 6	40 ± 1	810 ± 330	−760 ± 500	−80 ± 240
1.9	−16 ± 4	−77 ± 5	3 ± 1	78 ± 5	60 ± 5	40 ± 1	850 ± 330	−1060 ± 500	−120 ± 240
2.0	−10 ± 5	−85 ± 4	4 ± 1	83 ± 6	55 ± 6	39 ± 1	940 ± 340	30 ± 460	−90 ± 210
2.1	−22 ± 5	−78 ± 5	3 ± 1	80 ± 6	58 ± 6	38 ± 1	1060 ± 320	470 ± 470	−130 ± 220
2.2	−16 ± 5	−76 ± 7	−13 ± 1	77 ± 6	63 ± 7	42 ± 1	690 ± 330	0 ± 550	230 ± 260
2.3	−1 ± 6	−81 ± 5	−9 ± 1	81 ± 7	58 ± 6	40 ± 1	750 ± 340	−500 ± 470	10 ± 230
2.4	−12 ± 5	−85 ± 6	−12 ± 1	80 ± 6	59 ± 6	40 ± 1	410 ± 330	−710 ± 480	180 ± 240
2.5	−18 ± 6	−85 ± 6	−2 ± 1	78 ± 7	63 ± 7	42 ± 1	630 ± 330	−1870 ± 600	−40 ± 270
2.6	−25 ± 6	−90 ± 6	1 ± 1	71 ± 7	66 ± 7	42 ± 1	530 ± 310	−2390 ± 660	−460 ± 290
2.7	−23 ± 6	−98 ± 6	7 ± 1	72 ± 7	62 ± 7	39 ± 1	870 ± 290	−1150 ± 520	−1360 ± 300
2.8	−29 ± 6	−95 ± 7	−2 ± 1	76 ± 7	62 ± 7	41 ± 1	880 ± 320	−1740 ± 580	−1310 ± 300
2.9	−17 ± 6	−115 ± 7	5 ± 1	83 ± 7	68 ± 7	40 ± 1	840 ± 340	−820 ± 590	−1230 ± 300
3.0	−7 ± 6	−126 ± 7	8 ± 1	90 ± 7	67 ± 8	42 ± 1	1440 ± 410	−50 ± 610	−460 ± 280
3.1	−6 ± 7	−129 ± 8	1 ± 1	101 ± 8	67 ± 8	43 ± 1	1610 ± 450	−1150 ± 690	−650 ± 290
3.2	2 ± 6	−131 ± 8	−1 ± 1	99 ± 8	63 ± 8	42 ± 1	1470 ± 440	−1820 ± 650	−440 ± 280
3.3	8 ± 7	−140 ± 8	0 ± 1	101 ± 8	66 ± 8	44 ± 2	2180 ± 490	−2610 ± 740	−310 ± 300
3.4	12 ± 7	−140 ± 8	1 ± 1	98 ± 8	63 ± 9	43 ± 2	2100 ± 480	−1880 ± 660	260 ± 280
3.5	18 ± 7	−135 ± 9	5 ± 1	95 ± 9	63 ± 9	44 ± 2	1950 ± 460	−3390 ± 740	260 ± 280
3.6	29 ± 7	−137 ± 8	9 ± 1	101 ± 9	64 ± 9	44 ± 2	2290 ± 510	−2170 ± 700	360 ± 290
3.7	7 ± 7	−132 ± 9	0 ± 2	91 ± 9	61 ± 10	44 ± 2	2180 ± 460	−2150 ± 610	610 ± 280
3.8	15 ± 7	−134 ± 10	1 ± 2	92 ± 9	68 ± 10	43 ± 2	1990 ± 450	−2120 ± 700	−30 ± 300
3.9	5 ± 8	−139 ± 9	−6 ± 2	94 ± 10	66 ± 10	43 ± 2	1540 ± 430	−520 ± 620	210 ± 280
4.0	7 ± 8	−143 ± 10	−1 ± 2	93 ± 9	66 ± 11	41 ± 2	1910 ± 420	−100 ± 610	10 ± 280
4.1	17 ± 8	−145 ± 10	1 ± 2	92 ± 10	68 ± 10	44 ± 2	1570 ± 440	720 ± 640	−30 ± 300
4.2	22 ± 8	−142 ± 11	2 ± 2	95 ± 10	72 ± 11	45 ± 2	1480 ± 460	880 ± 710	10 ± 320
4.3	12 ± 9	−148 ± 10	−2 ± 2	94 ± 10	69 ± 11	48 ± 2	2320 ± 510	1170 ± 680	180 ± 330
4.4	17 ± 9	−152 ± 11	0 ± 2	96 ± 10	72 ± 11	46 ± 2	2330 ± 500	1370 ± 710	−50 ± 330
4.5	19 ± 10	−158 ± 11	−2 ± 2	93 ± 11	76 ± 11	46 ± 2	2180 ± 480	3960 ± 930	240 ± 360

Download table as:  ASCIITypeset image

In each bin, the mean value and dispersion of each velocity component were measured by means of the analysis of the corresponding probability plot (Hamaker 1978; Lutz & Hanson 1992). This powerful tool is highly insensitive to outliers, and it can be employed even in the case of poorly populated bins. The data are first ordered with increasing value of the velocity. Each point i is then assigned a value σ, equal to the deviation from the mean expected for the ith point of a normal distribution, in units of the standard deviation. When the underlying distribution is Gaussian, the data follow a linear relation in the σ–velocity plane, whose intercept and slope are given by the mean value and the dispersion, respectively. The probability plots were fitted with a straight line, weighting each data point by its error; two examples are shown in Figure 3. The intrinsic velocity dispersions were then obtained quadratically subtracting the mean velocity error in the bin. The wings of the distribution showing any deviation from linearity were excluded from the fit, suspected to be affected by objects with wrong distance or problematic proper motion, or by residual halo members.

Zoom In Zoom Out Reset image size

The analysis of artificial data samples, generated through Monte Carlo simulations and analyzed as the real data, indicated that the velocity dispersion is underestimated by 1–2 km s⁻¹ when the thin disk contamination approaches 10%. The contamination was therefore neglected for $\overline{\vert Z\vert }\ge 2.5$ kpc, where it is expected to be lower than 5% (i.e., two stars per bin; see Figure 2). Following Bochanski et al. (2007), in the presence of a non-negligible thin disk contamination ($\overline{\vert Z\vert }\le 2.5$ kpc) we derived the mean velocity and the dispersion of thick disk stars fitting only the wings of the probability plot, and only the negative one for V. The cut excluded the ranges |U| ⩽ 50 km s⁻¹, |W| ⩽ 30 km s⁻¹, and V ⩾ −60 km s⁻¹, approximately 1.5 times the expected thin disk dispersion, thus removing about 90% of the contaminants. The two wings of the probability plot should have the same slope but a different intercept, and the mean velocity was obtained averaging the two intercepts obtained from the fit. Although the mean values thus derived agree well with the trends observed in the contamination-free bins at higher |Z|, we regard them as less reliable.

The formal errors of the least-squares procedure used to obtain the intercept and the slope in the probability plots mainly indicate the deviations from Gaussianity of the observed distribution. Hence, they are not a good estimate of the real uncertainties. The errors were therefore derived by means of Monte Carlo simulations. For each bin, we generated 1000 artificial samples of 50 stars, changing the proper motion, distance, and RV of each star assuming a Gaussian random noise with dispersion equal to the estimated errors (3 mas yr⁻¹, 20%, and 0.4–0.7 km s⁻¹, respectively). The measurements were repeated in each artificial bin as done for the real data, and the dispersion of these artificial measurements was assumed as an estimate of the observational errors. The final errors are given in Table 3.

The non-diagonal terms of the dispersion matrix, σ²_UV, σ²_UW, and σ²_VW, were calculated by means of Equation (2). The errors were estimated from Monte Carlo simulations, as done for the velocity dispersions. For each bin, 1000 artificial data sets of 50 values were simulated, drawn from a distribution characterized by the same dispersion matrix as the real data. In each simulation, a Gaussian random noise with a dispersion equal to the observational errors was added, and the non-diagonal terms were estimated by means of Equation (2). The errors were defined by the scatter of these simulated measurements. With the same simulations we also evaluated the systematic error introduced by the observational uncertainties (e.g., Siebert et al. 2008), comparing the results when the random noise was considered or not. We found that the expected systematic thus introduced is negligible, being about one order of magnitude lower than the uncertainties on the measured values. Therefore, we did not correct the observed non-diagonal terms for observational errors.

The trend of the mean velocities with distance from the plane is shown in Figure 4. $\overline{W}$ is consistent with no significant departure from zero in the whole range. Some deviations are observed for |Z| ⩽ 2.5 kpc, but in this range the results are affected by large errors and, as already commented, not very trustworthy. This result agrees with previous studies that found no net vertical motion up to various kpc from the Galactic plane (Bochanski et al. 2007; Bond et al. 2010; Casetti-Dinescu et al. 2011). On the contrary, $\overline{U}(Z)$ has a more complex behavior: between |Z| = 1.5 and 3 kpc we detect a non-null mean velocity directed toward the Galactic anticenter, with an average value $\overline{U}=-19\,{\pm}$ 3 km s⁻¹. Beyond 3 kpc from the plane $\overline{U}$ abruptly increases and changes sign, and the net velocity between 3.5 and 4.5 kpc is 15 ± 2 km s⁻¹ toward the Galactic center. Another possible interpretation of the results is that $\overline{U}$ linearly increases with |Z|, and an S-shaped feature is superimposed to this trend between 2.5 and 3.5 kpc, as observed for σ_U and other kinematical quantities (see Section 5.5). After the exclusion of this range, the fit returns a vertical increase of 15 ± 2 km s⁻¹ kpc⁻¹. Adopting the new values of Schönrich et al. (2010) for the solar peculiar motion does not cause a relevant change in the results, as $\overline{U}$ would be higher by about 1 km s⁻¹ only.

Zoom In Zoom Out Reset image size

Previous studies have shown evidence that stars more distant than 1 kpc from the plane show a non-null net radial motion of the order of ∼9 km s⁻¹ toward the Galactic anticenter (Smith et al. 2009; Rybka & Yatsenko 2010; Casetti-Dinescu et al. 2011). Our results between |Z| = 1.5 and 3 kpc agree with their conclusion, although our mean value is larger by about a factor of two. The sign flip observed at larger |Z|, on the contrary, has never been reported in the literature. Nevertheless, Smith et al. (2009) analyze halo stars, Rybka & Yatsenko (2010) do not reach these heights above the plane, and Casetti-Dinescu et al. (2011) have too few stars in this range. Very interestingly, Bond et al. (2010) detect a small positive mean radial component out to |Z| ≈ 3 kpc, and a negative value in their more distant bins, although the authors conclude that these results are consistent, within errors, with a null net motion. Casetti-Dinescu et al. (2011) proposed that an inward motion of the local standard of rest (LSR) of the order of 10 km s⁻¹ is responsible for the observed non-null mean value of $\overline{U}$, as suggested by the fact that this is observed among both disk and halo stars more distant than 1 kpc from the Sun. This hypothesis cannot account for any vertical trend of $\overline{U}$ other than a constant non-zero value at any |Z|. Thus, while it is not contradicted by our results, at least another effect must be invoked. For example, an outward motion of ∼25 km s⁻¹ of the stars between 3.5 and 4.5 kpc, coupled with the mentioned LSR motion, could explain our results. This clumpy kinematical distribution would not be surprising, because perturbations produced by the bar and the spiral arms are expected to form groups of stars with velocity offset as large as 50 km s⁻¹ (Quillen et al. 2011).

Girard et al. (2006) studied the proper motions of our same stellar sample, and their results for $\overline{U}(Z)$ are overplotted to ours in the upper panel of Figure 4, after correcting for the solar peculiar motion and changing the sign of U to account for the different direction of the reference axis. The agreement with our results is poor: while their nearly flat profile roughly coincides with our mean value (−7.7 km s⁻¹), they did not detect any change of sign or a steep positive gradient. The different approach to the same data must have caused the different results. For example, the features observed by us could have been masked in Girard et al. (2006) by their smoothed, statistical distance determination or by dwarf stars and halo contaminants, removed in our work.

In the middle panel of Figure 4 we show the measured trend of $\overline{V}(Z)$, compared with other results from the literature. The mean rotational velocity of thick disk stars decreases with distance from the plane, owing to the well-known vertical shear first detected by Majewski (1992). The fit of the data, after the exclusion of the less reliable bins at |Z| ⩽ 2.5 kpc, yields $\overline{V}(Z)=(-25\pm 12)-(30\pm 3)\cdot \vert Z\vert$ (km s⁻¹). As can be seen in Figure 4, this solution excellently matches the results of Girard et al. (2006), and even the data points of Chiba & Beers (2000), at z ⩽ 2 kpc, are well aligned with the derived linear relation. Both these investigations measure a vertical shear of −30 km s⁻¹ kpc⁻¹, as also recently found by Brown et al. (2008) and Bond et al. (2010). Had we fitted all the data points down to |Z| = 1.5 kpc, we would have found a steeper slope (−35.1 ± 1.8 km s⁻¹ kpc⁻¹), at the upper limit of the range spanned by literature values, which vary from −16  ± 4 (Allende Prieto et al. 2006) to −36  ± 1 km s⁻¹ (Carollo et al. 2010). The revised values for the solar peculiar motion proposed by Schönrich et al. (2010) cause an upward revision of the results by 7.1 km s⁻¹.

Casetti-Dinescu et al. (2011) and Spagna et al. (2010) measured a shallower slope (−25.2  ±  2.1 and −19  ±  2 km s⁻¹ kpc⁻¹, respectively) between 0.7 and 2.8 kpc. While our results are compatible with theirs in our nearest bins, we find a lower mean rotational velocity beyond |Z| = 2 kpc. Nevertheless, the different vertical rotational gradient found by these studies is not necessarily a disagreement, because they sample a different |Z|-range, and the underlying shear is not required to be strictly linear. In fact, Ivezic et al. (2008) have proposed the nonlinear relation $\overline{V}(Z)=-20.1-19.2\cdot \vert Z\vert ^{1.25}$ km s⁻¹ from the analysis of Sloan Digital Sky Survey (SDSS) data. The combined data points of Chiba & Beers (2000), Spagna et al. (2010), Casetti-Dinescu et al. (2011), and the present work closely follow this equation (see Figure 5), and a fit in the form $\overline{V}(Z)=\alpha +\beta \cdot \vert Z\vert ^{\gamma }$ km s⁻¹ returns a very similar solution, with α = −22.5 km s⁻¹, β = −22.2 km s⁻¹ kpc⁻¹, and γ = 1.23. Very noticeably, the results of four surveys finding a different linear relation are all well reproduced by a single nonlinear curve proposed by an independent work. The underlying vertical trend of the thick disk rotational velocity is therefore most probably close to but not exactly linear and better represented by a low-order power law.

Zoom In Zoom Out Reset image size

The vertical profile of the velocity dispersions is shown in Figure 6, compared with other works from the literature. We excluded from comparison the results of Bond et al. (2010), because they do not distinguish between different disk components, thus finding steeper gradients as a result of the mix of thin and thick disk stars. Caution must also be taken when comparing the results obtained in different ranges of |Z|, because the underlying trend of the dispersions is not necessarily linear and the gradient can assume different values. For example, the models of Girard et al. (2006) indicate that the vertical profile should be progressively flatter at increasing |Z|.

Zoom In Zoom Out Reset image size

The linear fit of the trends with |Z| yields

$$\begin{equation} \sigma _\mathrm{U}(Z)=(82.9\pm 3.2)+(6.3\pm 1.1)\cdot (\vert Z \vert -2.5) \,\mathrm{km\,s^{-1}}, \end{equation} \tag{ 3 }$$

$$\begin{equation} \sigma _\mathrm{V}(Z)=(62.2\pm 3.1)+(4.1\pm 1.0)\cdot (\vert Z \vert -2.5) \,\mathrm{km \,s^{-1}}, \end{equation} \tag{ 4 }$$

$$\begin{equation} \sigma _\mathrm{W}(Z)=(40.6\pm 0.8)+(2.7\pm 0.3)\cdot (\vert Z \vert -2.5) \,\mathrm{km \,s^{-1}}, \end{equation} \tag{ 5 }$$

where |Z| is in kpc. The quoted errors were obtained restricting the fit to a subset of eight non-overlapping bins, to avoid the underestimate induced by the correlation between the data points. The results are identical to those of Moni Bidin et al. (2010), despite the different binning scheme, except for σ_V, for which we derive a slightly smaller gradient, but the difference is not significant (0.4 km s⁻¹ kpc⁻¹). The gradients are small, and the three dispersions increase by only ∼7% between 2.5 and 3.5 kpc. This explains why the change of the thick disk kinematics with distance from the plane has not been clearly detected for about two decades after its discovery.

The results for σ_W are the most precise and reliable, because ∼90% of the information on W comes from RVs, whose errors are an order of magnitude smaller than those in proper motions. σ_W(Z) shows a small but clear and constant increase, with small scatter around the derived linear relation. The vertical gradient is smaller than the one found by Casetti-Dinescu et al. (2011) by about a factor of two, but the results are consistent at the 1σ level. Yoss et al. (1987) propose an even steeper profile (∼10 km s⁻¹ kpc⁻¹) in the range |Z| = 0–2 kpc, but the same authors suspect that this could be due to contamination by halo stars increasing with Galactic height.

The dispersion of the rotational velocity component, σ_V, also shows a clear increase with |Z|, although the data points are affected by larger errors and are more scattered. The vertical gradient is smaller than the results of both Casetti-Dinescu et al. (2011) and Girard et al. (2006), by a factor of three and two, respectively, and the difference is at the 2σ level in both cases.

The vertical profile of σ_U shows a very peculiar behavior, with large deviations from linearity between 2.5 and 3.5 kpc. This feature will be discussed in Section 5.5. Outside this range, the data points closely follow a linear relation, and the solution given in Equation (3) was obtained excluding this interval from the fit. The derived vertical gradient of σ_U is very similar to that found by Girard et al. (2006), but their solution is offset toward higher values by about 10 km s⁻¹, while the gradient measured by Casetti-Dinescu et al. (2011) is 4σ times higher than ours.

In conclusion, we confirm that the velocity dispersions of the Galactic thick disk increase with distance from the plane, as suggested by previous investigations (Yoss et al. 1987; Yoss & Griffin 1997; Soubiran et al. 2003; Girard et al. 2006; Ivezic et al. 2008; Casetti-Dinescu et al. 2011), but we derive vertical gradients that are in general smaller than other studies. The differences with Casetti-Dinescu et al. (2011) can be at least in part due to the aforementioned expected decrease of the gradient with distance from the plane, because the results for σ_W are consistent, and the solutions proposed for the other two components overlap in the Z-range in common (1.5–2.2 kpc). On the contrary, Girard et al. (2006) studied our same sample, and the different data analysis must be the cause of the disagreement. Since they measured both higher dispersions and steeper vertical gradients than us, their results could have been affected by the dwarf stars and halo contaminators removed by us, and/or the thin disk contamination, that they did not take into account.

The three dispersions increase with distance from the plane proportionally, and the anisotropy is approximatively constant at all |Z|. In Figure 7 the vertical profile of the ratios σ_U/σ_W and σ_U/σ_V is shown, where the errors come from propagation of the uncertainties on the dispersions.

Zoom In Zoom Out Reset image size

The mean value of σ_U/σ_W is 2.08 ± 0.06, where the error was estimated from the statistical error on the mean for eight uncorrelated bins, as before. The linear fit of the data points returns a negligible gradient (0.09 ± 0.07 kpc⁻¹), and the data are consistent with a flat profile. Literature values for σ_U/σ_W span a wide range from 1.1 (Veltz et al. 2008) to 1.9 (e.g., Vallenari et al. 2006), and our results are at the upper end of this distribution. As shown in the upper panel of Figure 7, the results of Casetti-Dinescu et al. (2011) agree with our measurements in the Z-range in common, but they deviate from our measurements if extrapolated to higher distance from the plane.

The results for σ_U/σ_V also show no significant gradient (0.06 ± 0.05 kpc⁻¹) and a mean value of 1.36 ± 0.05. This is lower than the value predicted by the epicyclic approximation (1.49; cf. Equation (3-76) of Binney & Merrifield 1998), indicating that the population under analysis cannot be assumed in nearly circular orbits. Our results are well within the range spanned by literature values, which vary from ≈1 (e.g., Chiba & Beers 2000; Carollo et al. 2010) to ≈1.6 (Soubiran et al. 2003). As in the case of σ_U/σ_W, the results of Casetti-Dinescu et al. (2011) overlap with our data points in the range in common, but their extrapolation to higher |Z| does not. The analytical expressions given by Girard et al. (2006) for σ_U and σ_V, on the contrary, return a value of their ratio much lower than our data points at any Galactic height. This is probably due to their much higher σ_V, as shown in the middle panel of Figure 6.

The non-diagonal term of the dispersion matrix σ²_ij gives a measurement of the correlation between the ith and jth velocity components. In fact, if two principal axes of the dispersion ellipsoid are aligned with the ith and jth axes, we have σ²_ij = 0, and the two velocity components are uncorrelated. The orientation of the dispersion ellipsoid of an old, dynamically relaxed population is related to the shape of the Galactic potential (Lynden-Bell 1962; Ollongren 1962; Hori & Liu 1963) and can indicate the presence of non-axisymmetric structures in the disk (Kuijken & Tremaine 1991). Moreover, Bienaymé (2009) and Siebert et al. (2008) have shown that the vertical tilt of the ellipsoid allows one to estimate the flattening of the dark halo. However, calculating the expectation of a specific halo model through integration of orbits is beyond the scope of the paper, and here we will only present the observational results, for use in future works.

Our results for the three cross-terms are shown in Figure 8. The profile of σ²_UV closely follows a decreasing linear relation up to |Z| = 4 kpc but abruptly jumps to positive values at larger heights. This is probably due to the sign flip of $\overline{U}$ analyzed in Section 5.5. In any case, the results indicate that the U and V velocity components are correlated, and the velocity ellipsoid is titled in the radial-longitudinal plane. The rotation angle, i.e., the vertex deviation ψ, was calculated in each bin by means of the relation

$$\begin{equation} \psi =-\frac{1}{2}{\rm arctg}\left(\frac{2\sigma ^2_\mathrm{UV}}{\sigma ^2_\mathrm{U}-\sigma ^2_\mathrm{V}}\right) \end{equation} \tag{ 6 }$$

(Amendt & Cuddeford 1991). The results are shown in the upper panel of Figure 9. We measure a non-negligible vertex deviation, increasing from nearly zero at |Z| = 1.5 kpc to ∼20° at 3.5 kpc. The linear fit in the range 1.5–4 kpc yields the relation ψ = −1.0 + 5.4 · |Z|. Our results agree with Casetti-Dinescu et al. (2011), who measured ψ = 82 ± 32 at Z = 1.1 kpc. Previous investigations showed evidence that the vertex deviation decreases from ≈20° for young populations to nearly zero for old, metal-poor disk stars (Bienaymé 1999; Dehnen & Binney 1998; Soubiran et al. 2003; Fuchs et al. 2009). Our results are not at variance with this conclusion if the vertex deviation increases with |Z|, as suggested by our observations, because previous studies were limited to small Galactic heights.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

σ²_VW shows no significant deviation from zero in the whole range of |Z|, except between 2.5 and 3.5 kpc, as discussed in Section 5.5. This indicates that the V and W velocity components are not correlated. On the contrary, σ²_UW is significantly different from zero and steadily increases with |Z|. Proposing a reliable linear expression for σ²_UW(Z) is not straightforward, because it shows irregular fluctuations. As expected by Binney & Spergel (1983) and Binney & Merrifield (1998), σ²_UW is always bound between zero and $\sigma ^2_\mathrm{UW,max}=(\overline{U^2}-\overline{W^2})\cdot (Z/R)$, the value assumed when the velocity ellipsoid is aligned with the spherical coordinate system, because we find that σ²_{UW, max} ⩾ 4500 km² s⁻² at any Galactic height. We also find that the expression proposed by Kuijken & Gilmore (1989), obtained under the assumption that the ellipsoid points toward the Galactic center and its axis ratio is constant in the spherical coordinate system, is a relatively good approximation of the measured value up to about 3 kpc from the plane, but the predicted slope is too shallow and the agreement with observations degrades with |Z| (see Figure 8).

The tilt angle in the U–W plane can be calculated, analogously to the vertex deviation, from the equation

$$\begin{equation} \alpha =-\frac{1}{2}{\rm arctg}\left(\frac{2\sigma ^2_\mathrm{UW}}{\sigma ^2_\mathrm{U}-\sigma ^2_\mathrm{W}}\right), \end{equation} \tag{ 7 }$$

and the results are shown in Figure 9. We measure a slight increase of α with |Z|, and the linear fit of the data, after the exclusion of two deviating points, yields the relation α(|Z|) = 96 + 24 · [(|Z|/kpc) − 2)]. The results agree with previous works at |Z| = 1–1.5 kpc (Siebert et al. 2008; Carollo et al. 2010; Casetti-Dinescu et al. 2011), but not with the measurements of Fuchs et al. (2009) at |Z| ⩽ 1 kpc, whose sample is most probably dominated by thin disk stars. The dashed curve in Figure 9 indicates the value of α when the dispersion ellipsoid points toward the Galactic center, assuming R_☉ = 8 kpc. The tilt angle is constantly lower than this in the range 1.5–4.5 kpc, and the ellipsoid thus is directed toward a point located behind the Galactic center, at a Galactocentric distance increasing with |Z| from R_o = 2.3 to 9 kpc. Orbit integration studies indeed predict this result (Kuijken & Gilmore 1989; Binney & Spergel 1983; Kent & de Zeew 1991; Shapiro et al. 2003), but the inclination is noticeably higher than the expectations, because theoretical calculations return R_o = 5–10 kpc at |Z| = 1.1 kpc. This could indicate that the Galactic potential used in these studies needs to be refined to match the observations. On the contrary, the increase of α(|Z|) modeled by Bond et al. (2010) is too steep, and their expectation α(3.5 kpc) = 26° is at variance with our results.

Bienaymé et al. (2009) showed that the vertical trend of the tilt angle is an excellent observational signature of the underlying gravity law, and it can be used as a test for the MOdified Newtonian Dynamics (MOND) theory (Milgrom 1983), because its expectation diverges from that of the Newtonian dynamics with distance from the plane. Unfortunately, the calculations of Bienaymé et al. (2009) are limited to |Z| ⩽ 2 kpc, and the range of overlap with our data is very narrow. The value of α expected by Bienaymé et al. (2009) at |Z| = 2 kpc is 12° for Newtonian dynamics and 10° for MOND. We find α(2 kpc) = 95 ± 36, and the mean of the five measurements in the range 1.8–2.2 kpc is 100 ± 05, in excellent agreement with MOND expectations. The linear fit presented in Section 5.3 also implies α(2 kpc) = 96. However, we derive a vertical gradient (24 kpc⁻¹) much shallower than the predictions of Bienaymé et al. (2009) for both MOND (5° kpc⁻¹) and Newtonian dynamics (6° kpc⁻¹). In conclusion, our observations agree better with the models that Bienaymé et al. (2009) derived for MOND than with those from Newtonian dynamics, but the extension of Bienaymé et al. (2009) calculations to higher |Z|, and to specific initial conditions for the thick disk, is required to perform a reliable test of the gravitational law.

The vertical profile of σ_U(Z) shows a puzzling behavior between 2.5 and 3.5 kpc from the Galactic plane (top panel, Figure 6), where the dispersion deviates from the linear trend. A similar feature is clearly observable in the trend of other quantities related to the RV component U, like $\overline{U}$ and σ²_UW (as already noted by Moni Bidin et al. 2010), but deviations from the linear trend in the range |Z| = 2.5–3 kpc could be present even in other profiles, such as $\overline{V}(Z)$, $\overline{W}(Z)$, σ_W(Z), and σ²_VW(Z). The origin of this behavior is unclear. It is possible that a group of comoving stars, forming a sub-structure in the Galactic thick disk, is affecting the measured kinematics between 2.5 and 3.5 kpc from the plane. In particular, the similarity of the profile of two totally unrelated quantities such as σ²_VW(Z) and σ_U(Z) is instructive, because it excludes the possibility that this behavior is only due to some bad measurements. In any case, the existence of a kinematical substructure among our stars cannot be claimed on the basis of these results only.

Our targets are distributed in a very narrow range of R, and the data do not provide any direct information about the radial behavior of the kinematical quantities. Nevertheless, some indirect indication can be derived. Manipulating the radial Jeans equation in cylindrical coordinates in steady state, with the radial component of the force expressed as F_R = −v²_c/R, we obtain

$$\begin{equation} \frac{\partial \sigma ^2_\mathrm{U}}{\partial R}\,{=}\,\frac{\overline{UW}}{h_{Z,\rho }}\,{-}\,\frac{\partial \overline{UW}}{\partial Z} -\frac{v^2_c}{R}\,{+}\,\frac{\sigma ^2_\mathrm{V}+\overline{V}^2}{R}\,{-}\,\sigma ^2_\mathrm{U}\left(\frac{1}{R}\,{-}\,\frac{1}{h_{R,\rho }}\right), \end{equation} \tag{ 8 }$$

where we also assumed that $\overline{\mathrm{U}}^2\ll \sigma ^2_\mathrm{U}$, $\partial \overline{\mathrm{U}}/\partial R$ = 0, and the density decays exponentially with both R and Z, with scale length h_{R, ρ} and scale height h_{Z, ρ}. Inserting the observed quantities in the right-hand side of Equation (8), we can thus estimate, as an exercise, the radial behavior of σ_U. The large uncertainties involved prevent a precise measurement, but the comparison with the expectations provides a consistency check for the kinematical data presented here, because an unphysical result could indicate a problem with them. The results are shown in Figure 10, where we assumed h_{R, ρ} = 3.6 kpc and h_{Z, ρ} = 0.9 kpc (Jurić et al. 2008), R_☉ = 8 kpc, and v_c = 220 km s⁻¹. The data points are scattered around the mean value −10 ± 1 km s⁻¹ kpc⁻¹, following the behavior of σ_U(Z) discussed in Section 5.5, with no clear vertical trend. In the nearest bins the results are partially consistent with, but more negative than, those found by Casetti-Dinescu et al. (2011) and Neese & Yoss (1988, −3.8 ± 0.6 km s⁻¹ kpc⁻¹), although Neese & Yoss refer to the thin disk. The observations show (van der Kruit & Searle 1981, 1982) that σ_W should exponentially decay in the radial direction, with the same scale length of the mass density. This is often considered valid even for σ_U, under the assumption of a constant anisotropy, and Cuddeford & Amendt (1992) found that this should be the best approximation at the solar position. The results also roughly agree with the expectation of this model, as shown in Figure 10. In conclusion, the presented kinematics implies a radial behavior of σ_U consistent with both previous observations and the expectations of theoretical predictions.

Zoom In Zoom Out Reset image size

The presence of kinematical gradients in the Galactic thick disk is a powerful diagnostic to discriminate between the various models of its formation (Majewski 1993), but unfortunately the model expectations concerning the vertical gradient have not been investigated so far. It would be interesting to compare, in the near future, the predictions of the various models with the observed trend with Z of the velocity dispersions. On the contrary, the radial gradient of the dispersions has been modeled in the context of thick disk formation through disk heating by the merging of minor satellites (Villalobos & Helmi 2008; Bekki & Tsujimoto 2011). The observations of Casetti-Dinescu et al. (2011) agree with the results of these simulations, provided that a low inclination orbit of the merging satellite is assumed. A similar conclusion is drawn comparing our measurements of σ_U/σ_W with the theoretical expectations of the merging scenario. In fact, the simulations showed that this ratio is strongly linked to the inclination angle of the orbit (Villalobos & Helmi 2008; Purcell et al. 2009), and our result (σ_U/σ_W = 2.08) favors a small inclination angle, i ≈ 0°–30° (see Figure 15 of Villalobos et al. 2010). A low inclination orbit of the infalling satellite is also required to reproduce the rather large vertical gradient of the rotational velocity, found in this work as in previous investigations (Villalobos & Helmi 2008; Bekki & Tsujimoto 2011). On the contrary, the radial migration model predicts a much shallower gradient (−17 km s⁻¹ kpc⁻¹; Loebman et al. 2011), incompatible with the observations.

In conclusion, the models of thick disk formation have so far provided only fragmentary predictions about its kinematical properties, and our observations cannot still be used to fully discriminate between them. However, we find that our results are consistent with the scenario where the thick disk formed through dynamical heating of a pre-existing Galactic disk, induced by the merging of a minor satellite. Moreover, all the kinematical evidence shows that, if this is the correct model, a low-latitude (⩽30°) merging event is strongly preferred.