The wave surveyor technique for fast plasma wave detection in multi-spacecraft data

Multi-satellite missions like Cluster allow to study the full spatio-temporal variability of plasma processes in near-Earth space, and both the frequency and the wave vector dependence of dispersion relations can be reconstructed. Existing wave analysis methods include high-resolution beamformers like the wave telescope or k-filtering technique, and the phase differencing approach that combines the correlations measured at pairs of sensors of the spacecraft array. In this paper, we make use of the eigendecomposition of the cross spectral density matrix to construct a direct wave identification method that we choose to call the wave surveyor technique. The analysis scheme extracts only the dominant wave mode but is much faster to apply than existing techniques, hence it is expected to ease survey-type detection of waves in large data sets. The wave surveyor technique is demonstrated by means of synthetic data, and is also applied to Cluster magnetometer measurements.


Introduction
Near-Earth space is a dynamic plasma environment that creates and supports wave activity on a broad range of temporal and spatial scales.The inherent ambiguity of singlespacecraft data makes it difficult to identify waves as, e.g.Doppler shifts may significantly affect the frequency determination.Multi-satellite missions can overcome this problem.
Estimation of wave vectors k from such multipoint measurements, however, is not as straightforward as a Fourier transformation because of the small number of sensors in the Correspondence to: J. Vogt (j.vogt@jacobs-university.de) spacecraft array.In the context of the Cluster mission, several approaches to the problem have been presented.Dunlop et al. (1988) and Neubauer and Glassmeier (1990) introduced the term wave telescope for a method based on a linear filter bank approach, and quantified the spatial aliasing condition in terms of the reciprocal lattice of the spacecraft tetrahedron.The k-filtering technique constructed by Pincon and co-workers (e.g.Pincon andLefeuvre, 1991, 1992;Pincon and Motschmann, 1998) by means of a minimization principle is based on an estimator for the spatio-temporal power spectrum P (ω, k).Sensor weights are chosen such that the contribution of plane waves with wave vectors k outside a small spectral window around k to the resulting spectral energy density estimator is minimum.Such techniques were originally developed for seismic arrays (e.g.Capon et al., 1967;Capon, 1969;Cox, 1973), and are commonly referred to as Capon estimators, minimum variance estimators, high resolution beamformers or simply beamformers.Minimum variance estimators have been used to identify MHD waves in the magnetosheath and the foreshock region (Glassmeier et al., 2001;Narita et al., 2003;Narita and Glassmeier, 2005) and to study turbulence (Sahraoui et al., 2003(Sahraoui et al., , 2006;;Narita et al., 2006).Furthermore, Constantinescu et al. (2007) constructed a wave detection scheme based on spherical waves instead of plane waves to identify not only wave vectors but also the location of the wave source.
The term phase differencing approach refers to a class of wave analysis techniques where projections of the wave vector k onto the spacecraft separation vectors are estimated from phase differences of the signal measured between the corresponding pairs of sensors.If four or more point measurements are available, the full wave vector can be reconstructed from the projections.In the case of three sensors or less, physical constraints such as ∇•B=0 can be taken into account to partially make up for the missing information.The phase differences can be estimated, e.g. using spectral correlation measures based on Fourier transformations of Published by Copernicus Publications on behalf of the European Geosciences Union.the observed signals.In the preparation phase of the Cluster mission, the phase differencing approach was presented by Balikhin and Gedalin (1993).Dudok de Wit et al. (1995) introduced a technique based on the Morlet wavelet transform, and used AMPTE-UKS and AMPTE-IRM data to demonstrate that several waves at the same frequency can be identified.Matsui et al. (2007) applied a phase differencing technique to Cluster observations to study broadband ULF waves near the dayside polar cap boundary.For a detailed comparison of the phase differencing approach with the k-filtering technique using Cluster STAFF and EFW measurements, the reader is referred to Walker et al. (2004).For a review of the wave distribution determination problem, see Storey (1999).
A key analysis step in the application of the minimum variance estimators mentioned above is a peak search in threedimensional k space for each wave that may be present at a given frequency.Phase differencing schemes require peak finding in spectral cross correlation measures between various pairs of sensors in the array.Such search procedures can be quite time-consuming and also ambiguous in some cases.In this paper we propose a wave detection method that we choose to call the "wave surveyor technique".It allows to compute the wave vector and the polarization vector as a function of frequency directly from the data.At a given frequency, the method works only for the dominant wave mode whereas minimum variance estimators and phase differencing techniques can in principle identify a number of different modes.As the name suggests, the wave surveyor should be a useful tool for survey-type screening of large data sets for waves and wave parameters.
The wave surveyor technique makes use of the eigendecomposition of the cross spectral density (CSD) matrix that also plays a key role in the so-called multiple signal classification (MUSIC) scheme (Schmidt, 1979(Schmidt, , 1981)).For a concise introduction to the subject of array signal processing, the reader is referred to Pillai (1989) where minimum variance estimators as well as the MUSIC scheme and several other approaches are discussed.In the space physics context, methods based on the eigendecomposition of the CSD matrix have been widely used by Samson and co-workers (e.g. Samson and Olsen, 1980;Samson, 1983;Samson et al., 1990), e.g. to evaluate the significance of analysis results, and to yield general polarization measures.Santolík et al. (2003) carried out singular value decompositions (SVDs) of magnetic and electromagnetic spectral matrices to identify and analyze plasma waves in the auroral region, and provided a more physical interpretation of the eigenstructure of the CSD matrix.
This paper is organized as follows.The terminology, key variables and identities are introduced in Sect. 2. In Sect.3, the wave surveyor technique is constructed, and demonstrated by means of synthetic signals in Sect. 4. The new technique is applied to Cluster magnetometer data in Sect. 5. Advantages and limitations are discussed in Sect.6.We conclude in Sect.7.

Notation and key variables
In this paper vectors a, b, c, . . .are always understood as column vectors.Unit vectors are indicated by •, for example, â or b.Superscripts t, * , † denote transpose, complex conjugate, and hermitian adjoint, respectively.Accordingly, a t and a † are row vectors, the dot product of a and b is a•b=a t b, and the hermitian product is a † b.Matrices are typeset in sans serif font.The symbol E is used to denote identity matrices (of various dimensions).• • • stands for mathematical expectation which in practice is approximated through an averaging procedure.

Data representation and cross spectral density matrix
We consider vector time series B σ (t) with J components B j σ (t), j =1, . . ., J , measured at S points in space r σ , σ =1, . . ., S. If we consider CLUSTER magnetometer data, then J =3 and S=4.Let b σ (ω) denote the respective Fourier transforms which for continuous functions are defined through (1) In the more relevant case of observations taken at discrete times and over a finite measurement interval, we write where the constant depends on the chosen implementation of the Fourier transform (for details see Eriksson, 1998).Complex data vectors with L=J •S components can be formed through The matrix comprises all possible covariances of the Fourier transforms, and is at the heart of the wave analysis techniques considered here.For brevity, we refer to C as the FSC (Fourier Space Covariance) matrix.The matrix is hermitian and can be diagonalized.The eigenvalues γ , =1, . . ., L are real and non-negative, and the (normalized) eigenvectors are ĉ .Just as the FSC matrix, its eigenvalues and eigenvectors depend on ω but not on k.The eigenvalues are assumed to be in descending order.Of particular importance are the largest eigenvalue γ 1 and the associated eigenvector ĉ1 .Since they appear in many formulas below, we will most often drop the subscript and write γ for γ 1 as well as ĉ for ĉ1 .The FSC matrix C differes from the cross spectral density matrix M only by a constant scalar factor M 0 : which implies that both matrices share the same eigenvectors, and the eigenvalues µ of M are related to the eigenvalues of C through µ =M 0 γ .The constant factor M 0 depends on the implementation of the Fourier transform.For notational convenience, we choose to develop the wave surveyor formalism on the basis of the FSC matrix, and express the results also in terms of the eigenvalues µ of the CSD matrix when required.

The FSC matrix of the plane wave model
We intend to construct a direct technique to detect a plane wave in multi-spacecraft data, and to estimate the wave parameters such as the wave vector k and the polarization vector a as functions of (angular) frequency ω: k=k(ω), and thus also a(ω, k)=a(ω, k(ω))=a(ω).
In general, an individual Fourier component gives rises to a model signal that varies in time t and space r as which means that the Fourier transform of the model signal with respect to time only can be written as The signal is measured in space at r σ , σ =1, . . ., S to give The second term on the right-hand side is the mismatch of the model and the data, and is modeled as isotropic and homogeneous noise of variance η 2 .Forming the complex data vector b as described in the previous Sect.2.1 yields where the wave vector k=k(ω) is understood as a unique function of the (angular) frequency ω as explained above, hence b can be written as a function of ω only.The (L×J ) matrix encodes the array geometry, and E is the (J ×J ) identity matrix.
The FSC matrix of this model is given by where N=η 2 E for isotropic noise.
Since all other eigenvalues are simply η 2 and thus smaller, the first eigenvector ĉ1 ≡ ĉ is proportional to Ha, or, more precisely, and the other eigenvectors are orthogonal to Ha.

Scalar data and projection operators
The individual components of vector time series are scalar time series.In Sect.3, the wave surveyor technique is constructed first for the scalar case, and then formulated for the general case of vector-valued time series.The correspondence of the scalar and the vector technique can be conveniently quantified using the operators j : C L →C S (projection, note that L=J •S) and I j : C S →C L (injection) defined below.
Scalar time series measured at S points in space r σ , σ =1, . . ., S are written as B σ (t).The Fourier transforms b σ (ω) can be assembled into a complex data vector with S components: As before, the FSC matrix is defined through and differs from the cross spectral density matrix M only by a constant factor M 0 also in the scalar case.The complex vector function encodes the geometry of the array.Normalization yields ĥ(k)=h(k)/ √ S.

J. Vogt et al.: Wave surveyor technique
The projection operator j is defined through i.e. in matrix notation, Here 0 and êj are the zero vector and the j -th unit base vector in C J , respectively.The injection operator I j is the transpose of j , i.e.
It is easy to verify that for all b∈C L we have Hence J j =1 I j j ≡ J j =1 j † j is the identity E on C L .Furthermore, and thus for all eigenvectors ĉ .

The wave surveyor technique
In this section we derive the wave surveyor technique.As explained already, the wave surveyor is a direct wave identification and dispersion analysis technique in the sense that the wave vector is computed directly as a function of the angular frequency, i.e. k=k(ω), and a peak search in discretized three-dimensional wave vector space is not required.We first look at the case of scalar data, and then generalize the ideas to vector-valued time series.

The wave surveyor technique for scalar data
The construction of the wave surveyor technique is guided by the properties of the single plane wave model presented in Sect.2.2.In the case of scalar data, J =1, H=h, and the FSC matrix reads where a=a(ω) and h=h(k).The largest eigenvalue of the FSC matrix is γ =S |a(ω)| 2 +η 2 , and the first eigenvector is given by ĉ= ĥ(k)≡h(k)/ √ S. Hence the signal amplitude |a| can be determined from where the noise parameter η 2 can be estimated from the remaining eigenvalues γ , ≥2.Alternatively, if the eigenvalues of the CSD matrix M are to be used, the signal amplitude can be expressed as where µ=µ 1 =M 0 γ is the largest eigenvalue of M, and η 2 M =M 0 η 2 is estimated from the smaller eigenvalues µ , ≥2.
Since the eigenvector ĉ is proportional to the vector h(k) evaluated at the actual wave vector k of the signal, and k is part of the arguments of the complex exponentials in h, we expect that the wave vector can be estimated from the phases θ σ =θ σ (ω) of the eigenvector components ĉ1,σ =| ĉ1,σ | exp(iθ σ ).A component-wise comparison of the eigenvectors and the vector h(k) suggests that the phases θ σ should deviate from the expressions k•r σ by a constant phase delay φ only, and thus should minimize the cost function with respect to k and φ.
In the Appendix it is shown that the solution for k can be written as follows: Here the positions r σ are relative to the center of the sensor array which implies that σ r σ =0.
If S=4 as is the case for the Cluster mission, the solution can be explicitly given in terms of the reciprocal vectors κ σ of the spacecraft tetrahedron (for details of the reciprocal vector concept see, e.g.Chanteur, 1998).As demonstrated in the Appendix, the wave vector can be expressed in terms of the eigenvector phases and the reciprocal vectors as follows: Since the θ σ are determined directly from the FSC matrix C(ω), and the reciprocal vector κ σ are functions of the array geometry only, Eq. ( 28) can be directly evaluated to yield the wave vector in survey-type wave analyses of large scalar data sets.

The wave surveyor technique for vector data
In the case of vector data, the FSC matrix of the single plane wave model presented in Sect.2.2 is given by Note that now, both, the amplitude (polarization) vector a and the wave vector k enter the eigenvector ĉ ∝ H(k)a.
We apply the projection operators j (introduced in Sect.2.3) to Eq. ( 13) and note that j H=h êj † to write As in the scalar case, the vector h(k) (i.e.evaluated at the actual wave vector) can be written in terms of the first eigenvector.In fact, we now have a total of J such relationships that lead to J scalar cost functions, and we can combine them to estimate the wave vector k from the phases of the components of the first eigenvector.The weights of the partial cost functions are chosen to be proportional to |a j | 2 , i.e. to the square of the j -th component of the amplitude vector a, and this component is proportional to j ĉ as can be seen from the relationship given above.Hence the (total) cost function can be written as where θ j σ denotes the phase of the σ component of the projection j ĉ, and α j =| j ĉ| 2 /| ĉ| 2 , hence j α j =1.Minimizing the cost function works as for scalar data, and for the special case of the Cluster tetrahedron and FGM data (S=4, J =3) we finally obtain the wave surveyor estimate of the wave vector as where κ σ are the reciprocal vectors of the tetrahedron as before.For the general case of S sensors, we can write Equation ( 13) also allows to construct an estimator for the polarization vector a.We apply the operator H † to Eq. ( 13) and note that H † H=SE to obtain Since the CSD matrix M has the same eigenvectors as the FSC matrix, and the eigenvalues are related through µ =M 0 γ , the amplitude vector may also be expressed as where µ=µ 1 =M 0 γ is the largest eigenvalue of M, and η 2 M =M 0 η 2 is estimated from the smaller eigenvalues µ , ≥2.Equations ( 32), (34), and (35) allow to compute the wave vector k and the polarization vector a directly from the eigendecomposition of the FSC or the CSD matrix.If measurements from more than four sensors are available, Eq. ( 33) can be used instead of Eq. ( 32).The wave surveyor techniques does not require a peak search in the three-dimensional wave vector space.

Demonstration of the wave surveyor technique
We now demonstrate the wave surveyor technique by means of a synthetic model signal composed of two plane waves and isotropic noise: The time lag t n is a function of position, the term N represents white noise with zero mean and unit variance, and the coherent part of the signal consists of two harmonic (cosine) wave trains in a Gaussian envelope: Note that the amplitude spectrum of the coherent part is determined completely by the amplitudes A n and the model signals W T n ,τ n (t), and does not depend on the position r.
The model parameters are the periods T n of the two plane waves, the slowness vectors u n , the amplitude vectors A n , the widths τ n of the Gaussian envelope function, and the noise amplitude ν.The wave parameter values used here are summarized in Table 1.The noise amplitude is set to the value ν=0.2.
The synthetic signal is assumed to be sampled in space at four locations, namely, at the origin of the cartesian coordinate system, and at three points on the coordinate axes, each one at a distance of 200 km from the origin.The upper panel of Fig. 1 shows the generated signal at the sampling point in the origin.In the lower panel the amplitude spectrum of the coherent (noise-free) part of the model signal is displayed.As noted above, this amplitude spectrum does not depend on r, and is thus identical at all sensor locations.The two coherent contributions to the model signal are located in two frequency bands that are well separated from each other.www.ann-geophys.net/26/1699/2008/The square modulus of the Fourier amplitude, i.e. |b(ω)| 2 in our case, is a measure of the total signal power in the frequency domain.In minimum variance estimators like the wave telescope or the k-filtering method, such a power spectrum estimate is used to identify the frequency bands with sufficient power to support waves.Since this approach is equivalent to using the trace of the FSC matrix for inspecting the frequency domain.In the case of the CSD matrix M=M 0 C, its trace gives the total power spectral density.
In its principal axes system, the eigenvalues µ =µ (ω) of the CSD matrix reside on the diagonal, hence trace M= µ .As explained in Sect.2.2 by means of the plane wave model, the wave signature shows up in the first mode, whereas the noisy part contributes equally to all eigenvalues.Therefore, the eigenvalues of the CSD matrix effectively decompose the signal power into a number of modes of decreasing significance.
For the synthetic signal considered here, the eigenvalues and the trace of the CSD matrix are shown as functions of the frequency f =ω/2π in Fig. 2. The peaks associated with the waves can be seen in both the trace of the CSD matrix and in its first eigenvalue (the remaining eigenvalues collect the contribution of the noisy part of the signal), however, the peaks in the first eigenvalue stand above the noise background much more clearly than the peaks in the trace.In this sense, eigenstructure based methods like the wave surveyor technique can yield a better separation of the Fourier modes and the noise background in the frequency domain.
After peaks in the frequency domain are identified, existing multi-spacecraft wave analysis techniques have to discretize the three-dimensional k-space, compute a power spectrum estimator P (ω, k) on the resulting grid of wave vectors at least for the frequencies ω of interest, and then carry out a peak search in k-space.The wave surveyor technique allows to work out the wave parameters in a much more direct way by means of the explicit formulas ( 32) and (34): the slowness vector u, the wave vector k=ωu, and the Ann.Geophys., 26,[1699][1700][1701][1702][1703][1704][1705][1706][1707][1708][1709][1710]2008 www.ann-geophys.net/26/1699/2008/amplitude (polarization) vector a are computed directly as functions of frequency.
For the synthetic signal considered here, the components of the estimated slowness vector as functions of frequency f =ω/2π are shown in Fig. 3.At each frequency, the area of the diamond represents the power contained in the respective dominant mode as given by the largest eigenvalue.The sizes of the symbols are meant to serve as significance measures of the slowness vector estimates in the following way.In the frequency range outside the two bands supported by the plane waves where only the noise term contributes to the signal, there is relatively little power in the dominant modes, and the wave vector estimates cannot be considered meaningful.In the two frequency bands with significant wave power, the symbol sizes are larger, there is little scatter, and the results compare nicely with the parameters of the synthetic signals.

Application to Cluster FGM observations of foreshock waves
In this section we show results of an application of the wave surveyor technique to magnetic field fluctuations recorded by the fluxgate magnetometer on board the four Cluster spacecraft (Balogh et al., 2001).These observations were analysed already by Narita et al. (2007) using the well-established and thoroughly tested wave telescope analysis method (e.g.Neubauer and Glassmeier, 1990;Motschmann et al., 1996;Glassmeier et al., 2001;Narita et al., 2003) to study the dispersion of foreshock waves.We may thus validate the wave surveyor approach by comparing our results with the findings of Narita et al. (2007).
The analysis example makes use only of one magnetic field component, namely, the B z (northward) component in the GSE coordinate system.We thus follow the procedure outlined in Sect.3.1 for scalar data, see Eq. ( 28).The time interval of interest is 16 February 2002, 07:00-07:45 UT, and it comprises Cluster observations of a representative case of foreshock waves.Narita et al. (2007) identified the whistler wave dispersion branch and demonstrated how it becomes Alfvén wave dispersion (ω=kV A ) at small wave numbers.Background plasma and magnetic field values were as follows: the mean magnetic field was pointing away from the sun (B x =−5.6 nT, B y =−1.4 nT, B z =−1.4 nT in the GSE coordinate system), the plasma bulk velocity was almost 300 km/s (V x =−300.7 km/s, V y =24.3 km/s, V z =2.8 km/s), and the ion density had the value n=5.9 cm −3 .The plasma velocity and density were provided by the ion measurements of the Cluster CIS-HIA instrument (Rème et al., 2001).
The determination of the wave vectors further allows to transform the wave frequencies from the spacecraft frame (ω sc ) into the plasma rest frame (hereafter, the rest frame, ω re ), a frame which is co-moving with the plasma bulk ve-Fig.3. Demonstration of the wave surveyor technique.The components of the slowness vector u are given as functions of the frequency f =ω/2π .For each frequency, the area occupied by the plotting symbol is a measure of the signal power given by the largest eigenvalue.Only the frequency bands associated with the coherent part of the model signal yield significant power.The smaller diamonds that show much scatter are associated with the contribution of the noise term to the model signal.
locity.This transformation is carried out using the Doppler relation: where V =(V x , V y , V z ) t denotes the plasma bulk velocity given above.Figure 4 displays the dispersion relation, ω re =ω re (|k|), and the propagation angle with respect to the mean magnetic field direction, θ kb , derived by the wave surveyor technique.The dispersion relation exhibits a phase speed close to the Alfvén speed, (ω re kV A ) for wave numbers smaller than the ion inertial wave number, k in = i /V A =0.011 rad/km (here i =0.64 rad/s is the ion cyclotron frequency and V A =59.7 km/s is the Alfvén speed).However, for larger wave numbers it starts to deviate from the Alfvén wave branch toward higher frequencies (ω re >kV A ).This is characteristic to the low frequency part of the whistler mode dispersion.The propagation direction is almost anti-parallel to the mean magnetic field, therefore the waves propagate intrinsically away from the bow shock.
In conclusion, the results obtained through the wave surveyor technique are fully consistent with the findings of Narita et al. (2007) using the wave telescope analysis.

Discussion
The wave surveyor technique is a direct method to estimate the parameters of a dominant plane wave in multipoint measurements.Isotropic white noise has no influence on the eigenvectors of the CSD matrix and hence does not change the estimated wave vectors or amplitude vectors.The presence of other waves at the same frequency, however, may limit the applicability of the analysis method.Since their contributions to the total variance affect the eigenvalue distribution of the CSD matrix, the eigenvalue ratios may be Fig. 5. Eigenvalues of the CSD matrix for the analysis example presented in Sect. 5.The first eigenvalue is clearly much larger than the other eigenvalues throughout the whole frequency range which confirms that we are dealing with a single dominant mode in this case.
used to check the validity of the model assumptions.The single (dominant) plane wave model is expected to provide an appropriate characterization of the measured signal if the first eigenvalue proves to be much larger than the remaining ones.
The distribution of eigenvalues with frequency for the foreshock wave analysis event of Sect. 5 is shown in Fig. 5. Throughout the entire frequency range, the first eigenvalue is about three orders of magnitude larger than the other ones.Hence it is indeed quite safe to assume that the event is well characterized by the single (dominant) plane wave model.The practical significance of the eigenvalue ratios for the robustness of the parameter estimation is shown also in Fig. 3 for the case of synthetic data: the frequency ranges where the first eigenvalues are large (corresponding to large plotting symbols) yield stable parameter estimates whereas the frequency ranges where noise dominates (small plotting symbols) exhibit a lot of scatter.
A different kind of quality indicator for the wave surveyor parameter estimation is suggested by the construction principles discussed in Sects.3.1 and 3.2.The single plane wave model presented in Sect.2.2 implies that the (estimated) vector ĥ(k) coincides with the (observed) first eigenvector ĉ in the scalar case, or with properly normalized versions of the vectors j ĉ, j =1, . . ., J in the case of vector data.In fact, the cost functions ( 26) and ( 31) are quadratic measures of the angular mismatch between these sets of vectors, corrected for a possible constant phase offset φ, and may thus serve as quality indicators for wave parameter estimates.Since this paper was meant to introduce the wave surveyor technique and to provide a proof of concept, we did not try to quantify threshold values neither for the angular mismatch quality indicator, nor for the eigenvalue ratios discussed above, and leave this issue for future studies.
Wave analysis techniques based on the cross spectral density matrix implicitly concentrate on second order moments.In order to address more complex associations in multipoint measurements, the wave surveyor technique could be generalized by means of a singular value decomposition (SVD) applied to the L×N data matrix b(ω) to yield Here N is the number of ensembles (subintervals in time), diag( √ γ ) denotes the diagonal matrix with elements √ γ , and the columns of U are identical with the eigenvectors of the FSC matrix.The matrix V, however, provides new information: it allows to address the variability in time and to test for stationarity.
We conclude this section with a few comments on how the eigenstructure decomposition of the FSC matrix can help to address selected aspects of the two main classes of wave identification methods discussed in the introduction.Provided that we are dealing with a signal that can be described by the single (dominant) plane wave model, the linear parts of the sensor pair correlations that constitute the basis of the phase differencing approach are implicitly encoded in the phases of the first eigenvector.This follows from the relations h(k)∝ ĉ1 for the scalar case (see Sect. 3.1), h(k)∝ j ĉ1 for the case of vector data (see Sect. 3.2), and the definition of the vector function h(k).As demonstrated by Dudok de Wit et al. (1995) for the case of two-point measurements, the phase differencing approach is suited to identify several waves at the same frequency whereas the wave surveyor technique extracts the dominant wave only.Using four sensors instead of two allows to improve the effective signal-to-noise ratio, and in this case the inversion of the position tensor and thus the wave vector estimation in the wave surveyor technique can be carried out directly (see the Appendix) and with little effort.In the case of four-point phase differencing method, the improved signal-to-noise ratio goes along with a more involved reconstruction scheme (e.g.Matsui et al., 2007).
To gain additional insight into the performance of minimum variance estimators, it is instructive to rewrite them using the eigendecomposition of the FSC matrix: Since we obtain which means that the ĉ , =1, . . ., L, are eigenvectors also of the inverse matrix C −1 , and the corresponding eigenvalues are γ −1 .Therefore, we can write For brevity, we consider the scalar case only.Hence L→S, →σ , and the minimum variance estimator for the power spectral density estimate is given by P =(h † C −1 h) −1 (e.g.Motschmann et al., 1995).The eigenstructure representation of C −1 allows to rewrite P =P (ω, k) as follows: therefore, For the single plane wave model, h(k)= √ S ĉ1 which implies that h(k) ⊥ ĉσ for σ =1, or, equivalently, h † (k) ĉσ =0 for the remaining eigenvectors ĉσ , σ =2, . . ., S. Furthermore, γ −1 1 =(S|a| 2 +η 2 ) −1 , and the result is thus This has to be compared with the minimum value of P in the case when h(k ) lies in the noise subspace, i.e. the subspace spanned by the eigenvectors ĉσ , σ =2, . . ., S.Here we find P (ω, k )=η 2 /S and thus Hence the resolving power of the scalar minimum variance estimator measured by this analytical expression increases (quadratically) with the signal-to-noise ratio (as expected).
In practice, however, the numerical inversion of the CSD matrix C may cause problems if C is near singular which happens, e.g. in the case of a very large signal-to-noise ratio (|a| 2 η 2 /S) in the single plane wave model.In such a case one might be tempted to perform an inversion in the singular value sense and disregard the contributions of the smallest eigenvalues to obtain and P SV =γ 1 |h † ĉ1 | −2 .However, this would make the method completely useless.Although P SV (ω, k) gives the correct (and finite) result if k is the actual wave vector of the single plane wave model and thus h(k)=S ĉ1 , the values P SV (ω, k )

Fig. 1 .
Fig. 1.Demonstration of the wave surveyor technique.Upper panel: Synthetic time series sampled at one of the four points in space.Lower panel: Amplitude spectrum of the coherent part of the model signal as a function of frequency f =ω/2π .

Fig. 2 .
Fig. 2. Demonstration of the wave surveyor technique.Eigenvalues and trace of the CSD matrix as functions of frequency f =ω/2π .First eigenvalue: solid line.Remaining eigenvalues: dotted lines.Trace of the CSD matrix: dashed line.

Fig. 4 .
Fig. 4. Experimental dispersion relation of the foreshock waves (top) and propagation angles from the mean magnetic field direction identified by the wave surveyor technique using Cluster magnetic field data.The frequencies are represented in the plasma rest frame.The ion inertial scale and the ion cyclotron frequency are k in =0.011 rad/km and i =0.64 rad/s, respectively.

Table 1 .
Wave parameter values for the synthetic signal consisting of two planes waves in an isotropic noise background, see Eq. (36).Model parameters are the wave periods T n , the slowness vectors u n , the amplitude vectors A n , and the widths τ n of the Gaussian envelope function.The wave frequencies f n =1/T n and the wavelengths λ n =T n /|u n | are added for convenience.