© ESO, 2012

1. Introduction

The differences between the positions of a set of common sources in two astrometric catalogues are conveniently described mathematically by a vector field on a sphere. Each vector materialises the difference between the two unit vectors giving the direction of the common sources in each catalogue. This feature extends to the differences between the proper motions of each source which also generates a spherical vector field. Typically when connecting two positional catalogues with n common sources, one has the coordinates for the ith source in the first catalogue and for the same source in the second catalogue. Provided the two catalogues are close to each other, the differences can be mapped as a vector field with components in the local tangent plane given by (1)for each common source between the two catalogues. We wish to analyse these fields in order to summarise their largest overall or local features by means of a small set of base functions, which is much smaller in any case than the number of sources.

Several papers in the astronomical literature have applied this overall idea with either scalar or vectorial functions. As far as we have been able to trace the relevant publications regarding the use of vector spherical harmonics in this context, this has been initiated by Mignard & Morando (1990) and published in a proceedings paper that is not widely accessible. One of the motivations of this paper is therefore to update and expand on this earlier publication and to provide more technical details on the methodology.

The general idea of the decomposition has its root in the classical expansion of a scalar function defined on the unit sphere with a set of mutually orthogonal spherical harmonics functions. The most common case in astronomy and geodesy is the expansion of the gravitational potential of celestial bodies, in particular that of the Earth. The generalisation of such expansions to vectors, tensors, or even arbitrary fields was introduced in mathematics and theoretical physics decades ago (Gelfand et al. 1963). For example, these generalisations are widely used in gravitational physics (Thorne 1980; Suen 1986; Blanchet & Damour 1986). Restricting it to the case of vectors fields, which is the primary purpose of this paper, we look for a set of base functions that would allow representing any vector field on the unit sphere as an infinite sum of fields, so that the angular resolution would increase with the degree of the representation (smaller details being described by the higher harmonics). One would obviously like for the lowest degrees to represent the most conspicuous features seen in the field, such as a rotation about any axis or systematic distortion on a large scale. It is important to stress at this point that expanding a vector field on this basis function is not at all the same as expanding the two scalar fields formed by the components of the vector field: the latter would depend very much on the coordinate system used, while a direct representation on a vectorial basis function reveals the true geometric properties of the field, regardless of the coordinate system, in the same way as the vector field can be plotted on the sphere without reference to a particular frame without a graticule drawn in the background.

Several related methods using spherical functions to model the errors in astrometric catalogues, the angular distribution of proper motions, or more simply as a means of isolating a rotation have been published and sometimes with powerful algorithms. Brosche (1966, 1970) was probably the first to suggest the use of orthogonal functions to characterise the errors in all-sky astrometric catalogue, and his method was later improved by Schwan (1977) to allow for a magnitude dependence. But in both cases the analysis was carried out separately for each component Δαcosδ,Δδ of the error, by expanding two scalar functions. Therefore, this was not strictly an analysis of a vector field, but that of its components on a particular reference frame. The results therefore did not describe the intrinsic geometric properties of the field, which should reveal properties not connected to a particular frame (see also Appendix C). Similarly, a powerful method for exhibiting primarily the rotation has been developed in a series of papers of Vityazev (2010, and references therein). The vector spherical harmonics (hereafter abridged in VSH) have been in particular used in several analyses of systematic effects in VLBI catalogues or in comparison of reference frames, such as in Arias et al. (2000), Titov & Malkin (2009), or Gwinn et al. (1997), the galactic velocity field (Makarov & Murphy 2007), the analysis of zonal errors in space astrometry (Makarov et al. 2012), but none of these papers deals with the fundamental principles, the very valuable invariance properties and the relevant numerical methods. Within the Gaia preparatory activities the VSH are used also to compare sphere solutions and in the search of systematic differences at different scales by Bucciarelli et al (2011) and more is planned in the future to characterise and evaluate the global astrometric solution.

In this paper we provide the necessary theoretical background to introduce the VSHs and the practical formulas needed to explicitly compute the expansion of a vector field. The mathematical principles are given in Sect. 2 with a few illustrations to show the harmonics of lower degrees. In the next section, Sect. 3, we emphasise the transformation of the VSHs and that of the expansion of a vector field under a rotation of the reference frame, with applications to the most common astronomical frames. Then in Sect. 4 we discuss the physical interpretation of the harmonics of first degree as the way of representing the global effect shown by a vector field, such as the axial rotation and its orthogonal counterpart, for which we have coined the term glide, and show its relationship to the dipole axial acceleration resulting from the Galactic aberration. The practical implementation is taken up in Sect. 5, where we also discuss the statistical testing of the power found in each harmonic. The results of particular applications to the FK5 Catalogue and to the simulated QSO catalogue expected from Gaia are respectively given in Sects. 6 and 7. Appendix A provides the explicit expressions of the VSH up to degree l = 4, while Appendix B deals with the practical formulas needed for numerical evaluation of the VSH, and Appendix C focuses on the relationship between expansions of the components of a vector field on the scalar spherical harmonics and the expansion of the same field on the VSHs.

2. Mathematical principles

2.1. Definitions

In this section we give the main definitions needed to use the VSHs and the theoretical background from which practical numerical methods are established and discussed later in the paper. The VSHs form an orthogonal set of basis functions for a vector field on a sphere and split into two categories referred to as the toroidal T_lm and spheroidal S_lm functions (in the physical literature the former are also called “magnetic” or “stream”, and the latter “poloidal”, “potential”, or “electric”): where u = r/r, r = |r|, r is the radius-vector of the point, and ∇ denotes the gradient operator. Clearly, for points on the surface of a unit sphere r = 1. It is also obvious that at each point on the sphere and for each l and m one has u·S_lm = 0, u·T_lm = 0, and T_lm·S_lm = 0. Taking into account that (see Fig. 1) (so that |u| = 1, |e_α| = 1, |e_δ| = 1), one gets explicit formulas (7)for the toroidal functions, and (8)for the spheroidal functions.

Fig. 1 Local frame associated to the spherical coordinates α,δ, with the unit vectors along the longitude and latitude lines.

The Y_lm are the standard spherical functions defined here with the following sign convention (9)for m ≥ 0 and (10)for m < 0, where superscript “ ∗ ” denotes complex conjugation. The associated Legendre functions are defined as (11)Equations (9)–(10) agree with the well-known formula for the associated Legendre functions (12)Different sign conventions appear in the literature, in particular regarding the place of the (−1)^m, which is sometimes used in the definition of the associated Legendre functions instead (see, e.g., Chap. 8 of Abramowitz & Stegun 1972). However, these sign differences do not influence the vector spherical functions themselves. Also one sometimes considers instead of T_lm as the toroidal vector spherical functions.

From (10) and (7)–(8) one immediately has

2.2. Orthogonality properties

The spherical functions Y_lm form an orthonormal sequence of functions on the surface of a sphere since (15)which is also complete in the Hilbert space of the square-integrable functions on a sphere (making it an orthonormal basis of in the L² norm). As in Fourier expansions, the completeness property is hard to establish and is related to the fact that spherical harmonics are eigenfunctions of a special kind of differential equations (see, for example, Morse & Feshbach 1954). This is accepted and not further discussed in this paper directed towards astronomical applications. Here and below one has dΩ = cosδ   dδ   dα, and the integration is taken over the surface of the unit sphere: 0 ≤ α ≤ 2π, −π/2 ≤ δ ≤ π/2, and δ^ij is the Kronecker symbol: δ^ij = 1 for i = j and δ^ij = 0, otherwise.

Thanks to the completeness property a (square-integrable) complex-valued scalar function defined on a unit sphere f(α,δ) can be uniquely projected on the Y_lm: (16)where the Fourier coefficients f_lm age given by (17)The equality in (16) strictly means that the right-hand side series converge (not necessarily pointwise at discontinuity points or singularities but at least with the L² norm) to some function f, and given the completeness and the orthonormal basis, the coefficients f_lm are related to f by (17). A truncated form of (16) with l ≤ l_max < ∞ is an approximate expansion for which the equality in (16) does not strictly hold (in general).

Similarly, the VSH form a complete set of orthonormal vector functions on the surface of a sphere (with the inner product of the L² space): As we noted after (3) above, one also has the orthogonality between two vectors in the usual Euclidean space, (21)which holds for any point on the sphere.

Any (square-integrable) complex-valued vector field V(α,δ) defined on the surface of a sphere and orthogonal to u (radial direction), (22)can be expanded in a unique linear combination of the VSH, (23)where the coefficients t_lm and s_lm can again be computed by projecting the field on the base functions with

2.3. Numerical computation of the vector spherical harmonics

Analytical expressions of the VSH are useful for the lower orders to better understand what fields are represented with large-scale features and also to experiment by hand on simple fields in an analytical form. This helps develop insight into their properties and behaviour, which pays off at higher orders when one must rely only on numerical approaches. This is also a good way to test a computer implementation by comparing the numerical output to the expected results derived from the analytical expressions. Analytical expressions for the vector spherical harmonics of orders l ≤ 4 are given explicitly as a function of α and δ in Appendix A.

To go to higher degrees, one needs to resort to numerical methods. To numerically compute the two components of T_lm(α,δ) and S_lm(α,δ) given by (7)–(8), one first needs to have a procedure for the scalar spherical harmonics. Given the form of the Y_lm(α,δ) in (9), the derivative with respect to α is trivial, and for δ only, the derivative of the Legendre associated functions needs special care. There are several recurrence relations allowing the Legendre functions to be computed, but not all of them are stable for high degrees. Other difficulties appear around the poles when δ ± π/2 where care must be exercised to avoid singularities. Non-singular expressions and numerically stable algorithms for computing the VSH components are available in the literature, and the expressions we have implemented and tested are detailed in Appendix B.

2.4. Expansions of real functions

Describing real-valued vector fields with complex-valued VSHs leads to unnecessary redundancy in the number of basis functions and fitting coefficients. In general, the coefficients in (23) are complex even for a real vector field. But for a real field V(α,δ), the expansion must be real. Given the conjugation properties of the VSH, it is clear that for a real-values function one gets Then from the decomposition of the coefficients into their real and imaginary parts as one gets at the end by rearranging the summation on m ≥ 0 (30)which is real. It is obvious that and, therefore, and .

It is useful to note that the use of real and imaginary parts of the VSHs as in (30) is analogous to the use of “sin” and “cos” spherical functions suggested e.g. by Brosche (1966). However, we prefer to retain notations here that are directly related to the complex-values VSHs and the corresponding expansion coefficients.

Equation (30) provides the basic model (with only real numbers) to compute the coefficients by a least squares fitting of a field given on a finite number of points, which are not necessarily regularly distributed. In this case the discretised form of the integrals in (18), (19), and (20) is never exactly 0 or 1, and the system of VSH on this set of points is not fully orthonormal. Then one cannot compute the coefficients by a direct projection. However, one can solve the linear model (30) on a finite set of coefficients t_lm and s_lm. Provided the errors are given by a Gaussian noise, the solution should produce unbiased estimates of the true coefficients. This is discussed further in Sect. 5.

2.5. Relation between the expansion in vector and scalar spherical harmonics

Once we have fitted a vector field on the model (30), we have the expansion in the form (23) with l ≤ l_max. The components V^α and V^δ of the vector field V in the local basis e_α and e_δ as in (22) are expressed with the same set of coefficients t_lm and s_lm on their respective components of the VSH. But since the V^α and V^δ are scalar functions on a sphere, they could also have been expanded independently in terms of scalar spherical functions Y_lm as and , providing another representation, which has been used as mentioned in the Introduction to analyse stellar catalogues in the same spirit as with the VSHs. It is therefore important to relate the two different sets of coefficients for the purpose of comparing and discussing the major difference between the two approaches (see Sect. 3). This issue is not central to this paper so the mathematical details are deferred to in Appendix C. One must, however, note that the relations between coefficients and , on the one hand, and t_lm and s_lm, on the other, are rather complicated and involve infinite linear combinations. It means, for example, that the information contained in a single coefficient t_lm is spread over infinite number of coefficients and (see, e.g. Vityazev 2010, for a particular case l = 1).

3. Transformation under rotation

3.1. Overview

A very important mathematical and practical feature of the expansion of a vector field on the basis of the VSHs is how they transform under a rotation of the reference axes. In short, when a vector field is given in a frame and decomposed over the VSHs in this frame, one gets the set of components t_lm,s_lm in this frame. Rotating the reference frame (e.g. transforming from equatorial coordinates to galactic coordinates) to the frame , one gets the transformed vector field that can be fitted again on a set of VSHs defined in the frame . This results in a new set of components and , representing the vector field in this second frame. Significant additional work may be required to carry out the whole transformation of the initial field and perform the fit anew in the rotated frame. Fortunately due to the narrow relationship between VSHs and group representation, there is in fact a dramatic shortcut to this heavy procedure that adds considerably to the interest of using the VSH expansions in astronomy. This is illustrated in the self-explanatory accompanying commutative diagram, showing that to a space rotation (in the usual three-dimensional space) R corresponds an operator ℛ acting in the space of the VSHs and allows transforming t_lm and s_lm given in the initial frame into an equivalent set in the rotated frame:

The amount of computation is negligible since this is a linear transformation between the coefficients of a given degree l. Overall, the set of t_lm and s_lm is transformed into itself without mixing coefficients of different degrees.

This operator ℛ has been introduced by Wigner in 1927 as the D-matrices. A convenient reference regarding its definition and properties is the book of Varshalovich et al. (1988). Section 7.3 of that book is specifically devoted to this topic. The Wigner matrices have been used for decades in quantum mechanics, but probably they are not so well known in fundamental astronomy. We confine ourselves here to a quick introduction of the main formulas and conventions used in this paper. All the mathematical formulas have been implemented in computer programs (independently in FORTRAN and Mathematica), which can be requested from the authors.

For practical astronomical applications, there is an unexpected benefit in this transformation under space rotation, which is not shared by the separate analysis on spherical harmonics of the components like Δαcosδ and Δδ, or their equivalent for the proper motions. In general the two components for a given source obtained from observations do not have the same accuracy and are correlated. In the case of space astrometry missions like Hipparcos or Gaia, the correlations are generally smaller in the ecliptic frame due to the symmetry of the scanning in this frame compared to the correlations in the equatorial frame, where the solution is naturally available. Performing a least squares adjustment to compute t_lm and s_lm with correlated observations is a complication and standard pieces of software often allow only for diagonal weight matrix. One can easily get round this difficulty by carrying out the fitting in a frame where the correlations can be neglected (assuming such a frame does exist), and then rotate the coefficients with the Wigner matrix into the equatorial frame. This is equivalent to rigorously performing the fit in the equatorial coordinates with correlated observations. On the other hand, by fitting the components with scalar spherical harmonics, it is cumbersome to take the correlations between the components properly into account, and the coefficients in one frame do not transform simply into their equivalent into a rotated frame.

More generally, if we allow for the correlations in the frame where the least squares expansion over the VSHs is carried out, one can still apply the Wigner rotation to the results and obtain the expansion in a second frame. The results would be exactly the same as if a new least squares fit had been performed in this frame by rigorously propagating the non-diagonal covariance matrix of the observations to this second frame to generate the new weight matrix. Using the Wigner matrix in this context is conceptually much simpler and in keeping with the underlying group properties of the VSHs, even though there might be no definite advantage for numerical efficiency (but no disadvantage as well).

3.2. Mathematical details

Consider two rectangular Cartesian coordinate systems and of the same handedness. Two coordinate systems are related by a rotation parametrised by three Euler angles a, b, and c. For a given vector with components x in , the components in read as (31)where R_z and R_y are the usual matrix of (passive) rotation, respectively, about axes z and ywhere the order of arguments in the rotational matrices in (31) must be exactly followed. The inverse transformation reads (34)Here the convention z − y − z is used for the sequence of rotations instead of the usual z − x − z that is more common in dynamics and celestial mechanics. The former choice has the great advantage of giving a real matrix (see below) for the intermediate rotation, and this convention is universally used in quantum mechanics and group representation. The values of a, b, and c for transformations between the three usual frames in astronomy are listed in Table 1. The relation between the z − y − z and z − x − z conventions is given by (35)where R_x are the usual matrix of rotation about axis x: (36)Then the transformations between the scalar and vector spherical functions defined in and are similar to each other and read as where and are angles defined as right ascension and declination (generally speaking, the longitude and latitude of the spherical coordinates) from the components , while α and δ are those derived from the components x. In both cases Eq. (4) is used. Here , |m| ≤ l, |n| ≤ l are the Wigner D-matrices (generalised spherical functions) defined as Many efficient ways to evaluate and can be found in relevant textbooks. For example, Sect. 4.21 of Varshalovich et al. (1988) gives the expression of in the form of a simple Fourier polynomial (42)where can be precomputed (using the obvious simplification of (41)) as constants and used in further calculations for a given rotational angle b. However, even a straight implementation of (40)–(41) does the job without any problem for low to medium degrees.

The inverse transformations of the spherical functions read as where the superscript “∗” denotes complex conjugation.

Now what matters for our purpose is the transformation of the VSH coefficients t_lm and s_lm as in (23) under a rotation of the coordinate system, and, more generally, that of coefficients f_lm in the case of an expansion of a scalar field in the scalar spherical harmonics as given by (16).

Given the completeness of the scalar and vector basis functions, we have defined the following expansions for scalar field f and vector field V on sphere (cf. (16) and (23)): Substituting (37)–(39) or (43)–(45) into (46)–(47) and changing the order of summations, one gets the important transformation rules for the coefficients and the inverse transformations Once the expansion has been computed in one frame, it is very easy to compute the expansion in any other frame related by a rotation from the initial frame. Numerical checks have been performed with the analysis of a vector field in two frames, followed by the transformations (direct and inverse) of the coefficients t_lm and s_lm) with (48)–(50) and (51)–(53).

Fig. 2 Examples of spheroidal harmonics vector fields: S1,1, S2,1, S5,3, S10,5, where one sees the change in the angular resolution with higher degrees. The S1,1 is a glide flow in the direction of the y-axis.

These transformation rules allow us also to prove the invariance of the Euclidean norm of the set of coefficients of a given degree: This can be easily demonstrated by substituting (48)–(50) into corresponding equation and using the well-known relation (57)which is a special case of the general addition theorem for the Wigner D-matrices (Varshalovich et al. 1988, Sect. 4.7). The invariance of will be used in Sect. 5.2 to assess the level of significance of the coefficients.

4. Relation of the first-order VSH to global features

The multipole (VSH) representation of a vector field has some similarity with a Fourier or a wavelet decomposition since as the degree l increases, one sees smaller details in the structure of the field. In the case of the comparison between two catalogues, increasing the degree reveals the zonal errors of the catalogue, while the very first degrees, say l = 1 and l = 2, show features with the longer “wavelengths”. In particular, the first degree both in the toroidal and spheroidal harmonics is linked to very global features, such as rotation between the two catalogues or a systematic dipole displacement, like a flow from a source to a sink located at the two poles of an axis. It must be stressed that these global features are found in a blind way, without searching for them. They come out naturally in the decomposition as the features having the least angular resolution, and they are interpreted as a rotation or a dipole glide a posteriori.

4.1. Rotation

Consider a infinitesimal rotation on a sphere given by an the rotation vector R whose components in a particular frame are (R₁,R₂,R₃). This rotation generates the vector field (58)where u is the unit radial vector given by (4). Equation (58) can be projected on the local tangent plane and leads to From the explicit expressions of the toroidal harmonics given in Table A.1, one sees that any infinitesimal rotation field has non-zero projection only on the T_1m harmonics, and is orthogonal to any other set of toroidal or spheroidal harmonics. Therefore (58) can be written as a linear combination of the T_1m. One has explicitly the equivalence (61)where allowing the three parameters R of the rotation to be extracted from the harmonic expansion of degree l = 1. A change of reference frame will lead to the same rotation vector expressed in the new reference frame. Given the importance of a rotation between two catalogues or in the error distribution of a full sky astrometric catalogue, having the rotational vector field as one of the base functions is a very appealing feature of the VSH. This contrasts with the analysis of each components of the vector on scalar spherical harmonics, in which a plain rotation projects on harmonics of any degree as shown in Vityazev (1997, 2010) or in Appendix C. However, the analysis of each component as scalar fields retains its interest if one has good reasons to assume that the field has been generated by a process in which the two spherical coordinates have been handled more or less separately. In this case there is even no reason to change the frame for another one since the two components may have very different statistical properties.

4.2. Glide

The rotation field is the most common fully global signature on a spherical vector field, and often it even seems to be the only one deserving the qualification of global. However, if we consider a rotational field G × u, one may draw at each point on the sphere a vector perpendicular to the rotation field and lying in the tangent plane to the sphere. In this way we generate another field V^G with axial symmetry, fully orthogonal to the rotation field G × u. This new field has no projection on the T_1m harmonics and, in fact, on no other T_lm. Its components are given by where the vector G = (G₁,G₂,G₃) determines the orientation and the magnitude of the field, as R did for the rotational field. Going to Table A.1 it is easy to see that it is equivalent to the S_1m field with the decomposition (66)where showing its close relationship with a rotation field. A corresponding picture can be seen in Fig. 2 for an axis along the y-axis.

Fig. 3 Examples of toroidal harmonics vector fields: T1,1, T2,1, T5,3, T10,5. One sees the change in the angular resolution with higher degrees. The T1,1 is a simple rotation about the y-axis.

The same can be stated somewhat differently: we have introduced the glide¹ vector field V^G, which is orthogonal to both radial vector u and the rotation field G × u, or equivalently one has (69)which shows clearly that the associated vector field is just the projection of G on the surface of the unit sphere (we have above G minus its radial component), that is to say the projection of an axial, or dipole, field on the sphere.

As a pattern, a glide field is as global a feature as a rotation, and can be seen as a regular flow between a source and a sink diametrically opposed. From the astronomical point of view this is a field associated to a motion of the observer toward an apex, with all the stars showing a kinematical stream in the opposite direction. For extragalactic sources like quasars, this also characterises the effect induced on QSO systematic proper motions resulting from the acceleration of the observers with respect to the frame where the QSOs are at rest on average (Fanselow 1983; Sovers et al. 1998; Kovalevsky 2003; Kopeikin & Makarov 2006; Titov et al. 2011). The main source of this acceleration is thought to be the centripetal acceleration from the galactic rotation. Its unambiguous determination should be achieved with Gaia. It is important to recognise that in the general analysis of catalogues or spherical vector fields, rotation is not more fundamental a feature than glide, and both must be searched for simultaneously. Otherwise, for most of the distributions of the data on the sphere, when the discrete orthogonality conditions between the VSHs are not fully satisfied, there will be a projection of the element not included in the model in the subspace generated by the other one. Thus, fitting only the rotation, the rotation vector will be biased if there is a glide not accounted for in the data model. And reciprocally for the glide without rotation properly fitted.

Fig. 4 Rotation and glide with the z-axis shown here with a Mercator projection to emphasise the close relationship between the two vector fields. The constant displacement in longitude from the rotation has its counterpart with a constant glide displacement in the reduced latitude (Mercator latitude).

4.3. Duality of the rotation and glide effects

Finally the duality between the infinitesimal rotation and glide can be made striking if the fields are plotted in the Mercator projection, as illustrated in Fig. 4. Both R and G are along the z-axis and have the same magnitude. Thus they are rotation and glide with the same axis. The effect of a small rotation is seen as a uniform translation in longitude, identical to all latitudes, as expected for a cylindrical projection. For the glide, the translation is uniform in the reduced latitude given by the Mercator transformation (70)and does not depend on the longitude.

There is a clear similarity between the two fields (both represent a mapping of the sphere onto itself with two fixed points: the poles), but also an important difference seen with the Mercator projection: for the rotations, the mapping of a flow line (the continuous line tangent to the vectors of the field) is bound between the two extreme boundaries −π and  + π, which must be identified as the same line on the sphere. The rotation axis is the natural source of a cylindrical symmetry shown by the Mercator projection (or any other cylindrical projection). For the glide, it goes very differently since the two end points of a flow line are sent to infinity, in both directions, and can only be identified in the projective space, a much more complex manifold than the Euclidean cylinder. Could this important topological difference be the underlying reason why only the rotation is usually considered as the only obvious global structure? We do not have the answer and leave this issue open for the moment. As long as infinitesimal transformations are concerned, which is exactly what matters for the analysis of astrometric catalogues, the two effects are in fact very similar and must be taken into account together as a whole.

5. Practical implementation

In principle, once the difference between two catalogues has been formed or when a vector field on a sphere (e.g. a set of proper motions) is available, it is normally trivial to project the underlined vector field on the VSH, by computing the integrals (24) with an appropriate numerical methods. This will work as long as the distribution of the sources on the unit sphere preserves the orthogonality conditions between the different base functions. What is precisely meant is the following: using the same numerical method to evaluate the Fourier integrals, one should check that Eqs. (18)–(20) hold with a level of accuracy compatible with the noise up to some degree l_max. Unless one has many sources almost evenly distributed, this condition is rarely fulfilled. However, the fact that the coefficients can be expressed by integrals, that is to say by direct projection on the base functions, is a simple consequence of the application of a least squares criterion leading to a diagonal normal matrix when the orthogonality condition holds. Starting one step back to the underlying least squares principle we can directly fit the coefficients to the observed components of the field with the model (30), by applying a weighting scheme based on the noise in the data. In case the orthogonality conditions hold, the solution will be equivalent to the evaluation of the integrals. When this assumption is not true, meaning with the discrete sampling and/or uneven weight distribution at the source coordinates, the harmonics are not numerically orthogonal, but the least squares fitting provides an unbiased estimate of the coefficients t_lm and s_lm, together with their covariance matrix.

5.1. Least squares solution

The fitting model is given by (30) up to the degree l = l_max. It is fitted to a set of N data points sampled at (α_k, δ_k), k = 1,2,···   ,N, where at each point the two components of the vector field Δ (e.g., difference between two catalogues, or proper motion components) are given as (71)together with the estimates of their associated random noise and . For the sake of simplicity we assume here that the covariance matrix between the observations is diagonal, meaning the observations are realisations of independent random variables. The observations are fitted to the model (30) with the least squares criterion (72)From (30) and (22) we see that with N observations one has 2N conditional equations to determine the 2l_max(l_max + 2) coefficients , m = 1,2,..., l, l = 1,2,..., l_max. The size of the design matrix is 2N    ×    2l_max(l_max + 2), which can be very large in analysis of astrometric surveys, with N above 10⁵−10⁶ and a sensitivity on a small enough angular scale to represent the regional errors leading to explore up to l_max ~ 15−20. The data storage could become a problem if one wishes to store and access the full design matrix. One solution to get round this problem is to build up the normal matrix on the fly, taking its symmetry into account. For a large data set (), this bookkeeping part will be the most demanding in terms of computation time. For fixed N the computational time for this part of calculations grows like . In general the problem is rather well conditioned, with small off-diagonal terms, and is not prone to generate numerical problems to get the least squares solution.

5.2. Significance level

The least squares solution ends up with an estimate for each of the coefficients included in the model. The model must go to a certain maximum value of l_max by including all the component harmonics of order m, for each l ≤ l_max. The total number of unknowns is then 2l_max(l_max + 2) and grows quadratically with l_max. Within a degree l, one can decide on the significance of each coefficient of order m, from its amplitude compared to the standard deviation. However, these individual amplitudes are not invariant by a change of reference frame, like the components of a vector in ordinary space, and the significance of a particular component is not preserved through a rotation of the frame. On the other hand, the whole set of t_lm or s_lm for a given l behaves like a vector under rotation, and the significance (or lack of) of the vector itself, that is to say its modulus, is intrinsic and remains true or false in any other rotated frame. Therefore it is more interesting and more useful in practice to investigate the significance level for each degree l rather than that of the individual coefficients. The reason behind this statement is deeply rooted in the invariance property of the Y_lm under a rotation in the three-dimensional Euclidean space and the fact that the set of spherical harmonics Y_lm of degree l is the basis of an irreducible representation of the rotation group of the usual three-dimensional space. Therefore, this set is globally invariant under a rotation, meaning there is linear transformation linking the Y_lm in one frame to their equivalent of same degree l in the other frame (see Sect. 3).

Therefore, one should consider the 2l + 1 elements of the set Y_lm, for m = −l, ...l like the components of a vector that are transformed into each other under a rotation, keeping the Euclidean norm unchanged in the process. Given the linearity of the derivative operator, this property extends naturally to the T_lm and S_lm. For the same reason, the coefficients t_lm and s_lm in the expansion of a vector field must be considered together within a degree l and not examined separately, unless one has good reasons to interpret physically one or several coefficients in a particular frame: the Oort coefficients in the Galactic frame are one of the possible examples, as shown in Vityazev (2010). More generally a wide range of parameters in Galactic kinematics have very natural representations on vector harmonics, such as the solar motion towards the apex or a more general description of the differential rotation than the simple Oort’s modelling (see for example Mignard 2000).

Therefore it is natural to look at the power of the vector field and at how this power is spread over the different degrees l. Let  be the power of a real vector field V(α,δ) on a sphere defined by the surface integral², (73)or its discrete form used hereafter with numerical data, (74)where n is the number of samples on the sphere.

When projected on the VSH one finds easily with (18), (19), and (23) that (75)For a real vector field having the symmetries (26)–(27) one gets (76)The quantities and are the powers of toroidal and spheroidal components of the degree l of the field. Both quantities are scalars invariant under rotation of the coordinate system (see Sect. 3). Thus and represent an intrinsic property of the field, similar to the magnitude of a vector in an ordinary Euclidean vector space.

We can now take up the important issue of the significance level of the subset of the coefficients s_lm and t_lm for each degree l. In harmonic or spectral analysis, this is always a difficult subject, because it needs to scale the power against some typical value that would be produced by a pure noisy signal. In general one can construct a well-defined test, with secure theoretical foundations, but with assumptions about the time or space sampling, which are rarely found in the real world. Therefore the theory provides a good guideline, but not necessarily a fully safe solution applicable in all situations. The key issue here is to construct a criterion that is well understood from the statistical view point and robust enough to be usable for mild departure from the underlying assumptions, like regular time sampling in time series. We offer below three possibilities for testing whether the power  and  is significant and therefore that there is a true signature of the VSH of degree l in the vector field. In the following we drop superscripts “s” and “t” and write simply . Relevant formulas below are valid separately for and or for . The involved coefficients will be denoted as r_lm to represent either t_lm or s_lm or both.

Start from the null hypothesis H₀ that the signal is just made of noise, producing some values of in the analysis. One needs to build a test to tell us that above a certain threshold, the value of  is significant and would only be produced by pure noise with a probability p ≪ 1. To build the test properly, the probability density function (PDF) of is required. It can be obtained with the additional assumptions that the sources are rather evenly distributed on the sphere. The least squares fitting also gives us the variance-covariance matrix of the unknowns s_lm and t_lm. Given the orthogonality properties of the VSH, the normal matrix is almost diagonal if the distribution of sources over the sky and the associated errors of the quantity representing the vector field are homogeneous. In this case the covariance matrix is also almost diagonal.

Consider the coefficients r_l0, , and in (30) with a fixed degree l. Each diagonal term of the relevant part of the normal matrix (unweighted here) takes the form where the sum extends over all the sources on the sphere. Here R_lm stands for either T_lm or S_lm. The first term with m = 0 does not depend on the right ascension α, while the others have a part in cos²mα and sin²mα, whose average is 1/2 for α ranging from 0 to 2π and a relatively regular distribution in α. Given the normalisation of the VSHs, we now see that the diagonal terms with m ≠ 0 are just two times the term with m = 0. The covariance matrix being the inverse of the normal matrix, one eventually has the variances of the unknowns as where is the variance of r_l0. Then dividing in (76) by leads to the normalised power: (83)With (81), (82) and the assumption that the signal is pure white Gaussian noise, this is the sum of 2l + 1 squares of standard normal random variables with zero mean and unit variance. It is well known that W_l follows a χ² distribution with n = 2l + 1 degrees of freedom. With this normalisation, W_l is strictly proportional to the power .

Fig. 5 Differences in positions and proper motions between the FK5 catalogue and Hipparcos in 1991.25. The typical arrow sizes in the plots are respectively 100 mas and 2.5 mas/yr.

Under the H₀ hypothesis that the signal is just white noise, one can build a one-sided test to decide whether a harmonic of degree l is significant with the probability level γ: (84)where Γ(a,x) and Γ(a) are the incomplete and complete gamma functions, respectively. This quantity can be easily computed. In some cases it may be easier to use the transformation of Wilson & Hilferty (1931), (85)where n is the number of degrees of freedom. It is well known that Z approximately follows, even for small n, a standard normal distribution of zero mean and unit variance. A test risk of γ = 0.01 is achieved if Z > 2.33 or at the level of γ = 0.025 is Z > 1.96. Experiments on simulated data have shown that the key elements in the above derivation, namely that is practically achieved when sources are rather regularly distributed in longitudes, even with irregularities in latitudes. The reason is basically that this factor follows from the average of cos²mα and sin²mα, while the latitude effects are shared in a similar way in all the r_lm.

Alternatively, instead of scaling the power by only and taking advantage of the asymptotic properties (80)–(82), it would have been more natural to define a reduced power with (86)which is not strictly proportional to the power. But what matters in practice is whether one can construct an indicator, related to the power, with well defined statistical properties, so one may relax the constraint of having exactly the power, provided that a reliable test of significance can be built. For least squares solutions with Gaussian noise, and small correlations between the estimates of the unknowns, it is known that each coefficient follows a normal distribution of zero mean whose variances can be estimated with the inverse of the normal matrix. Therefore, is also a sum of squares of standard normal variables and its PDF is that of a χ² distribution with 2l + 1 degrees of freedom. For regular distributions of the data points, the two reduced powers W_l and are equivalent. So either test should give the same result. In (86) the multiplying factor two before the sum in (76) has been absorbed by (i.e. by and being approximately a factor smaller than σ_l0).

Finally, we can also go back to the first test in (83), which is the power scaled by the variance of the component s_l0. Since the variance is in practice an estimate of a random variable based on the data, it may happen that the value produced by the least squares is unusually small, which will make W_l too large and trigger a false detection. Therefore, while this test has a good statistical foundation for regular distribution of sources, it is not very robust to departure from this assumption. We suggest improving this robustness, a very desirable feature in the real data processing, by replacing the scaling factor by the average of the set , m = 1,...,l whose mathematical expectation is the true value of , but with a smaller scatter than . Let be this average, we have then instead of (83) a new expression for the reduced power (87)which is directly proportional to the power and follows a χ² distribution with 2l + 1 degrees of freedom.

From a theoretical standpoint, the three indicators are equivalent for a large number of data points evenly distributed on a sphere. In practical situations, with some irregularities in the angular distribution and/or a limited number of sources, the significance testing with should by construction be more robust against random scatter affecting the estimated standard deviations. Preliminary experiments have confirmed this feature, although a wider Monte Carlo simulation is needed to quantify the advantage.

6. Application to FK5 Catalogue

The comparison between the FK5 and Hipparcos was published soon after the release of the Hipparcos catalogue by Mignard & Froeschlé (2000) primarily provides the orientation and the rotation (spin) of the FK5 system with respect to the Hipparcos frame. The differences in position and proper motions are shown in Fig. 5 for the whole set of stars used in the comparison. In the positional plot one clearly sees a general shift toward the south celestial pole and few regional effects of smaller scales.

The global analysis carried out at that time was restricted to fitting these six parameters from the positional and proper motion differences from Hipparcos to enable the transformation of astronomical data from the FK5 frame into the Hipparcos frame or, more or less equivalently, to the ICRF. The remaining residuals were large and considered as zonal errors, not reducible to a simple representation and presented in the form of plots of Δαcosδ and Δδ as a function of the position of the stars. Large zonal errors were clearly shown by these plots, but no attempt was made to model these residuals.

Using the VSH now, with the software developed during this work, it is possible to draw refined conclusions and to show that the terms of low degree beyond the rotation are truly significant. In particular we can see in Table 4 that the positional glide term from the spheroidal harmonic of degree l = 1 is larger than the orientation term. Terms of degree l = 2 are also large, both in position and proper motion, a feature that had not been noticed as clearly with the earlier analysis based on a priori model.

6.1. Global differences

This concerns the rotation and the glide between the two systems, or equivalently using the VSH terminology, the toroidal and spheroidal harmonics with l = 1. The results for the rotation shown in Table 2 are fully compatible with the earlier analysis of Mignard & Froeschlé (2000) and confirm a systematic difference between the FK5 inertial frame and the ICRS materialised by Hipparcos. The interpretation of the systematic spin is also discussed in the same paper and seems to be related to the precession constant. These new values also agree with a similar analysis performed by Schwan (2001), provided the epoch is corrected from 1949.4 (mean epoch of the FK5 observations) to 1991.25 (Hipparcos Catalogue reference epoch). A new result is given in Table 3 for the glide, but is only significant in the positions. It describes a systematic difference of about 60   cosδ  mas in the declination system between the two frames. This is a very large systematic effect given the claimed accuracy of the FK5 below 50 mas. Nothing similar is seen for the proper motions. Any transformation of a catalogue tied by its construction to the FK5 system and later on aligned to the HCRF (Hipparcos frame aligned to the ICRF) must absolutely take this large systematic difference into account in addition to the rotation. The results just shown correspond to the harmonics l = 1 from a fit of the position and proper motion difference to l_max = 10. They are slightly different with another choice of the maximum degree, but within the formal error. In both tables the standard deviations quoted are derived from the covariance matrix of the least squares solution scaled by the rms of the post-fit residuals.

6.2. Regional differences

The analysis at higher degree, up to l = 10, clearly shows significant power in most degrees, indicating a definite structure in the space distribution of the differences. Given the absence of zonal error in the Hipparcos catalogue, at the level seen in this analysis, they must come entirely from the FK5. A detailed analysis of these regional variations was performed by Schwan (2001) by projecting separately each components (difference in right ascension and declination) on a set of scalar orthogonal polynomials, instead of using the vectorial nature of the differences. Given the numerous positional instruments used to produce the FK5, some measuring only the right ascension, others only the declination, others giving both, analysing the components in equatorial coordinates separately makes sense and may reveal structure scattered in many degrees with the VSHs. Here we simply want to illustrate the relevance of the VSH to analyse the difference between two real catalogues. A detailed discussion of the application to the FK5 system is deferred to an independent paper.

6.2.1. Regional differences in position

The amplitudes, defined as the square root of the unweighted power, in each degree for the positional differences between the FK5 and Hipparcos in 1991.25 are given in Table 4, for the toroidal (left part of the table) and spheroidal (right part of the table) harmonics. The level of significance computed with Eq. (87) and the normal approximation of Eq. (85) is shown in the columns labelled Z. With l = 1 one recovers the very large significance of the rotation and the glide. Harmonics l = 2 are also very large, explaining together with the glide a significant part of the systematic differences between the two catalogues. Then the powers remain statistically very significant up to l = 7, although with small amplitudes, and then decline. With the values of the t_lm and s_lm (not given here), one could produce an analytical transformation relating the two systems, in the same spirit as in Schwan (2001). With the last column one sees that the expansion to l = 10 explains about 58% of the variance seen in the data for the position and 32% for the proper motion.

6.2.2. Regional differences in proper motions

The same analysis has been done for the proper motions. The powers and their significance level expressed with a standard normal variable are listed in Table 5. Again the signal is strong in the harmonics of the first two degrees and up to l = 7 primarily for the toroidal components.

7. Analysis of the expected Gaia results

Another important application of the VSH technique is an analysis of mathematical properties (e.g. systematic errors) of an astrometric catalogue, as well as extraction of physical information from the catalogue. Here we concentrate on the QSO catalogue as part of the expected Gaia catalogue. Our goal is to investigate (i) the anticipated accuracy of the determination of the rotational state of the reference frame and the acceleration of the solar system with respect to the QSOs, and (ii) the expected estimate of the energy flux of the primordial (ultra-low-frequency) gravity waves. From the mathematical point of view, this amounts to checking the accuracy of determining the VSH coefficients of orders l = 1 and l = 2. The toroidal coefficients of order 1 describe the rotational state of the reference frame with respect to the QSOs, while the spheroidal coefficients of order 1 give the acceleration of the solar system with respect to the QSOs, which shows up as a glide (Fanselow 1983; Sovers et al. 1998; Kovalevsky 2003; Kopeikin & Makarov 2006; Titov et al. 2011). The VSH coefficients of order 2 can be related to the energy flux of the ultra-low-frequency gravity waves (Gwinn et al. 1997).

7.1. Simulated QSO catalogue from Gaia

We simulated the Gaia QSO catalogue with all the details that we know at the time of writing. The total number of QSOs observable by Gaia – 700 000 – was chosen according to the analysis of Mignard (2012). The realistic distribution of the QSOs in the Gaia integrated magnitude G was taken from Slezak & Mignard (2007; see also Robin et al. 2012). The distribution of the Gaia-observed QSOs over the sky was chosen according to the standard Gaia extinction model (Robin et al. 2012). The model follows the realistic distribution of the dust in the Galaxy, so that the QSO distribution used in this section is not homogeneous on the sky and depends on both angular coordinates with the most pronounced feature being the avoidance area close to the Galactic plane. The proper motions of QSOs are generated from a randomly selected rotation of the reference frame and an acceleration of the solar system’s barycentre with respect to the QSOs with magnitude and direction expected from the circular motion of the solar system with respect to the galactic centre (i.e. the proper motions are generated from the full set of VSH coefficients of order 1). This resulted in our simulated “true” (noise-free) catalogue of QSOs.

The expected accuracies of astrometric parameters (positions, proper motions, and parallaxes) as functions of the G magnitude and the position on the sky are taken from the standard Gaia science performance model (ESA 2011; de Bruijne 2012). Then the QSO simulated solution with Gaia is obtained by adding Gaussian noise with the corresponding standard deviation to each parameter of each source.

One more feature used in our simulations is that we take the correlations between the estimates of the two components of proper motion into account. As a model for the correlations, the empirical distribution is calculated from the corresponding histogram for the Hipparcos catalogue as shown in Fig. 3.2.66 of ESA (1997, vol.1, Sect. 3.2). We neglect herewith the complicated distribution of the correlation on the sky as shown in Figs. 3.2.60 and 3.2.61 of ESA (1997). The generated random correlation is then used to generate correlated Gaussian noise for the components of the proper motion of each source.

7.2. Main results of the analysis of the simulated QSO catalogue

From the resulting realistic QSO catalogue expected from Gaia, we fitted the VSH coefficients for maximal order l_max ranging from 1 to 15. The correlations between the components of the proper motion for the same star are taken into account in the fit by using block diagonal weight matrix with 2 × 2 blocks corresponding to each source. Accounting for these correlations in the fit generally modifies the estimates at the level of one tenth of the corresponding standard deviations. We repeated these simulations (including generation of the QSO catalogue) tens of times. For a typical simulation, the resulting biases of the estimates of the rotation of the reference frame and the acceleration of the solar system, as well as the standard errors of those estimates are shown in Fig. 6. The biases are computed as differences between the estimated values and true ones (those put in the simulated “noise-free” catalogue). For the standard errors of the estimates, we took the formal standard deviations from the fit multiplied by the square root of the reduced χ² of the fit. Comparing of the estimates for different values of l_max as shown in Fig. 6 is a useful indicator of the reliability of the estimate. Our analysis generally confirms the assessments of the anticipated Gaia accuracy for the reference frame and acceleration published previously: both the rotational state of the reference frame and the acceleration of the solar system could be assessed with an accuracy of about 0.2 μas/yr. However, this accuracy can only be reached under the assumption that (1) Gaia will be successfully calibrated down to that level of accuracy (that is, the errors of astrometric parameters remain Gaussian down to about 0.2 μas); and (2) real proper motions of the QSOs (transverse motions of the photocentres due to changing QSO structure, etc.) are not too large and do not influence the results. Whether these assumptions are correct will not be known before Gaia flies.

Fig. 6 True errors and formal errors of the estimates of the absolute values of the rotation and acceleration from the fits with 1 ≤ lmax ≤ 15 for the simulated QSO catalogue matching expected Gaia accuracy.

Finally, we discuss the expected upper estimate of the energy flux of the ultra-low-frequency gravity waves that Gaia could measure. We use the computational recipe described by Gwinn et al. (1997, see Eq. (11) and the discussion in the first paragraph of their Sect. 4). Our analysis shows that Gaia can be expected to give an estimate at the level of Ω_GW < 0.00008   h^-2 for the integrated flux with frequencies ν < 6.4 × 10^-9   Hz. Here, h = H/(100   km   s^-1/Mpc) is the normalised Hubble constant. This is to be compared with the VLBI estimate Ω_GW < 0.11   h^-2 for ν < 2 × 10^-9   Hz of Gwinn et al. (1997). The same approach applied to the VLBI catalogue used by Titov et al. (2011) gives Ω_GW < 0.009   h^-2 for ν < 1.5 × 10^-9   Hz. Therefore, one can expect that Gaia will improve current estimates from VLBI data by two orders of magnitude while covering a larger interval of frequencies. However, this conclusion relies on the ability to see large-scale features of small amplitude in the Gaia proper motion thanks to the improvement of the detection in . This can be hindered by correlated noise in the intermediates results or more likely the various sorts of unmodelled systematic errors, not known today, although every effort is made in the instrument design and manufacturing to keep them low and compatible with the above results. These systematic errors can be expected to influence lower degree VSH coefficients to a greater extent than the higher degree ones. Therefore, for analysing the gravity waves, it is advantageous to use not only the quadrupole harmonics as in Gwinn et al. (1997) and Titov et al. (2011), but also an ensemble of harmonics that includes higher order harmonics (see e.g. Book & Flanagan 2011). Given the number of QSOs in the Gaia catalogue, such a refined approach is feasible and will help establish more reliable limits on the energy flux resulting from primordial gravitational waves.

7.3. Correlations between the VHS coefficients

It is interesting to briefly discuss the correlations between the VSH harmonics resulting from the least squares fits with the simulated catalogue of 700 000 QSOs described above. These correlations were computed form the inverse of the weighted normal matrix from the VHS fit. Depending on the maximal order l_max, the root mean square value of the correlations is 0.04, while a few correlations may attain the level of 0.4. Having such a large number of sources one might expect that the correlations should be extremely small. Indeed, for a catalogue of 700 000 sources homogeneously distributed on the sky, the rms of the correlations between the VSH coefficients is typically 0.0006 and the maximum is at the level of 0.006. We see that the correlations for the realistic QSO catalogue are about 70 times larger. One can find two reasons for the larger correlations: (1) the inhomogeneous distribution of the QSOs on the sky (this can be called the “kinematical inhomogeneity” of the catalogue) and (2) the fact that the accuracies of the QSOs’ proper motions, which are used in the weight matrix, also depend on the position on the sky (this can be called the “dynamical inhomogeneity” of the catalogue). For the QSO catalogue used in this section the dynamical inhomogeneity alone gives correlations with an rms of 0.01 and the maximum of 0.2, while the kinematical inhomogeneity alone leads to the rms of 0.03 and the maximum of 0.35. Generally speaking, correlations lead to higher errors in the fitted parameters. We finally note that the dynamical inhomogeneity will be slightly greater for the real Gaia catalogue than we see in our simulation. This is caused by the fact that the standard Gaia science performance model (ESA 2011; de Bruijne 2012) is smoothed in the galactic longitude. However, our calculations show that this additional inhomogeneity does not significantly change the results given above.

8. Conclusion

In this paper we have shown the relevance of the vector spherical harmonics (VSH) to decomposing, more or less blindly, spherical vector fields frequently encountered in astronomical context. This happens very naturally in comparing of stellar positional or proper motion catalogues or the different solutions produced during the data analysis in global astrometry. We have provided the mathematical background and considered several practical issues related to the actual computation of the VSH.

We have in particular:

• shown how to perform the decomposition with a least squaresminimisation, enabling processing of irregularly distributed setsof points on the sphere;

• provided a way to test the significance level of any degree in the expansion and then set a criterion for where to stop the expansion;

• stressed the importance of the invariance properties against rotation of the expansion and given the necessary formulas that express the decompositions in the most common astronomical frames;

• provided applications with a reanalysis of the FK5 catalogue compared to Hipparcos, and an example of how this tool should be used during the Gaia data processing to prepare the construction of the inertial frame and derive important physical parameters;

• compiled an extensive set of practical expressions and properties in the appendices.

Acknowledgments

F.M. thanks the laboratory Cassiopée for supporting his stay at the Lohrmann Observatory, Dresden, during the early phase of this work. S.K. was partially supported by the BMWi grant 50 QG 0901 awarded by the Deutsches Zentrum für Luft- und Raumfahrt e.V. (DLR).

References

Appendix A: Explicit formulas for the vector spherical functions with l  ≤  4

Appendix B: Practical numerical algorithm for the scalar and vector spherical harmonics

B.1. Scalar spherical harmonics

Scalar spherical functions Y_lm can be computed directly using definitions (9)–(10). The Legendre functions can be evaluated numerically with the following stable recurrence relation on (B.1)starting with The algorithm based on (B.1)–(B.3) is discussed in Sect. 6.8 of Press et al. (1992). The algorithm described there aims at computing a single value of P_lm(x) for given values of l ≥ 0, 0 ≤ m ≤ l and |x| ≤ 1. It is trivial to generalise it to compute and store all the values of P_lm(x) for l ≤ l_max, l_max ≥ 0, and a given x. The algorithm takes the form \arraycolsep1.75ptHere l_max ≥ 0 is the maximal value of l to be used in the computation. The notation “a = ...,  b = ...   ” denotes outer cycle for b and inner cycle for a with the specified boundaries. As a result one gets a table of values of P_lm(x) for a specified x (|x| ≤ 1) and for all l and m such that 0 ≤ m ≤ l and 0 ≤ l ≤ l_max. An extension of this and related algorithms useful for very high orders l is given by Fukushima (2012).

B.2. Vector spherical harmonics

To compute vector spherical functions as defined by (7), (8) one also needs to compute derivatives of the associated Legendre functions. A closed-form expression for the derivatives reads as (B.8)Special care should be taken for the computations with δ close to  ± π/2. Indeed, for δ =  ± π/2, one has x = sinδ =  ± 1, and both the factors 1/cosδ in (7), (8) and the derivative go to infinity. Therefore, numerical computations (in particular, Eq. (B.8)) become unstable for δ close to  ± π/2. To avoid this degeneracy the definitions of T_lm and S_lm can be rewritten as One can easily see that (e.g. from (11)) that A_lm(x) and B_lm(x) remain regular for any 0 ≤ m ≤ l and |x| ≤ 1. It is easy to see that numerically stable algorithm for A_lm and B_lm read as (B.19)\arraycolsep1.75ptAs a result tables of A_lm(x) and B_lm(x) are computed for a given |x| ≤ 1, 0 ≤ m ≤ l and 1 ≤ l ≤ l_max. Explicit expressions for the vector spherical harmonics for 1 ≤ l ≤ 4 are given in Appendix A.

Appendix C: VSH expansion of a vector field vs. the scalar expansions of its components

As explained earlier with the VSH, we can expand a vector field defined on a sphere on vectorial basis functions preserving the vectorial nature of the field and behaving very nicely under space rotations. However, based on the usual practice in geodesy and gravitation, following Brosche (1966) and Brosche (1970) astronomers have also used, for many years, the expansions of each component of the vector field, namely V^α = V·e_α and V^δ = V·e_δ in the scalar spherical harmonics. It is interesting to show the relationship between the two expansions and how their respective coefficients are related to each other.

C.1. Definition of the expansions

Let us consider a vector field V on the surface of a sphere. Its VSH expansion reads as (C.1)On the other hand, the components of this vector field V^α,V^δ are also expandable independently of each other in terms of the usual scalar spherical functions Y_lm: or for the vector field itself: (C.4)As with all our infinite sums, the equalities (C.1)–(C.4) only hold in the sense of convergence in the L² norm as already stated in Sect. 2.2. Pointwise convergence is not guaranteed in the general case.

Therefore, one has the identity (C.5)

C.2. Formal relations between the coefficients

Multiplying both sides of (C.5) by or (again, superscript “ ∗ ” denotes complex conjugation) and integrating over the surface of the sphere, one gets where where as usual dΩ = cosδ   dδ   dα, and the integration is taken over the surface of the unit sphere: 0 ≤ α ≤ 2π, −π/2 ≤ δ ≤ π/2.

On the other hand, multiplying both sides of (C.5) by and integrating over the surface of the sphere one gets It remains to compute A_l′m′lm and B_l′m′lm in a convenient way. But formally we have obtained the two-way correspondence between the coefficients (t_lm,s_lm) and .

C.3. Explicit formulas for A_lml′m′ and B_lml′m′

It is straightforward to show that where δ^ij is the Kronecker symbol (δ^ij = 1 for i = j and δ^ij = 0 otherwise) and k ∈ Z is arbitrary integer. Both α_lmk and β_lmk are positive rational numbers provided that the indices l, m, and k are selected in such a way that the Legendre polynomials under the integral are different from zero: |m| ≤ l and |m| ≤ l + 2k + 1 for α_lmk, and |m| ≤ l and |m| ≤ l + 2k for β_lmk. Otherwise α_lmk and β_lmk are zero. Numbers α_lmk and β_lmk can be computed directly or by using recurrence formulas that can be obtained from the well-known recurrence formulas for the associated Legendre polynomials. The following relations allow one to consider only the case of β_lmk with positive indices: Finally, a number of equivalent formulas for β_lmk valid for any l ≥ 0, 0 ≤ m ≤ l and k ≥ 0 can be derived. Thus, using that any associated Legendre polynomial can be represented as a finite hypergeometric polynomial and integrating term by term, one gets the following formula for β_lmk valid for any l ≥ 0, 0 ≤ m ≤ l and k ≥ 0: We note that a number of explicit formulas for β_lmk can be found for l, m, and k satisfying some specific conditions (e.g. m = l). For the general case, an alternative formula for β_lmk can be derived using the Gaunt formula (Gaunt 1929, Appendix, pp. 192−196). However, Eq. (C.22) is sufficient for practical computations. These formulas can be easily implemented numerically. Equation (C.22) involves a double sum of terms of alternating signs. This means that computational instabilities can appear if floating point arithmetic is used.

C.4. Final relations between the coefficients

Since A_lml′m′ and B_lml′m′ vanish unless certain relations between the indices hold, one can significantly simplify the transformations (C.6), (C.7), and (C.10), (C.11) can be significantly simplified and represented as a single sum over integer k: In the above formulas we simplified the notations of the lowest values of k in the sums and assume now that , , t_0m = 0, and s_0m = 0 and that the terms containing them vanish identically.

Thus, we have proved that the information of each particular t_lm and s_lm is distributed over infinite number of and and vice versa.

C.5. Relation to space rotations

It is important here to draw attention to two major differences between the expansions (C.1) and (C.2), (C.3).

• A rotation between two catalogues is very simply represented in the vectorial expansion with the harmonic T_1m and the three coefficients t_1,0,   t_1, −1,   t_1,1, while from (C.10), (C.11), this will require an infinite number of coefficients in the component representations. One may rightly argue that the opposite is equally true: a simple scalar decomposition only over, say, l = 1 will be much more complex in the vectorial expansion. This is true, but physically the really meaningful global effects are precisely the rotation and the glide, which generate very simple vector fields and project on VSHs of first degree only and requires many degrees if done with the scalar components. We know of nothing equivalent to the scalar representation, unless, obviously, one builds an ad-hoc field for the purpose of illustration.

• The behaviour of the VSHs or of the scalar spherical harmonics Y_lm under space rotation are very similar and show the same global invariance within a given degree l. Their transformations are fully defined with the Wigner matrix as shown in Sect. 3. Therefore it seems at first glance that the transformations under space rotation of the t_lm and s_lm in (C.1) are not simpler than that of the and in (C.2), (C.3). There is, however, a very important difference that makes the use of the VSH so valuable. By applying the Wigner matrix associated to a space rotation to (C.1), the new coefficients correspond to the expansion of the vector field with its components given in the rotated frame. Now in (C.10), (C.11), the two scalar fields are considered independently of each other, and the new coefficients deduced from the application of the Wigner operator correspond to the expansion of the initial scalar fields (V^α,V^δ) expressed in the rotated fields, but these components are not those of the field V projected on the rotated coordinates, and they are still the components in the initial frame, since a scalar field is invariant by rotation.

All Tables

All Figures

	Fig. 1 Local frame associated to the spherical coordinates α,δ, with the unit vectors along the longitude and latitude lines.
In the text

	Fig. 2 Examples of spheroidal harmonics vector fields: S1,1, S2,1, S5,3, S10,5, where one sees the change in the angular resolution with higher degrees. The S1,1 is a glide flow in the direction of the y-axis.
In the text

	Fig. 3 Examples of toroidal harmonics vector fields: T1,1, T2,1, T5,3, T10,5. One sees the change in the angular resolution with higher degrees. The T1,1 is a simple rotation about the y-axis.
In the text

	Fig. 4 Rotation and glide with the z-axis shown here with a Mercator projection to emphasise the close relationship between the two vector fields. The constant displacement in longitude from the rotation has its counterpart with a constant glide displacement in the reduced latitude (Mercator latitude).
In the text

	Fig. 5 Differences in positions and proper motions between the FK5 catalogue and Hipparcos in 1991.25. The typical arrow sizes in the plots are respectively 100 mas and 2.5 mas/yr.
In the text

	Fig. 6 True errors and formal errors of the estimates of the absolute values of the rotation and acceleration from the fits with 1 ≤ lmax ≤ 15 for the simulated QSO catalogue matching expected Gaia accuracy.
In the text