4. Spatial Correlation Analysis

Correlation analysis studies the relationship between measurements from different observation points. The observers can be either separated spacecraft in space or different ground stations. When a disturbance or a wave travels from one observer to another, it is observed with a timing difference in the time domain or a phase difference in the frequency domain between the two observers plus any change of the disturbance during the travel time period. From the timing or phase difference and the separation between the two observers, one is able to determine the wavenumber of the wave. The change in the signal while propagating reduces the level of the correlation between the two observed signals and thus the time or phase difference may be a function of frequency and possibly time.

4.1 Background.

4.1.1 Correlation analysis in the time domain.

If two observers measure the same quantity at different locations and the two time series data are Q(0, t) and , where we have assumed that the location of the first observer as the origin and the separation between the two observers is , the cross-correlation coefficient is defined as

where t1 and t2 are the start and stop times of the selected interval, and K is the normalization factor and is

The time lag with the highest cross-correlation coefficient can be used as the timing difference, , between the two observers. The integration is realized through the summation over the time window. In general, a narrower width of the window provides a greater maximum of the correlation coefficient because the noise within the window width may be likely weaker. For discontinuity analysis, one should select a narrower window because for an ideal discontinuity, the cross-correlation coefficient decreases as . When the window is much greater than the time shift, the correlation coefficient appears flat. For wave analysis, a wider window, which includes several wave cycles will be statistically more significant.

If a wave is monochromatic, the correlation coefficient will periodically increase to unity and hence there is an ambiguity in determining the time lag. Fortunately, the waves in space are not purely monochromatic and one usually has no difficulty to find a peak in the correlation coefficient with the greatest value. One still has to be careful since in some situations it is possible that the greatest correlation occurs associated with observations of different wave fronts.

What limits the resolution of the timing difference? If a wave is purely monochromatic, the resolution of the timing difference may be infinitely high, and depends on the accuracy of the clocks rather than the resolution of the data. However, because the waves in space are not monochromatic, as just discussed above, and also because a wave may change slightly as it propagates from one observer to another, or be slightly different at some distance along the wavefront where the second observation is made, there is uncertainty in the timing difference. This uncertainty is anti-correlated with the value of the highest correlation coefficient and equal to or smaller than the resolution of the data. For example, since a discontinuity departs most from an ideal monochromatic wave, the uncertainty in the timing difference is the highest and equals half of the time resolution of the data. To minimize the uncertainty in timing discontinuities, one should use the data with the highest resolution and filter the noise before an analysis. Otherwise, in common practice cross-correlation analysis is replaced by simply timing the delay of the discontinuity between the two observers. When multiple waves are present, each of them may have a different time delay. One needs to band-pass filter the data before the correlation analysis. Otherwise, the time delay will be determined by the wave of greatest power and the correlation coefficient will be smaller.

4.1.2 Correlation analysis in the frequency domain.

In the frequency domain, correlation analysis studies the phase difference between observations at two locations. From this phase difference, one may be able to determine the wavenumber.

If two time series data are Q(0, t) and and their Fourier spectra are and , the cross-correlation coefficient is

where the asterisk stands for the conjugate and <> indicates an ensemble average. The ensemble average is performed in time domain and equals the Fourier spectrum of the corresponding quantity in equation (4.1) [e.g. Jenkins and Watts, 1968]. Note here that equation (4.2) is different from that given in most theoretical studies. In theoretical studies, the ensemble average is performed over wavevector. While the wavevectors are obviously known in theory, observationally we believe that with a very limited number of observers we are still far from being able to resolve multiple wave modes at a single frequency.

Since the correlation coefficient in the frequency domain is complex, it can be written as

where and are coherence and phase lag, respectively. The coherence gives a measure of wave versus noise. The noise can be either due to random fluctuations or due to waves with a finite coherence length compared with the separation between the two observers. The coherence length along a wavefront is the spatial scale of the wave. The coherence length normal to the wavefront is roughly the traveling distance of the wave, including both convection and propagation, within its decay time. For example, since the convection speed of the solar wind is much greater than the phase velocity, the waves in the solar wind often have larger coherence length along the flow than transverse to the flow. The difference between random fluctuations and waves with a coherence length shorter than the separation is that the power spectra are broad for the former and at least one of them is peaked for the latter.

If the noise is due to random fluctuations, the phase difference of a wave between two observers is

where , ...... Thus, the wave number can be determined, but with the so-called ambiguity. Here the direction of the wave vector is determined from the methods discussed in sections 3.1 and 3.2. Since the frequency in equation (4.4) is in the spacecraft frame, further determination of the wave properties needs to consider the Doppler shift, equation (3.18). The time lag can be calculated from , which is a function of frequency.

In summary, correlation analysis can determine either the propagation speed of a discontinuity or the wavelength of a wave. For a discontinuity, the correlation analysis has to be performed in the time domain. For a wave, the analysis can be performed in either the time domain or the frequency domain. The analysis in the frequency domain contains the information about the phase as functions of frequency and hence can provide a better overall view of the waves at different frequencies. The result of an analysis in the frequency domain is in principle the same as in the time domain with prefiltering. These analyses are restricted to waves of a single source.

4.1.3. Relationship between the quantities in satellite frame and plasma frame

In section 1.1, we discussed the relationship between a quantity in the satellite frame and in the plasma frame for a single spacecraft. When comparing observations from two satellites at different times, the variations other than simple convection from one satellite to the other will affect the data analysis. These variations include spatial (between the two satellites) and temporal (between the two times) variations due to the physical evolution of the phenomenon being studied, and changes in the velocity and spacecraft separation.

The difference between the quantity Q measured by satellite i at location at time t and by satellite j separated by a distance from satellite i at time is, see Fig. 4.1 for the geometry,

where should be understood as . A time domain correlation analysis is to provide time shift with a minimal difference, i.e. goes to zero. Taylor expanding equation (4.5) and retaining the terms up to the second order, the difference operator is

We should point out that the validity of the Taylor expansion relies on to be much smaller than Lc and to be much smaller than tc, where Lc and tc are the coherence length and coherence time, respectively. Whereas the measured quantity can contain discontinuities. Namely, is not necessarily much smaller than the scale of the discontinuities (in the normal direction). As long as the correlation analysis provides a meaningful correlation, the Taylor expansion is valid. The terms in the first bracket are due to the linear effects. The main differences observed by the two spacecraft are caused by the temporal variations in the plasma frame and the spatial non-uniformity of the plasma along the motion direction and along the separation direction of the two spacecraft. The function of the correlation analysis is to provide a that minimizes the value in this bracket. Namely, the correlation analysis converts linearly the temporal variations in the spacecraft frame into spatial variations in the plasma frame of reference. It is less capable of providing information about spatially static structures in the spacecraft frame. In general the nonlinearities in the temporal and spatial variations will result in a residue for these terms. This residue reduces the correlation coefficient of the analysis. The terms in the second bracket are due to the nonlinear effects of temporal variations. The terms in the third bracket are due to the nonlinear effects of nonsteady state, nonuniform plasma motion. Similarly, the terms in the fourth bracket are due to the temporal and spatial changes in the separation between the two spacecraft.

 Fig. 4.1. Geometry of two separated spacecraft relative to a wave front. Satellite A is located at point A and satellite B at point B. The plasma element that is measured by satellite A at time t and the plasma element at point A' is measured at .Similarly, the plasma element that is measured by satellite B at is from point B'. A correlation analysis compares measurements at different locations and/or different times.

In most cases, the variations in the separation vector (the fourth bracket) is negligibly small. The velocity variation (the third bracket) could become significant near an oscillating boundary and for a circularly polarized wave the velocity of which could be twisted by the wave perturbations, in particular when the background velocity is small.

For a single spacecraft, L = 0, the first order terms recovers equation (1.1). Note that the velocity in equation (4.6) is not necessarily the plasma velocity but could be the velocity of any frame relative to the spacecraft. In a shock frame or a wave frame, = 0. The first order term is . In the ideal 1-D case, where is defined as the time delay with the maximum correlation coefficient, and it yields equations (2.20) and (3.20) by noting that the speed of the wave frame relative to the spacecraft is the sum of the phase velocity and the plasma frame speed.

4.1.4 Data analysis of observations from a cluster of satellites.

A cluster'' of spacecraft refers to four closely separated satellites. The relative locations of the satellites can be represented in general by an irregular tetrahedron. Regulated by orbital dynamics, the shape of the tetrahedron changes with time. If the change in their relative distance is much less than in the relative distance itself during the interval of the event being studied, the change of the tetrahedron may be neglected. The uncertainties in the determination of the relative locations and in the intracalibration of the instruments among the satellites determine the lower limit of the uncertainty of any analysis.

In principle, when the cluster is not coplanar, the propagation direction of a 1-D discontinuity or a 1-D steady state wave can be determined without the assumption that the propagation direction is along the minimum variance direction. With the approximation that the propagation direction is along the minimum variance direction, either the unsteadiness of the wave (in the wave frame), or a 2-D steady state structure or wave can be resolved under certain conditions. A moving flux tube with a constant shape is an example of a 2-D steady state structure. A purely periodic surface wave in the wave frame is an example of the 2-D steady state wave. For a given perturbation, even with a noncoplanar cluster of satellites, there is no general conclusive interpretation. One has to make approximations about each of the extension in three dimensions (bounded, periodic, or uniform) and its time dependence (pulsating, periodic, or steady). Differences in these approximations lead to different interpretations, which has been the major cause of existing controversies in the science community. Each of these approximations, however, can be tested both using either the information from other instruments or techniques that are not included in this article as standard techniques, and using the constraints provided by statistics and physical consequences.

In the following sections, we discuss the application of correlation analysis to multiple spacecraft observations. Most discussion is based on the authors' experience from previous missions and recent results from other groups. We emphasize that the problems with massive multiple spacecraft data analysis have not yet been fully exposed. Significant modification may be needed when applying the following principles.

4.2 One Dimensional Discontinuity Analyses.

As discussed in section 2.1.1, to completely determine observationally the properties of a 1-D steady state discontinuity, more than one satellite is needed. The central issue is to determine the direction and magnitude of the velocity of the shock frame relative to the spacecraft without imposing the conservation laws which should be verified by the results. The timing difference between spacecraft i and j is

Note that, as discussed at the end of section 4.1.3, the above expression is accurate to the first order and only for an ideal steady state 1-D discontinuity. In comparison, equation (1.2) does not require the steadiness approximation. With three spacecraft, can be uniquely determined if the plane containing the spacecraft is not parallel to the discontinuity. With four noncoplanar spacecraft, equation (4.7) can become overdetermined and can be determined using best-fit, for example, or average over the several obtained from each subset of the three time differences [Russell et al., 1983]. In practice, one may obtain several time delays crossing a discontinuity for each pair of spacecraft, for example, one at the leading edge, one in the middle and one at the trailing edge, when using high resolution data. If one assumes that the difference between the two time delays is not due to temporal variation of the thickness, averaging these sets will improve the statistical significance of the results.

Using the resultant and equation (2.7), one can derive the quantities in the R-H relations. The accuracy of the time delay is limited by the temporal resolution of the data when the fluctuations within the discontinuity cannot be resolved. A greater time delay compared with the temporal resolution implies a smaller uncertainty in the timing, whereas during this longer interval, there is a greater chance for the discontinuity to change its properties. If the resolution is sufficiently high, the thickness of the discontinuity then can be determined from the duration of the crossing.

A major uncertainty arises from the determination of the relative locations of the spacecraft.

If the type of the discontinuity is known or the techniques discussed in section 2 are applicable, combination of the timing method with those in section 2 will improve the statistical significance of the results. In addition, one may be able to find intervals when the discontinuity lies between the spacecraft and hence some of them provide the upstream conditions while others the downstream ones. Quantities in equations (2.10) and (2.11) can be substituted with those simultaneously measured by different satellites instead of those measured at different times from the same satellite.

Since each method has its own uncertainty, when averaging results from different methods, how to weigh each of them becomes important. Russell et al. [1983] documented several cases and provided comprehensive discussions on the uncertainties.

Nonsteadiness of a discontinuity is observed as the difference in the duration of the discontinuity crossing registered among different satellites and the difference in the time shifts. It can be treated as the variation around the average. The nonsteadiness includes the temporal variations in the discontinuity's orientation (i.e. its normal direction) and the thickness, and the changes in the shock velocity relative to the satellites.

4.3 Two Dimensional Structures

4.3.1. Surface Wave Analyses.

Surface waves in magnetospheric physics refer to the wavy motion of a discontinuity. To describe such a motion, the most important parameters are the wavelength and the displacement, as well as the period. The wavelength is measured along the discontinuity surface in the direction of the surface wave propagation and the displacement is normal to the average surface. Simple boundary oscillations are a 1-D approximation when the wavelength is infinity. With a finite wavelength, the discontinuity is bent. As expected, the boundary normal varies. If the thickness of the boundary is much smaller than the curvature of the boundary the discontinuity analysis methods discussed in section 2 are applicable. The normal direction of the boundary can be determined by each of the crossings. The variation in the normal direction can then be measured by the satellites, see Figure 4.2 for example. With multiple satellites, one is able to determine the characteristics of a surface wave [Song et al., 1989]. Mottez and Chanteur [1994] discussed a similar technique to analyze a general 3-D thin surface using a cluster of satellites.

 Fig. 4.2. As a surface wave passes by two spacecraft, a few timing differences and normal directioncan be obtained in order to determine the parameters of the surface wave.

A different approach has been recently developed using the remote sensing'' technique [Walthour et al., 1993]. As a boundary oscillates, the plasma and field parameters in the regions adjacent to the boundary vary. For example, the mathematical description of a surface wave is an exponential decay in the amplitude of the perturbation away from the boundary [ e.g., Southwood, 1968]. From these variations, a satellite can remotely sense the boundary and its motion through the course when it approach to and departs from the boundary. This new technique requires high quality measurements. It has been tested by using two satellites data [Walthour et al., 1994].

4.3.2. Flux Tube.

When the plasma beta is not extremely large, the anisotropy in a plasma introduced by the field often results in such structures that have similar sizes across the field and their greatest dimension along the field. Usually such a structure has properties distinct from its surrounding plasmas. Russell et al. [1990] provided a collection of comprehensive discussions on flux tubes from various regions in space. The majority of previous observations are based on single spacecraft measurements. The geometry of a structure, for example, a flux transfer event [Russell and Elphic, 1978] or a plasmoid, is often deduced from time series data (which provides the variation in one dimension), the region where it is observed on a statistical basis (which provides the scale of the second dimension), and the assumption that the structure is very long in the third dimension. Sometimes, observations from multiple satellites are available. Examination of these events confirms the general geometry people have developed for these phenomena.

A cluster of spacecraft can provide much more detailed information about these flux tubes. A general approach to analyzing these structures is similar to that of surface waves.

It is worth mentioning that since a flux tube does not satisfy the 1-D approximation, the minimum variance direction across the event may not have any definite meaning. For example, unless there are other indications, one should not take the minimum variance direction as the direction of the symmetric axis because the axial field may increase toward the axis due to a pinch process that arises from an axis-aligned current within the tube.

4.4. Determination of Electric Currents.

The electric current is a very important physical quantity in space physics. However it is extremely difficult to measure in space. Because the conductivity in a collisionless plasma is poorly known, Ohm's law cannot be used to evaluate the current. In theory, there are two ways to measure the current. They are

where subscript i denotes each charge species. As discussed in section 1.3.3, it is difficult to measure accurately the fluxes of the species that carry most of the current due to the intrinsic limitation of particle detectors and spacecraft charging. This is particularly true in collisionless plasmas where particle gyro motions are dominant. In most regions in the magnetosphere, the relative drift between ions and electrons is small. The limitation of instruments further reduces the ability of resolving such a small difference. For example, as shown in Equation (1.10), the effective cutoff velocity of a detector depends on the species. The spacecraft charge will have much greater effects on electron detection than on ion detection. When the spacecraft is strongly negatively charged, the difference in the errors in the velocity measurements will produce a large error in evaluating Equation (4.8) even if the detector has a similar small lower cutoff energies for electrons and ions. Although there are several reports on the determination of the electric current using equation (4.8), see Frank et al., [1981], for a partially successful test, we remain cautious about quantitative use of the current from direct particle measurements.

4.4.1. Single Satellite.

Equation (4.9) has been widely used to estimate the current. Almost all previous single satellite usage was based on a 1-D current sheet assumption. In 1-D

where D is the thickness of the current sheet across which is evaluated. To determine D from a time difference, one needs to know the velocity of the current layer relative to the spacecraft. Therefore all discussions given in section 2 on determination of the shock frame apply here.

4.4.2. Dual Satellites

When measurements from two satellites are available, there are two simple approaches to estimate the current. The first approach is to take the difference between the simultaneous field measurements from the two satellites assuming that the current sheet is one dimensional and lies between the two spacecraft. This corresponds to = 0 case in Equation (4.6). The first order terms gives the same expression of the current as in Equation (4.10), but the field differences are evaluated from simultaneous measurements at different locations. The current can then be determined as functions of time [Elphic, 1989]. The resultant time series is often highly fluctuating and spiky, see Fig. 4.3, for example. The causes for such fluctuations could be multiple: there may exist physically meaningful small scale structures within the current layer, the relative motion among sub-current sheets in the normal direction can result in varying time-space conversion, and the 1-D approximation to these structures may break down. When the 1-D approximation breaks down, the time series of the field differences cannot be simply interpreted as currents.

 Fig. 4.3. Field differences at two satellite locations which can be used to infer the current [Elphic, 1989]. The upper four panels show the magnetic field measurements from two magnetopause crossings from ISEE 1. The lower four panels show the differences between the ISEE 1 and 2 field measurements.

Taking an orthogonal approach, one may compare the timing of the same field measured by the two satellites. This corresponds to setting the left-hand-side of Equation (4.6) zero. Under the assumption that the structure of the current is stationary, but the motion of the current changes with time, retaining the first order terms, the velocity of the current layer can then be calculated by dividing the satellite separation along the normal of the current by the timing difference between the two satellites' measuring of the same magnetic field, as shown by equation (4.7). When the velocity of the current layer is known, from equation (4.10) one can determine the current according to the time series of each satellite, see Figure 4.4. In theory, the currents derived from the two satellites should be the same if the current is in steady state. In fact, the resultant currents are different, indicating that the current is not in a steady state. How to reconcile this inconsistency remains to be investigated.

 Fig. 4.4. Two-satellite determination of the current [McComas et al., 1986]. Under the 1-D steady state assumption of the current structure, the velocity of the current motion can be determined, middle panel, and then the current density can be derived, lower panel.

Often the two spacecraft are on the same side of the current sheet (after passing through or before the crossing) at different times. These data on the same side can be used to help test the 1-D time stationary assumption by combining the two approaches above.

4.4.3. Multiple Satellites.

The above techniques rely on the 1-D assumption of the current. The measurements from a cluster of satellites will enable us to determine the spatial derivatives and hence the current in equation (4.9) without the 1-D assumption. The spatial derivative method applied to a quasi-stationary current system is very different from the analysis of a fast moving current structure. In the former case, correlation analysis may not be able to provide a meaningful timing difference. Therefore the approach of setting equal to zero is used. In the latter case, can be determined and hence more information is available to measure the current and its motion velocity.

Quasi-stationary currents. A current system can be considered as quasi-stationary if where V is the velocity of the current system relative to the spacecraft, (note here that when V=0, goes to infinity). In other words, the measured temporal variation in the field is much smaller than the difference among satellites. This assumption is applicable for a cluster moving slowly through a current region that is thicker than the separation of the satellites. The electric current density can be determined by [Dunlop et al., 1988],

The direction of is outward along the normal direction of the surface formed by spacecraft i, j, and k. Four such can be obtained. The current density within the cluster can be calculated by averaging the four current densities. The direction of the current, is parallel to the vector summation of s, and its magnitude is

An important test of the result is to check whether the sum of is much smaller than each of them.

We should emphasize that the approach under the steady state approximation is not the same as the approach when the spacecraft separations are much smaller than the spatial scale of the current structure. When two spacecraft see basically different things, the correlation analysis does not provide anything meaningful.

Traveling current structures. Current systems in space often move quickly relative to the satellites. A boundary or surface, that is neither necessarily thinner than the spacecraft separation nor planar, can pass a cluster of spacecraft quickly without significant changes in itself. Under this circumstance, the correlation analysis can provide very important information. Note here that more than one time shift (here one may use the measured timing difference directly instead of performing a correlation analysis because to perform the correlation analysis at each of the edges is usually very difficult) can be determined for a single passage of the current for each pair of spacecraft, at its leading and trailing edges. Averaging the velocities obtained from equation (4.7) for all satellite pairs, one yields the velocity of the current structure. Here we recall that there is no assumption made on the shape or normal of the structure when using equation (4.7) and hence it is applicable to structures of any regular shape. The geometry of the current structure can be examined according to the timing difference of each pair from the average and their relative locations.

With the derived velocity , the current density at spacecraft i can be calculated by

where prefix denotes the difference measured from the same spacecraft at different times.

4.5. Wave Analyses.

In wave analysis, the time derivatives in Equation (4.6) are, in general, not zero in the spacecraft frame. The approach that sets will eliminate all the terms of frequency and hence cannot be used to analyze waves. With a non-coplanar cluster of four spacecraft, in principle, one is able to determine the frequency, wavelength and phase velocity of a planar wave if its wavelength is of the order of spacecraft separation. For waves with a wavelength that is either much greater or smaller than the spacecraft separation, there is no general method to derive all these three wave parameters.

For perturbations that can be approximated as 1-D, but not necessarily monochromatic, one can follow the general steps.

(1) Select the perturbations of interest and band-pass filter the time series data. Note that if the wave is not monochromatic, the spectrum will have finite width.

(2) Perform wave analysis to determine the direction of the wave propagation for each spacecraft, and take the average of these vectors as the direction of the wavevector.

(3) Perform correlation analysis for each component of the filtered data using each pair of spacecraft.

(4) Calculate the apparent velocity of the perturbation from the correlation analysis and compare its direction with the direction of the wave derived from the wave analysis.

(5) Analyze the deviation of each quantity from its corresponding average. In general, these deviations can be attributed to temporal variations that may not be directly related to the wave.

The biggest difficulty one may face here is that there is too much information and some of these pieces are in conflict. While one could use average and normalization to enhance the robustness of the results, we suggest one leave enough room for different possible interpretations because we think that at this moment we have very limited knowledge about the problems inherent in an analysis of such complexity.

Neubauer and Glassmeier [1990] proposed the concept of wave-telescope. As illustrated in Fig. 4.5, when a wave propagates across a cluster of satellites, each satellite observes a different phase according to the wave propagation direction and the relative location among them. From these phase shifts, in theory, one is able to derive the wave properties. There are established techniques in radio science to analyze such information. While this proposal is intriguing, we note that in space wave analysis the signal-noise ratio is often low. Multiple modes in space plasma often exist, see section 3.3 for more details, and are difficult to distinguish. Often when testing with a known wave, a scheme may work, but there may be multiple solutions when the input is unknown.

 Fig. 4.5. Many wave modes pass across a cluster of satellites [Neubauer and Glassmeier, 1990].