A&A 427, 131-144 (2004)
DOI: 10.1051/0004-6361:20041144

Cumulants and symmetries in a trivariate normal mixture

A qualitative study of the local velocity distribution

R. Cubarsi - S. Alcobé

Dept. Matemàtica Aplicada IV, Universitat Politècnica de Catalunya, 08034 Barcelona, Catalonia, Spain

Received 22 April 2004 / Accepted 21 July 2004

Abstract
An alternative parameterization of a trivariate normal mixture for the stellar velocity distribution provides a set of constraint equations between global cumulants that are used to estimate characteristic constants of the mixture and population parameters. Under particular distribution symmetries these relationships become simpler and easy to evaluate, and they are used to test some meaningful mixtures. This method for the analysis of trivariate normal mixtures is applied to local star catalogs with known velocity space in order to find Gaussian components and underlying distribution symmetries. A large local sample (13 678 stars, 300 pc) obtained from HIPPARCOS catalog allows us to identify two mixed velocity ellipsoids with parameters corresponding to thin and thick disk populations, similar to those of samples selected from non-kinematical criteria. In addition to the usual assumptions of distribution symmetry plane and non-significant differential movement in the radial direction, our analysis also detects a slight but clear loss of axial symmetry associated with the vertex deviation of both population ellipsoids.

Key words: Galaxy: kinematics and dynamics - stars: Population II - methods: analytical - methods: statistical - Galaxy: solar neighbourhood

1 Introduction

In galactic dynamics a stellar system, with an enough large number of stars, may be represented by a continuous mass distribution. From properties of the stellar distribution function it is possible to transform local knowledge of the velocity distribution into knowledge about densities and velocities at other points of the Galaxy. This information is condensed in some fundamental equations explicitly depending on the moments of the stellar velocity distribution.

More precisely, the phase space density function of the stars, $F(t, {\vec r}, {\vec V})$, is subject to a continuity condition in the six-dimensional phase space, defined by the three components of the star's position, ${\vec r}$, and the three components of its velocity, $ {\vec V}$, which is expressed by the Liouville equation. It has has linear properties with respect to the density function, so that the superposition principle is satisfied. In order to isolate information about the spatial properties of the stellar system it is useful to integrate the Liouville equation over the velocity space. The resulting equation contains moments of the velocity distribution, such as the mean velocity, or the velocity dispersions. Furthermore the Liouville equation may be multiplied through by any powers of the velocities before integrating, and each choice of powers leads to a different result which involves different velocity moments. These differential equations are the stellar hydrodynamic equations or moment equations (e.g. Gilmore et al. 1989, p. 118). Therefore it is essential to have good estimations of the velocity moments, and in particular to know whether some of them are intrinsically null, in order to solve these equations.

On the other hand, according to Jeans' theorem, the phase space density function is a function of the integrals of motion of the stars, that involve the potential function. Hence the nature of the phase space density function must be explained from the dynamics of large stellar groups, sharing a common potential, through integrals like the energy, the angular momentum, or more general quadratic integrals, rather than from the kinematics of specific groups of stars, such as those producing streaming motions (Kurth 1957). Thus, by taking into account a basic set of integrals of motion, and also for the sake of simplifying the calculations in solving the hydrodynamic equations, the phase space density function may be approximated in two ways. First by assuming some specific functional dependences, and second by considering some plausible symmetry hypotheses of the distribution.

The first approximation concerns the superposition principle. Stellar populations can be identified with galactic components, like thin disk, thick disk, stellar halo, etc. For some of these stellar populations the velocity distribution function, that is $F(t_0, {\vec r_0}, {\vec V})$ for fixed time t0 and position ${\vec r_0}$(although sometimes F is assumed time-independent), can be represented by a trivariate Gaussian function, depending on an ellipsoidal integral of motion, according to Oort's approach, or to the more general Chandrasekhar's approach (Chandrasekhar 1942). This fact can be interpreted, according to the theory of gas dynamics, that the population has reached the statistical equilibrium. Conversely, a Gaussian distribution can be used to define a pure statistical population. In this situation all the odd-order central moments vanish, and higher-order moments are explicitly dependent on the second ones. Of course this produces a dramatic simplification of moment equations, but actual samples do not show such a pure stellar population alone in the solar neighborhood. However a mixture of trivariate normal populations is sufficient to explain the most relevant local kinematic features.

In the present work we shall see that an alternative parameterization of the mixture, based on the symmetry of the velocity distribution around the direction containing both population means, introduces a set of constraint equations between the overall distribution cumulants, that are associated with characteristic constants of the mixture.

The second approximation concerns the geometry of the distribution. For example, the velocity distribution of most disk stellar samples shows a galactic plane of symmetry. Then all the partial and total central moments which involve an odd-number of times the velocity component perpendicular to this plane are null. Similarly, some stellar populations, like thick disk or halo, may present an axially symmetric distribution with a velocity ellipsoid whose major axis always points toward the galactic center (Soubiran et al. 2003). Then their velocity ellipsoid would have no vertex deviation on two orthogonal planes of the velocity space and, for these stellar populations, two of the non-diagonal central moments would be null. Thus, some perturbations of those velocity moments may be explained through the influence of a galactic bar (Mühlbauer & Dehnen 2003). For a normal population alone it is possible to test the consistency of the above or similar hypotheses by seeking the indices of the vanishing central moments (Erickson 1975; Vandervoort 1975), while for a mixture of two normal components it must be done by evaluating the cumulant constraint equations that have been obtained. Moreover under specific distribution symmetries only a subset of non-vanishing moments may be used, and then the constraint equations adopt a simpler form.

Therefore we shall study how the foregoing assumptions are transferred to the moments and cumulants, so that these statistics can serve for testing the consistency of actual samples versus the diverse hypotheses. In addition, the analysis based on the above approximations is converted into a numerical algorithm in order to evaluate the population parameters.

The method is applied to real stellar samples of the solar neighborhood, mainly obtained from two catalogs: The Third Catalog of Nearby Stars (CNS3, Gliese & Jahreiss 1991), that is used to compare the results with early works, and the more recent HIPPARCOS Catalog (ESA 1997), with a subsample up to 300 pc from the sun with known radial velocities, which gives an actual portrait of the local kinematics.

The comparison of samples shows that CNS3 sample slightly overestimates high velocity stars, due to high proper motion sampling. It contains a few stars (3%) belonging to a high velocity component, the thick disk, and a predominant old thin disk. Nevertheless the local kinematics is better described when the larger HIPPARCOS subsample is studied. In addition to verify that a two-component normal mixture is a good approximation of the local velocity distribution, we find that the cumulant constraints are consistent with the following features: (1) A symmetry plane of the distribution, as is commonly assumed. (2) A lag in rotation of thick disk stars behind the thin disk of $51 \pm 3$ km s-1. (3) No significant differential movement in the radial direction. (4) Although CNS3 sample shows a velocity distribution closer to the axisymmetry hypothesis, the more accurate data of HIPPARCOS sample, with 13 678 stars, reveal an incipient deviation from axisymmetry due to the vertex deviation of both, thin and thick disk, velocity ellipsoids.

The paper is structured as follows. In Sect. 2 we introduce the basic notation and statistical definitions. In Sect. 3 we summarize the classical parameterization of the superposition model and the set of equations to be solved, that is, the equations of the total cumulants depending on the mixture parameters. In Sect. 4 the new parameterization is presented, that leads to the cumulant constraint equations. In Sect. 5 the constraint equations are studied by assuming specific distribution symmetries. In Sect. 6 we explain the steps composing the segregation algorithm. The application to local samples is carried out in Sect. 7, where the population parameters are estimated. In Sect. 8 kinematical features of HIPPARCOS sample are discussed. Finally, in Sect. 9, the symmetry hypotheses are tested and we discuss and compare our results.

2 Notation and basic definitions

Let us assume that the kinematic behavior of the stellar system is described from a conservative linear dynamic system through the Liouville equation. Then the superposition principle is satisfied and a composite velocity distribution function may be assumed. Both ideal stellar populations are associated with Gaussian components. Notice that, laying aside astronomical considerations, from a Bayesian criterion the Gaussian distribution is the less informative one with known means and covariances (e.g. Koch 1990), and hence it is the usual and less restrictive approach for the components of a mixture model without any other prior information. With these basic assumptions the notation is hereafter introduced, so that the following concepts and definitions may be applied to a stellar population alone in the present section, as well as to the total mixture, in the following section.

For fixed time t0 and position ${\vec r_0}$ let be $f({\vec V}) = F(t_0, {\vec r_0}, {\vec V})$ the velocity density function. The mean velocity, or velocity of the centroid, is noted as ${\vec v}$. The n-rank symmetric tensor of the n-order central moments is defined according to the following expected value, depending on the peculiar velocity ${\vec u}={\vec V}-{\vec v}$, that is the velocity referred to the centroid,

 \begin{displaymath}{\vec M}_{n} = E\left[({\vec u})^{n}\right] = \int_{\vec V}({\vec V}-{\vec v})^{n} ~ f({\vec V}) ~ {\rm d}{\vec V} \end{displaymath} (1)

where $(\cdot)^{n}$ denotes the n-tensor power. The symmetric tensor  ${\vec M}_{n} $ of the n-order central moments has $\left( \begin{array}{c} n+2 \\ 2 \end{array} \right)$ different elements according to the following expression,

 \begin{displaymath}\mu_{\alpha_{1}\alpha_{2}\ldots\alpha_{n}} = E\left[ u_{\alpha_{1}}u_{\alpha_{2}}\ldots u_{\alpha_{n}}\right] \end{displaymath} (2)

where the indices, depending on the velocity components, belong to the set $\{1,2,3\}$.

The velocity density function of an ideal stellar population, of normal type, may be written according to the expression

 \begin{displaymath}f({\vec V}) \equiv \psi(Q) \propto {\rm e}^{- \frac {1}{2}Q} \end{displaymath} (3)

being $Q = {\vec u} ^{T} \cdot {\vec A}_{2} \cdot {\vec u}$ a semi-positive definite quadratic form, and ${\vec A}_{2} $a second-rank symmetric tensor. Note that $\psi$ is here representing the Gaussian function. In this case the relationship between ${\vec A}_{2} $ and the second central moments ${\vec M}_{2}$, i.e. the covariance matrix, is ${\vec M}_{2} = {\vec A}_{2}^{-1}$, and the quadric equation Q=1 defines the velocity ellipsoid. In general the velocity ellipsoid is not oriented according to the coordinate axes. A rotation of angle $\epsilon_{ij}$ in the plane $u_{i}u_{j}, (i \ne j)$ around their orthogonal k-axis is measured from the following relationship involving the second moments

 \begin{displaymath}\tan 2\epsilon_{ij}=\frac{2\mu_{ij}}{\mu_{ii}-\mu_{jj}}\cdot \end{displaymath} (4)

When this rotation is referred to the major semi-axis it is called vertex deviation. Furthermore, due to the symmetry properties of Gaussian distributions, all the odd-order central moments are zero and the even-order central moments may be computed from the second ones. The fourth central moments satisfy, by components, the following relationship,

 \begin{displaymath}\mu_{ijkl} = A^{-1}_{ij} A^{-1}_{kl}+ A^{-1}_{ik} A^{-1}_{jl}+ A^{-1}_{il} A^{-1}_{jk}; \; i,j,k \in \{1,2,3\} \end{displaymath} (5)

that can also be expressed in terms of the second moments according to the equality

 \begin{displaymath}\mu_{ijkl} = \mu_{ij} \mu_{kl}+ \mu_{ik} \mu_{jl}+ \mu_{il} \mu_{jk}; \; i,j,k \in \{1,2,3\}. \end{displaymath} (6)

These set of equations are explicitly written in the Appendix A, and they may be used in order to identify a trivariate normal distribution alone.

The foregoing relationships can be expressed in a more compact notation (Cubarsi 1992) in order to simplify the algebraic notation of the following sections. If ${\vec A}_{m}$ and ${\vec B}_{n}$ are two m- and n-rank symmetric tensors, we define the tensor ${\vec A}_{m}\star {\vec B}_{n}$ as the obtained by symmetrizing the tensor product ${\vec A}_{m} \otimes {\vec B}_{n}$, and by normalizing with respect to the number of summation terms. The result is a (m+n)-rank symmetric tensor, whose components are

 \begin{displaymath}{\vec A}_{m} \star {\vec B}_{n}\vert _{i_{1}i_{2} \ldots i_{m... ...\ldots \alpha i_{m}} B_{\alpha i_{m+1} \ldots \alpha i_{m+n}} \end{displaymath} (7)

where $\alpha$ belongs to the symmetric group ${\cal S} (m + n)$. Then Eq. (6) simply becomes

 \begin{displaymath}{\vec M}_{4} = 3 \; {\vec M}_{2} \star {\vec M}_{2}. \end{displaymath} (8)

On the other hand the cumulants, namely ${\vec K}_{n}$, (their tensor elements will be noted with the greek letter $\kappa$) may also be used for describing the velocity distribution. Originally introduced by Fisher in 1928, they present several attractive properties coming from their symmetry. Up to fourth order, and with the notation of Eq. (7), the cumulants can be written from the corresponding central moments according to the following relationships (e.g. Stuart & Ord 1987). For first, second and third order, ${\vec K}_{1}={\vec v}, {\vec K}_{2}={\vec M}_{2}$, and ${\vec K}_{3}={\vec M}_{3}$. The fourth cumulants satisfy:

 \begin{displaymath}\begin{array}{l} {\vec K}_{4} = {\vec M}_{4} - 3 \; {\vec M}_{2}\star {\vec M}_{2}. \end{array}\end{displaymath} (9)

In particular, for a multivariate normal distribution, $ {\vec K}_{n} = {\vec 0}$ is fulfilled, with $ n \geq 3$ (obviously the case n=4 is a consequence of Eq. (8)). Hence for a Gaussian population the only non-vanishing cumulants are the second ones, which, in addition to the mean, completely characterize the distribution.

It is known that the sample moments, namely ${\vec m}_{n}$, are biased estimators of the population moments ${\vec M}_{n} $, conversely, under the assumption of homogeneous observational errors, the population cumulants have as unbiased estimators the corresponding k-statistics, namely ${\vec k}_{n}$, or sample cumulants. Let us remark that whereas the sample moments ${\vec m}_{n}$ are the same function of the sample values as the population moments ${\vec M}_{n} $, the same relation does not hold for ${\vec k}_{n}$ and ${\vec K}_{n}$. Basically the k-statistics are sums of products of the sample moments, which - like the cumulants and the central moments - are invariant under change of the origin, except for the first order. The tensor forms for the k-statistics of a multivariate distribution were published by Kaplan (1952). The first k-statistic is equal to the mean, and up to fourth-order they may be obtained depending on the sample moments according to

 \begin{displaymath}\begin{array}{l} \displaystyle {\vec k}_{2}=\frac{N}{N-1}{\ve... ...(N+1)} ~ 3 ~ {\vec m}_{2}\star {\vec m}_{2} \right] \end{array}\end{displaymath} (10)

where N is the size of the sample.

3 Cumulants of the mixture

The overall density function is now obtained from the superposition of two normal density functions according to Eq. (3), each one associated with the corresponding stellar component, (') or ('') for the first or second population,

 \begin{displaymath}f({\vec V}) = n'\psi(Q') + n''\psi(Q'') \end{displaymath} (11)

where n' and n'' represents the respective mixing proportions, obviously satisfying n'+ n''= 1. The quantities ${\vec v}$ and ${\vec M}_{n} $ for the total velocity distribution - defined in the above section - may easily be deduced starting from those of the partial ones. Hence the total mean velocity, expressed from the population means, satisfies $ {\vec v} = n'{\vec v}'+ n''{\vec v}''$.

Let us review a way of computing the mixture moments and cumulants. The total central moments are written, taking into account Eqs. (1) and (11), by using the centroid differential velocity,

 \begin{displaymath}{\vec w} = {\vec v}'- {\vec v}''. \end{displaymath} (12)

By using the product defined in Eq. (7), the tensor of the total n-order central moments, expressed from the partial ones, has the following form:

 \begin{displaymath}{\vec M}_{n} = \sum_{k=0}^{n} \left( \begin{array}{c}n \\ k... ...ight)^{k} n'' {\vec M}''_{n-k} \right\} \star ({\vec w})^{k}. \end{displaymath} (13)

We explicitly write the total central moments up to fourth-order. Obviously for n=0 we get ${\vec M}_{0}={\vec M}'_{0}={\vec M}''_{0}=1$, and for n=1, ${\vec M}_{1}={\vec M}'_{1}={\vec M}''_{1}= {\vec 0}$. The symmetric tensor  ${\vec M}_{2}$ of the second-order moments, with six different elements, becomes

 \begin{displaymath}{\vec M}_{2}= n'{\vec M}'_{2}+ n''{\vec M}''_{2}+ n'n''({\vec w})^{2}. \end{displaymath} (14)

The tensor ${\vec M}_{3}$ of the third moments, with ten different elements, being ${\vec M'}_{3}={\vec M''}_{3}={\vec 0}$, is written as

 \begin{displaymath}{\vec M}_{3}= 3n'n''({\vec M}'_{2}-{\vec M}''_{2}) \star {\vec w} + n'n''(n'' - n')({\vec w})^{3} . \end{displaymath} (15)

And the fourth moments ${\vec M}_{4}$, with fifteen different elements, can be expressed, by using Eq. (8), as follows
 
                      $\displaystyle {\vec M}_{4}$ = $\displaystyle 3n'{\vec M}'_{2} \star {\vec M}'_{2} + 3n''{\vec M}''_{2} \star {\vec M}''_{2}$  
    $\displaystyle + 6n'n'' (n''{\vec M}'_{2} + n'{\vec M}''_{2}) \star ({\vec w})^{2}$  
    $\displaystyle + n'n''(1-3n'n'') ({\vec w})^{4}.$ (16)

However, if the population cumulants are used (Eq. (9)), it is possible to write the above relationships in a shorter form. This is done by introducing the following new variables,

 \begin{displaymath}{\vec D} = \sqrt{n'n''} ~ {\vec w} ; \qquad q = \sqrt{n'/n''} - \sqrt{n''/n'} \end{displaymath} (17)

(we can appoint the populations so that $n'\geq n''$, then q is non-negative) and the following second-rank tensors,

 \begin{displaymath}{\vec a}_{2} = n'{\vec M}'_{2} + n''{\vec M}''_{2}; ~~ {\vec ... ...{2}+ 4}} ({\vec M}'_{2} - {\vec M}''_{2}) - q ({\vec D})^{2}\ \end{displaymath} (18)

(all the above square roots are taken as positive values). With the foregoing definitions, Eqs. (14)-(16) can be rewritten explicitly depending on the overall mixture cumulants in a shorter form, as follows:

 \begin{displaymath}\begin{array}{l} {\vec K}_{2} = {\vec a}_{2} + ({\vec D})^{2... ...2} \star {\vec C}_{2} - 2(q^{2}+1) ({\vec D})^{4}. \end{array}\end{displaymath} (19)

From the total cumulants of the mixture the parameters of the partial distributions have to be determined by inverting the above non-linear system of equations. These unknowns are the partial moments - six components for ${\vec M}'_{2}$ and six for ${\vec M}''_{2}$ -, the population fraction n', and the three components of the centroid differential velocity ${\vec w}$. Sixteen unknowns in total. Since we have a set of thirty-one scalar equations, involved in Eq. (19), we must also find a set of fifteen constraint equations, which can be used as a test in order to verify whether a given sample is consistent with a two normal mixture.

In the following section we describe how the total cumulants of the mixture are related. A set of equations that generalize Eq. (6) is obtained, which provide some characteristic mixture constants, such as a vector ${\vec d}$ in the direction along both sub-centroids, and two constants, ${\cal A}$ and ${\cal B}$, that can be linearly estimated from total cumulants, with useful information about the geometry of the mixture.

4 Constraint equations

We study the general case where the difference between population means, ${\vec w}$, (and hence the vector ${\vec D}$ of Eq. (17)) is not null. Let us assume the vector component $D_{3} \neq 0$ (in order to minimize the error propagation this component may be chosen to be $\max_{i}\vert D_{i}\vert$), and let us define a normalized vector ${\vec d} = {\vec D}/D_{3}$ in the direction containing both sub-centroids C1 and C2 (Fig. 1). Since every normal distribution is symmetric with respect to its centroid, then the total velocity distribution will be symmetric in whatever direction within a plane $\Phi$ orthogonal to the vector ${\vec d}$, and in particular the one containing the global centroid Ct. Thus, in order to take profit of this symmetry, it is convenient to work with a transformed vector ${\vec U}$ instead of the peculiar velocity ${\vec u}$, whose components are three non-orthogonal projections of the peculiar velocity ${\vec u}$ on the directions ${\vec W_1}=(d_{3},0,-d_{1})^{t}$ and ${\vec W_2}=(0, d_{3},-d_{2})^{t}$, on the plane $\Phi$, and another independent direction, for example ${\vec W_3}=(0, 0, d_{3})^{t}$.


  \begin{figure} \par\includegraphics[width=8.8cm]{1144fig1.ps}\end{figure} Figure 1: Directions ${\vec W_1}, {\vec W_2}, {\vec W_3}$ where the peculiar velocity ${\vec u}$ is projected.
Open with DEXTER

The transformed peculiar velocity ${\vec U}$ can be expressed from the following isomorphic transformation of the vector ${\vec u}$,

 \begin{displaymath}{\vec U} = {\vec H}_{2} \cdot {\vec u} ; \; {\vec H}_{2} = \... ... 0 & d_{3} & -d_{2} \\ 0 & 0 & d_{3} \end{array} \right)\cdot \end{displaymath} (20)

Note that $\det({\vec H}_{2})=d_{3}^3$, and d3=1, but in case of permuting indices it is helpful to maintain this notation. Also note that ${\vec U}$, as well as ${\vec u}$, have null mean. In fact U3=u3 is an invariant component. From the definition of ${\vec H}_{2}$ the following equality is deduced,

 \begin{displaymath}\sum_{j} H_{ij}D_{j}=D_{3}\delta_{3i};\; i, j \in \{1,2,3\} \end{displaymath} (21)

where $\delta$ is the Kronecker delta.

If the third and fourth moments of ${\vec U}$ are calculated in function of the central velocity moments $\mu_{ijk}$ and $\mu_{ijkl}$, the following equalities are obtained,

 
    $\displaystyle E[ U_{\alpha} U_{\beta} U_{\gamma}] = \sum_{ijk} {H}_{\alpha i} {H}_{\beta j} {H}_{\gamma k} \mu _{ijk};$  
    $\displaystyle ~~~~~~~~~~~~~~~~~~~~~~~~~~\alpha, \beta, \gamma, \delta, i, j, k, l \in \{1,2,3\}$  
    $\displaystyle E[ U_{\alpha} U_{\beta} U_{\gamma} U_{\delta}] = \sum_{ijk} {H}_{\alpha i} {H}_{\beta j} {H}_{\gamma k} {H}_{\delta l} \mu_{ijkl}.$ (22)

And also the the third and fourth cumulants of the transformed peculiar velocity ${\vec U}$ can be computed in function of the corresponding cumulants of ${\vec u}$. With the indices $\alpha,\beta,\gamma,\delta \in \{1,2\}$ and i,j, $k,l \in \{1,2,3\}$, the ten components of the third ${\vec U}$-cumulants are,

 \begin{displaymath}\begin{array}{l} o_{\alpha \beta \gamma} \equiv \sum_{ijk} {... ...um_{i} {H}_{\alpha i} \kappa_{i33}\\ \kappa_{333} \end{array}\end{displaymath} (23)

and the fifteen components of the fourth ${\vec U}$-cumulants are,

 \begin{displaymath}\begin{array}{l} X_{\alpha \beta \gamma \delta} \equiv \sum_... ...H}_{\alpha i} \kappa_{i333} \\ [1mm] \kappa_{3333}. \end{array}\end{displaymath} (24)

Thus the ${\vec U}$-cumulants can be grouped according to the above two-dimensional tensors ${\vec o}_{3}$, ${\vec p}_{2}$, ${\vec s}$, ${\vec X}_{4}$, ${\vec Y}_{3}$, ${\vec Z}_{2}$, and ${\vec T}$, depending on the total cumulants of the distribution. These tensor components are explicitly written in the Appendices B and C.

Let us remark that all the above quantities are explicitly dependent on the velocity component that remains invariant under the transformation of Eq. (20). Hence they would have to be noted, for example, with a super-index (3) indicating that component, since the described procedure is also valid under permutation of indices of the velocity components. However, in order to simplify the notation, this super-index has been here omitted, although it appears in the Appendices.

Hereafter the main properties and the steps in order to obtain the cumulant constraints, the mixture constants, and the population parameters are summarized.

By substitution in Eq. (23) of ${\vec K}_{3}$, from Eq. (19), and taking into account Eq. (21), four independent vanishing linear combinations of the third ${\vec U}$-moments are obtained:

 \begin{displaymath}o_{\alpha \beta \gamma} = 0; \; \alpha,\beta,\gamma \in \{1,2\}. \end{displaymath} (25)

These are the third cumulants of the vector components U1 and U2, that vanish as consequence of having defined the vectors  ${\vec W_1}$ and ${\vec W_2}$ on the plane $\Phi$, with symmetrical properties of the distribution. Similar equalities are satisfied for higher odd-order cumulants. From the above equations (see Appendix B) the values d1 and d2 can be computed. The solution of this system of nonlinear equations provides us with the first characteristic constant of the mixture, ${\vec d}$, that is the orientation of both sub-centroids. Now it is already possible to calculate the elements of the matrix ${\vec H}_{2}$, in Eq. (20), and the tensor elements of Eqs. (23) and (24).

Let us remark that, while the method of moments for an univariate normal mixture requires to solve a fundamental nonic equation (Cohen 1967), originally derived by Pearson, for a trivariate mixture only three-degree polynomials have to be solved, and the moments method loses, or substantially reduces, the consideration of ill-conditioned problem.

For the elements of tensors ${\vec p}_{2}$ and ${\vec s}$ in Eq. (23), by substitution of the third cumulants ${\vec K}_{3}$ from Eq. (19), and taking into account Eq. (21), the following equivalences are obtained:

 \begin{displaymath}\begin{array}{l} p_{\alpha \beta} = D_{3} \sum_{ij} {H}_{\al... ...s_{\alpha} = D_{3} \sum_{i} {H}_{\alpha i} C_{i3} . \end{array}\end{displaymath} (26)

These relationships are used to compute the elements of ${\vec C}_{2}$ other than C33.

Similarly, for the elements of tensors ${\vec X}_{4}$, ${\vec Y}_{3}$, ${\vec Z}_{2}$, and ${\vec T}$ in Eq. (24), by substitution of the fourth cumulants ${\vec K}_{4}$ from Eq. (19), and taking into account Eq. (26), the following set of constraint equations is obtained,

 \begin{displaymath}\begin{array}{l} {\vec X}_{4}=3 ~ {\cal A} ~ {\vec p}_{2} \s... ... B} ~ {\vec p}_{2}\\ {\vec T}=3 {\cal B}~ {\vec s} \end{array}\end{displaymath} (27)

where ${\cal A} = D_{3}^{-2}$ and ${\cal B} = C_{33}D_{3}^{-1}$. These fourteen scalar relationships are explicitly written in the Appendix D.

The set of relationships in Eq. (27) represents an overdeterminate linear system, which can be solved by means of weighted least squares in order to find optimal values for the mixture constants ${\cal A}$ and ${\cal B}$. Note that this step provides the absolute values of C33 and D3.

The mixing proportions are evaluated from the parameter q, so that the following two relationships are fulfilled,

 \begin{displaymath}\begin{array}{l} \kappa_{333} = 3C_{33} D_{3} + 2 q D_{3}^3 ... ...ppa_{3333} = 3C_{33}^{2} - 2(q^{2} + 1) D_{3}^{4}. \end{array}\end{displaymath} (28)

Hence, by elimination of q, a new constraint equation is hold.

The remaining five unknowns of the tensor ${\vec C}_{2}$ may be evaluated from Eq. (26), and finally, from Eqs. (17) and (18), the population parameters n', ${\vec w}$, ${\vec M}'_{2}$, and ${\vec M}''_{2}$, can be determined.

5 Distribution symmetries

The constraint equations of the trivariate normal mixture can be used in order to test underlying distribution symmetries by working only from the total k-statistics of the sample. This is a usual procedure for a normal population alone, where this information can be deduced, for example, from the orientation of the velocity ellipsoid. Now this idea is generalized for the case of a two population mixture.

In the first section we have explained that it is useful to characterize the geometry of the stellar velocity distribution in order to assume specific symmetries that simplify the set of distribution moments. We shall pay attention to the following assumptions: (a) existence of a distribution symmetry plane; (b) non differential movement in any specific direction, that is, the same mean for a given velocity component; and (c) velocity distribution with axial symmetry. These particular situations are explicitly analyzed below, where the index 1 will apply to the radial velocity, the index 2 to the rotation velocity - which is now taken as the invariant component in the transformation given by Eq. (20) -, and the index 3 to the velocity perpendicular to the galactic plane.

a)
Symmetry plane. Let us assume that the velocity component V3 has a symmetrical distribution around their mean for both population components, that is v'3=v''3=0. This is actually satisfied for the velocity component perpendicular to the galactic plane of most disk stellar samples. Hence there is a symmetry plane u3=0 in the space of the peculiar velocities, so that for each population, as well as for the total distribution, the probability density function has the form f=f(u1, u2, u32). Thus all the partial and total central moments and cumulants involving an odd-number of times the index 3 are null on this plane. This is similarly satisfied by the elements of ${\vec D}, {\vec a_2}$, and ${\vec C}_{2}$, defined in Eqs. (17) and (18). In particular it implies that each population component has an independent distribution of the variate u3. Then, from Eq. (25), the following relationship is deduced,

 \begin{displaymath}\frac{D_1}{D_2}=\frac{\kappa_{133}}{\kappa_{233}}\cdot \end{displaymath} (29)

Moreover in Eqs. (26) and (27) there is a set of expressions that does not explicitly depend on the elements of vector ${\vec d}$, since they can be computed directly from the mixture cumulants, with a minimum error estimation, as follows,

 \begin{displaymath}\begin{array}{l} \displaystyle \frac{1}{C_{33}^2}=\frac{3}{\k... ...ac{D_1}{\kappa_{133}}=\frac{D_2}{\kappa_{233}}\cdot \end{array}\end{displaymath} (30)

b)
Non differential movement. A more general situation occurs when the population components have the same mean in one direction, for example v'1=v''1. Older stellar populations usually satisfies this condition for the radial velocity component (e.g., Chiba & Beers 2000). Notice that (a) is a particular case of (b). Hence D1=d1=0 and, by Eq. (25), $\kappa_{111}$ must vanish. Similarly to the previous case, d1=D1/D2 is computed from Eq. (29).

If the ratios

 \begin{displaymath}\alpha_2 = \frac{\kappa_{122}}{2\kappa_{112}}, \quad \alpha_3 = \frac{\kappa_{133}}{2\kappa_{113}} \end{displaymath} (31)

are defined, the following relationships can be computed directly from the total cumulants, without depending explicitly on d1,
$\displaystyle % \begin{array}{l} \displaystyle \frac{1}{C_{11}^2}=\frac{3}{\kap... ...frac{3 \alpha_2}{\kappa_{1112}}=\frac{3 \alpha_3}{\kappa_{1113}}\\ \end{array}$      
 
                                   $\displaystyle \frac{1}{C_{11}}$ = $\displaystyle \frac{C_{22}}{\kappa_{1122}-\frac{2}{3}\alpha_2^2 \kappa_{1111}}= \frac{C_{33}}{{\kappa_{1133}-\frac{2}{3}\alpha_3^2 \kappa_{1111}}}$  
  = $\displaystyle \frac{\alpha_2}{C_{12}}= \frac{\alpha_3}{C_{13}}=\frac{D_2}{\kapp... ...}}= \frac{D_3}{\kappa_{113}}=\frac{3\alpha_2 \kappa_{112}C_{22}}{\kappa_{1222}}$  
  = $\displaystyle \frac{3\alpha_3 \kappa_{113}C_{33}}{\kappa_{1333}}\cdot$ (32)

The foregoing equations represents a particular case of Eqs. (26) and (27) under the assumption of non differential movement in the direction of the velocity component V1.

c)
Axial symmetry. In stellar dynamics one of the more usual assumptions is the axial symmetry hypothesis of mass and velocity distributions. Then one axis of the velocity ellipsoid points toward the axis of rotation of the galaxy, and another axis remains in the direction of rotation (e.g. Gilmore et al. 1989, p.132). Thus each population velocity ellipsoid has no vertex deviation on two planes of the peculiar variates, that is u1u2 and u2u3. Hence the population components satisfy $ \kappa'_{12}=\kappa''_{12}= \kappa'_{23}=\kappa''_{23}=0$. Therefore, for each population component, but not for the total velocity distribution, the only variates that are correlated are u1 and u3, while u2 has an independent distribution from the others. From Eqs. (17)-(19) it is easy to deduce that the difference between population means is the only contribution to the total non-diagonal second moments,

\begin{displaymath}\kappa_{12}=D_1 D_2; \; \kappa_{23}=D_2 D_3. \end{displaymath} (33)

Two of the non-diagonal elements of the tensor ${\vec C}_{2}$ are proportional to second cumulants, $C_{12}=-q \kappa_{12}$ and $ C_{23}=-q \kappa_{23}$. Then by defining the quantities

 \begin{displaymath}\begin{array}{l} \kappa^{*}_{1122}=\kappa_{1122}+2 \kappa^2_{... ...mm] \kappa^{*}_{2233}=\kappa_{2233}+2 \kappa^2_{23} \end{array}\end{displaymath} (34)

the following relationships, involving the total cumulants, can be obtained in a similar way as in previous cases,

 \begin{displaymath}% \begin{array}{l} \displaystyle \frac{1}{C_{22}^2}=\frac{1}{... ...{1223}}=\frac{\kappa_{233}}{\kappa^{*}_{2233}}\cdot \end{array}\end{displaymath} (35)

Notice that the relationships obtained for the different symmetry hypotheses introduce a worthy simplification of the equations in Appendices B, C, and D, so that sometimes the involved parameters can be evaluated nearly by hand.

Three combinations of above hypotheses are particularly interesting: (a+b) for samples in the galactic plane with differential rotation alone, (a+c) for samples with axial and galactic plane symmetries, and (a+b+c) for the complete set of hypotheses.

a+b)
In this case it is easy to see that $\alpha_3=0$ in Eq. (31). Hence the relationships in Eq. (32) are hold only if they do not contain D3 and $\alpha_3$. The vanishing third and fourth cumulants are, like in the case of symmetry plane, those having an odd-number of times the index 3. In addition, the following third moments are also null:

 \begin{displaymath}\kappa_{111}=\kappa_{133}=0. \end{displaymath} (36)

Therefore, in general, there are four third cumulants and nine fourth cumulants which are non-null.

a+c)
In order to estimate the ratio between differential velocity components, Eq. (29) may also be used. The combination of both hypotheses leads to the following subset of equations:

 \begin{displaymath}\begin{array}{l} \displaystyle\frac{1}{C_{22}^2}=\frac{1}{3} ... ...*}_{1122}}= \frac{\kappa_{233}}{\kappa_{2233}}\cdot \end{array}\end{displaymath} (37)

Obviously, for the overall sample, the only non-vanishing non-diagonal second moment is $\kappa_{12}=D_1 D_2$, and the cumulants with an odd-number of times the index 3 are also zero.

a+b+c)
This is a quite simple case, where all the non-diagonal second moments are null, $\kappa_{12}=\kappa_{13}=\kappa_{23}=0$. In addition to the six third moments of case (a+b), which vanish, the moment $\kappa_{122}$ is also null. Hence the only non-vanishing third moments are $\kappa_{112}$, $\kappa_{222}$, and $\kappa_{233}$. In general there are only six non-null fourth cumulants: $\kappa_{1111}$, $\kappa_{2222}$, $\kappa_{3333}$, $\kappa_{1122}$, $\kappa_{1133}$, and $\kappa_{2233}$.
When applying the method to real stellar samples, we shall test whether the complete set of constraint equations can be reduced to any of the foregoing symmetry cases.

6 Numerical algorithm

The method is converted into a numerical procedure (Alcobé 2001) paying special attention to the error analysis. In the earlier work (Cubarsi 1992) a partial application of constraint equations for axisymmetric stellar systems was carried out, using interval arithmetic (Moore 1966). It leads however to low precision results. Now the general set of cumulant constraints is taken into account, since actual velocity samples, in particular the HIPPARCOS one, have more accurate data. In each computational step, statistical propagation of errors and weighted least squares estimation have been adopted in order to get minimum variance estimates. In addition, the goodness of the approximation is evaluated by means of a $\chi ^2$test, that is commented at the end of the section. The main steps which compose the complete numerical procedure are hereafter described:

i)
Computation of the sample moments and evaluation of the k-statistics in order to estimate the mixture cumulants. The standard errors are evaluated through higher-order sampling cumulants. The algorithm allows to neglect those values which are smaller than a certain number of times their standard error.

Since the asymmetry of the distribution around the mean is provided by the third moments (by the odd-order central moments in general) it is convenient to take the velocity component that remains invariant under the transformation of Eq. (20) as the one having more accurate non-vanishing third moment. In general this is the rotation velocity.

ii)
Calculation of the orientation vector ${\vec d}$ by using the equations in the Appendix B. The two first equations, for o111 and o222, are solved separately as one variable cubic equation. A first couple of values, d1 and d2, is found, which are taken as initial values for a Newton-Raphson iterative process in order to solve the complete system of four linear equations with two unknowns by means of an iterative least squares estimation.

iii)
The elements of tensors ${\vec p}_2,{\vec s},{\vec X}_4,{\vec Y}_3, {\vec Z}_2$, and ${\vec T}$ of Appendix C are computed, as well as the respective standard errors. The relationships between these quantities are shown in Eq. (27) and Appendix D, where a linear least squares system with 14 equations and 2 unknowns is solved. The equations are weighted according to the inverse of the error covariance matrix. The solutions for ${\cal A}$ and ${\cal B}$, if they are computed with a reasonable certainty, are assumed to be constants of the mixture, otherwise the superposition hypothesis should be rejected. The vector ${\vec D}$ and C33 are also computed.

iv)
The remaining elements of ${\vec C}_{2}$ are also computed from Eq. (26). Its components are calculated with their definite sign. From the $\chi ^2$ test it is found that the formulas involving the third cumulants lead to better results than those involving the fourth ones, since, in general, the higher order moments, the higher sampling errors.

v)
The population factor q is computed according to Eq. (28). It may be evaluated from the third or fourth order cumulants. Once again the third cumulants provide a better $\chi ^2$ estimation. Finally the population second moments, ${\vec M}_2^{\prime }$ and  ${\vec M}_2^{\prime \prime }$, and the centroid differential velocity ${\vec w}$are computed by using Eqs. (17) and (18).
The $\chi ^2$ test is used in order to estimate the goodness of the approximation of the total cumulants through the two-normal mixture. The cumulants of the mixture up to fourth order, that we can represent by a vector ${\vec y}$, with R components, are fitted from the corresponding mixture parameters, namely ${\vec x}$, so that a set of relationships ${\vec y}={\vec g}({\vec x})$, according to Eq. (19), is satisfied. Thus, a number of R=31 equations are involved, and an amount of $\upsilon=16$ parameters have to be determined. If ${\vec e}$ denotes the error vector of ${\vec y}$, a calibration of the approximation can be done from the weighted mean of the squared errors,

 \begin{displaymath}\chi^2 = \sum_{i=1}^{31} \frac{1}{e_{i}^{2}} \left\vert y_{i}- g_{i}({\vec x})\right\vert^2. \end{displaymath} (38)

In order to apply the test we assume that the total error has a 31-Gaussian probability distribution, so that the random variate expressed in Eq. (38) is expected to be a $\chi ^2$ distribution with $n=R-\upsilon=15$ degrees of freedom. Notice that if the errors are given not by a Gaussian but by a Poisson distribution, Eq. (38) obeys the $\chi ^2$ distribution in the large limit n anyway. It is known that, if $P(\chi ^2;n)$ denotes its probability distribution, the relevant quantity to make decisions about goodness of the fit is the $\chi ^2$ probability given by the integral

 \begin{displaymath}\pi(\chi^2;n)=\int_{\chi^2}^{\infty} P\left(\zeta;n\right) {\rm d}\zeta \end{displaymath} (39)

which gives the probability that a function which describes a set of n data points would give a value of $\chi ^2$ as large, or larger, than the one we already have. Nevertheless, when studying actual samples, more than the values $\chi ^2$ and $\pi$, what is more significant is the increase or decrease of these quantities, since they may become distorted due to over or under estimation of observational errors, or to an undesired error distribution.

7 Stellar samples

Two catalogs are used in order to apply the algorithm to local stellar samples, where only the velocity space is taken into account. Hence the resulting population components, distribution symmetries, etc. will reflect strictly kinematic data. Since, in general, the stars with known radial velocity produce kinematically biased samples (Binney et al. 1997), we shall compare our kinematic-based segregation with those obtained from non-kinematically biased samples. Nevertheless Skuljan et al. (1999) show that their samples with radial velocities, similarly to those of Figueras et al. (1997) and Asiain et al. (1999), demonstrate the same basic features than Dehnen (1998), although the latter distribution was obtained without radial velocities. They also conclude that the kinematic bias does not significantly affect the inner parts of the velocity ellipsoid. In the next section we shall analyze in detail the possible kinematic bias of our main sample, by selecting nested subsamples containing higher velocity stars.

The first catalog, Third Catalog of Nearby Stars (CNS3, Gliese & Jahreiss 1991), is used to test our improved method in order to compare the results with previous works. It is composed of all known stars within a distance of 25 pc from the Sun. It was the most statistically complete stellar sample available with known space velocity in the galactic disk (Jahreiss & Gliese 1993), although it overestimates the high velocity stars due to high proper motion sampling. It contains 1946 stars with known velocity space. As subdwarfs are not considered to belong to the Galactic disk (Erickson 1975), six stars among them which are so described have been also rejected. Anyway the obtained results are not significantly different to those derived if such stars are included, but as errors are slightly smaller if we exclude these stars, the criteria has been maintained. The catalog is a natural extension of CNS2, used by Erickson (1975).

Table 1: Means, central moments and cumulants of CNS3 sample selected by $\left\vert {\vec V}\right\vert<183.5$ km s-1, with 1916 stars.

For the sake of working with a homogeneous and a nearly complete sample the catalog is filtered by using the $\chi ^2$ test together with the following condition. Let $\left\vert {\vec V}\right\vert$ be the module of the star velocity, and let v0 be its mean and $\sigma_0 $ its standard deviation for the overall sample. Then, in order to avoid stars with an extreme kinematic behavior the selection of the working sample is made by taking the one having the minimum $\chi ^2$within the subsamples containing stars with $\left\vert {\vec V}\right\vert < v_0 + k \sigma_0 $, for $k \geq 2$.

The sample with the minimum fitting error is the corresponding to $\left\vert {\vec V}\right\vert<183.5$ km s-1 leading to an acceptable Gaussian mixture. Let us remark that the first populations of the subsamples selected by $\left\vert {\vec V}\right\vert< l$, with $l \geq 145$ km s-1, have nearly the same moments, while second components have increasing moments. This clearly indicates that the first population has been included in all the subsamples, while the entering stars are continuously merged to the second component. The final sample is composed of 1916 stars, with moments and cumulants listed in Table 1. For all the actual samples the velocities are expressed in the heliocentric galactic coordinate system, where V1 is the radial velocity toward the galactic center, V2 is the component in the direction of the galactic rotation, and V3 is the component in the direction of the north galactic pole.

The algorithm segregates two main normal populations with kinematic parameters that may be associated with thin disk (97%), noted as Pop-I, and with a very short and extreme thick disk (57 stars), noted as Pop-II (Table 2). In km s-1, the respective velocity dispersions are $(\sigma'_1 ,\sigma'_2 ,\sigma'_3)= (37\pm 1, 24 \pm 2, 19 \pm 1)$ and $(\sigma''_1 ,\sigma''_2 ,\sigma''_3)= (80\pm 7, 52 \pm 23, 52 \pm 6)$. The differential movement is $(w_1 ,w_2 ,w_3)= (2\pm 6, 65\pm 7, -1 \pm 4)$. We get a clearly non-null vertex deviation for the thin disk component, $\mu'_{12}= 145 \pm 33$ $(\epsilon' \sim 9^{\circ})$, while for the thick disk stars we obtain $ \mu''_{12}= 895\pm 423$ $(\epsilon'' \sim 13^{\circ})$, which is non-null within a $2\sigma$ confidence level.

Table 2: Means, central moments and population fraction for CNS3 subsamples, selected by $\left\vert {\vec V}\right\vert<183.5$ km s-1.

We get more accurate results than in Cubarsi (1990, 1992), where the sample of Erickson (1975) was used. Now the vertex deviation of both disk components is asserted. Dispersions are consistent with those published by Sandage (1987), Wyse & Gilmore (1995), and Bienaymé & Séchaud (1997) for similar solar samples.

HIPPARCOS Catalog: This more recent catalog is composed of a large number of stars with known velocity space. The stellar sample that we are using has been obtained by Figueras (2000), basically by crossing HIPPARCOS Catalog (ESA 1997) with radial velocities coming from Hipparcos Input Catalog HINCA (ESA 1992). We assume that, similarly as in Figueras et al. (1997) and Asiain et al. (1999), it is not significantly biased, as we shall see in the next section. In order to work with a representative local sample, the overall sample has been limited up to a distance of 300 pc, since up to this distance the computation of moments is very stable. This distance corresponds to a sphere inside the local thin disk component (e.g. Majewski 1993; Ojha et al. 1999; Chiba & Beers 2000), so that the fractions of thin and thick disk might be representative of the solar neighborhood. The resulting sample is composed of 13 678 stars. Similarly to the CNS3 sample, in order to avoid some non-representative high-velocity stars we select a subsample with minimum fitting error, and $\left\vert {\vec V}\right\vert < v_0 + k \sigma_0 $, for $k \geq 2$, which now corresponds to $\left\vert {\vec V}\right\vert < 210$ km s-1, with 13 531 stars, leading a nearly perfect normal mixture. Also for this catalog, the first segregated populations of the subsamples selected by $\left\vert {\vec V}\right\vert< l$, with $l \geq 145$ km s-1, have nearly the same moments, as it is shown in Table 5 of the following section, while second populations have increasing moments, so that the entering stars are continuously merged to it. The moments and cumulants are displayed in Table 3.

Table 3: Means, central moments and cumulants of HIPPARCOS sample selected by $\left\vert {\vec V}\right\vert < 210$ km s-1, with 13 531 stars.

Table 4: Means, central moments and population fraction for HIPPARCOS subsamples, selected by $\left\vert {\vec V}\right\vert < 210$ km s-1.

The algorithm segregates two main populations, with kinematic parameters that can be associated with the thin disk (91%), noted as Pop-I, and thick disk (9%), noted as Pop-II (Table 4). The respective velocity dispersions are $(\sigma'_1 ,\sigma'_2 ,\sigma'_3)$ = ($28\pm 1$, $16 \pm 2$, $12 \pm 1$) and $(\sigma''_1 ,\sigma''_2 ,\sigma''_3)= (66\pm 2, 40 \pm 9, 41 \pm 2)$. The differential movement is $(w_1 ,w_2 ,w_3)= ( 3\pm 2, 51\pm 3, 1\pm 4)$. For this catalog we get a clearly non-null vertex deviation for both populations, $\mu'_{12}= 93 \pm 15$ $(\epsilon' \sim 9^{\circ})$, and $ \mu''_{12}=314 \pm 112$ ( $\epsilon'' \sim 7^{\circ}$).

Therefore we find that our local HIPPARCOS subsample is composed of two main Gaussian populations, according to thin and thick disk, both with non-vanishing vertex deviation in the galactic plane.

8 Features of HIPPARCOS sample

It is clear (Binney et al. 1997) that the solar velocity obtained from the local standard of rest of solar samples is very sensible to star selection, specially if these stars have high-proper motions, as in the case of stars with known radial velocities. The mean velocity of the sample, in particular the rotation component, significantly varies if the sample contains young disk stars, old disk stars, or thick disk stars. In general the higher velocity stars included in the global sample, the greater biased estimation of the solar motion. However this not implies a biased velocity distribution of population components, since the solar motion might be only deduced from a local stellar group, which the sun belongs to, and probably not from the total thin disk.

In order to study how the high velocity stars modifies the velocity distribution, or the mixture parameters, we select several subsamples containing the whole thin disk component, and also an increasing and faster part of the thick disk. Following a similar criterion than in the previous section, the thin disk component (noted as Pop I in Table 4), which now is well fitted by a normal distribution, can be approximately drawn from the overall sample by taking into account the module of its mean velocity, v, and its total standard deviation, $\sigma$, so that $\left\vert {\vec V }\right\vert= v + 2 \sigma$. This value corresponds to 125 km s-1, and we select subsamples from this value on. Of course all the subsamples will contain an increasing amount of high velocity stars, perhaps a tail of the thick disk, where the second component is not exactly normally distributed. However in this range of velocities the method provides a quite good approximation of mixture components.

Table 5: Components of HIPPARCOS sample selected by star velocity lesser than $\left\vert {\vec V}\right\vert$.

Some segregations for values of $\left\vert {\vec V}\right\vert$ that indicate some kind of discontinuity in the merging process are shown in Table 5. In a future work we shall study in more detail how to detect this kind of discontinuities. We can see that the thin disk component is well isolated in cases of samples selected higher than 145 km s-1, so that the high velocity stars do not contaminate the characteristic thin disk parameters, while these stars are added to a second population. Thus, the second component of the subsample labeled as 125 km s-1 has dispersion parameters $(\sigma''_1 ,\sigma''_2 ,\sigma''_3)= (51 \pm 2, 29 \pm 5, 29 \pm 1)$, and differential centroid velocity $w_2=35\pm2$ km s-1, similar to those of old disk populations, like in Sandage (1987), Wyse & Gilmore (1995), or Bienaymé & Séchaud (1997). The second component of samples from 145 km s-1 to 190 km s-1 have increasing dispersions, although the lag in rotation between populations remains constant at $w_2 \approx 44\pm3$ km s-1. Finally the sample of 210 km s-1, obtained in the previous section accordingly to the minimum fitting error, has already merged all the thick disk population, and has increased the rotational lag and the dispersions.


  \begin{figure} \par\includegraphics[angle=-90,width=8.5cm,clip]{1144fig2.ps}\end{figure} Figure 2: Dependence of the mean heliocentric rotational velocity v2 (km s-1) on the total dispersion S2 (km2 s-2) for the stellar components in HIPPARCOS subsamples. The dashed line corresponds to the thin disk ( $\Delta =89\pm 15$ km s-1) and the continuous line to the thick disk ( $\Delta =116\pm 22$ km s-1).
Open with DEXTER

In order to compare the kinematic features of the foregoing stellar components with others of non-kinematically biased samples we shall also analyze their asymmetric drift. Similarly to Dehnen & Binney (1998), by using Strömberg's quadratic relation for the mean heliocentric rotational velocity, v2, and the total dispersion, S2, of each stellar component, we can fit the linear relationship

 \begin{displaymath}-v_2=\frac{S^2}{\Delta}\cdot \end{displaymath} (40)

The plot we obtain in Fig. 2 is completely consistent with two differentiated stellar galactic components. The dashed line corresponds to the asymmetric drift relation within the thin disk. We obtain a constant $\Delta =89\pm 15$ km s-1 in total agreement with the calculated by Dehnen & Binney (1998). Moreover the thick disk, continuous line, may also be fitted leading to $\Delta =116\pm 22$ km s-1. It also points to the origin of the coordinate axes, since we are representing heliocentric velocities. Between both population components, at $v_2 \approx -58$ km s-1, we find the points representing the incomplete thick disk component of the intermediate and partially mixed subsamples. Hence both galactic components are clearly distinguishable. Thus we can now confirm that the sample of 125 km s-1 still has rotational properties similar to the thin disk component.

Therefore HIPPARCOS sample show a homogeneous and consistent kinematic behavior concerning the velocity distribution of each mixture component. Similar trends are found in the CNS3 sample, but with less accuracy.

9 Discussion

Some of the tensor elements involved in Eqs. (23) and (24) are roughly zero, since their standard errors are similar or greater than their values. In general we are working within a $2\sigma$-confidence level. For example, in both samples the only non-vanishing non-diagonal second moment is $\mu_{12}$. The third moments that may be assumed zero are $\kappa_{111}, \kappa_{113}, \kappa_{123}, \kappa_{223}, \kappa_{133}$, and $\kappa_{333}$, although for the CNS3 sample the moment $\kappa_{122}$ is very low. The vanishing four cumulants are those having the index 3 an odd-number of times. On the other hand the cumulants $\kappa_{1112}, \kappa_{1222}$, and $\kappa_{1233}$ are non-null, but much smaller than the other cumulants, and with a higher relative error. For the CNS3 sample the moment $\kappa_{1233}$ is also null. Hence the equations of the Appendix D may be simplified and reduced to some particular cases of symmetry.

In general the samples from both catalogs are qualitatively very close, and they have total central moments that are consistent with the assumption (a) of a symmetry plane u3=0. The hypothesis of case (b) is also fulfilled for both samples, according to Eq. (36) of case (a+b). Thus the differential centroid velocity can be considered null in the radial direction. However the constraint equations of the axisymmetric hypothesis (c) are not satisfied. The partial moments $\mu'_{12}$ and $\mu''_{12}$ are non-null, and the total moment $\mu_{12}$ can not be explained from a superposition of two axisymmetric normal distributions without vertex deviation. The set of non-null but small fourth cumulants, that have been outlined above, are indicators that the samples are not far from the axisymmetry hypothesis, according to case (a+b+c), but definitively the distribution has lost the axial symmetry.

The segregation algorithm has led to population parameters that are completely consistent with the actual portrait of the local kinematics. Although CNS3 subsample shows slightly high moment values, HIPPARCOS subsample is in a total agreement with population parameters obtained from non-kinematically biased samples. Thus, the proportion of 9% thick disk population is similar to the obtained by authors like Soubiran (1994), Ojha et al. (1996), or Chiba & Beers (2000). Also thin and thick disk velocity dispersions, $(\sigma'_1 ,\sigma'_2 ,\sigma'_3)= (28\pm 1, 16 \pm 2, 12 \pm 1)$ and $(\sigma''_1 ,\sigma''_2 ,\sigma''_3)= (66\pm 2, 40 \pm 9, 41 \pm 2)$, are similar to the values obtained, for example, by Soubiran (1993, 1994), Beers & Sommer-Larsen (1995), Ojha et al. (1996, 1999), Soubiran et al. (2003), as well as the lag in rotation between thick and thin disk stars, $51 \pm 3$ km s-1. However we must point out that the thin disk component, which is a complex mixture of stellar moving groups (e.g. Figueras et al. 1997; Dehnen 1998; Chereul et al. 1998; Skuljan et al. 1999) is globally well fitted by an ellipsoidal distribution.

The non-null vertex deviation of the thin disk component, $\mu'_{12}= 93 \pm 15$ $(\epsilon' \sim 9^{\circ})$, is consistent with the obtained by Dehnen & Binney (1998) and Muhlbauer & Dehnen (2003). But we must also remark the similar non-null vertex deviation of the thick disk population, $ \mu''_{12}=314 \pm 112$ ( $\epsilon'' \sim 7^{\circ}$), what is suggesting a non-axisymmetric distribution also for this component, according to Dehnen (1998).

Therefore the constraint equations provide a way for testing the geometry of the mixture distribution only through evaluation of the total cumulants, without computing the population parameters of the stellar components. In addition those parameters can be estimated from the segregation algorithm that has been deduced from the constraint equations, so that plausible values and error bars have been obtained. The use of other statistical techniques such as SEM (Stochastic, Expectation, Maximization) (Celeux & Diebolt 1985; Soubiran et al. 1990; Ojha et al. 1996), EM (Expectation, Maximization), other maximum likelihood-based methods (Robin et al. 1996; Ratnatunga & Upgren 1997; Dehnen 1998), other related multivariate sampling algorithms (Bougeard & Arenou 1990), and specially the more recent wavelets-based algorithms (e.g. Figueras et al. 1997; Chereul et al. 1998; Skuljan et al. 1999), are nowadays common in astronomy. They are very efficient in segregating particular stellar groups with common properties. Generally a global multivariate analysis from kinematic, photometric, spectroscopic, and all the available star attributes is carried out, and a lot of clusters, generally with few stars each one, are isolated. In some cases the background velocity distribution is not relevant for their analysis, although sometimes spherical or bivariate Gaussian distributions are assumed (Soubiran 1993; Ojha et al. 1996, 1999). Some of these works discuss kinematic features of HIPPARCOS Catalog, and are mainly devoted to the detection of stellar moving groups (e.g. Dehnen & Binney 1998; Asiain et al. 1999). As in our work, some of them conclude that the axisymmetric hypothesis is not completely fulfilled, although, in order to assert it, they need sometimes considerations about the mass distribution and the gravitational potential. Thus, our work must be considered a qualitative approach to the study of the velocity distribution complementary to the foregoing techniques. In the future this work may be continued by generalizing the method to n-component mixtures. Also other symmetry cases, and samples far from the galactic plane may be studied. In order to segregate normal populations it seems feasible a recursive application of the algorithm, focusing in the selection of subsamples that minimize the fitting error.

In conclusion the application of our mixture model and symmetry analysis to the subsample drawn from HIPPARCOS catalog gives a good approximation of the local kinematics, starting from statistics involving only the velocity space. In particular, the cumulants of the mixture provide meaningful information of the velocity distribution leading to the following main results: (i) The solar neighborhood can be described from two populations with normal velocity distribution, that are associated with thin and thick disk components. (ii) The cumulant constraints are consistent with the hypotheses of symmetry plane and non-differential motion in the radial direction. (iii) The velocity distribution shows a deviation of the axisymmetry hypothesis and both population velocity ellipsoids have vertex deviation in the galactic plane.

Acknowledgements
The authors wish to thank Dra. Francesca Figueras for providing the subsample from HIPPARCOS Catalog.

Appendix A


\begin{eqnarray*}1 &=& \frac {\mu_{1111}} {3\mu_{11}^{2}} = \frac {\mu_{2222}} {... ...\frac {\mu_{1233}} {\mu_{33}\mu_{12} + 2\mu_{13}\mu_{23}}\cdot \end{eqnarray*}


Appendix B


\begin{displaymath}\begin{array}{l} o^{(3)}_{111} = - \kappa_{333} d_{1}^{3} + 3... ...223} d_{3}^{2} d_{2} + \kappa_{222} d_{3}^{3} = 0 \end{array} \end{displaymath}


\begin{eqnarray*}o^{(3)}_{112} &=& - \kappa_{333} d_{1}^{2} d_{2} + \kappa_{233}... ...& - \kappa_{223} d_{1} d_{3}^{2} + \kappa_{122} d_{3}^{3} = 0. \end{eqnarray*}


Appendix C


\begin{displaymath}\begin{array}{l} p^{(3)}_{11} = \kappa_{333} d_{1}^{2} - 2\ka... ...kappa_{333} d_{2} + \kappa_{233} d_{3}\right) \\ \end{array} \end{displaymath}


\begin{eqnarray*}X^{(3)}_{1111} &=& \kappa_{3333} d_{1}^{4} - 4\kappa_{1333} d_{... ... T^{(3)}_{2} &=& - \kappa_{3333} d_{2} + \kappa_{2333} d_{3} . \end{eqnarray*}


Appendix D


\begin{eqnarray*}D_{3}^{2} &=& \frac {3\left[p^{(3)}_{11}\right]^{2}} {X^{(3)}_... ...}_{1}} {3s^{(3)}_{1}} = \frac {T^{(3)}_{2}} {3s^{(3)}_{2}} \cdot \end{eqnarray*}


References

 
Copyright ESO 2004