Annales Geophysicae High-resolution wave number spectrum using multi-point measurements in space – the Multi-point Signal Resonator ( MSR ) technique

A new analysis method is presented that provides a high-resolution power spectrum in a broad wave number domain based on multi-point measurements. The analysis technique is referred to as the Multi-point Signal Resonator (MSR) and it benefits from Capon’s minimum variance method for obtaining the proper power spectral density of the signal as well as the MUSIC algorithm (Multiple Signal Classification) for considerably reducing the noise part in the spectrum. The mathematical foundation of the analysis method is presented and it is applied to synthetic data as well as Cluster observations of the interplanetary magnetic field. Using the MSR technique for Cluster data we find a wave in the solar wind propagating parallel to the mean magnetic field with relatively small amplitude, which is not identified by the Capon spectrum. The Cluster data analysis shows the potential of the MSR technique for studying waves and turbulence using multi-point measurements.


Introduction
Recent developments of multi-spacecraft measurements have considerably advanced our understanding of space plasma dynamics, as they provide the three-dimensional spatial resolution in the measurements.The wave telescope technique (also named k-filtering) developed by Neubauer and Glassmeier (1990); Pinc ¸on and Lefeuvre (1991); Motschmann Correspondence to: Y. Narita (y.narita@tu-bs.de)et al. (1996), and Glassmeier et al. (2001) is one of the most important applications of multi-spacecraft measurements: it enables us to analyze the data directly in the threedimensional wave vector domain.Based on this technique, various analysis methods have been developed and used for Cluster data, e.g., dispersion relations (Narita et al., 2003;Sahraoui et al., 2003;Vogt et al., 2008;Sahraoui et al., 2010b), phase speed diagrams (Schäfer et al., 2005), energy spectra (Narita et al., 2006;Sahraoui et al., 2006Sahraoui et al., , 2010b;;Narita et al., 2010), magnetic helicity density (Narita et al., 2009), and higher order moments such as bispectra (Narita et al., 2008).In short, various wave properties can be determined directly in the wave vector domain.
The wave telescope technique is based on the minimum variance projection proposed by Capon (1969).While this technique is now widely used for multi-point measurements in space, it also has several limitations about which one should keep in mind.In particular, the wave telescope technique cannot resolve waves properly when their wavelengths are very close to each other; and it gives a moderately high background noise level in the spectrum.Here we present a new analysis method which is a generalization of the wave telescope technique and provides a high-resolution energy spectrum in the wave vector domain.It has the ability to reduce the background level considerably to enhance the quality of the spectrum (signal-to-noise ratio) while maintaining the energy of the signal.Therefore, it has the ability to resolve waves even though they have similar wavelengths.We provide a mathematical description of the new analysis method in Sect.2; it is applied to synthetic data and to real spacecraft data in Sects.3 and 4, respectively.Section 5 concludes the manuscript.

Array signal processing
The new analysis method is called the MSR technique (Multi-point Signal Resonator) and is essentially a combination of Capon's minimum variance estimator (Capon, 1969) and the MUSIC algorithm (Multiple Signal Classification) developed by Schmidt (1986).These two estimators are briefly reviewed, and then the MSR technique is presented.

Capon's minimum variance method
Suppose we have an array of L sensors and each sensor obtains a time series of signal such as a set of waves.While transforming the data from the temporal domain to the frequency domain is performed by Fourier transform such as FFT (Fast Fourier Transform), it is not an easy task to transform the data from the spatial coordinates into a broad wave number domain if the number of sampling points is too small in the spatial domain.Capon's method provides an estimator of the power spectrum in the frequency and wave vector domain using a limited number of measurement points.In Capon's method we first establish the measurement vector or the state vector S defined as where S i (ω) (i = 1,2,••• ,L) denotes the Fourier transformed data into the frequency domain ω at the i-th sensor.For simplicity we consider a measurement of a scalar field such as density, temperature, or magnitude or any components of a vector.The cross spectral density (CSD) matrix R is constructed from the measurement vector as Here the angular bracket denotes the statistical average and the dagger is the operation of Hermitian conjugate.This is an L × L matrix and the task is to reduce the matrix into a scalar power using a suitable projection or weight vector w such that the projected quantity provides an estimator of the wave power in the frequency and wave vector domain, The essence of Capon's method is to use the Lagrangian multiplier method to minimize the wave power (that is to minimize the background level in the spectrum) while maintaining the gain along the scanning direction or the wave vector (unit gain constraint): Here h denotes the steering vector which has the information on the wave vector, The constraint optimization problem can be solved analytically (Capon, 1969;Haykin, 1991), and the weight vector is obtained as It is worthwhile to note that the weight vector is a function of the frequency ω and the wave vector k and it is determined by the measurement itself through the state vector S(ω).The estimator of the wave power in Capon's method is hence given as As an example, Fig. 1 displays synthetic magnetic field data (represented as time series plots) at four distinct points in an irregular, one-dimensional array (located at 100, 180, 290, and 375 km on the array axis).The data exhibits a plane wave with period 30 s and wavelength 300 km with smallamplitude, random, isotropic fluctuations (as noise) added to the wave.The Capon's method is applied to the four-point synthetic data set and the wave number spectrum is obtained (Fig. 2, dotted curve).The Capon's spectrum displays a peak at the signal wave number (0.021 rad km −1 ).The noise field is barely visible in the time series as ripples on top of the plane wave, but this causes a higher background level in the Capon spectrum.All the synthetic data used in the test of the analysis method consist of plane waves as signal and smallamplitude, random and isotropic fluctuating field as noise.In Capon's method the wave power is optimized to give the squared amplitude of the signal and it does not require a priori the knowledge of the number of signals.Therefore the method is of wide use.On the other hand, the estimated power still exhibits a higher background level compared to that of the MSR technique (solid curve, discussed later).The reason for the relatively higher background level in the Capon's spectrum stems from the fact that smallest eigenvalues (which may be referred to as the noise floor) in the CSD matrix contribute to the projection, i.e., due to the finite signal-to-noise ratio in the data.Tjulin et al. (2005) also pointed out that the analytic solutions of the Capon method still show that even if the signal consists of a single plane wave, the resulting wave number spectrum is broadened.The eigenvector-based projection such as the MUSIC algorithm makes use of the orthogonality of eigenvectors in order to enhance the signal-to-noise ratio in the spectrum.Relationship between the Capon's method and the eigenvector-based spectra are discussed in great detail by Haykin (1991).

Eigenvector analyses and the MUSIC algorithm
It was discussed by Samson (1983aSamson ( ,b, 1995) ) and Schmidt (1986) that eigenvalues and eigenvectors of the CSD matrix can be used in analyzing multi-channel data.The MU-SIC algorithm (Schmidt, 1986) suppresses the noise contribution considerably, providing a sharp contrast between the signal power and the noise level in the spectrum by investigating the eigenvector structure of the CSD matrix.Samson et al. (1990) applied the MUSIC method to obtain the twodimensional wave number spectra for ground-based radar observations of gravity waves.The MUSIC algorithm assumes that the data contains signal and noise such that the CSD matrix can be decomposed into these two parts.The state vector is therefore interpreted as a combination of the signal term and the noise term: which we write symbolically as In the MUSIC algorithm the noise is assumed to be isotropic, random fluctuating field.The symbol Q represents the M true signals in the data and it is transmitted to the measurement of S through the matrix A; N represents the noise at each sensor.It can be shown that the CSD matrix becomes decomposed into the signal term and the noise term, too, as The eigenvalues and eigenvectors of R are denoted by respectively.For the noise part the eigenvalues are given as We split the eigenvectors of the CSD matrix R into the signal subspace E s = [e 1 ,e 2 ,••• ,e M ] and the noise subspace The power estimation in the MUSIC method is given as Equation ( 12) represents the noise reduction or filtering procedure performed in the both frequency and wave vector domain.
The MUSIC spectrum is also expressed as where F is the eigenvector matrix of R sorted after the magnitude of eigenvalues in the descending order: The matrix F is an arrangement of the eigenvectors of the CSD matrix, sorting the signal-associated eigenvectors on the left side in the matrix and the noise-associated eigenvectors on the right side.The matrix L is a diagonal matrix and is defined as The MUSIC algorithm is based on finding the eigenvectors associated with noise that are orthogonal to the steering vector with the signal wave vector.The spectrum estimated by the MUSIC algorithm uses the product of the noise eigenvectors and the steering vectors and therefore the method gives the spectrum in the dimensionless unit.It should also be noted that the MUSIC algorithm requires that the number of signals must be known in the analysis to extract the set of the eigenvectors associated with noise.One method to determine the number of signals is to investigate the noise floor of eigenvalues of the CSD matrix (Haykin, 1991).We mean the number of signals by the number of wave vectors associated with signal.This has to be known at each frequency in the analysis.In the context of multi-spacecraft data analysis, the number of wave vectors or the wave modes must be known at each frequency.The problem of the number of signal sources to be known in the MUSIC algorithm has been addressed by Choi et al. (1993), who replaced the diagonal matrix L by −n with Here the power −n is an adjustable parameter in the analysis that controls the asymptotic behavior of the matrix −n that becomes L in the limit n → ∞.In other words, replacing the L matrix by the −n matrix automatically selects the noise subspace of the CSD matrix R. It should be noted that the procedure of the matrix replacement by −n does not stem from a mathematical theory guaranteeing the better functionality of the technique, but it represents an intuitive picture of generalization of the L matrix to soften its sharp transition in the diagonal elements from zero to unity.Therefore, other extensions or generalizations are possible for the MUSIC algorithm.Choi et al. (1993) found that even a small number of n such as n = 2 can successfully reproduce the MUSIC spectrum without knowing the number of signal sources.We use n = 2 in the applications presented in Sects.3 and 4. The spectrum using the extended MUSIC algorithm is given as

Multi-point Signal Resonator technique
The MSR technique makes use of the Capon estimator as well as the extended MUSIC estimator.The notion of the MSR technique is as follows.We use Capon's estimator and obtain the power spectrum that provides the proper value of the spectrum at the signal wave number, and we use additionally the extended MUSIC spectrum with normalization as a dimensionless filter to make the signal-to-noise contrast of the Capon spectrum sharper while keeping the spectral peak values of the Capon spectrum.The power spectrum in the MSR technique is therefore given as Here the factor P EM0 denotes the normalization factor.It is determined by the maximum value of the spectrum P EM .
In the MSR technique P EM (ω,k)/P EM0 serves as a filter that returns the value of unity at the signal wave number and almost zero otherwise (Fig. 2 bottom).An example of the MSR spectrum with comparison to the Capon spectrum is displayed in Fig. 2 top panel.It is worthwhile to note that the eigenvalues are based on estimates in real data, and we have therefore bias and variance: These parameters can be dependent on the number of degrees of freedom in the estimators.
The spectral matrix should be statistically significant based on a larger number of degrees of freedom so that the eigenvalues are well determined.In practice (shown below), we split the time interval by 100 sub-intervals and average the spectral matrix over them.Assuming that the eigenvalues determined at each sub-interval follow the normal distribution, the confidence interval is determined by the degree of freedom (Jenkins and Watts, 1968).The MSR technique assumes that the measured fluctuations represent a set of propagating waves.Therefore, a careful analysis is needed when analyzing wave properties almost at frequency zero.Otherwise the spectrum yields a spurious peak.The MSR technique can also be used for a measurement of vectors such as magnetic field.Furthermore, it is also possible to set the divergence-free condition as an additional constraint.The MSR technique for the vectorial magnetic field measurement can be written in the following fashion: with The trace of the matrix P MSR (ω,k) gives the total fluctuation power:

Numerical tests
To examine the quality of the MSR technique, three kinds of tests are presented using synthetic data.First, we examine if two waves can be resolved even though they have wavelengths close to each other and cannot be resolved in the Capon spectrum.Second, effect of uncertainty of sensor positions is investigated.Third, proportionality of energy of two waves in the spectrum is studied.

Separation of two waves
The MSR technique is applied to synthetic data and compared with Capon's technique as follows.We generate synthetic data of scalar magnetic field exhibiting two waves that have almost the same frequency but different wave numbers (the frequencies are 0.0333 Hz and 0.0335 Hz, and the wave numbers are 0.021 rad km −1 and 0.012 rad km −1 for the wave 1 and the wave 2, respectively).The synthetic waves are sampled at four spatial points with the same array configuration as presented in the previous section.The wave number spectra are evaluated from the four time series data sets using Capon's method (Fig. 3 dotted line) and the MSR technique (solid line), respectively.Both methods show peaks in the spectra at the signal wave numbers, but the background level is much reduced in the MSR spectrum.
We then modify the values of the wave numbers in the synthetic data.The wave 2 has now the wave number 0.020 rad km −1 , which is close to that of the wave 1. Figure 4 displays the wave number spectra for this data set.The MSR spectrum can resolve two peaks, whereas Capon's spectrum exhibits a peak with a larger broadening.In Fig. 4, the wave number resolution is k = 1.0 × 10 −5 rad km −1 , while the difference between the two signal wave numbers is 1.0 × 10 −3 rad km −1 .There are enough data points between the two peaks, and identifying the two signals is successful using the MSR technique.It is also worthwhile to note that the MSR technique gives the same spectral peak values as that of the Capon technique.

Effect of uncertainty of sensor positions
Effect of uncertainty of sensor positions is investigated, as well.We use the data set of the numerical test 2 and determine the spectra for shifted sensor positions (that represents the uncertainty of sensor positions).Figure 5 displays the peak wave numbers identified in the spectra as a function of the relative uncertainty of the sensor positions.When the positions are known accurately enough ( r/r ≤ 1%), then  the spectra exhibit peaks at almost the right wave numbers.For increasing uncertainty of the sensor position, the spectra exhibit a shift of the peaks toward lower wave numbers although the separation between the two wave numbers is still possible.At the position uncertainty 10% or larger, the peak shift becomes significant and the spectra do not exhibit the right wave numbers any more.Therefore, it is ideal that the sensor positions are accurately known and do not change during the measurement.

Proportionality in the spectrum
We investigate the proportionality of the MSR spectrum.Here, we mean the proportionality by the proper order of peak values in the energy spectra, i.e., the spectrum should exhibit larger peak values for larger fluctuation amplitudes, otherwise a wrong order of the spectral peaks causes a distorted shape of the spectrum.We change the amplitude of the wave component 1 in the synthetic data for Fig. 4: 1/2, 1/5, Fig. 7. Time series of the magnetic field magnitude, the ion bulk speed and the ion number density in the solar wind measured by Cluster-1 spacecraft.and 1/10 of that of the wave component 2. In the ideal case, the spectra should exhibit peaks with the ratios 1/4, 1/25, and 1/100, respectively.The MSR spectra are displayed in Fig. 6.We find that the order of the spectral peaks is the same as the inputs (i.e., larger peak values for larger fluctuation amplitudes) but on the other hand the ratios of the spectral peaks tend to be emphasized roughly by the square of the ideal energy ratios (O(10 −1 ), O(10 −3 ), and O(10 −5 ), respectively).Therefore, our technique can qualitatively reconstruct the shape of the wave number spectrum, but a further correction method is needed for quantitative analysis.

Application to spacecraft data
The MSR technique is applied to the observation of magnetic field fluctuations in the solar wind performed by four Cluster spacecraft (Escoubet et al., 2001;Balogh et al., 2001) in the time interval from 14:15 to 14:45 UT on 12 February 2002 (Fig. 7).This time exhibits a good example in which two waves are identified by the MSR technique while they cannot be seen by the Capon method.We chose the interval for several reasons: (1) a good quality of tetrahedron with spacecraft separation about 100 km in order to minimize the distortion of energy distribution caused by an irregular tetrahedral configuration of Cluster (Sahraoui et al., 2010a), (2) stationary fluctuation such that the mean fields are regarded as almost constant, and (3) small variation of flow velocity to improve the accuracy of the Doppler correction.The mean magnetic field strength is about 6.3 nT and the mean flow speed and density are 512 km s −1 and 3.3 cm −3 , respectively.
Fluctuations in the solar wind are believed to be in a fully developed turbulent state, as energy spectra of both magnetic field and flow velocity fluctuations typically exhibit an power-law curve in the frequency domain (Coleman, 1968;Matthaeus and Goldstein, 1982;Matthaeus et al., 1982;Marsch and Tu, 1990;Podesta et al., 2007), which is reminiscent of the inertial-range spectrum of turbulence in Kolmogorov theory with the power-law index close to −5/3.At frequency near around 0.1 Hz to 1 Hz, measurements further show a spectral breakpoint.The spectrum becomes a steeper power-law above these frequencies (Leamon et al., 1998;Smith et al., 2006;Alexandrova et al., 2008) and is characterized by the index steeper than −2 (Behannon, 1978;Denskat et al., 1983;Goldstein et al., 1994;Leamon et al., 1998;Bale et al., 2005).Recent analyses of data from the multi-spacecraft Cluster mission have shown that magnetic spectra in the range between 0.5 Hz and 20 Hz scale with the index −2 to −3 (Sahraoui et al., 2009;Kiyani et al., 2009;Alexandrova et al., 2009), and suggest that there is a second breakpoint with still more steeply decreasing spectra at frequencies beyond 20 Hz.
The frequency spectrum of the magnetic field in the analyzed time interval (Fig. 8) exhibits the inertial-range spectrum at lower frequencies and the steeper spectral curve above the spectral break at about 0.3 Hz.We focus on the wave field at 0.15 Hz, representing the high-frequency limit of the inertial range, for the wave analysis using the Capon and MSR techniques.It is shown that magnetic field fluctuations at this frequency exhibit two distinct wave vectors that can be resolved using the MSR technique while they cannot be separated clearly by the Capon method.The total length of time interval used in the analysis is 30 min with the time resolution 0.2 s.We split the total interval into 100 sub-intervals with the length 4096 data points for statistical significance in the analysis.The sub-interval is 819.2 s long, while we investigated the fluctuation at the period 6.7 s (at 0.15 Hz).The ratio of the wave period to the length of the sub-interval is about 1/123.
Figure 9 displays two energy distributions in the wave vector domain at this frequency.The spectra are determined by the Capon method (left panel) and the MSR technique (right panel), respectively.The 3-D wave-vector spectra are projected into three planes (k x -k y , and k x -k z , and k y -k z ) in the MFA (mean field aligned) coordinate system (spanned by the mean magnetic field in the z-direction and the flow velocity direction in the x-z-plane).For projection we summed the spectrum over the wave numbers perpendicular to the respective planes.The Capon spectrum exhibits a peak at the wave vector about (0.0022, 0.0005, 0.0000) rad km −1 (wave 1) as well as two extended structures, one along the y-axis and the other nearly in the positive z-direction with slight inclination to the negative x-direction.The MSR spectrum exhibits a peak around this wave vector, but it also shows another peak at the wave vector about (−0.0010, −0.0010, 0.0058) rad km −1 (wave 2).The overall shape of the spectrum looks similar to each other, but the existence of the secondary peak can be clearly identified by the MSR technique.It is difficult to identify the secondary peak in the Capon spectrum.It can also be seen that the noise level is much reduced and the spectral shape is sharper in the MSR technique.
Using the two wave vectors obtained by the MSR technique, frequencies and phase speeds of the waves are estimated in the plasma rest frame by correcting for the Doppler shift.Table 1 summarized the wave numbers, propagation angles from the mean magnetic field direction, the rest-frame frequencies and phase speeds, respectively.We find that the wave 1, characterized by quasi-perpendicular wave vector to the mean magnetic field, has a very small frequency and phase speed in the plasma rest frame.However, the uncertainty of the wave vector determination and the large flow velocity causes a large uncertainty in computing the rest-frame properties and the accurate frequency and phase speed cannot be obtained in the analysis.Therefore it is not clear if the wave 1 really represents a convected spatial structure or a propagating wave.The wave 2, on the other hand, exhibits a parallel wave vector and the rest-frame frequency and phase speed were determined as 0.45 rad s −1 and 76 km s −1 within the accuracy about 40-50%.It is interesting that the phase speed of the wave 2 is close to the Alfvén speed at about 75 km s −1 .Considering the factor cosθ kB for oblique propagation, the Alfvén wave phase speed is estimated as about 72 km s −1 .Though the uncertainty of the phase speed determined in the analysis is still relatively large, the agreement between the two speeds is suggestive of the the existence of the Alfvén wave in the solar wind.Fig. 9. 3-D energy distributions in the wave vector domain in the mean-field-aligned coordinate system at frequency 0.15 Hz, determined by Capon's method (left) and the MSR technique (right).In this coordinate system the k z -axis (vertical axis) points to the direction of the mean magnetic field, and the k x -k z plane is spanned by the mean field and the flow velocity direction.The 3-D spectrum is represented using projection into 3 planes: it is summed over the wave vector component perpendicular to the respective planes.The vertical bar at color scale represent the 95% confidence interval.
Table 1.Wave properties at spacecraft-frame frequency 1.5 Hz: wave number k, angle between the wave vector and the mean magnetic field θ kB , frequency in the plasma rest frame ω rest , and phase speed in the rest frame v ph .k θ kB ω rest v ph rad km −1 degree rad s −1 km s −1 Wave 1 0.0023 ± 0.0003 90 ± 13 0.00 ± 0.17 2 ± 48 Wave 2 0.0060 ± 0.0003 14 ± 3 0.45 ± 0.17 76 ± 33 MUSIC algorithm that the number of signals must be known can be solved by using the eigenvalue ratios as proposed by Choi et al. (1993).While Capon's technique is already widely applied to Cluster data and is certainly a useful wave analysis method, the MSR technique offers further improvements in two points: significantly reduced background level of the spectra, and the capability of identifying wave components that have similar wavelengths.The MSR technique was applied to the synthetic magnetic field data as well as the Cluster observation of magnetic field fluctuations in the solar wind.We have found a wave component in the solar wind that is propagating parallel to the mean magnetic field at about the Alfvén speed.This wave was not found in the Capon spectrum.The MSR technique offers further improvements in performing dispersion relation analysis and spectral analysis in a broad wave vector domain.Finally, this method can be applied to characterizing fluctuations in various regions of space plasma.Magnetosheath, for example, would be a good place to use our mehod, where mirror mode fluctuations are known to present and separion between signal of interest and noise (unwanted fluctuations) is important.

Fig. 1 .
Fig. 1.Synthetic scalar magnetic field data at four points.The data set is used for Capon's method and the MSR technique displayed in Fig. 2.

Fig. 2 .
Fig. 2. Top: Capon-spectrum in the wave number domain determined by four-point measurement of synthetic scalar magnetic field data shown in Fig. 1 (dotted curve); the MSR spectrum for the same data set (solid curve).Bottom: Spectrum determined by the normalized, extended-MUSIC algorithm for the same data set.The vertical bar represents the confidence interval.

Fig. 4 .
Fig. 4. Another example of the Capon and MSR spectra for the synthetic data.

Fig. 5 .
Fig. 5. Identified peak wave numbers for the two wave components for shifted sensor positions for the data set used in Fig. 4. The xaxis is the relative uncertainty of the sensor position in percentage.

Fig. 8 .
Fig. 8. Frequency spectrum of the magnetic field fluctuation for the time interval shown in Fig. 7.