© ESO, 2012

1. Introduction

Stellar radial velocity surveys such as the RAdial Velocity Experiment (RAVE, Steinmetz et al. 2006; Zwitter et al. 2008; Siebert et al. 2011b) can provide, for all the observed stars, three of the six components of the Galactic stellar phase space – two directional components and one velocity component. The lack of the other three components makes it difficult to fully characterize the dynamical behaviour of our Galaxy and its present phase-space description from such surveys.

Despite these limitations, RAVE has the unique merit of systematically measuring stellar radial velocities with a precision never previously realized for such a large sample of stars, providing us with a statistically significant amount of data that continues to grow, recently surpassing 500 000 spectra. RAVE is a southern-hemisphere survey obtaining spectra in the near-infrared CaII triplet region of magnitude selected stars (9 < I < 13). The usefulness of this vast dataset has already been demonstrated in the context of kinematics in order to deduce the characteristics of the stellar velocity distributions (e.g., Veltz et al. 2008; Casetti-Dinescu et al. 2011; Siebert et al. 2008, 2011a; Hahn et al. 2011; Coşkunoǧlu et al. 2011; Wilson et al. 2011; Karataş & Klement 2012; Bilir et al. 2012), to identify the stellar streams (e.g., Seabroke et al. 2008; Klement et al. 2008; Williams et al. 2011; Kiss et al. 2011), and to study the chemistry of MW components (e.g., Ruchti et al. 2010, 2011b; Fulbright et al. 2010; Boeche et al. 2011; Coşkunoǧlu et al. 2012). Here, we carry out an analysis on topics recently addressed by Coşkunoǧlu et al. (2011) and Casetti-Dinescu et al. (2011), leaving aside for the moment the interpretative tools based on the orbit integration as developed by Wilson et al. (2011). In other words, although we will make use of a Galactic potential model, the evolution with time t will not be the direct object of our analysis. The magnitude range of RAVE (9 < I < 13) (see e.g., Steinmetz et al. 2006; Zwitter et al. 2008; Siebert et al. 2011a) implies that most of the stars targeted by RAVE belong to the thin or the thick disk of the MW. Moreover, working with RAVE data to investigate the kinematics of our Galaxy in a statistical sense, basically means deriving the time-independent single-component disk-like distribution function (DF) f⁽x,v⁾, of the stars in the phase space mapped with coordinates {x,v}, or an overlapping set of DFs ${\mathit{f}}^{\mathrm{tot}}{}^{\mathrm{\right(}}{x}\mathit{,}{v}^{\mathrm{\left)}}\mathrm{=}\sum {\mathit{f}}_{\mathit{j}}{}^{\mathrm{\right(}}{x}\mathit{,}{v}^{\mathrm{\left)}}$ for each stellar population j that we are able to disentangle. In terms of direct observations of the phase space, and within reasonable errors, RAVE provides the direction (line of sight, l.o.s.) and the velocity along the l.o.s. for each star. The radial velocity is the observed projection of the true heliocentric observed velocity vector v_hel of the star along the line of sight, (with ⟨...,...⟩ denoting the standard inner product). The direction is just the heliocentric position vector divided by the unknown heliocentric distance of each star, ∥r_hel∥ = ∥x − x_⊙∥ ≡ r_hel (with components in a suitable reference system, see Sect. 3).

The analysis of the kinematic properties of the thin and thick disks is of paramount importance for the comprehension of the origin and formation of our Galaxy as well as of every disk galaxy. The object of our analysis in the present paper is mainly the thick disk. A separate paper is dealing with the thin disk analysis (Pasetto et al. 2012, hereafter Paper II).

In this paper we investigate the kinematic properties of the Galactic disk components by applying a methodology based on singular value decomposition (SVD). SVD permits one to find the solution, in a least-squares sense, of an inhomogeneous system of N linear equations in the case of singular matrices (where N is the number of stars we select). A number of studies have used this method in the context of proper motion survey analyses: see, e.g., a fully analytical exercise developed for Hipparcos proper motion and parallax data in Dehnen & Binney (1998) or a recent work by Fuchs et al. (2009) combining data release 7 of the Sloan Digital Sky Survey (Abazajian et al. 2009) with astrometric USNO-B data supplemented by photometric distances.

Finally, several authors devised methods to derive photometric distances for RAVE (e.g., Klement et al. 2008, 2011; Zwitter et al. 2010; Breddels et al. 2010; Burnett et al. 2011). Here we will adopt the internal release of the catalogue of Zwitter et al. (2010) which includes different distance determinations depending on different sets of isochrones for about 260 000 stars (see Zwitter et al. 2010, for an extended discussion). We will use the Yonsei-Yale distances based on isochrones of Demarque et al. (2004).

The structure of this paper is as follows: in Sect. 2 we will prepare the data for the analysis by cleaning the radial velocities from the effects of Galactic rotation. In Sect. 3 we briefly review the inversion techniques we are going to apply for the analysis of the data. In Sect. 4 we present the results and discuss them in Sect. 5.

2. Preparing the data

Our procedure requires the following steps. We need to

• 1.

  correct the RAVE survey data for large scale effects (thisSect. 2),

• 2.

  isolate the thick disk (see Sect. 3),

• 3.

  evaluate the thick disk kinematics (see Sect. 4).

As an additional result of this procedure, we will be able to estimate two of the three components of the solar motion relative to the LSR and the velocity lag of the thick disk component.

In this way, we can determine some of the moments of the underlying multi-component distribution function f^tot(x,v⁾. While strictly speaking only the infinite series of the moments is equivalent to the original distribution function, only the first few moments will be the subject of our exercise (for a review on the DF functional form and the number of its independent moments see e.g., Cubarsi 2010). Thus, our study is intended to be of an exploratory nature rather than exhaustive, for a set of problems that are far from being fully solved.

2.1. Radial velocity component from the Galaxy differential rotation

For a given position of a star in the Galactic disk, we want to estimate the influence of the differential rotation of the Galaxy in the radial velocity component. This is necessary to minimize the influence of the radial component due to the Galactic rotation on the observed radial velocity.

To compute this first step we note that, for a star at a given position x in the Milky Way’s disk, the following general vector relation holds: ${{v̅}}_{\mathrm{LSR}}\left({{x}}_{\mathrm{\odot }}\right)\mathrm{+}{{v}}_{\mathrm{\odot }}\mathrm{+}{{v}}_{\mathrm{hel}}\mathrm{=}{{v̅}}_{\mathrm{c}}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{+}{{v}}_{\mathrm{p}}\mathit{,}$(1)where is the mean rotational velocity at the given position x in the Galactic disk reference system centred on the barycentre of the Milky Way (MW) (considered at rest or in rectilinear unperturbed motion). is the solar local standard of rest (LSR), which differs from the LSR speed at any other location in the Galaxy. v_⊙ is the Solar peculiar velocity relative to the Solar LSR and v_p is the peculiar velocity of the star relative to its own mean Galactic rotational speed.

It is also convenient to introduce here the orthonormal standard system of reference in the velocity space ${\mathrm{(}}^{}\mathit{O}\mathrm{;}\mathit{U,V,}{\mathit{W}}^{\mathrm{)}}$. This reference system is centred on the velocity of the Solar LSR, O, with U aligned with the reference system in the configuration space pointing to the Galactic centre, V aligned with the rotation of the Galaxy and W pointing to the north Galactic pole (NGP). We will call the generic velocity components in this reference system v = {v_U,v_V,v_W}. From Eq. (1), if we express the component of the radial velocity due to Galactic rotation as a function of the star’s longitude and latitude ${\mathrm{(}}^{}\mathit{l,}{\mathit{b}}^{\mathrm{)}}$, we get $\hspace{0.17em}\begin{array}{}\begin{array}{}\end{array}{\mathit{v}}_{\mathit{r}}^{\mathit{G}}\mathrm{=}\left(\parallel {{v̅}}_{\mathrm{c}}\mathrm{(}{{x}}^{\mathrm{)}}\parallel \frac{{\mathit{R}}_{\mathrm{\odot }}}{\mathit{R}}\mathrm{-}{\mathit{V}}_{\mathrm{LSR}}\left({\mathit{R}}_{\mathrm{\odot }}\right)\right)\cos \mathit{b}\sin \mathit{l}\mathit{.}\end{array}$(2)For simplicity we assume cylindrical symmetry (O;R,φ,z) in the configuration space for the Galactic model, v_LSR(x_⊙) = {0,V_LSR(R_⊙),0}. v_⊙ = {v_U, ⊙ ,v_V, ⊙ ,v_W, ⊙ } are the components of the peculiar motion of the Sun relative to the solar LSR. R is the Galactocentric radial distance in cylindrical coordinates for the position of a given star. R_⊙ is the position of the Sun on the plane of symmetry of the MW at an azimuthal position φ_⊙ = 0. The height of the Sun relative to the plane of symmetry of the Galaxy is hereafter neglected, z_⊙ ≅ 0. The orthonormal rotation matrix R^TR = 1, i.e. $\mathbf{R}\mathrm{\in }\mathit{SO}{}^{\mathrm{\right(}}\mathrm{3}^{\mathrm{\left)}}$ with det =  + 1 is adopted here to transform a velocity vector (given in terms of radial velocity v_r and motion along the Galactic coordinates ${}^{\mathrm{(}}\mathit{v}_{\mathit{l}}\mathit{,}{{\mathit{v}}_{\mathit{b}}}^{\mathrm{)}}$) into a vector in the velocity space reference system vector ${\mathrm{(}}^{}\mathit{O}\mathrm{;}\mathit{U,V,}{\mathit{W}}^{\mathrm{)}}$. This orthonormal rotation matrix is defined as $\mathbf{R}\mathrm{=}\left(\begin{array}{c}\\ \cos \mathit{b}\cos \mathit{l}& \cos \mathit{b}\sin \mathit{l}& \sin \mathit{b}\\ \mathrm{-}\sin \mathit{l}& \cos \mathit{l}& \mathrm{0}\\ \mathrm{-}\sin \mathit{b}\cos \mathit{l}& \mathrm{-}\sin \mathit{b}\sin \mathit{l}& \cos \mathit{b}\end{array}\right)\mathit{.}$(3)

2.2. The influence of distance

We show now how the influence of the Galactic rotation on the radial velocity is strongly dependent on the viewing direction, i.e., the direction ${\mathrm{(}}^{}\mathit{l,}{\mathit{b}}^{\mathrm{)}}$ for a specific star of the survey, and on the mean streaming velocity $\parallel {{v̅}}_{\mathrm{c}}{\mathrm{(}}^{}{x}{\mathrm{)}}^{}\parallel$ that will be our free tuning parameter. On the other hand, we can choose suitable directions (l,b) for which the influence of the Galactic rotation on the radial velocity is weakly dependent on the distance r_hel of the star along (l,b). This result holds for the range of distances from the Sun in which RAVE dwarf stars are mainly located, say r_hel ≤ 1.0   kpc as shown in Appendix C.

In Eq. (2)the Galactic radial distance of a template star can be expressed as a function of the heliocentric distance, r_hel. In particular the first term in the sum of Eq. (2)reads $\parallel {{v̅}}_{\mathrm{c}}\mathrm{(}{{x}}^{\mathrm{)}}\parallel \frac{{\mathit{R}}_{\mathrm{\odot }}}{\mathit{R}}\cos \mathit{b}\sin \mathit{l}\mathrm{=}\frac{\parallel {{v̅}}_{\mathrm{c}}\mathrm{(}{{x}}^{\mathrm{)}}\parallel \cos \mathit{b}\sin \mathit{l}}{\sqrt{\mathrm{1}\mathrm{+}\frac{{\mathit{r}}_{\mathrm{hel}}}{{\mathit{R}}_{\mathrm{\odot }}}\cos \mathit{b}\mathrm{(}\frac{{\mathit{r}}_{\mathrm{hel}}}{{\mathit{R}}_{\mathrm{\odot }}}\cos \mathit{b}\mathrm{-}\mathrm{2}\cos \mathit{l}\mathrm{)}}}\mathrm{·}$(4)We will treat the terms $\parallel {{v̅}}_{\mathrm{c}}{\mathrm{(}}^{}{x}{\mathrm{)}}^{}\parallel$ independently later. For simplicity, we set $\mathit{\epsilon }\mathrm{=}\frac{{\mathit{r}}_{\mathrm{hel}}}{{\mathit{R}}_{\mathrm{\odot }}}$. For our approximation we only consider the kinematics in the solar neighbourhood, defined, e.g., as the zone where $\mathit{\epsilon }\mathrm{=}\frac{{\mathit{r}}_{\mathrm{hel}}}{{\mathit{R}}_{\mathrm{\odot }}}\cong \frac{\mathrm{1}}{\mathrm{8}}$. Then it follows from Eqs. (2)and (4)that in this neighbourhood, by expanding the previous equation in a Maclaurin series in ε, we obtain $\begin{array}{ccc}\frac{\cos \mathit{b}\sin \mathit{l}}{\sqrt{\mathrm{1}\mathrm{+}\mathit{\epsilon }\cos \mathit{b}\mathrm{(}\mathit{\epsilon }\cos \mathit{b}\mathrm{-}\mathrm{2}\cos \mathit{l}\mathrm{)}}}& \mathrm{\simeq }& \cos \mathit{b}\sin \mathit{l}\\ & & \mathrm{+}\mathit{\epsilon }{\cos }^{\mathrm{2}}\mathit{b}\cos \mathit{l}\sin \mathit{l}\mathrm{+}\mathit{o}\mathrm{(}{{\mathit{\epsilon }}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{·}\end{array}$Thus this expansion suggests that ${\mathit{v}}_{\mathit{r}}^{\mathrm{G}}$ is primarily dominated by the direction of observation ${\mathrm{(}}^{}\mathit{l,}{\mathit{b}}^{\mathrm{)}}$ and the value (which differs from star to star and from the thin disk to the thick disk), but we expect no strong dependence on the heliocentric distance of the star within our neighbourhood. This is supporting the idea that the error affecting the photometrically determined distances will not be dramatically relevant for the thick disk parameter determinations (see also Casertano et al. 1990). At present, we must consider the term as a free parameter. We postpone its treatment and its dependence on the heliocentric distance ${{v̅}}_{\mathrm{c}}{}^{\mathrm{(}}{x}{}^{\mathrm{)}}\mathrm{=}{{v̅}}_{\mathrm{c}}{}^{\mathrm{\right(}}{{r}}_{\mathrm{hel}}^{\mathrm{\left)}}$ until Sect. 4, because its behaviour will be treated in the context of the inversion technique. Different arguments hold for the thin disk analysis (see Paper II) where the photometric distance errors are the main source of uncertainties for the trend of the velocity ellipsoid in the meridional plane.

Fig. 1 Isocontour plot on an Aitoff projection of the differences in radial velocity due to Galactic rotation $\mathit{\delta }{\mathit{v}}_{\mathit{r}}^{\mathrm{G}}$ from Eq. (5)for test values of calculated using Eq. (24)and the values in Table 1. The gradients of the velocity dispersion tensor are treated as in Vallenari et al. (2006). The assumed values for the motion of the Sun relative to the LSR are from this paper. The yellow points show the RAVE data selection from Sect. 2.1. The black contour lines represent the limit of 14 km   s-1 adopted in our paper.

We can test the validity of the previous approximation without using the expansion into a series, but instead by numerically solving for the dependence ${\mathit{v}}_{\mathit{r}}^{\mathrm{G}}\mathrm{=}{\mathit{v}}_{\mathit{r}}^{\mathrm{G}}{}^{\mathrm{\right(}}{\mathit{r}}_{\mathrm{hel}}^{\mathrm{\left)}}$ at different values of ${\mathrm{(}}^{}\mathit{l,}{\mathit{b}}^{\mathrm{)}}$. In particular, we are interested in controlling the errors introduced by large scale effects such as differential rotation and photometric distances for selected distances. The dispersion of the thick disk velocity ellipsoid is generally assumed to be constant and thus should not introduce any adverse large scale effects on the small volume considered. Nevertheless, we will check this assumption “a posteriori” in Sect. 4.4.2. If we limit ourselves to consider the RAVE stars within a sphere of radius 1.0 kpc centred at the Sun’s location, from Eq. (4)we can calculate the value of the radial velocity due to the Galactic rotation ${\mathit{v}}_{\mathit{r}}^{\mathrm{G}}$ along the line of sight r_hel ∈ [0,1] kpc for each direction in the sky specified by (l,b). To evaluate this equation, we test our initial guess values for Eq. (4)with the help of a Galactic population model by Vallenari et al. (2006) where the MW Galactic potential and velocity dispersion profiles are consistently determined (see also our Eq. (24)and Sect. 4.4.1).

This value of the Galactic component of the radial velocity ranges from a maximum, $\max {}^{\mathrm{\right(}}{\mathit{v}}_{\mathit{r}}^{\mathrm{G}}^{\mathrm{\left)}}$, and a minimum value, $\min {}^{\mathrm{\right(}}{\mathit{v}}_{\mathit{r}}^{\mathrm{G}}^{\mathrm{\left)}}$, whose difference we will call $\mathit{\delta }{\mathit{v}}_{\mathit{r}}^{\mathrm{G}}\mathrm{\equiv }{\left|\begin{array}{}\max \end{array}{\mathit{v}}_{\mathit{r}}^{\mathrm{G}}\mathrm{-}\begin{array}{}\min \end{array}{\mathit{v}}_{\mathit{r}}^{\mathrm{G}}\right|}_{{\mathit{r}}_{\mathrm{hel}}\hspace{0.17em}\mathrm{\in }\hspace{0.17em}\left[\mathrm{0}\mathit{,}\mathrm{1}\right]\mathit{.}}$(5)This difference $\mathit{\delta }{\mathit{v}}_{\mathit{r}}^{\mathrm{G}}$ can be evaluated at different values of ${\mathrm{(}}^{}\mathit{l,}{\mathit{b}}^{\mathrm{)}}$ and for an optimized value of . In this way we obtain the background in Fig. 1 that shows the Aitoff sky projection contour plot obtained from Eq. (5). As we can see from this contour plot, the highest $\mathit{\delta }{\mathit{v}}_{\mathit{r}}^{\mathrm{G}}$ (yellowish zones) occur at low latitude and the smallest in the green zones. For b = 0, neglecting the motion of the Sun relative to the Solar LSR, Eq. (5)is just the terminal velocity relation (e.g., Binney & Merrifield 1998). The highest induced error is $\mathrm{\left|}\mathit{\delta }{\mathit{v}}_{\mathit{r}}^{\mathrm{G}}\mathrm{\right|}\cong \mathrm{16} {\mathrm{kms}}^{-1}$ at a distance of 1.0   kpc. This may seem quite large if compared with the mean radial velocity error of RAVE data but this error occurs at low latitude regions that are not of interest for RAVE because these are not covered by the survey (see Fig. 1). Also the dependence of on the range of interest of r_hel is again more relevant at low latitude.

2.3. RAVE data selection criteria

It becomes evident from the previous considerations that when we introduce an indetermination in the distance d_i of the ith-star due to the error in the photometrically determined distance Δd_i, we are also introducing an error in the correction for the radial velocity contribution due to the Galactic rotation. If we proceed with a Monte Carlo approach by randomly generating a distance d_i within Δd_i for the ith-star, the correction ${\mathit{v}}_{\mathit{r,i}}^{\mathrm{G}}$ for this star will assume only values between the two extremes $\underset{{\mathit{r}}_{\mathrm{hel}}\mathrm{\in }\left[\mathrm{0}\mathit{,}\mathrm{1}\right]}{\max }{}^{\mathrm{\right(}}{\mathit{v}}_{\mathit{r}}^{\mathrm{G}}^{\mathrm{\left)}}$ and $\underset{{\mathit{r}}_{\mathrm{hel}}\mathrm{\in }\left[\mathrm{0}\mathit{,}\mathrm{1}\right]}{\min }{}^{\mathrm{\right(}}{\mathit{v}}_{\mathit{r}}^{\mathrm{G}}^{\mathrm{\left)}}$ in Eq. (5). If this correction is smaller than the intrinsic error on the radial velocity as derived by RAVE, $\mathrm{\Delta }{\mathit{v}}_{\mathit{r}}^{\mathrm{RAVE}}{}^{\mathrm{(}}\mathit{l}_{\mathit{i}}\mathit{,}{{\mathit{b}}_{\mathit{i}}}^{\mathrm{)}}$, i.e., $\mathrm{\Delta }{\mathit{v}}_{\mathit{r}}^{\mathrm{RAVE}}{}^{\mathrm{(}}\mathit{l}_{\mathit{i}}\mathit{,}{{\mathit{b}}_{\mathit{i}}}^{\mathrm{)}}\mathit{>}\mathit{\delta }{\mathit{v}}_{\mathit{r}}^{\mathrm{G}}{}^{\mathrm{(}}\mathit{l}_{\mathit{i}}\mathit{,}{{\mathit{b}}_{\mathit{i}}}^{\mathrm{)}}$, we can either apply this correction or we can safely neglect it. This is because its contribution to ${\mathit{v}}_{\mathit{r,i}}^{\mathrm{G}}$ is too small to bias our consideration, i.e. its effect cannot be accounted for within the error $\mathrm{\Delta }{\mathit{v}}_{\mathit{r}}^{\mathrm{RAVE}}{}^{\mathrm{(}}\mathit{l}_{\mathit{i}}\mathit{,}{{\mathit{b}}_{\mathit{i}}}^{\mathrm{)}}$ of the observed stellar radial velocity. Equivalently, by shifting the ith-star by its distance uncertainty, it does not affect the radial velocity contribution ${\mathit{v}}_{\mathit{r}}^{\mathrm{G}}$ due to the Galactic rotation because it is negligible compared to the intrinsic error $\mathrm{\Delta }{\mathit{v}}_{\mathit{r}}^{\mathrm{RAVE}}{}^{\mathrm{(}}\mathit{l}_{\mathit{i}}\mathit{,}{{\mathit{b}}_{\mathit{i}}}^{\mathrm{)}}$. We point out that the observational errors, Δv_r, for which we use capital delta (Δ), are of a different nature than the errors $\mathit{\delta }{\mathit{v}}_{\mathit{r}}^{\mathrm{G}}$ described by Eq. (5), which have different systematic trends and an origin that can be accounted for with a complete Galactic dynamics model (see Vallenari et al. 2006, their appendix), and our Eq. (24)).

In order to adopt an even more restrictive condition, we will use all those stars from the RAVE catalogue that have an intrinsic error $\mathrm{\Delta }{\mathit{v}}_{\mathit{r}}^{\mathrm{RAVE}}{}^{\mathrm{(}}\mathit{l}_{\mathit{i}}\mathit{,}{{\mathit{b}}_{\mathit{i}}}^{\mathrm{)}}$ greater than twice the error introduced by neglecting the stars’ distance dependence on the velocity due to the Galactic rotation. In other words, for every star i, we require: $\begin{array}{ccc}\mathrm{\Delta }{\mathit{v}}_{\mathit{r}}^{\mathrm{RAVE}}\mathrm{(}{\mathit{l}}_{\mathit{i}}\mathit{,}{{\mathit{b}}_{\mathit{i}}}^{\mathrm{)}}\mathit{>}\mathrm{2.0}\mathit{\delta }{\mathit{v}}_{\mathit{r}}^{\mathrm{G}}\mathrm{(}{\mathit{l}}_{\mathit{i}}\mathit{,}{{\mathit{b}}_{\mathit{i}}}^{\mathrm{)}}\mathit{.}& & \end{array}$(6)A test of the validity of all these arguments to recover the correct results has been performed on a mock catalogue as explained in Sect. 4.3.

From the sample of stars that survive the selection defined by Eqs. (5)and (6)shown in Fig. 1, we proceed further with a few extra selection criteria. We consider stars with signal-to-noise ratio S/N > 20. We require the near-infrared colours of the RAVE stars taken from 2MASS observations to be in the range of J − K ∈ [0.2,1.1] in order to clean the data set from extremely young stars. Furthermore, we require ${\mathit{v}}_{\mathit{r}}^{\mathrm{RAVE}}\mathit{<}\mathrm{300}{\mathrm{kms}}^{-1}$ to avoid stars dynamically not representative of the thick disk within 2 kpc from the Sun’s position, $\mathrm{\Delta }{\mathit{v}}_{\mathit{r}}^{\mathrm{RAVE}}\mathit{<}\mathrm{14}{\mathrm{kms}}^{-1}$ to reduce the propagation error while retaining a considerable amount of stars (see the velocity contour highlighted in black in Fig. 1), |b| > 10° because we are not interested in thin disk stars, and log ₁₀g > 3.5 to avoid (with g surface gravity) the influence of giant stars which may sample a much more distant part of the Galaxy. Together with Eq. (6)the last criterion represents the most severe cut on the total number of stars that we can use. Where possible a further cut [Fe/H] >  − 1 dex was applied adopting the metallicity determination described by Zwitter et al. (2008).

Starting with an initial number of roughly 260 000 stars in the data release of Zwitter et al. (2010), the final, remaining number of stars for which we can perform the analysis is N_tot ≃ 38 805 which is orders of magnitude larger than the samples currently available in the literature (see Fig. 3 for a plot of the error distribution).

3. Inversion techniques for radial velocities

In the following Sect. 3.1 we will present the technique to determine the principal moments of a DF for a mixture of thin and thick disk stars from a radial velocity survey. With a few extra assumptions we will be able to disentangle the two components and focus our attention on the thick disk as will be shown in Sect. 3.2.

3.1. The mixture of the thin and thick disk distribution of moments

Once we have selected the sample of RAVE stars that satisfy the conditions laid out in Sect. 2.3 (see yellow dots in Fig. 1), we can deduce the first moments of the underlying mixture distribution function with simple algebra based on a very popular technique for dealing with sets of linear equations that can be written in the matrix form Ax = b, for the unknown vector x. In the following we recall the basics applied to our specific case of a radial velocity survey. For more details we refer to Appendix A.

Based on Eq. (1)we can express the components of the radial velocity vector, say v _∥ , by ${{v}}_{\mathrm{\parallel }}\mathrm{=}\left(\begin{array}{c}\\ {\mathit{v}}_{\mathit{r}}\cos \mathit{b}\cos \mathit{l}\\ {\mathit{v}}_{\mathit{r}}\cos \mathit{b}\sin \mathit{l}\\ {\mathit{v}}_{\mathit{r}}\sin \mathit{b}\end{array}\right)\mathit{,}$(7)where with p being the idempotent matrix of the projection operator along the line of sight. Proceeding component by component we obtain: ${{v}}_{\mathrm{\parallel }}\mathrm{=}\mathbf{p}{{v}}_{\mathrm{hel}}\mathrm{=}{{r̂}}^{\mathrm{\otimes }\mathrm{2}}{{v}}_{\mathrm{hel}}\mathrm{=}\left(\begin{array}{c}\\ \mathit{x̂}\mathrm{2}& \mathit{x̂}\mathit{ŷ}& \mathit{x̂}\mathit{ẑ}\\ \mathit{x̂}\mathit{ŷ}& \mathit{ŷ}\mathrm{2}& \mathit{ŷ}\mathit{ẑ}\\ \mathit{x̂}\mathit{ẑ}& \mathit{ŷ}\mathit{ẑ}& \mathit{ẑ}\mathrm{2}\end{array}\right)\mathit{.}\left(\begin{array}{c}\\ {\mathit{v}}_{\mathit{U}}\\ {\mathit{v}}_{\mathit{V}}\\ {\mathit{v}}_{\mathit{W}}\end{array}\right)\mathit{,}$(8)where ${}^{\mathrm{\right(}}\mathit{...}^{\mathrm{\left)}}\mathrm{·}{}^{\mathrm{\right(}}\mathit{...}^{\mathrm{\left)}}$ is the ordinary matrix product, a ^⊗ n is the standard tensor power of the generic vector a (and n ∈ N any non-negative integer) and is again the unitary vector of a star whose components are in the configuration space collinear with the velocity space ${\mathrm{(}}^{}\mathit{O}\mathrm{;}\mathit{U,V,}{\mathit{W}}^{\mathrm{)}}$ and centred on O. Clearly the matrix p is singular (its ) thus not permitting us to determine its inverse.

Nevertheless, we can proceed by taking into consideration the whole subsample of the selected RAVE data in order to statistically invert the overdetermined system Ax = b that we obtain by defining A as the block diagonal matrix of all the projection operators p once the off-diagonal blocks are small (see Appendix A for further details). We call this (3Nx3) block matrix of all the projection operators simply P. b will be defined from the vectors of the observed radial velocities for the selected sample, which we call for simplicity again v _∥ . We obtain in this way the system ${\mathbf{P}\overline{{v}}}_{\mathrm{hel}}\mathrm{=}{{v}}_{\mathrm{\parallel }}$ that we multiply as usual by the transposed matrix P^T in order to obtain the square matrix product P^TP that can be inverted in order to give ${{v̅}}_{\mathrm{hel}}\mathrm{=}{{}^{\mathrm{\right(}}{\mathbf{P}}^{\mathit{T}}\mathbf{P}^{\mathrm{\left)}}}^{-1}{\mathbf{P}}^{\mathit{T}}{{v}}_{\mathrm{\parallel }}\mathrm{\equiv }{\mathbf{P}}^{\mathrm{+}}{{v}}_{\mathrm{\parallel }}$, where P ⁺  is the pseudo inverse of P (Penrose 1956, 1955). Probably the most widely used orthogonal decomposition suitable for the solution of our system is the Singular Value Decomposition due to its numerical stability (e.g., Press et al. 1986), which we will also adopt here (e.g., by writing ${\mathbf{P}}^{\mathrm{+}}\mathrm{\equiv }\frac{{\mathbf{WU}}^{\mathit{T}}}{\mathrm{diag}\mathbf{Q}}$ with W and U as orthogonal matrices and Q as a diagonal matrix (see e.g., Press et al. (1986) for an extended discussion and implementation techniques). Finally, in order to check the consistency of this approximation we will keep track of the sky distribution by weighting the obtained matrix P with the isotropic case laid out in Appendix A (see Sect. 4.3 and Eq. (A.3)).

There is an implicit hypothesis that is underlying our approach. When we are considering the system of equations , we need to implicitly assume that either the group of stars we are sampling is so extremely local that it is not (or only weakly) influenced by the spatial trend of the velocity dispersion tensor, or that the population we are sampling has an isothermal dispersion tensor within the distances sampled. In the latter case the sample does not need to be local. For the thick disk stars we can exploit the latter, while the first has been already exploited in the work of Dehnen & Binney (1998). A bi-dimensional schematic representation of this assumption can be seen in Fig. 2: the light blue arrows represent the unknown velocity vectors of which we observe the radial velocities (green arrows) along a given direction (grey arrows). The velocity dispersion described by the light blue arrows is found (see below) not to depend on the position inside or outside the solar radius R_⊙ throughout the volume covered by our sample. For comparison see Fig. 2 of Paper II where the trend of the thin disk velocity dispersion tensor along the meridional plane is considered.

Fig. 2 Schematic representation of the hypothesis of an anisotropic and isothermal distribution of an unbiased sample of thick disk velocities superimposed on a thin disk sample. The observer is located at the Sun’s position R⊙ at the centre. The Galactic centre (GC) is at the bottom of the figure. The directions on the celestial sphere to each star si are indicated with dashed arrows for thick disk stars and solid arrows for the thin disk stars. The green thick arrows represent the radial velocities of the thick disk, whose corresponding unknown true velocity vector is in light blue. The thin disk radial velocities are represented by the thin green arrows, and the corresponding unknown olive-green arrows represent the true velocity vectors. Note that the olive-green arrows are generally longer in the GC direction and shorter in the anti-Galactic centre direction as expected from the velocity dispersion trend.

This complex procedure provides our first moment of the composite distribution function, hereafter simply . With this approach the solar motion relative to a selected sample of stars is simply the mean motion of the stars relative to the Sun with a change in sign for the U and W directions, say ${\mathrm{\left(}{\mathit{v}}_{\mathit{U}}\mathit{,}{\mathit{v}}_{\mathit{W}}\mathrm{\right)}}_{\mathrm{\odot }}\mathrm{=}\mathrm{-}\mathrm{\left(}\overline{{\mathit{v}}_{\mathit{U}}}\mathit{,}\overline{{\mathit{v}}_{\mathit{W}}}\mathrm{\right)}$. In the V direction the lag of the component motion has to be taken into account with an extra v_lag whose contribution overlaps with v_V, ⊙  to form the mean value $\overline{{\mathit{v}}_{\mathit{V}}}$.

In the same way we can get the second, third, and fourth central moment, by defining and averaging, e.g., for the moment of order 2: $\overline{{{{v}}^{\mathbf{\prime }}}_{\mathrm{\parallel }}^{\mathrm{\otimes }\mathrm{2}}}\mathrm{=}\overline{{\left(\mathbf{P}{{v}}^{\mathbf{\prime }}\right)}^{\mathrm{\otimes }\mathrm{2}}}\mathrm{=}\overline{{\mathbf{P}}^{\mathrm{\otimes }\mathrm{2}}{{{v}}^{\mathbf{\prime }}}^{\mathrm{\otimes }\mathrm{2}}}\mathit{.}$(9)Hereafter we exploit the symmetries of the matrices: ${{\sigma }}_{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\overline{\mathrm{(}\mathbf{P}\mathrm{\otimes }{\mathbf{P}}^{\mathrm{)}}\mathrm{+}{\mathrm{(}\mathbf{P}\mathrm{\otimes }{\mathbf{P}}^{\mathrm{)}}}^{\mathit{T}}}\mathrm{·}{{\sigma }}^{\mathrm{2}}\mathrm{\equiv }\mathbf{M̅}\left[\mathrm{2}\right]\mathrm{·}{{\sigma }}^{\mathrm{2}}$(10)where ${{\sigma }}_{\mathrm{\parallel }}^{\mathrm{2}}$ is the dispersion tensor of the radial velocity and σ² is the ordinary pressure tensor. In the same way, we define the projected third moment along the l.o.s. as ${\mathbf{S}}_{\mathrm{\parallel }}\mathrm{\equiv }\overline{{{{v}}^{\mathbf{\prime }}}_{\mathrm{\parallel }}^{\mathrm{\otimes }\mathrm{3}}}\mathrm{=}\overline{{\mathbf{P}}^{\mathrm{\otimes }\mathrm{3}}{{{v}}^{\mathbf{\prime }}}^{\mathrm{\otimes }\mathrm{3}}}\mathrm{\equiv }\mathbf{M̅}\left[\mathrm{3}\right]\mathrm{·}\mathbf{S}\mathit{,}$(11)and the projected fourth moment as ${\mathbf{T}}_{\mathrm{\parallel }}\mathrm{=}\overline{{{{v}}^{\mathbf{\prime }}}_{\mathrm{\parallel }}^{\mathrm{\otimes }\mathrm{4}}}\mathrm{=}\overline{{\mathbf{P}}^{\mathrm{\otimes }\mathrm{4}}{{{v}}^{\mathbf{\prime }}}^{\mathrm{\otimes }\mathrm{4}}}\mathrm{\equiv }\mathbf{M̅}\left[\mathrm{4}\right]\mathrm{·}\mathbf{T}\mathit{,}$(12)where the suitable symmetries have been exploited. The interested reader can find the element by element characterization of the matrices in Appendix A.

The errors in these quantities are given by standard statistical tools coming from linear regression, taking into account that each star i contributes to v _∥  with its own error Δv _∥ i (referred to as $\mathrm{\Delta }{\mathit{v}}_{\mathit{r}}^{\mathrm{RAVE}}{}^{\mathrm{(}}\mathit{l}_{\mathit{i}}\mathit{,}{{\mathit{b}}_{\mathit{i}}}^{\mathrm{)}}$ in Eq. (6)). Hence, for instance, the error on the mean will be $\mathrm{\Delta }{v̅}{}^{\mathrm{2}}\mathrm{=}{{}^{\mathrm{\right(}}\overline{{\mathbf{P}}^{\mathrm{+}}}^{\mathrm{\left)}}}^{\mathrm{2}}\overline{\mathrm{\Delta }{{v}}^{\mathrm{2}}}\mathit{,}$ on the dispersion tensor it will be $\mathrm{\Delta }{{\sigma }}^{\mathrm{2}}\mathrm{=}\mathrm{2}{\left(\overline{{{}^{\mathrm{(}}\mathbf{P}\mathrm{\otimes }{\mathbf{P}}^{\mathrm{)}}}^{\mathrm{+}}}\right)}^{\mathrm{2}}\overline{{\left({{{v}}^{\mathbf{\prime }}}_{\mathrm{\parallel }}\mathrm{\otimes }\mathrm{\Delta }{{{v}}^{\mathbf{\prime }}}_{\mathrm{\parallel }}\right)}^{\mathrm{2}}}$, and so on (exactly as in Fuchs et al. 2009).

3.2. The scattering processes and the thick disk component disentanglement

So far, we have a procedure that, through the Eqs. (9) − (12), provides us with the first four true moments (i.e., ) of the distribution function of the thin and thick disk mixture from the projected moments along the l.o.s. (i.e., ) that we obtain directly from the observations.

To proceed further with the thick disk analysis, we want to disentangle the two distribution functions (DFs) and proceed by determining the thick disk velocity dispersion tensor alone.

We assume that the two single-particle DFs for the thin and thick disk, ${\mathit{f}}_{\mathrm{1}}{}^{\mathrm{\right(}}{x}\mathit{,}{v}^{\mathrm{\left)}}\mathrm{\equiv }{\mathit{f}}_{\mathrm{thin}}$ and ${\mathit{f}}_{\mathrm{2}}{}^{\mathrm{\right(}}{x}\mathit{,}{v}^{\mathrm{\left)}}\mathrm{\equiv }{\mathit{f}}_{\mathrm{thick}}$, can be added linearly to yield the overall distribution function f^tot(x,v) = f₁ + f₂. For the mixture distribution function f^tot in the previous Section we computed the first four moments. The first and the second moments correspond to the mean of the sample and to the dispersion tensor, respectively, ${{\sigma }}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{\mathit{N}}{\mathrm{\int }}_{}^{}{{{v}}^{\mathbf{\prime }}}^{\mathrm{\otimes }\mathrm{2}}{\mathit{f}}^{\mathrm{tot}}\left({x,v}\right)\mathrm{d}{v}\mathit{,}$(13)where the integral extends over the whole velocity space and N is the number of stars. The third and fourth moments are: $\mathbf{S}\mathrm{=}\frac{\mathrm{1}}{\mathit{N}}{\mathrm{\int }}_{}^{}{{{v}}^{\mathbf{\prime }}}^{\mathrm{\otimes }\mathrm{3}}{\mathit{f}}^{\mathrm{tot}}\left({x,v}\right)\mathrm{d}{v}\mathit{,}$(14)and$\mathbf{T}\mathrm{=}\frac{\mathrm{1}}{\mathit{N}}{\mathrm{\int }}_{}^{}{{{v}}^{\mathbf{\prime }}}^{\mathrm{\otimes }\mathrm{4}}{\mathit{f}}^{\mathrm{tot}}\left({x,v}\right)\mathrm{d}{v}\mathit{,}$(15)with the same definition interval for the integrals. We see that $\mathit{N}\mathrm{=}\mathit{N}{}^{\mathrm{\right(}}{x}\mathit{,t}^{\mathrm{\left)}}$, as well as ${{\sigma }}^{\mathrm{2}}\mathrm{=}{{\sigma }}^{\mathrm{2}}{}^{\mathrm{\right(}}{x}\mathit{,t}^{\mathrm{\left)}}$, $\mathbf{S}\mathrm{=}\mathbf{S}{}^{\mathrm{\right(}}{x}\mathit{,t}^{\mathrm{\left)}}$ and $\mathbf{T}\mathrm{=}\mathbf{T}{}^{\mathrm{\right(}}{x}\mathit{,t}^{\mathrm{\left)}}$ because, in general, ${\mathit{f}}^{\mathrm{tot}}\mathrm{=}{\mathit{f}}^{\mathrm{tot}}{}^{\mathrm{\right(}}{x}\mathit{,}{v}\mathit{,t}^{\mathrm{\left)}}$. The form of this distribution function is generally unknown but its sub-components, f_thin or f_thick, are often assumed to be quadratic in the peculiar velocity components, i.e., of a generalized Schwarzschild type.

In the time evolution of the total DF f^tot, we can argue that a very different role is played by f_thin or f_thick. A few processes of scattering are generally assumed to influence the time evolution of f_thin: mostly the scattering due to transient spiral arms and the encounters with giant molecular clouds (e.g., Binney & Tremaine 1987). Both processes lead to an increase of the velocity dispersion with increasing distance from the plane and provide insight into the origin of the age-velocity dispersion relation. Moreover, the Schwarzschild characterization of f_thin, on large scales in the velocity space, is not as good a description on smaller scales (e.g., Seabroke & Gilmore 2007). Sub-structures are mainly found for young stars (but not only for these). Such stars in a moving group are born at the same place and time, and then disperse into a stream that may intersect the solar neighbourhood. In this scenario, the stars that are moving in the same group should share the same age, metallicity, and azimuthal velocity (e.g., Dehnen 1998; Chereul et al. 1998). A different explanation for the sub-structures in the velocity DF is given by De Simone et al. (2004). These authors suggest that sub-structures arise naturally from the same spiral gravitational fluctuations that excite the growth of the velocity dispersion. In this picture, sub-structures are caused by homogeneous star formation in an irregular potential, as opposed to inhomogeneous star formation in a regular potential. Once the thin disk component is disentangled from the thick disk component, its analysis will proceed with the full exploitation of the proper motions (Paper II).

The situation is less clear for the thick disk where the role of the scattering processes is still undetermined. Its formation may be related to the influence of massive satellites that can either heat pre-existing disks (preserving a vertical metallicity gradient if present, (Katz et al. 2011) or get accreted (e.g., Quinn & Goodman 1986; Quinn et al. 1993). Other formation scenarios include gas-rich mergers (e.g., Brook et al. 2004), or internal radial migration processes (e.g., Schönrich & Binney 2009) or formation induced by satellites (e.g., Di Matteo et al. 2011). We will start working with the hypothesis that the thick disk DF can be approximated as ${\mathit{f}}_{\mathrm{thick}}\mathrm{\propto }{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\left({{{v}}^{\mathbf{\prime }}}^{\mathit{T}}{{C}}_{\mathrm{thick}}{{v}}^{\mathbf{\prime }}\right)}$ where the covariance matrix ${{C}}_{\mathrm{thick}}^{-1}$ is σ²_thick and the quadratic equation v^′TC_thickv^′ = 1 defines a velocity ellipsoid that characterizes the distribution function of the thick disk fully and that we want to determine. We will check “a posteriori” (see Sect. 4.4.2) if this approximation is consistent with our results for the group of stars selected as in Sect. 2.3. The thick disk is more diffuse in the solar neighbourhood than the thin disk and is characterized by a larger scale height (e.g., Binney & Merrifield 1998). We assume that the dependence of the velocity ellipsoid in the sample analysed with RAVE can be adequately considered to be independent from the Galactic plane distance, i.e., an isothermal but anisotropic DF, as commonly assumed in most of the thick disk velocity ellipsoid determinations available in the literature (e.g., Layden et al. 1996; Alcobé & Cubarsi 2005; Chiba & Beers 2000; Holmberg et al. 2007). We will return to the validation of this hypothesis in Sect. 4.4.

Under these assumptions, the decoupling between the thin and thick disk component is straightforward and analytically developed in, e.g., Cubarsi (1992, see also Cubarsi 2007). We are going to apply this inversion technique in its version based on the cumulants (Cubarsi & Alcobé (2004)) because this is numerically more stable and hence suitable to minimize the numerical error propagation.

We outline the methodology applied here because it differs slightly from the original work of Cubarsi & Alcobé (2004). The interested reader is referred to Appendix B, where the technique contained in Cubarsi & Alcobé (2004) is explicitly derived and adapted for our specific case. The derivation of the cumulants from the Eqs. (13) − (15)is straightforward: $\hspace{0.17em}\begin{array}{}\begin{array}{}\end{array}{\kappa }\left[{\sigma }\right]\mathrm{=}\frac{\mathit{N}}{\mathit{N}\mathrm{-}\mathrm{1}}{\sigma }\\ \begin{array}{}\end{array}{\kappa }\left[\mathbf{S}\right]\mathrm{=}\frac{{\mathit{N}}^{\mathrm{2}}}{\mathrm{(}{\mathit{N}\mathrm{-}\mathrm{1}}^{\mathrm{)}}\mathrm{(}{\mathit{N}\mathrm{-}\mathrm{2}}^{\mathrm{)}}}\mathbf{S}\\ \begin{array}{}\end{array}{\kappa }\left[\mathbf{T}\right]\mathrm{=}\frac{{\mathit{N}}^{\mathrm{2}}\mathrm{(}{\mathit{N}\mathrm{+}\mathrm{1}}^{\mathrm{)}}}{\mathrm{(}{\mathit{N}\mathrm{-}\mathrm{1}}^{\mathrm{)}}\mathrm{(}{\mathit{N}\mathrm{-}\mathrm{2}}^{\mathrm{)}}\mathrm{(}{\mathit{N}\mathrm{-}\mathrm{3}}^{\mathrm{)}}}\left(\mathbf{T}\mathrm{-}\frac{\mathit{N}\mathrm{-}\mathrm{1}}{\mathit{N}\mathrm{+}\mathrm{1}}{{\sigma }}^{\mathrm{\otimes }\mathrm{2}}\right)\mathit{,}\end{array}$(16)hereafter simply indicated as κ_ij, κ_ijk, κ_ijkl, respectively. If we let N_thin and N_thick be the unknown numbers of stars for the thin and thick disk components, we can define the parameter q as follows: $\mathit{q}\mathrm{\equiv }\sqrt{\frac{{\mathit{n}}_{\mathrm{thick}}}{{\mathit{n}}_{\mathrm{thin}}}}\mathrm{-}\sqrt{\frac{{\mathit{n}}_{\mathrm{thin}}}{{\mathit{n}}_{\mathrm{thick}}}}\mathit{,}$(17)where n_thin = N_thin/N and n_thick = N_thick/N are the unknown fractions of stars in the thin and thick disk. This parameter varies from 0 when 50% of the stars belong to each component to ∞ when 100% of the stars belong to the thick disk. Our goal is to maximize this parameter q in the subsample of stars selected from the RAVE catalogue as in Sect. 2.3.

Cubarsi & Alcobé (2004) show the great advantage of exploiting the symmetries of the normal distribution functions by working with a transformed peculiar vector. We follow this convention by defining a normalized vector as ${d̂}\mathrm{\equiv }\sqrt{{\mathit{n}}_{\mathrm{thin}}{\mathit{n}}_{\mathrm{thick}}}\frac{{w}}{{\mathit{w}}_{\mathrm{2}}}\mathit{,}$(18)where w = {w₁,w₂,w₃} is the velocity vector of the difference between the thin and thick disk centroids; . The reason for the normalization to the highest velocity difference component of w is to reduce the error propagation in the SVD solution using the best-determined differential velocity, i.e., w₂. We define the transformed peculiar velocity V that keeps the differential velocity in the azimuthal direction invariant: $\begin{array}{ccc}\hspace{0.17em}\begin{array}{}\begin{array}{}\end{array}{V}\mathrm{=}\mathbf{H}\mathrm{·}{v}\\ \begin{array}{}\end{array}\mathbf{H}\mathrm{=}\left(\begin{array}{c}\\ & & \\ & & \\ & & \end{array}\right)\mathrm{·}\end{array}& & \end{array}$This is an isomorphic transformation of the peculiar velocity vector v, with $\det {}^{\mathrm{\right(}}\mathbf{H}^{\mathrm{\left)}}\mathrm{=}\mathrm{1}$ due to the previously assumed normalization of Eq. (18). The new cumulant equation for the transformed DF function of V can be derived as a function of the old cumulants of Eq. (16)(see Appendix B). Once the new V-mixture cumulants have been computed, say K ^[3]  and K ^[4]  for the third and fourth cumulants respectively, we determine the parameter q of Eq. (17)by making again use of the SVD technique to determine the solution of the linear system of over-determined equations as already done before for the system of equations in Sect. 3.1. We present the description of the system in Appendix B to which we refer the interested reader.

Finally we extrapolate the values of n_thick, σ_thick, and w₂, i.e. the characterizing thick disk kinematics parameters we were searching for, with a simple numerical iterative procedure that maximizes the thick disk component alone (e.g., Simulated Annealing Methods, Press et al. 1986; Cubarsi et al. 2010). Within the preselected set of data satisfying the selection criteria of Sect. 2.3, we consider a Monte Carlo generation of the stellar parameters within their observational errors and we isolate the thick disk component alone with the methodology described above. The kinematic characterization of the thick disk is achieved with the same purely geometrical approach developed for the mixture but now on the thick disk selected sample alone. Once we apply this Monte Carlo approach to generate several catalogues, and then apply the procedure to each of them, we summarize our results after a simple statistical analysis.

4. Results

Here we present the results for the determination of the motion of the Sun relative to the LSR and the velocity dispersion tensor for the thick disk.

4.1. Solar motion relative to the LSR and thick disk rotation velocity

Computing the first moment of the velocity part of the distribution function, namely the mean, we are left with a free parameter . We assume that the dependence of on the model is rather weak for the two components, v_U, ⊙  and v_W, ⊙ , which we hence derive as $\begin{array}{c}\\ {\mathit{v}}_{\mathit{U,}\mathrm{\odot }}\mathrm{=}\mathrm{(}{\mathrm{9.87}\mathrm{\pm }\mathrm{0.37}}^{\mathrm{)}} {\mathrm{kms}}^{-1}\\ {\mathit{v}}_{\mathit{W,}\mathrm{\odot }}\mathrm{=}\mathrm{(}{\mathrm{8.01}\mathrm{\pm }\mathrm{0.29}}^{\mathrm{)}} {\mathrm{kms}}^{-1}\mathit{.}\end{array}$(19)The determination of the last value of : v_V, ⊙  is more complicated. This is mostly for the following two reasons:

• 1.

  The overlap of the rotational delay of the populationv_lag and the contribution of the peculiar motion of the Sun in the same direction v_V, ⊙ , both acting together on the mean streaming velocity $\parallel {\overline{{v}}}_{\mathrm{c}}\parallel$, prevents us from obtaining the relative velocity between the LSR and the standard of rest centred on the motion of the Sun from our first distribution moment, . Thus v_V, ⊙  remains undetermined.

• 2.

  A full determination of could in principle be achieved with this methodology by including all stars down to a low latitude, say, b < 15 deg. However, the survey data are missing this low-latitude sky coverage. Moreover, there is no mapping in the Galactic rotation direction (see Fig. 1).

Nevertheless, while we cannot carry out a clean determination of the first moment of the velocity section of the DF, by assuming the literature value from Francis & Anderson (2009), we gain a value for the velocity lag of our thick disk sample of stars, ${\mathit{v}}_{\mathrm{thick,lag}}\cong \mathrm{49}\mathrm{\pm }\mathrm{6}{\mathrm{kms}}^{-1}\mathit{.}$(20)As explained at the end of the previous section, the procedure works iteratively in order to select a single population of thick disk stars from which to deduce the kinematic parameters. This is done because we cannot determine the solar motion relative to the Solar LSR from a non-local mixture because the thin disk velocity dispersion trend on the meridional plane could bias the results. Thus we need to disentangle one single component from which to derive (v_U, ⊙ ,v_W, ⊙ ).

4.2. Thick disk velocity ellipsoid

The values determined for the velocity dispersion tensor, σ_RR = 59.2 ± 4.4  kms^-1, σ_φφ = 47.3 ± 7.5  kms^-1 and σ_zz = 35.9 ± 4.1  kms^-1 (see Table 1), are in good agreement with what was already presented in the literature for the thick disk component (e.g., Alcobé & Cubarsi 2005; Chiba & Beers 2000; Holmberg et al. 2007). The presence of the thin disk is expected to bias these values only slightly in the sense that the random selection of the stars used in order to maximize the parameter q in Eq. (17)is proven to be able to separate the two components as shown in Cubarsi & Alcobé (2004) and as confirmed by the simulated mock catalogues (see Sect. 4.3).

As can be seen, the vertical tilt of the velocity ellipsoid is deduced from the mixed term component σ_Rz = 9.3 ± 2.3  kms^-1, leading to a tilt angle of about $\mathit{\phi }\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\arctan \left(\frac{\mathrm{2}{\mathit{\sigma }}_{\mathit{Rz}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{RR}}^{\mathrm{2}}\mathrm{-}{\mathit{\sigma }}_{\mathit{zz}}^{\mathrm{2}}}\right)\cong \mathrm{3.07}\mathrm{\pm }\mathrm{1.1}\mathrm{deg}$. This is similar to the result of Siebert et al. (2008), who applied different selection criteria to the same data, but deviates from the recent determination by Fuchs et al. (2009). Regarding the comparison with these studies it is worth noting that we used Eq. (6), while in Fuchs et al. (2009) the tilt depends on the mixture of the two components (thin and thick disk) and it was considered without any technique capable of disentangling the relative statistical influence of the two, nor accounting for the radial dependence on R in σ_ij = σ_ij(R,z). In the meridional plane the thin disk’s radial component profile could partially influence the determination of the velocity ellipsoid tilt, if the binning is done only in the vertical direction (as proposed in their paper) but not in the radial direction.

4.3. Testing the results with the Padova Galaxy model

As indicated at the end of Sect. 3.2, in the spirit of a Monte Carlo approach to the error analysis, we proceed by generating a set of synthetic catalogues, numbered from j = 1,...,N_cat (with N_cat being a high number) by assigning to the ith star a random distance d_i within its photometric distance error Δd_i and a random radial velocity v_r,i within its radial velocity error Δv_r,i (no errors are assumed in the star’s coordinates (l,b)). Then, the procedure to disentangle the stellar thin and thick population is applied for the jth realization of the catalogue and for all the j = 1,...,N_cat. If N_cat is sufficiently high, the stability of the error is achieved with standard statistics tools on the N_cat realizations of the catalogue.

Nonetheless, to be sure our method is working correctly and to get more insight into its limits and potential, we apply the methodology just described to a fully analytical catalogue. Then, we apply the Padova Galaxy model (e.g., Ng et al. 1995, 2002; Vallenari et al. 2006) and Pasetto (2005) and references therein) to generate a full synthetic Galaxy where everything is under control to test the ability to recover the properties of the data generated by the new method developed here. In this procedure we proceed by producing mock RAVE catalogues with prescribed relations and parameters that we ultimately want to recover with the technique here developed. Initially three synthetic stellar populations, representative of the thin and thick disks and halo are generated. The stellar structure models are taken from Bertelli et al. (2009, 2008) from which we generate a stellar population by assuming an initial mass function (IMF) (e.g., Kroupa et al. 1993) and the star formation rate (SFR) is assumed as a free parameter. Once the synthetic Hertzsprung-Russel diagrams with stars of the stellar populations are generated and projected in the corresponding RAVE/2MASS colours and magnitude passbands, we can distribute the synthetic stars representing the mixture of three stellar populations in the phase space of the MW once a model of interstellar extinction is assumed (e.g., from Drimmel et al. 2003). We are particularly interested in the kinematic description the stellar populations of the disks. For them we assume double exponential spatial density profiles $\mathit{\rho }\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}\mathrm{(}{\mathit{t}}^{\mathrm{)}}\exp \left(\mathrm{-}\frac{\mathit{R}}{{\mathit{h}}_{\mathit{R}}\mathrm{(}{\mathit{t}}^{\mathrm{)}}}\mathrm{-}\frac{\left|\mathit{z}\right|}{{\mathit{h}}_{\mathit{z}}\mathrm{(}{\mathit{t}}^{\mathrm{)}}}\right)\mathit{,}$(21)where t is age of the stellar populations that we discretize for simplicity in four temporal ranges for the thin disk: i = 1 for t ∈ [0,3[Gyr with ρ_0,1 = 2.2 × 10⁸   M_⊙   kpc^-3, h_R,1 = 2.90 kpc, h_z,1 = 0.20 kpc, i = 2 for t ∈ [3,5[Gyr with ρ_0,2 = 1.0 × 10⁸   M_⊙   kpc^-3, h_R,2 = 2.90 kpc, h_z,2 = 0.25 kpc, i = 3 for t ∈ [5,7[Gyr with ρ_0,3 = 2.1 × 10⁸   M_⊙   kpc^-3, h_R,3 = 3.10 kpc, h_z,3 = 0.28 kpc, i = 4 for t ∈ [7,10[Gyr with ρ_0,4 = 8.0 × 10⁸   M_⊙   kpc^-3, h_R,4 = 3.10 kpc, h_z,4 = 0.35 kpc. The thick disk scale parameters adopted are ρ_0,thick = 1.1 × 10⁶   M_⊙   kpc^-3, h_R,thick = 2.50 kpc, h_z,thick = 1.31 kpc and interstellar medium parameters are ρ_0,ISM = 1.48 × 10⁸   M_⊙   kpc^-3, h_R,ISM = 4.54 kpc, h_z,ISM = 0.20 kpc. To complete the density profile from which we deduce the potential that is used in the Boltzman equation to implement the kinematic description (as, e.g., in Eq. (24)) we add a simple Hernquist bulge (Hernquist 1993) and logarithmic potential for the halo. We then tune a Poisson-solver to match the observational constraints on the rotation curve, i.e. the Oort functions, the local density ratio between the stellar populations, the terminal velocity for the inner Galaxy with respect to the solar radius, the total mass inside 100 kpc etc. (see Appendix A in Vallenari et al. 2006 for an extended description and Pasetto 2005).

It is beyond the goal of this paper to derive the structural parameters of the MW from the RAVE data. Thus, the previous values are frozen and have to be considered as guess values assumed to test our new methodology. Other literature scale parameters, shorter in the scale length of the thick disk, can be easily tested (e.g., Carollo et al. 2010; Bensby et al. 2011) as suggested in relation to the alpha-enhancements (e.g., Cheng et al. 2012; Brook et al. 2012). Once we have fixed these parameters, we fix the trend of the velocity dispersion tensor for the thin disk in the meridional plane ${\mathit{\sigma }}_{\mathrm{thin}}^{\mathrm{2}}$ (see Paper II) and arbitrary values for the velocity dispersion tensor of the thick disk ${\mathit{\sigma }}_{\mathrm{thick}}^{\mathrm{2}}$ with which to generate our mock catalogue and that we want to recover with our novel method. The Padova Galaxy modelling technique produces a mock catalogue with synthetic values for $\left\{\mathit{m̂}\mathit{,}\mathit{ĉ}\mathit{,}\log \mathit{ĝ}\mathit{,}{\mathit{r̂}}_{\mathrm{hel}}\mathit{,l,b,}{\mathit{v̂}}_{\mathit{r}}\mathit{,}{\mu ̂}\right\}$(22)where (and ĉ) are magnitude (and colour) in a desired passband, and log ₁₀ĝ is the logarithm of the stellar surface gravity. Considering that the cut adopted in the surface gravity in Sect. 2.3 lies around the turn-off of a thick disk population with t ∈ [10,12[Gyr, we have an almost bijective relation between stellar magnitude and distance once the Drimmel et al. (2003) 3D extinction map is considered, thus is a catalogue of photometrically determined distances for the dwarf stars of the RAVE catalogue¹. In Eq. (22), when , then l,b are exactly the observed RAVE stellar directions, (and ) are the radial velocities (and proper motions) for our synthetically generated catalogue.

We extensively tested the new technique against the ability to recover

• the correct thick disk velocity dispersion tensor with which wegenerate the RAVE mock catalogues;

• the influence of large-scale effects on the thick disk velocity ellipsoid (see also Paper II);

• the minimal sky coverage for which we can correctly recover the underlying thick disk kinematics;

• the possibility to correctly disentangle the v_V, ⊙ ;

• the role of the log ₁₀g to disentangle dwarf and giant stars in our sample.

But one of the most interesting results of our extended analysis on the stability and selection cut criteria is probably the exercise presented in Appendix C, where we show that the methodology developed here mathematically can retrieve the correct thick disk velocity ellipsoid even without the use of photometric distances.

It is possible to make further cuts, for example in the colour magnitude diagram (CMD), but as shown in Veltz et al. (2008) or Siebert et al. (2008) this drastically reduces the available number of stars. Large numbers of stars, however, are the means of our present method to reduce the error bars. Interesting different approaches have also been discussed in the literature, e.g. by Robin et al. (2003); Schuster et al. (1993) or Casertano et al. (1990); Ratnatunga et al. (1989). The applicability of their methods to RAVE data is also interesting but beyond the scope of the present paper.

4.4. Discussion of the results

In the following subsections we will first take a look at the assumptions and approximations trying to verify them a posteriori and then we will improve our results.

4.4.1. Looking to the past

Fig. 3 Plot of the errors in velocity and position: $\mathrm{\Delta }{\mathit{v}}_{\mathit{r}}^{\mathrm{RAVE}}$ is the error in the radial velocity as observed in the RAVE survey. Δrhel is the error in the photometric distances as inferred in Zwitter et al. (2010). The binning is 0.5 km s-1 and 50 pc with a scale colour normalized to the highest number of stars per bin (in light yellow) to 0 stars per bin (dark blue contour).

Up to now, we approached the problem of the determination of the thick disk mean velocity and its velocity dispersion tensor. This approach permitted us to determine the kinematic properties of the thick disk component in the framework of an isothermal velocity ellipsoid model (Table 1). Nevertheless, in the literature the ellipsoid is commonly described as a velocity distribution independent of Galactic position because of the difficulties in achieving a more accurate description. However, both these hypotheses are without a rigorous theoretical basis for the thick disk component of our Galaxy. Here, thanks to the RAVE survey, we can investigate this aspect further.

In Fig. 3 we plot the errors in the velocity and position for the sample of stars in the mixture of thin and thick disks as obtained from the catalogue (same stars as in Fig. 1). As already pointed out in Sect. 2.1, we work with the stars with larger errors in the velocity. Thus the bulk of the sample has radial velocity uncertainties up to 6kms^-1, slightly higher than the average error produced by the RAVE survey (e.g., Steinmetz et al. 2006). The photometric distance errors are typically of the order of 0.15 kpc for our selected sample (Zwitter et al. 2010).

A free parameter of our modelling approach is the mean azimuthal velocity of the RAVE sample we select. This quantity is related to the asymmetric drift of the stellar populations, the rotational lag of the thick disk and the circular velocity of rotation. The thick disk subsample of the RAVE catalogue analysed here appears to lag behind the Solar LSR by  ≅ 49   km   s^-1. The general behaviour of the asymmetric drift is not known a priori, especially in its vertical trend and for a mixed set of thin and thick disk populations. For a single stellar population in the disk, the mean circular streaming velocity can theoretically be predicted, starting from the Jeans equation in cylindrical symmetry. We can approximate the term ${{v̅}}_{\mathrm{c}}{}^{\mathrm{\right(}}{x}^{\mathrm{\left)}}$ in Eq. (1)as ${{v̅}}_{\mathrm{c}}\mathrm{(}{{x}}^{\mathrm{)}}\mathrm{=}{{v̅}}_{\mathrm{c}}\mathrm{(}{\mathit{R,z}}^{\mathrm{)}}\mathrm{=}\left(\begin{array}{c}\\ \left|{{v̅}}_{\mathrm{c}}\mathrm{(}{\mathit{R,z}}^{\mathrm{)}}\right|\frac{{\mathit{r}}_{\mathrm{hel}}\cos \mathit{b}\sin \mathit{l}}{\mathit{R}}\\ \left|{{v̅}}_{\mathrm{c}}\mathrm{(}{\mathit{R,z}}^{\mathrm{)}}\right|\frac{{\mathit{R}}_{\mathrm{\odot }}\mathrm{-}{\mathit{r}}_{\mathrm{hel}}\cos \mathit{b}\cos \mathit{l}}{\mathit{R}}\\ \mathrm{0}\end{array}\right)$(23)with $\begin{array}{c}\\ \left|{{v̅}}_{\mathrm{c}}\mathrm{(}{\mathit{R,z}}^{\mathrm{)}}\right|\mathrm{=}[{\mathit{V}}_{\mathrm{LSR}}^{\mathrm{2}}\mathrm{-}\frac{\mathit{\partial }\ln \mathit{\rho }}{\mathit{\partial }\ln \mathit{R}}\mathrm{(}{{\mathit{\sigma }}_{\mathrm{RR}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{Rz}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}\mathrm{(}{{\mathit{\sigma }}_{\mathrm{RR}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathit{\varphi \varphi }}^{\mathrm{2}}}^{\mathrm{)}}\\ \mathrm{+}\mathit{R}\left(\frac{\mathit{\partial }{\mathit{\sigma }}_{\mathrm{RR}}^{\mathrm{2}}}{\mathit{\partial R}}\mathrm{+}\frac{\mathit{\partial }{\mathit{\sigma }}_{\mathrm{Rz}}^{\mathrm{2}}}{\mathit{\partial z}}\mathrm{+}\frac{\mathit{\partial }{\mathrm{\Phi }}_{\mathrm{tot}}\left(\mathit{R,z}\right)}{\mathit{\partial R}}\right){]}^{\mathrm{1}\mathit{/}\mathrm{2}}\end{array}$(24)where Φ_Tot is the total potential of the MW. This indicates how the terms of the velocity dispersion profiles, the circular velocity and the density profile act together to play a crucial role in the prediction of the streaming circular velocity and in its axisymmetric approximation. Assuming that the selected data sample is mostly representative of the thick disk component, the previous Eq. (24)is reduced, to: $\begin{array}{ccc}\sqrt{\left(\frac{\mathit{\partial }\ln \mathit{\rho }\left({\mathit{R}}_{\mathrm{\odot }}\mathit{,}\mathrm{0}\right)}{\mathit{\partial }\ln \mathit{R}}\mathrm{+}\mathrm{1}\right){\mathit{\sigma }}_{\mathit{RR}}^{\mathrm{2}}\mathrm{-}{\mathit{\sigma }}_{\mathit{\phi \phi }}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathrm{c}}^{\mathrm{2}}}\cong \mathrm{201} {\mathrm{kms}}^{-1}& & \end{array}$where h_R ≅ 2.9   kpc is the adopted value for the scale length of a double exponential density profile, e.g., Vallenari et al. (2004). v_c is the circular velocity of the Galaxy in the solar neighbourhood as above, and ${\mathit{\sigma }}_{\mathit{RR}}^{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathit{\phi \phi }}^{\mathrm{2}}$ take the values derived in our paper. The expected value for a pure thick disk component is not so far from the best fit value derived in our study, namely . We point out here that this permits us to impose only an upper limit on the expected value of the streaming circular velocity for the mixed sample of thin and thick disks we analysed. This is thus not necessarily in disagreement with previous studies predicting lower values (e.g., Soubiran et al. 2003) because of the thin disk influence. A further cut in the latitude can in principle reconcile these values but unfortunately reduces the number of stars that we can retain with our iterative procedure and consequently increases greatly the error bars in the velocity dispersion tensor.

Thus we prefer to consider ∥v_c∥ as a free tuning parameter in our approach, remembering that if the thick disk results from a sudden heating due to an infalling satellite, it can be kinematically more decoupled from other components of the MW (e.g., Quinn & Goodman 1986; Walker et al. 1996; Spagna et al. 2010; Bekki & Tsujimoto 2011). The overlap between rotational velocity and peculiar motion of the Sun relative to the LSR can not be disentangled with the use of an isothermal, non-local sample of stars so it remains beyond the scope of the present work (see e.g., Schönrich et al. 2010).

4.4.2. Looking to the future

In principle, it is nowadays possible to disentangle the thick and thin disks using selection criteria based on the chemical properties of the stars (e.g., see the recent review by Nissen 2011). In brief, Neves et al. (2009) determined the trend of abundance ratios as a function of [Fe/H] from 451 high-resolution spectra of F, G, and K main-sequence stars in the solar neighbourhood confirming the long known bimodal distribution of the [α/Fe] ratio for the disk stars (with the thin disk stars less alpha-enhanced than the thick disk ones), see for instance Reddy et al. (2006) or Fuhrmann (2004) and references therein. Along the same line of thought, Casagrande et al. (2011) determined the ratios [Fe/H] and [α/Fe] for 1498 selected stars of the Geneva-Copenhagen Survey from Str$\begin{array}{}{\mathrm{¨}}_{}\\ \mathrm{o}\end{array}$mgren photometry, even if their data do not provide a clear bimodal distribution between thin and thick disk stars as far as the [α/Fe] ratio is concerned. Navarro et al. (2011) investigated the separation of thin and thick disks using a combined index of [Fe/H], [α/Fe] and the heavy element Eu for a sample of 306 stars. Ruchti et al. (2011a) studied the vertical and radial gradients in metallicity and alpha-elements for a sample of selected thick disk stars. Moni Bidin et al. (2012a,b) studied 412 red giant stars in the direction towards the South Galactic Pole trying to describe kinematics, chemistry, and content of dark matter in the MW disks (but see e.g., Bovy & Tremaine 2012). Finally, Afşar et al. (2012) conducting a detailed abundance analysis and atmospheric parameters of 76 stars in the thin and thick disks with very high-resolution spectra (ℜ ≃ 60   000) reaching conclusions similar to those found by other authors.

Common to all the above studies based on chemical properties (see Table 2), is the small number of stars that are considered and the limited spatial coverage in about the solar neighbourhood. In contrast, our method, owing to its statistical nature, first requires a high number of stars such as that produced by RAVE, and second does not need chemical information to work properly. It is worth recalling that data on radial velocities (from medium resolution spectroscopy from which the kinematics is derived) and photometry are easier to gather than good chemical abundances (from high resolution spectroscopy). Our approach can easily be applied to large sets of spectroscopic data such as the continuously growing RAVE survey (or Gaia in the near future). Finally, our method makes use only of kinematic data, leaving the coupling between kinematics and chemistry to be investigated in a forthcoming paper (Boeche et al. 2012).

Although in recent years we have seen continuous improvements in Galactic modelling, see e.g. the Besancon model (Robin et al. 2003), the Padova Galaxy model (used in this paper), DF based models (e.g., Binney 2012) or the Galaxia model (Sharma et al. 2011), we still lack a good theoretical framework coupling kinematics and dynamics with population synthesis and chemistry (e.g., Schönrich & Binney 2009; Just et al. 2011; Just & Jahreiß 2010).

Neglecting the radial and vertical dependence for thethick disk velocity ellipsoid is an almost universal assumption that indeed produces acceptable results for the vertical tilt in our case as well. Nevertheless, we can further explore this working hypothesis by making use of the available photometric distances. Although the proximity of the stars selected grants the validity of our approximations, it limits the exploration of the more distant zones of the thick disk. For every sample statistically representative of the thick disk population that we obtain with the previously outlined method, we split the data set into two regions: within and beyond the Sun’s Galactocentric position, and R > R_⊙ respectively (where R is the cylindrical Galactocentric radius). The values obtained for the velocity dispersion tensor are then statistically averaged, as done to obtain Table 1, and are listed in Tables 3 and 4.

The entries of Tables 3 and 4 hardly show any difference between the stars outside and inside the solar circle. Nonetheless, the suspicion arises that within the errors, there may be some evidence for a general increase in the velocity dispersion of the thick disk along the radial direction (the velocity dispersion is higher for stars inside the solar circle than outside). This implies that for the thick disk alone the classical isothermal picture needs to be improved. The thick disk velocity dispersion tensor seems to imply a dependence on the position in the meridional plane ${\mathit{\sigma }}_{\mathit{ij}}^{\mathrm{thick}}\mathrm{=}{\mathit{\sigma }}_{\mathit{ij}}^{\mathrm{thick}}{}^{\mathrm{\right(}}\mathit{R,z}^{\mathrm{\left)}}$. The effect is, however, small due to the small range of distances sampled with our data. We have higher error bars in the inner sample because of the cut we have applied on the radial velocity errors (see Eq. (6)) that retains a higher number of stars in the anti-centre direction (see Fig. 1). As a consequence, despite the stellar density decrease towards the outer regions of the Galaxy, the error bars of the inner and outer samples of stars are comparable in size.

This suggests a way of evaluating the effect of the distance if a gradient in the vertical dependence of the absolute value of the azimuthal velocity occurs along the thick disk (e.g., Chiba & Beers 2000). Equation (24), which holds for small distances from the plane (e.g., Vallenari et al. 2006, their Eqs. (31) to (33)) can be used to create mock catalogues with radial and vertical gradients in $\left|{{v̅}}_{\mathrm{c}}^{\mathit{i}}{}^{\mathrm{\right(}}\mathit{R,z}^{\mathrm{\left)}}\right|$ where i refers to the thin and thick disk once the populations are embedded in a common potential Φ_tot. In this way we can check the role of small gradients, if any, in ${\mathit{\sigma }}_{\mathit{ij}}^{\mathrm{thick}}\mathrm{=}{\mathit{\sigma }}_{\mathit{ij}}^{\mathrm{thick}}{}^{\mathrm{\right(}}\mathit{R,z}^{\mathrm{\left)}}$ on $\left|{{v̅}}_{\mathrm{c}}^{\mathrm{thick}}{}^{\mathrm{\right(}}\mathit{R,z}^{\mathrm{\left)}}\right|$ or we can simply artificially place a gradient in $\left|{{v̅}}_{\mathrm{c}}^{\mathrm{thick}}{}^{\mathrm{\right(}}\mathit{R,z}^{\mathrm{\left)}}\right|$. This approach has been tested with Eq. (24)and gradients up to 30    kms^-1   kpc^-1. The resulting values for the thick disk velocity ellipsoid are within the error bars presented in Table 1 (see also Bertelli et al. 2003). On the one hand, this is because of the limited vertical and radial extension reached by our sample, which reduces large scale effects. On the other hand, this happens due to the method itself not involving a selection based on chemical abundances. For example, a vertical gradient in mean azimuthal velocity of about 20    kms^-1   kpc^-1 was found by Spagna et al. (2010) even though their result refers to |z| ∈]1,3[  kpc, i.e., a region beyond the range spanned by our sub-sample of the RAVE catalogue.

The different result in the work of Veltz et al. (2008) is due to their data sample being based on fewer than 600 stars in the direction of the south Galactic pole (SGP). This leads to a higher value of σ_RR, on which the influence of halo stars is yet to be investigated (see, e.g., Chiba & Beers 2000; Smith et al. 2009). Particular attention should be paid to the interpretation of the components that are mixed with the azimuthal component: σ_Rφ = 29.4 ± 17.2   kms^-1, σ_φz = 5.8 ± 5.1     kms^-1. The first is reminiscent of the well-known vertex deviation. For the thin disk the vertex deviation can be predicted from σ_Rφ and amounts to an angle of  ≅ $\mathrm{21}{\mathrm{◦}}^{}$. Its large error is probably due to the thick disk components closer to the MW plane which could induce circular velocity gradients that we were unable to detect with our method. The σ_φz component is a further indication of the coupling of the vertical and azimuthal velocity components.

Finally we point out that the method we have developed can be extended to an arbitrary number of populations provided that first the formalism is expanded to include higher order cumulants, second each population is suitably sampled (sufficiently high number of stars), and third each population is characterized by at least one distinct kinematical parameter. In our case, we present evidence of the existence of at least two distinct populations whose characterizing parameters are the second order velocity dispersion tensors. To prove this statement we perform the following experiment. Suppose that the solar position is moved closer towards the Galactic centre, for instance at the position  kpc. Then we set up a mock catalogue centred on the new position of the Sun, $\mathit{R̂}{\mathrm{\odot }}_{}$. By construction, the velocity dispersions of the thin and thick disk stars are the same (no correction for different extinction is applied). Since the kinematical parameters of the populations are identical, the method fails to converge to a solution. This result can be better understood when looking at Fig. 2, where now we imagine that the Sun is located at $\mathit{R̂}{\mathrm{\odot }}_{}$ and the vectors representing the velocity distribution of thick disk stars (light blue arrows) have the same length of those for the thin disk stars (olive-green arrows). Therefore, there is no kinematic way of distinguishing two different populations from the radial velocity distribution of the mixture (the thick green arrows). However, this does imply that the two populations could be separated by considering other parameters such as chemistry, alpha-enhancements etc.

5. Discussion and conclusions

In the sample extracted from the RAVE catalogue that we analysed, two populations of stars with different kinematics are found that correspond to the thin and thick disks, even though our analysis actually focused only on the thick disk.

The major difference with respect to previous studies in the literature is the peculiar spatial volume covered by RAVE, which allows us to determine the kinematics of the thick disk not by extrapolating from a single small field-of-view to the whole thick disk, but by directly measuring the kinematics by studying the much larger volume of the thick disk observed by RAVE (see also Paper II). Over this extension, the key thick disk parameters are determined as:

• 1.

  two components of the solar motion relative to the solar LSR,namely as ${\mathit{v}}_{\mathit{U,}\mathrm{\odot }}\mathrm{=}{}^{\mathrm{\right(}}\mathrm{9.87}\mathrm{\pm }\mathrm{0.37}^{\mathrm{\left)}} {\mathrm{kms}}^{-1}$ and ${\mathit{v}}_{\mathit{W,}\mathrm{\odot }}\mathrm{=}{}^{\mathrm{\right(}}\mathrm{8.01}\mathrm{\pm }\mathrm{0.29}^{\mathrm{\left)}} {\mathrm{kms}}^{-1}$,

• 2.

  the rotational lag of the thick disk component relative to the LSR v_thick,lag ≅ 49 ± 6kms^-1,

• 3.

  the velocity dispersion tensor of the thick disk considered to be an isothermal population: ${\mathit{\sigma }}_{\mathit{RR}}^{}\mathrm{=}\mathrm{(}\mathrm{56.1}\mathrm{\pm }\mathrm{3.8}\mathrm{)} {\mathrm{kms}}^{-1}$, σ_Rϕ = (29.4 ± 17.2)   kms^-1, σ_Rz = (10.1 ± 3.3)   kms^-1, σ_ϕϕ = (46.1 ± 6.7)   kms^-1, σ_ϕz = (5.8 ± 5.1)   kms^-1, σ_zz = (35.1 ± 3.4)   kms^-1.

Moreover, we mention that the missing full determination of the Sun’s velocity vector relative to the LSR is just a choice. This does not mean at all the impossibility for RAVE to characterize completely the solar neighbourhood (see, e.g., Kiss et al. 2010; Coşkunoǧlu et al. 2011). We defer the study of the thin disk kinematics to a companion paper (Paper II).

The determination of the thick disk velocity dispersion tensor and its behaviour in the (O;R,z) plane is just a small step in the investigation of this Galactic component. The presence of a small gradient, especially in the σ_RR component, does not represent the failure of the classical picture of an isothermal description within the range of distances investigated but it is an example of the quality of data that a radial velocity survey like RAVE can provide to confirm, extend or investigate new ideas.

The thick disk is a prominent feature of our galaxy (Gilmore & Reid 1983; Yoshii 1982; Jurić et al. 2008; Yanny et al. 2009) and of external disk galaxies (e.g., van der Kruit & Freeman 2011; Yoachim & Dalcanton 2006). Several possible formation mechanisms have been suggested for the MW thick disk formation. For instance it may be related to the influence of massive satellites that can either heat pre-existing disks or contribute by being accreted (e.g., Quinn & Goodman 1986; Quinn et al. 1993; Abadi et al. 2003). Walker et al. (1996) showed in detail how low-mass satellites, while rapidly sinking into the potential well of a galaxy, could substantially heat a disk. Brook et al. (2004) investigated the influence of gas-rich mergers. Observational evidence of this process is presented by Robin et al. (1996b), Robin et al. (1996a), Gilmore et al. (2002) and Wyse et al. (2006). After the merger, it is plausible that the star formation stopped for a while until the gas assembled again in the thin disk (e.g., see extensive discussions in Freeman & Bland-Hawthorn 2002; Wyse 2002). Finally the MW can produce thick disk features in itself by radial migration processes (Roškar et al. 2008; Schönrich & Binney 2009) or from disruption of massive star clusters (e.g., Kroupa 2002).

In order to distinguish the role of these different scenarios the improved kinematics data expected from the forthcoming astrometric Gaia satellite are fundamental (e.g., Robin et al. 2012), especially in order to constrain time-evolving self-consistent dynamical and chemical models (e.g., Rahimi et al. 2011; Loebman et al. 2011; House et al. 2011; Moni Bidin et al. 2012a). If the thick disk is formed from accreted stars, e.g. during a merger event, then no vertical gradient, $\left[\mathrm{Fe}\mathit{/}\mathrm{H}\right]{}^{\mathrm{\right(}}\mathit{z}^{\mathrm{\left)}}$ is expected for the thick disk (but see Katz et al. 2011). The radial mixing for a MW in isolation (e.g. Schönrich & Binney 2009) is probably not an efficient mechanism to remove the metallicity (or alpha elements) radial gradients (if any, see e.g., Ruchti et al. 2011a) beyond R > 9 kpc because of the steeply decreasing probability of radial migration suggested by Bensby et al. (2011). Radial migration remains still a mechanism to be theoretically defined beyond the mere N-body numerical experiment (Solway et al. 2012; Minchev et al. 2012), while mergers can sensibly flatten the migration probability or also anti-correlate it with the radial density profile (Bird et al. 2012). A correlation between and [Fe/H] would probably disfavour a migration scenario and slow heating mechanisms (e.g., Navarro et al. 2011).

Online material

Appendix A: The whole sky symmetry

Here we outline the procedure to pass from a radial velocity set of data to the true velocity first order moments in the case of an all-sky-survey coverage. This is an application of what is shown in Sect. 3.1 that we used as a reference case. In this case the matrix in Eq. (10)can be computed directly using spherical coordinates . The generic element of the matrix, , ${\left(\mathbf{M̅}\left[\mathrm{2}\right]\right)}_{\mathit{i,j}}\mathrm{=}\left(\begin{array}{c}\\ \mathit{M̅}\left[\mathrm{2}\right]\mathrm{1}\mathit{,}\mathrm{1}& \mathit{...}& \mathit{M̅}\left[\mathrm{2}\right]\mathrm{1}\mathit{,}\mathrm{6}\\ \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}& {\begin{array}{}\mathit{.}\end{array}}^{\mathit{.}}\mathit{.}& \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ \mathit{M̅}\left[\mathrm{2}\right]\mathrm{1}\mathit{,}\mathrm{6}& \mathrm{·}\mathrm{·}\mathrm{·}\hspace{0.17em}& \mathit{M̅}\left[\mathrm{2}\right]\mathrm{6}\mathit{,}\mathrm{6}\end{array}\right)$(A.1)can be computed as $\mathbf{M̅}\left[\mathrm{2}\right]\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}{\mathrm{\int }}_{{\mathit{S}}^{\mathrm{2}}}^{}\mathit{M̅}\left[\mathrm{2}\right]\mathit{i,j}\mathrm{d\Omega }$(A.2)where S² = [0,2π[× [− π/2,π/2[ and the solid angle $\mathrm{d\Omega }\mathrm{=}\mathrm{d}\mathit{ld}{}^{\mathrm{\right(}}\cos \mathit{b}^{\mathrm{\left)}}$. Or, element by element: $\left(\begin{array}{c}\\ \mathit{M̅}\left[\mathrm{2}\right]\mathrm{1}\mathit{,}\mathrm{1}& \mathit{...}& \mathit{M̅}\left[\mathrm{2}\right]\mathrm{1}\mathit{,}\mathrm{6}\\ \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}& {\begin{array}{}\mathit{.}\end{array}}^{\mathit{.}}\mathit{.}& \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ \mathit{M̅}\left[\mathrm{2}\right]\mathrm{1}\mathit{,}\mathrm{6}& \mathrm{·}\mathrm{·}\mathrm{·}\hspace{0.17em}& \mathit{M̅}\left[\mathrm{2}\right]\mathrm{6}\mathit{,}\mathrm{6}\end{array}\right)\mathrm{=}\frac{\mathrm{1}}{\mathrm{5}}\mathit{.}\left(\begin{array}{c}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{3}}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{3}}\\ \mathrm{0}& \frac{\mathrm{2}}{\mathrm{3}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \frac{\mathrm{2}}{\mathrm{3}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \frac{\mathrm{1}}{\mathrm{3}}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{3}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{2}}{\mathrm{3}}& \mathrm{0}\\ \frac{\mathrm{1}}{\mathrm{3}}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{3}}& \mathrm{0}& \mathrm{1}\end{array}\right)\mathit{.}$(A.3)From its inverse we can easily obtain the elements of the velocity ellipsoid as: ${{\sigma }}^{\mathrm{2}}\mathrm{\equiv }\left(\begin{array}{c}\\ {\mathit{\sigma }}_{\mathit{RR}}^{\mathrm{2}}\\ {\mathit{\sigma }}_{\mathit{R\phi }}^{\mathrm{2}}\\ {\mathit{\sigma }}_{\mathit{Rz}}^{\mathrm{2}}\\ {\mathit{\sigma }}_{\mathit{\phi \phi }}^{\mathrm{2}}\\ {\mathit{\sigma }}_{\mathit{\phi z}}^{\mathrm{2}}\\ {\mathit{\sigma }}_{\mathit{zz}}^{\mathrm{2}}\end{array}\right)\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}\left(\begin{array}{c}\\ \mathrm{4}{\mathit{\sigma }}_{\mathrm{\parallel }\mathit{RR}}^{\mathrm{2}}\mathrm{-}{\mathit{\sigma }}_{\mathrm{\parallel }\mathit{\phi \phi }}^{\mathrm{2}}\mathrm{-}{\mathit{\sigma }}_{\mathrm{\parallel }\mathit{zz}}^{\mathrm{2}}\\ \mathrm{5}{\mathit{\sigma }}_{\mathrm{\parallel }\mathit{R\phi }}^{\mathrm{2}}\\ \mathrm{5}{\mathit{\sigma }}_{\mathrm{\parallel }\mathit{Rz}}^{\mathrm{2}}\\ \mathrm{4}{\mathit{\sigma }}_{\mathrm{\parallel }\mathit{\phi \phi }}^{\mathrm{2}}\mathrm{-}{\mathit{\sigma }}_{\mathrm{\parallel }\mathit{RR}}^{\mathrm{2}}\mathrm{-}{\mathit{\sigma }}_{\mathrm{\parallel }\mathit{zz}}^{\mathrm{2}}\\ \mathrm{5}{\mathit{\sigma }}_{\mathrm{\parallel }\mathit{\phi z}}^{\mathrm{2}}\\ \mathrm{4}{\mathit{\sigma }}_{\mathrm{\parallel }\mathit{zz}}^{\mathrm{2}}\mathrm{-}{\mathit{\sigma }}_{\mathrm{\parallel }\mathit{RR}}^{\mathrm{2}}\mathrm{-}{\mathit{\sigma }}_{\mathrm{\parallel }\mathit{\phi \phi }}^{\mathrm{2}}\end{array}\right)\mathit{.}$(A.4)In general, all this is valid when the data have a spherically symmetric distribution. With real data, owing to the partial sky coverage, the matrix (A.3)may substantially differ from the symmetric case. To take this into account, we use the matrix (A.3)as a mask, i.e. as a constraint on the relative weight that the generic matrix element $\mathit{M̅}{\mathrm{\left[}\mathrm{2}\mathrm{\right]}\mathit{i,j}}_{}$ has with respect to another element $\mathit{M̅}{\mathrm{\left[}\mathrm{2}\mathrm{\right]}\mathit{k,l}}_{}$ when trying to maximize the parameter q of Eq. (17), i.e. to control the coupling of the off-diagonal blocks of the operator matrix p. The generic matrix element to be determined requires integrals of many complex trigonometric functions, thus implying long tedious calculations. Fortunately, many elements are null by symmetry. For instance, because there exist only three linearly independent isotropic fourth-rank tensors related to the Kronecker delta tensor by ℐ_ijkl = δ_ilδ_ik, and , the generic symmetric fourth-rank isotropic tensor X _[4]  can be expressed as a linear combination of these (e.g., classically written with coefficient X^s,X^a,X^Tr): ${\mathit{X}}_{\mathrm{\left[}\mathrm{4}\mathrm{\right]}}\mathrm{=}{\mathit{X}}^{\mathit{s}}\left(\frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{ℐ}\mathrm{+}\mathrm{𝒥}\right)\mathrm{-}\frac{\mathrm{1}}{\mathrm{3}}\mathrm{𝒦}\right)\mathrm{+}{\mathit{X}}^{\mathit{a}}\left(\mathrm{ℐ}\mathrm{-}\mathrm{𝒥}\right)\mathrm{+}{\mathit{X}}^{\mathit{Tr}}\frac{\mathrm{1}}{\mathrm{3}}\mathrm{𝒦}\mathit{,}$(A.5)so that only the different terms can be easily singled out from Eq. (A.5). Comparing the shape of the matrix in Eq. (A.3)with the true case of RAVE data from Fig. A.1 shows that in most cases nearly symmetric conditions apply. Indeed, matrices in Eq. (10), , Eq. (11), , and (12), closely resemble the case of spherical symmetry (see, e.g., Fig. A.1 for the matrix of Eq. (10)). For the third moment we have $\mathbf{M̅}\left[\mathrm{3}\right]\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}{\mathrm{\int }}_{{\mathit{S}}^{\mathrm{2}}}^{}\mathit{M̅}\left[\mathrm{3}\right]\mathit{i,j}\mathrm{d\Omega }$(A.6)to get $\left(\begin{array}{c}\\ \mathit{M̅}\left[\mathrm{3}\right]\mathrm{1}\mathit{,}\mathrm{1}& \mathit{...}& \mathit{M̅}\left[\mathrm{3}\right]\mathrm{1}\mathit{,}\mathrm{10}\\ \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}& {\begin{array}{}\mathit{.}\end{array}}^{\mathit{.}}\mathit{.}& \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ \mathit{M̅}\left[\mathrm{3}\right]\mathrm{1}\mathit{,}\mathrm{10}& \mathrm{·}\mathrm{·}\mathrm{·}\hspace{0.17em}& \mathit{M̅}\left[\mathrm{3}\right]\mathrm{10}\mathit{,}\mathrm{10}\end{array}\right)\mathrm{=}\frac{\mathrm{1}}{\mathrm{7}}\mathrm{·}\left(\begin{array}{c}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{2}}{\mathrm{5}}& \mathrm{0}& \frac{\mathrm{2}}{\mathrm{5}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \frac{\mathrm{3}}{\mathrm{5}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{2}}{\mathrm{5}}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{5}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \frac{\mathrm{3}}{\mathrm{5}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{5}}& \mathrm{0}& \frac{\mathrm{2}}{\mathrm{5}}\\ \frac{\mathrm{2}}{\mathrm{5}}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{3}}{\mathrm{5}}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{5}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right)\mathit{,}$(A.7)

Fig. A.1 The inverse matrix of the Eq. (10)for the best fit model. The matrix closely resembles the fully analytical case except for the numerical part. Different shades of colour are applied to visualize the symmetries.

and from its inverse $\mathbf{S}\mathrm{\equiv }\left(\begin{array}{c}\\ {\mathit{S}}_{\mathit{RRR}}\\ {\mathit{S}}_{\mathit{RR\phi }}\\ {\mathit{S}}_{\mathit{RRz}}\\ {\mathit{S}}_{\mathit{R\phi \phi }}\\ {\mathit{S}}_{\mathit{R\phi z}}\\ \hspace{0.17em}\begin{array}{}\begin{array}{}\end{array}{\mathit{S}}_{\mathit{Rzz}}\end{array}\end{array}\right)\mathrm{=}\frac{\mathrm{35}}{\mathrm{24}}\mathit{.}\left(\begin{array}{c}\\ \mathrm{4}\mathrm{(}{\mathrm{2}{\mathit{S}}_{\mathrm{\parallel }\mathit{RRR}}\mathrm{-}{\mathit{S}}_{\mathrm{\parallel }\mathit{R\phi \phi }}\mathrm{-}{\mathit{S}}_{\mathrm{\parallel }\mathit{Rzz}}}^{\mathrm{)}}\\ \mathrm{11}{\mathit{S}}_{\mathrm{\parallel }\mathit{RR\phi }}\mathrm{-}\mathrm{4}{\mathit{S}}_{\mathrm{\parallel }\mathit{\phi \phi \phi }}\mathrm{-}{\mathit{S}}_{\mathrm{\parallel }\mathit{\phi zz}}\\ \mathrm{11}{\mathit{S}}_{\mathrm{\parallel }\mathit{RRz}}\mathrm{-}{\mathit{S}}_{\mathrm{\parallel }\mathit{\phi \phi z}}\mathrm{-}\mathrm{4}{\mathit{S}}_{\mathrm{\parallel }\mathit{zzz}}\\ \mathrm{11}{\mathit{S}}_{\mathrm{\parallel }\mathit{R\phi \phi }}\mathrm{-}{\mathit{S}}_{\mathrm{\parallel }\mathit{Rzz}}\mathrm{-}\mathrm{4}{\mathit{S}}_{\mathrm{\parallel }\mathit{RRR}}\end{array}\right)\mathit{.}$(A.8)From $\mathbf{M̅}\left[\mathrm{4}\right]\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}{\mathrm{\int }}_{{\mathit{S}}^{\mathrm{2}}}^{}\mathit{M̅}\left[\mathrm{4}\right]\mathit{i,j}\mathrm{d\Omega }\mathit{,}$(A.9)we get $\left(\begin{array}{c}\\ \mathit{M̅}\left[\mathrm{4}\right]\mathrm{1}\mathit{,}\mathrm{1}& \mathit{...}& \mathit{M̅}\left[\mathrm{4}\right]\mathrm{1}\mathit{,}\mathrm{15}\\ \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}& {\begin{array}{}\mathit{.}\end{array}}^{\mathit{.}}\mathit{.}& \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ \mathit{M̅}\left[\mathrm{4}\right]\mathrm{1}\mathit{,}\mathrm{15}& \mathrm{·}\mathrm{·}\mathrm{·}\hspace{0.17em}& \mathit{M̅}\left[\mathrm{4}\right]\mathrm{15}\mathit{,}\mathrm{15}\end{array}\right)\mathrm{=}\frac{\mathrm{1}}{\mathrm{9}}\mathrm{·}\left(\begin{array}{c}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{2}}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{3}}{\mathrm{35}}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{10}}& \mathrm{0}& \frac{\mathrm{3}}{\mathrm{35}}\\ \mathrm{0}& \frac{\mathrm{4}}{\mathrm{7}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{12}}{\mathrm{35}}& \mathrm{0}& \frac{\mathrm{8}}{\mathrm{35}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \frac{\mathrm{4}}{\mathrm{7}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{8}}{\mathrm{35}}& \mathrm{0}& \frac{\mathrm{12}}{\mathrm{35}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right)\mathit{,}$(A.10)and from its inverse $\begin{array}{ccc}& & \mathbf{T}\mathrm{\equiv }{\mathrm{(}{{\mathit{T}}_{\mathit{RRRR}}\mathit{,}{\mathit{T}}_{\mathit{RRR\phi }}\mathit{,}{\mathit{T}}_{\mathit{RRRz}}\mathit{,}{\mathit{T}}_{\mathit{RR\phi \phi }}\mathit{,}{\mathit{T}}_{\mathit{RR\phi z}}\mathit{,}{\mathit{T}}_{\mathit{RRzz}}\mathit{,}{\mathit{T}}_{\mathit{R\phi \phi \phi }}\mathit{,}{\mathit{T}}_{\mathit{R\phi \phi z}}\mathit{,}{\mathit{T}}_{\mathit{R\phi zz}}\mathit{,}{\mathit{T}}_{\mathit{Rzzz}}\mathit{,}{\mathit{T}}_{\mathit{\phi \phi \phi \phi }}\mathit{,}{\mathit{T}}_{\mathit{\phi \phi \phi z}}\mathit{,}{\mathit{T}}_{\mathit{\phi \phi zz}}\mathit{,}{\mathit{T}}_{\mathit{\phi zzz}}\mathit{,}{\mathit{T}}_{\mathit{zzzz}}}^{\mathrm{)}}}^{\mathit{T}}\mathrm{=}\\ & & \end{array}$(A.11)

Appendix B: The V-cumulants mixture distribution

To disentangle the first and second cumulants of the thick disk alone, we need up to the fourth order cumulants of the mixture. The V-cumulants mixture distribution of third and fourth order are introduced here. These moments exhibit a symmetry along the azimuthal velocity vector of the mixture. Hence they differ from what is laid out in the Appendices B, C, and D of Cubarsi & Alcobé (2004). Also the notation differs in order to be consistent with the notation used in our study. The third cumulants from which we can easily compute the two components of the normalized vector defined in Eq. (17), are $\hspace{0.17em}\begin{array}{}\begin{array}{}\end{array}{\mathrm{K}}_{\mathit{RRR}}\mathrm{=}{\mathit{\kappa }}_{\mathit{zzz}}\mathrm{-}\mathit{d̂}\begin{array}{}\mathrm{3}\\ \mathrm{3}\end{array}{\mathit{\kappa }}_{\mathit{\phi \phi \phi }}\mathrm{+}\mathrm{3}\mathit{d̂}\begin{array}{}\mathrm{2}\\ \mathrm{3}\end{array}{\mathit{\kappa }}_{\mathit{\phi \phi z}}\mathrm{-}\mathrm{3}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{\phi zz}}\mathrm{=}\mathrm{0}\mathit{,}\\ \begin{array}{}\end{array}{\mathrm{K}}_{\mathit{RRz}}\mathrm{=}\mathit{d̂}\mathrm{1}{\mathit{\kappa }}_{\mathit{\phi zz}}\mathrm{-}{\mathit{\kappa }}_{\mathit{Rzz}}\mathrm{+}\mathit{d̂}\mathrm{3}\mathrm{(}{\mathrm{-}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{R\phi \phi }}\mathrm{+}\mathrm{2}{\mathit{\kappa }}_{\mathit{R\phi z}}\mathrm{+}\mathit{d̂}\mathrm{1}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{\phi \phi \phi }}\mathrm{-}\mathrm{2}\mathit{d̂}\mathrm{1}{\mathit{\kappa }}_{\mathit{\phi \phi z}}}^{\mathrm{)}}\mathrm{=}\mathrm{0}\mathit{,}\\ \begin{array}{}\end{array}{\mathrm{K}}_{\mathit{Rzz}}\mathrm{=}{\mathit{\kappa }}_{\mathit{RRz}}\mathrm{-}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{RR\phi }}\mathrm{+}\mathit{d̂}\mathrm{1}\mathrm{(}{\mathrm{2}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{R\phi \phi }}\mathrm{-}\mathrm{2}{\mathit{\kappa }}_{\mathit{R\phi z}}\mathrm{-}\mathit{d̂}\mathrm{1}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{\phi \phi \phi }}\mathrm{+}\mathit{d̂}\mathrm{1}{\mathit{\kappa }}_{\mathit{\phi \phi z}}}^{\mathrm{)}}\mathrm{=}\mathrm{0}\mathit{,}\\ \begin{array}{}\end{array}{\mathrm{K}}_{\mathit{zzz}}\mathrm{=}\mathrm{-}{\mathit{\kappa }}_{\mathit{RRR}}\mathrm{+}\mathit{d̂}\mathrm{1}\mathrm{(}{\mathrm{3}{\mathit{\kappa }}_{\mathit{RR\phi }}\mathrm{+}\mathit{d̂}\mathrm{1}\mathrm{(}{\mathrm{-}\mathrm{3}{\mathit{\kappa }}_{\mathit{R\phi \phi }}\mathrm{+}\mathit{d̂}\mathrm{1}{\mathit{\kappa }}_{\mathit{\phi \phi \phi }}}^{\mathrm{)}}}^{\mathrm{)}}\mathrm{=}\mathrm{0.}\end{array}$(B.1)Once the values of and are computed as an SVD solution of the previous overdetermined system, the remaining cumulants can be calculated via $\hspace{0.17em}\begin{array}{}\begin{array}{}\end{array}{\mathrm{K}}_{\mathit{RR\phi }}\mathrm{=}{\mathit{\kappa }}_{\mathit{\phi zz}}\mathrm{-}\mathrm{2}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{\phi \phi z}}\mathrm{+}\mathit{d̂}\begin{array}{}\mathrm{2}\\ \mathrm{3}\end{array}{\mathit{\kappa }}_{\mathit{\phi \phi \phi }}\mathit{,}\\ \begin{array}{}\end{array}{\mathrm{K}}_{\mathit{R\phi \phi }}\mathrm{=}{\mathit{\kappa }}_{\mathit{\phi \phi z}}\mathrm{-}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{\phi \phi \phi }}\mathit{,}\\ \begin{array}{}\end{array}{\mathrm{K}}_{\mathit{R\phi z}}\mathrm{=}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{R\phi \phi }}\mathrm{-}{\mathit{\kappa }}_{\mathit{R\phi z}}\mathrm{+}\mathit{d̂}\mathrm{1}\mathrm{(}{{\mathit{\kappa }}_{\mathit{\phi \phi z}}\mathrm{-}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{\phi \phi \phi }}}^{\mathrm{)}}\mathit{,}\\ \begin{array}{}\end{array}{\mathrm{K}}_{\mathit{\phi \phi \phi }}\mathrm{=}{\mathit{\kappa }}_{\mathit{\phi \phi \phi }}\mathit{,}\\ \begin{array}{}\end{array}{\mathrm{K}}_{\mathit{\phi \phi z}}\mathrm{=}\mathit{d̂}\mathrm{1}{\mathit{\kappa }}_{\mathit{\phi \phi \phi }}\mathrm{-}{\mathit{\kappa }}_{\mathit{R\phi \phi }}\mathit{,}\\ \begin{array}{}\end{array}{\mathrm{K}}_{\mathit{\phi zz}}\mathrm{=}{\mathit{\kappa }}_{\mathit{RR\phi }}\mathrm{+}\mathit{d̂}\mathrm{1}\mathrm{(}{\mathrm{-}\mathrm{2}{\mathit{\kappa }}_{\mathit{R\phi \phi }}\mathrm{+}\mathit{d̂}\mathrm{1}{\mathit{\kappa }}_{\mathit{\phi \phi \phi }}}^{\mathrm{)}}\mathit{.}\end{array}$(B.2)The fourth order equation can be computed as $\hspace{0.17em}\begin{array}{}\begin{array}{}\end{array}{\mathrm{K}}_{\mathit{RRRR}}\mathrm{=}{\mathit{\kappa }}_{\mathit{RRRR}}\mathrm{+}\mathit{d̂}\mathrm{3}\mathrm{(}{\mathrm{-}\mathrm{4}{\mathit{\kappa }}_{\mathit{\phi zzz}}\mathrm{+}\mathit{d̂}\mathrm{3}\mathrm{(}{\mathrm{6}{\mathit{\kappa }}_{\mathit{\phi \phi zz}}\mathrm{-}\mathrm{4}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{\phi \phi \phi z}}\mathrm{+}\mathit{d̂}\begin{array}{}\mathrm{2}\\ \mathrm{3}\end{array}{\mathit{\kappa }}_{\mathit{\phi \phi \phi \phi }}}^{\mathrm{)}}}^{\mathrm{)}}\mathit{,}\\ \begin{array}{}\end{array}{\mathrm{K}}_{\mathit{RRRz}}\mathrm{=}\mathrm{-}{\mathit{\kappa }}_{\mathit{Rzzz}}\mathrm{+}\mathit{d̂}\mathrm{1}{\mathit{\kappa }}_{\mathit{\phi zzz}}\mathrm{+}\\ \begin{array}{}\end{array}\mathit{d̂}\mathrm{3}\mathrm{(}{\mathrm{3}{\mathit{\kappa }}_{\mathit{R\phi zz}}\mathrm{-}\mathrm{3}\mathit{d̂}\mathrm{1}{\mathit{\kappa }}_{\mathit{\phi \phi zz}}\mathrm{+}\mathit{d̂}\mathrm{3}\mathrm{(}{\mathrm{-}\mathrm{3}{\mathit{\kappa }}_{\mathit{R\phi \phi z}}\mathrm{+}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{R\phi \phi \phi }}\mathrm{+}\mathrm{3}\mathit{d̂}\mathrm{1}{\mathit{\kappa }}_{\mathit{\phi \phi \phi z}}\mathrm{-}\mathit{d̂}\mathrm{1}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{\phi \phi \phi \phi }}}^{\mathrm{)}}}^{\mathrm{)}}\mathit{,}\\ \begin{array}{}\end{array}{\mathrm{K}}_{\mathit{RRzz}}\mathrm{=}{\mathit{\kappa }}_{\mathit{RRzz}}\mathrm{-}\mathrm{2}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{RR\phi z}}\mathrm{+}\mathit{d̂}\begin{array}{}\mathrm{2}\\ \mathrm{3}\end{array}{\mathit{\kappa }}_{\mathit{RR\phi \phi }}\mathrm{+}\\ \begin{array}{}\end{array}\mathit{d̂}\mathrm{1}\mathrm{(}{\mathrm{-}\mathrm{2}{\mathit{\kappa }}_{\mathit{R\phi zz}}\mathrm{+}\mathit{d̂}\mathrm{1}{\mathit{\kappa }}_{\mathit{\phi \phi zz}}\mathrm{+}\mathit{d̂}\mathrm{3}\mathrm{(}{\mathrm{4}{\mathit{\kappa }}_{\mathit{R\phi \phi z}}\mathrm{-}\mathrm{2}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{R\phi \phi \phi }}\mathrm{-}\mathrm{2}\mathit{d̂}\mathrm{1}{\mathit{\kappa }}_{\mathit{\phi \phi \phi z}}\mathrm{+}\mathit{d̂}\mathrm{1}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{\phi \phi \phi \phi }}}^{\mathrm{)}}}^{\mathrm{)}}\mathit{,}\\ \begin{array}{}\end{array}{\mathrm{K}}_{\mathit{Rzzz}}\mathrm{=}\mathrm{-}{\mathit{\kappa }}_{\mathit{RRRz}}\mathrm{+}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{RRR\phi }}\mathrm{+}\\ \begin{array}{}\end{array}\mathit{d̂}\mathrm{1}\mathrm{(}{\mathrm{3}{\mathit{\kappa }}_{\mathit{RR\phi z}}\mathrm{-}\mathrm{3}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{RR\phi \phi }}\mathrm{+}\mathit{d̂}\mathrm{1}\mathrm{(}{\mathrm{-}\mathrm{3}{\mathit{\kappa }}_{\mathit{R\phi \phi z}}\mathrm{+}\mathrm{3}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{R\phi \phi \phi }}\mathrm{+}\mathit{d̂}\mathrm{1}{\mathit{\kappa }}_{\mathit{\phi \phi \phi z}}\mathrm{-}\mathit{d̂}\mathrm{1}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{\phi \phi \phi \phi }}}^{\mathrm{)}}}^{\mathrm{)}}\mathit{,}\\ \begin{array}{}\end{array}{\mathrm{K}}_{\mathit{zzzz}}\mathrm{=}{\mathit{\kappa }}_{\mathit{RRRR}}\mathrm{+}\mathit{d̂}\mathrm{1}\mathrm{(}{\mathrm{-}\mathrm{4}{\mathit{\kappa }}_{\mathit{RRR\phi }}\mathrm{+}\mathit{d̂}\mathrm{1}\mathrm{(}{\mathrm{6}{\mathit{\kappa }}_{\mathit{RR\phi \phi }}\mathrm{-}\mathrm{4}\mathit{d̂}\mathrm{1}{\mathit{\kappa }}_{\mathit{R\phi \phi \phi }}\mathrm{+}\mathit{d̂}\begin{array}{}\mathrm{2}\\ \mathrm{1}\end{array}{\mathit{\kappa }}_{\mathit{\phi \phi \phi \phi }}}^{\mathrm{)}}}^{\mathrm{)}}\mathit{,}\\ \begin{array}{}\end{array}{\mathrm{K}}_{\mathit{RRR\phi }}\mathrm{=}{\mathit{\kappa }}_{\mathit{\phi zzz}}\mathrm{-}\mathit{d̂}\mathrm{3}\mathrm{(}{\mathrm{3}{\mathit{\kappa }}_{\mathit{\phi \phi zz}}\mathrm{+}\mathit{d̂}\mathrm{3}\mathrm{(}{\mathrm{-}\mathrm{3}{\mathit{\kappa }}_{\mathit{\phi \phi \phi z}}\mathrm{+}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{\phi \phi \phi \phi }}}^{\mathrm{)}}}^{\mathrm{)}}\mathit{,}\\ \begin{array}{}\end{array}{\mathrm{K}}_{\mathit{RR\phi z}}\mathrm{=}\mathrm{-}{\mathit{\kappa }}_{\mathit{R\phi zz}}\mathrm{+}\mathit{d̂}\mathrm{1}{\mathit{\kappa }}_{\mathit{\phi \phi zz}}\mathrm{+}\mathit{d̂}\mathrm{3}\mathrm{(}{\mathrm{2}{\mathit{\kappa }}_{\mathit{R\phi \phi z}}\mathrm{-}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{R\phi \phi \phi }}\mathrm{-}\mathrm{2}\mathit{d̂}\mathrm{1}{\mathit{\kappa }}_{\mathit{\phi \phi \phi z}}\mathrm{+}\mathit{d̂}\mathrm{1}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{\phi \phi \phi \phi }}}^{\mathrm{)}}\mathit{,}\\ \begin{array}{}\end{array}{\mathrm{K}}_{\mathit{R\phi zz}}\mathrm{=}{\mathit{\kappa }}_{\mathit{RR\phi z}}\mathrm{-}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{RR\phi \phi }}\mathrm{+}\mathit{d̂}\mathrm{1}\mathrm{(}{\mathrm{-}\mathrm{2}{\mathit{\kappa }}_{\mathit{R\phi \phi z}}\mathrm{+}\mathrm{2}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{R\phi \phi \phi }}\mathrm{+}\mathit{d̂}\mathrm{1}{\mathit{\kappa }}_{\mathit{\phi \phi \phi z}}\mathrm{-}\mathit{d̂}\mathrm{1}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{\phi \phi \phi \phi }}}^{\mathrm{)}}\mathit{,}\\ \begin{array}{}\end{array}{\mathrm{K}}_{\mathit{\phi zzz}}\mathrm{=}\mathrm{-}{\mathit{\kappa }}_{\mathit{RRR\phi }}\mathrm{+}\mathit{d̂}\mathrm{1}\mathrm{(}{\mathrm{3}{\mathit{\kappa }}_{\mathit{RR\phi \phi }}\mathrm{+}\mathit{d̂}\mathrm{1}\mathrm{(}{\mathrm{-}\mathrm{3}{\mathit{\kappa }}_{\mathit{R\phi \phi \phi }}\mathrm{+}\mathit{d̂}\mathrm{1}{\mathit{\kappa }}_{\mathit{\phi \phi \phi \phi }}}^{\mathrm{)}}}^{\mathrm{)}}\mathit{,}\\ \begin{array}{}\end{array}{\mathrm{K}}_{\mathit{RR\phi \phi }}\mathrm{=}{\mathit{\kappa }}_{\mathit{\phi \phi zz}}\mathrm{+}\mathit{d̂}\mathrm{3}\mathrm{(}{\mathrm{-}\mathrm{2}{\mathit{\kappa }}_{\mathit{\phi \phi \phi z}}\mathrm{+}{\mathit{d̂}}_{\mathrm{3}}{\mathit{\kappa }}_{\mathit{\phi \phi \phi \phi }}}^{\mathrm{)}}\mathit{,}\\ \begin{array}{}\end{array}{\mathrm{K}}_{\mathit{R\phi \phi z}}\mathrm{=}\mathrm{-}{\mathit{\kappa }}_{\mathit{R\phi \phi z}}\mathrm{+}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{R\phi \phi \phi }}\mathrm{+}\mathit{d̂}\mathrm{1}\mathrm{(}{{\mathit{\kappa }}_{\mathit{\phi \phi \phi z}}\mathrm{-}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{\phi \phi \phi \phi }}}^{\mathrm{)}}\mathit{,}\\ \begin{array}{}\end{array}{\mathrm{K}}_{\mathit{\phi \phi zz}}\mathrm{=}{\mathit{\kappa }}_{\mathit{RR\phi \phi }}\mathrm{+}\mathit{d̂}\mathrm{1}\mathrm{(}{\mathrm{-}\mathrm{2}{\mathit{\kappa }}_{\mathit{R\phi \phi \phi }}\mathrm{+}\mathit{d̂}\mathrm{1}{\mathit{\kappa }}_{\mathit{\phi \phi \phi \phi }}}^{\mathrm{)}}\mathit{,}\\ \begin{array}{}\end{array}{\mathrm{K}}_{\mathit{R\phi \phi \phi }}\mathrm{=}{\mathit{\kappa }}_{\mathit{\phi \phi \phi z}}\mathrm{-}\mathit{d̂}\mathrm{3}{\mathit{\kappa }}_{\mathit{\phi \phi \phi \phi }}\mathit{,}\\ \begin{array}{}\end{array}{\mathrm{K}}_{\mathit{\phi \phi \phi z}}\mathrm{=}\mathrm{-}{\mathit{\kappa }}_{\mathit{R\phi \phi \phi }}\mathrm{+}\mathit{d̂}\mathrm{1}{\mathit{\kappa }}_{\mathit{\phi \phi \phi \phi }}\mathit{,}\\ \begin{array}{}\end{array}{\mathrm{K}}_{\mathit{\phi \phi \phi \phi }}\mathrm{=}{\mathit{\kappa }}_{\mathit{\phi \phi \phi \phi }}\mathit{.}\end{array}$(B.3)Once we have the whole set of the cumulants we define the tensor ${C}\mathrm{\equiv }\frac{\mathrm{1}}{\sqrt{{\mathit{q}}^{\mathrm{2}}\mathrm{+}\mathrm{4}}}\mathrm{(}{{{\sigma }}_{\mathrm{thick}}\mathrm{-}{{\sigma }}_{\mathrm{thin}}}^{\mathrm{)}}\mathrm{-}\mathit{q}{{d}}^{\mathrm{\otimes }\mathrm{2}}\mathit{,}$(B.4)where . The constraining equations derived in Cubarsi & Alcobé (2004) can be reduced to the following set of fourteen scalar relations $\hspace{0.17em}\begin{array}{}\begin{array}{}\end{array}\frac{\mathrm{3}{\mathrm{K}}_{\mathit{RR\phi }}^{\mathrm{2}}}{{\mathrm{K}}_{\mathit{RRRR}}}\mathrm{-}\mathit{d̂}\begin{array}{}\mathrm{2}\\ \mathrm{2}\end{array}\mathrm{=}\frac{\mathrm{3}{\mathrm{K}}_{\mathit{RR\phi }}{\mathrm{K}}_{\mathit{R\phi z}}}{{\mathrm{K}}_{\mathit{RRRz}}}\mathrm{-}\mathit{d̂}\begin{array}{}\mathrm{2}\\ \mathrm{2}\end{array}\mathrm{=}\frac{{\mathrm{K}}_{\mathit{RR\phi }}{\mathrm{K}}_{\mathit{zz\phi }}\mathrm{+}\mathrm{2}{\mathrm{K}}_{\mathit{R\phi z}}^{\mathrm{2}}}{{\mathrm{K}}_{\mathit{RRzz}}}\mathrm{-}\mathit{d̂}\begin{array}{}\mathrm{2}\\ \mathrm{2}\end{array}\mathrm{=}\frac{\mathrm{3}{\mathrm{K}}_{\mathit{R\phi z}}{\mathrm{K}}_{\mathit{\phi zz}}}{{\mathrm{K}}_{\mathit{Rzzz}}}\mathrm{-}\mathit{d̂}\begin{array}{}\mathrm{2}\\ \mathrm{2}\end{array}\mathrm{=}\\ \begin{array}{}\end{array}\frac{\mathrm{3}{\mathrm{K}}_{\mathit{\phi zz}}^{\mathrm{2}}}{{\mathrm{K}}_{\mathit{zzzz}}}\mathrm{-}\mathit{d̂}\begin{array}{}\mathrm{2}\\ \mathrm{2}\end{array}\mathrm{=}\frac{\mathrm{3}{\mathrm{K}}_{\mathit{RR\phi }}{\mathrm{K}}_{\mathit{R\phi \phi }}}{{\mathrm{K}}_{\mathit{RRR\phi }}}\mathrm{-}\mathit{d̂}\begin{array}{}\mathrm{2}\\ \mathrm{2}\end{array}\mathrm{=}\frac{{\mathrm{K}}_{\mathit{RR\phi }}{\mathrm{K}}_{\mathit{z\phi \phi }}\mathrm{+}\mathrm{2}{\mathrm{K}}_{\mathit{R\phi z}}{\mathrm{K}}_{\mathit{R\phi \phi }}}{{\mathrm{K}}_{\mathit{RR\phi z}}}\mathrm{-}\mathit{d̂}\begin{array}{}\mathrm{2}\\ \mathrm{2}\end{array}\mathrm{=}\\ \begin{array}{}\end{array}\frac{{\mathrm{K}}_{\mathit{zz\phi }}{\mathrm{K}}_{\mathit{R\phi \phi }}\mathrm{+}\mathrm{2}{\mathrm{K}}_{\mathit{R\phi z}}{\mathrm{K}}_{\mathit{z\phi \phi }}}{{\mathrm{K}}_{\mathit{R\phi zz}}}\mathrm{-}\mathit{d̂}\begin{array}{}\mathrm{2}\\ \mathrm{2}\end{array}\mathrm{=}\frac{\mathrm{3}{\mathrm{K}}_{\mathit{zz\phi }}{\mathrm{K}}_{\mathit{\phi \phi z}}}{{\mathrm{K}}_{\mathit{\phi zzz}}}\mathrm{-}\mathit{d̂}\begin{array}{}\mathrm{2}\\ \mathrm{2}\end{array}\mathrm{=}\frac{\mathrm{1}}{{\mathrm{K}}_{\mathit{RR\phi }}}\left({\mathrm{K}}_{\mathit{RR\phi \phi }}\mathrm{-}\frac{\mathrm{2}{\mathrm{K}}_{\mathit{R\phi \phi }}^{\mathrm{2}}}{\mathit{d̂}\begin{array}{}\mathrm{2}\\ \mathrm{2}\end{array}}\right)\frac{\mathit{d̂}\begin{array}{}\mathrm{3}\\ \mathrm{2}\end{array}}{{\mathit{c}}_{\mathit{\phi \phi }}}\mathrm{-}\mathit{d̂}\begin{array}{}\mathrm{2}\\ \mathrm{2}\end{array}\mathrm{=}\\ \begin{array}{}\end{array}\frac{\mathrm{1}}{{\mathrm{K}}_{\mathit{R\phi z}}}\left({\mathrm{K}}_{\mathit{Rz\phi \phi }}\mathrm{-}\frac{\mathrm{2}{\mathrm{K}}_{\mathit{R\phi \phi }}{\mathrm{K}}_{\mathit{\phi \phi z}}}{{\mathit{d̂}}_{\mathrm{2}}^{\mathrm{2}}}\right)\frac{\mathit{d̂}\begin{array}{}\mathrm{3}\\ \mathrm{2}\end{array}}{{\mathit{c}}_{\mathit{\phi \phi }}}\mathrm{-}\mathit{d̂}\begin{array}{}\mathrm{2}\\ \mathrm{2}\end{array}\mathrm{=}\frac{\mathrm{1}}{{\mathrm{K}}_{\mathit{zz\phi }}}\left({\mathrm{K}}_{\mathit{zz\phi \phi }}\mathrm{-}\frac{\mathrm{2}{\mathrm{K}}_{\mathit{z\phi \phi }}^{\mathrm{3}}}{{\mathit{d̂}}_{\mathrm{2}}^{\mathrm{2}}}\right)\frac{\mathit{d̂}\begin{array}{}\mathrm{3}\\ \mathrm{2}\end{array}}{{\mathit{c}}_{\mathit{\phi \phi }}}\mathrm{-}\mathit{d̂}\begin{array}{}\mathrm{2}\\ \mathrm{2}\end{array}\mathrm{=}\\ \begin{array}{}\end{array} \frac{{\mathrm{K}}_{\mathit{R\phi \phi \phi }}}{\mathrm{3}{\mathrm{K}}_{\mathit{R\phi \phi }}}\frac{\mathit{d̂}\begin{array}{}\mathrm{3}\\ \mathrm{2}\end{array}}{{\mathit{c}}_{\mathit{\phi \phi }}}\mathrm{-}\mathit{d̂}\begin{array}{}\mathrm{2}\\ \mathrm{2}\end{array}\mathrm{=}\frac{{\mathrm{K}}_{\mathit{z\phi \phi \phi }}}{\mathrm{3}{\mathrm{K}}_{\mathit{z\phi \phi }}}\frac{\mathit{d̂}\begin{array}{}\mathrm{3}\\ \mathrm{2}\end{array}}{{\mathit{c}}_{\mathit{\phi \phi }}}\mathrm{-}\mathit{d̂}\begin{array}{}\mathrm{2}\\ \mathrm{2}\end{array}\mathit{,}\end{array}$(B.5)that we solve in a least-squares sense with respect to the elements c_φφ and . The final step to calculate the desired results for q is to take the derived values for c_φφ and and to work out the last constraining equation from the relations $\hspace{0.17em}\begin{array}{}\begin{array}{}\end{array}{\mathrm{K}}_{\mathit{zzz}}\mathrm{=}{\mathit{\kappa }}_{\mathit{zzz}}\mathrm{=}\mathrm{3}{\mathit{c}}_{\mathit{\phi \phi }}\mathit{d̂}\mathrm{2}\mathrm{+}\mathrm{2}\mathit{q}\mathit{d̂}\begin{array}{}\mathrm{3}\\ \mathrm{2}\end{array}\mathit{,}\\ \begin{array}{}\end{array}{\mathrm{K}}_{\mathit{zzzz}}\mathrm{=}{\mathit{\kappa }}_{\mathit{zzzz}}\mathrm{=}\mathrm{3}{\mathit{c}}_{\mathit{\phi \phi }}^{\mathrm{2}}\mathrm{-}\mathrm{2}\mathrm{(}{{\mathit{q}}^{\mathrm{2}}\mathrm{+}\mathrm{1}}^{\mathrm{)}}\mathit{d̂}\begin{array}{}\mathrm{4}\\ \mathrm{2}\end{array}\mathit{.}\end{array}$(B.6)