Annales Geophysicae On the error estimation of multi-spacecraft timing method

The multi-spacecraft timing method, a data analysis technique based on four-point measurements to obtain the normal vector and velocity of an observed boundary, has been widely applied to various discontinuities in the solar wind and the magnetosphere studies. In this paper, we perform simulations to analyze the errors of the timing method by specifying the error sources to the uncertainties in the determination of the time delays between each spacecraft pair. It is shown that the timing method may have large errors if either the spacecraft tetrahedron is largely elongated and/or flattened, or the discontinuity moves much slower than the constellation itself. The results, therefore, suggest that some of the applications of the timing method require reexamination with special caution, in particular for the studies of the slow-moving discontinuities associated with, for example, the plasmaspheric plumes.


Introduction
Multi-spacecraft missions, such as the Cluster (Escoubet et al., 1997), THEMIS (Angelopoulos, 2008), and the future MMS mission, were designed with the potential to unambiguously separate the spatial and temporal variabilities.Since the idea of the Cluster mission with four identical satellites making simultaneous measurements at different locations was proposed in the 1980s, many data analysis techniques have been developed, based on time series of fourpoint observations, to enable the three-dimensional examination of the observed structures in space plasma.Among these techniques, the one most widely used is probably the Correspondence to: X.-Z.Zhou (xzhou@igpp.ucla.edu)multi-spacecraft timing method, sometimes referred to as the triangulation method and/or the time-delay method, in the studies of the orientation and motion of planar discontinuities (Russell et al., 1983;Harvey, 1998;Sonnerup et al., 2008).
Assuming that a one-dimensional discontinuity moves with a constant velocity (which registers a pure spatial variation as a set of time-series variations), the multi-spacecraft timing method (hereinafter referred to as the timing method) makes use of the measured time differences between the passage of the discontinuity over satellites, along with the relative positions between the crossing locations, to calculate the normal unit vector N and the normal velocity V. Taking the four-spacecraft Cluster mission as the example, if the discontinuity can be identified unambiguously, i.e., all of the discontinuity passage times t α are well determined, we have: where spacecraft 1 is arbitrarily taken as the reference, and R α1 represents the relative position between the discontinuity crossing locations observed by spacecraft α and spacecraft 1.
As long as the four satellites are not coplanar, the equation yields a unique solution of the normal vector and the normal velocity.
It should be kept in mind that the position vectors of different satellites are not simultaneous as the crossings occur at different times.Therefore, the instantaneous interspacecraft distances D α1 , which can be directly obtained (from the European Space Operations Centre for the Cluster mission), cannot simply be used to replace the relative position vectors R α1 in Eq. (1).Instead, the position changes of the spacecraft during these time intervals have to be taken into account, with the terms of (t α − t 1 )V C being added to D α1 to obtain the corresponding R α1 vectors, where V C represents the velocity of the Cluster constellation.
In practice, however, the observed discontinuity can hardly be unambiguously time-stamped on each spacecraft, either because there are generally other local structures superposed on the discontinuity, the discontinuity may not be a perfect Published by Copernicus Publications on behalf of the European Geosciences Union.
one-dimensional structure, or the discontinuity itself may evolve in time.One way to improve the results is to determine the optimal relative time delay for each spacecraft pair, which can be calculated either by maximizing the value of the cross-correlation (e.g., Song and Russell, 1999) between these data streams, or by the application of the Local Wavelet Correlation (LWC) technique (Soucek et al., 2004).A total of six time delays t αβ (1<α≤4, 1≤β<α) can be calculated, and the six corresponding relative position vectors R αβ are also obtained, by adding the t αβ V C term to the interspacecraft distances D αβ .Thus we can obtain an over-determined equation set containing three unknowns and six independent equations.By minimizing the function (2) a least-squares approach can be used to determine the optimal solution, i.e., the normal vector and the normal velocity of the observed discontinuity (Harvey, 1998).Equation ( 2) further highlights the different roles R αβ and D αβ play in the application of the timing method, namely, if one select the instantaneous interspacecraft distances D αβ instead of using the relative position vectors R αβ in the timing method, the resulting velocity would be in the frame moving with the Cluster constellation.As we know that in most cases the velocity of the discontinuity should be cast in the Earth's frame for physical significance and usefulness, it is important to make sure that the relative position vectors R αβ are used (instead of D αβ ) when performing the timing analysis.
Since the launch of the Cluster constellation, the timing method has been applied to various types of discontinuities in the near-Earth space, e.g., the interplanetary discontinuities (Knetter et al., 2003(Knetter et al., , 2004)), the terrestrial bow shock (Bale et al., 2003;Maksimovic et al., 2003), the magnetopause (Haaland et al., 2004;Owen et al., 2004), the plasma sheet (Dewhurst et al., 2004;Runov et al., 2005), the cusp region (Taylor et al., 2004;Zong et al., 2004), and the plasmapause/plume (Darrouzet et al., 2004(Darrouzet et al., , 2006)).The timing method has also been modified (Zhou et al., 2006a;Zhou et al., 2006b) to determine the orientations and the velocities of 2-D structures, such as flux ropes.Furthermore, it should be noted that there are modified versions of the timing method to consider the acceleration term and/or the curvature term (Chanteur, 1998;Dunlop et al., 2002;Haaland et al., 2004), although in this paper, we will focus only on the classical version (with a planar discontinuity moving at a constant speed) of the method described above.
In order to avoid the improper usage of the timing method, the estimation of the errors both in the normal direction and in the velocity becomes very important.However, only a few attempts have been made on this issue (Chanteur, 1998;Cornilleau-Wehrlin et al., 2003;Vogt et al., 2008).One of the most extensive studies was made by Knetter (2005), who applied the timing method on fast solar wind discontinuities and analyzed the parameters affecting the error of the timing method.As was noted in their works, the detected instantaneous satellite distances D αβ contain very minor uncertainties, and the error of the timing method mostly arises from the uncertainties in the determination of the six time delays t αβ .
In this paper, we examine the error introduced in the time delay determination, which is assumed to be normally distributed with zero mean and a standard deviation of 1 s (in comparison with the typical time-resolution of 4 s which is the Cluster spin period).This random error, being treated as the only error source, is embedded in the performance of the timing method on a set of simulated Cluster crossings over a discontinuity moving at a constant velocity, to compare the resulting normal directions and velocities with the modeled ones.We find that the error depends not only on the shape of the Cluster tetrahedron, but also on the speed of the discontinuity compared with that of the Cluster constellation.For those cases in which the discontinuity moves at a lower speed, the error would be significantly enlarged, which requires the application of timing method with special caution, to the slow-moving discontinuities such as the plasmaspheric plumes.

Accuracy of the timing method
For the Cluster mission, the tetrahedron composed by the four spacecraft has the configurations varying continuously along its orbit around the Earth.Typically, the tetrahedron appears to be fairly regular at the apogee of the Cluster constellation, but it is significantly stretched at the perigee, which can be also seen in Fig. 16.1 of Robert et al. (1998a).
The best way to describe the geometry of the Cluster tetrahedron is to use a 2-D geometric factor, as the combination of the elongation parameter E and the planarity parameter P (Robert et al., 1998b).Both of the parameters are derived from the eigenvalues of the tetrahedron volumetric tensor, which can be also determined by the six interspacecraft distance vectors D αβ (for more detail, see Robert et al. (1998b) and Harvey (1998)).The values of E and P vary from 0 to 1: for a regular tetrahedron, E and P are both zero; if the tetrahedron is stretched, the value of E will increase, until E=1 when the satellites are on a straight line; if the tetrahedron is squashed, the value of P will increase, and the satellites will lie in a plane if P equals 1.
To study the relationship between the tetrahedron geometry and the accuracy of the timing method, we first construct a "homogeneous tetrahedron reservoir," which contains 1600 tetrahedra with a wide variety of configurations.In generating the reservoir, the procedures described in Robert et al. (1998b) are applied to ensure that the E and P values of these tetrahedra cover the whole E-P plane homogeneously.
As the second step, we simulate a planar discontinuity moving along its normal direction (the −x-direction) with the speed of 50 km/s, and we have all of the tetrahedra moving in the +x-direction, with a relatively smaller speed of 5 km/s.The time delays between each spacecraft pair when they observe the discontinuity, along with their relative position vectors, can thus be unambiguously determined.Here we further normalize the size of each modeled tetrahedron to have the averaged absolute value of the six corresponding time delays equal to 30 s, which is the typical value and of practical interest in the timing method application to the real Cluster data: in those cases with much longer time delays, any real discontinuities would evolve significantly which might invalidate the timing method, while those cases with much shorter delays would suffer from the quantization effects owing to the limited time resolution of the data (typically 4 s which is the spin period of the Cluster satellites).
As we discussed before, there are further error sources in the determination of the six time delays, e.g., the nonplanar properties of the discontinuity as well as the superposed structures and/or fluctuations.In fact, the classical way to obtain the time delays requires a process of maximizing the cross-correlation function, and the uncertainties in the resulting time delays, according to Eq. (1.7) of Sonnerup et al. (2008), depend on several factors: the number of data points, the obtained maximum correlation coefficient, the average slope square of the signal, and the average magnitude square of the signal deviation from its average value.Furthermore, the cross-correlation process is somewhat arbitrary: the data window selection and the interpolation method selection may both produce errors (Song and Russell, 1999).To reduce these errors, a more sophisticated technique, called the Local Wavelet Correlation (Soucek et al., 2004), was also presented.However, as was admitted by Soucek et al. (2004), the error cannot be fully eliminated unless the signals are identical up to a shift in time.
Here, as an example, we assume a normally distributed random error with zero mean and a standard deviation of 1 s, for each spacecraft pair, to adjust the six "theoretical" time delays.As the interspacecraft distance vectors D αβ remain unchanged, the relative position vectors R αβ between each two satellites can be accordingly adjusted based on the resulting time delays and the modeled D αβ values, which is important as the uncertainties in the determination of the time delay would also produce uncertainties in the relative positions of the corresponding spacecraft pair.The timing analysis can thus be performed through Eq. (2) to calculate the normal direction and the normal velocity of the certain discontinuity.
The last step is to compare the results obtained by the timing method with the modeled ones (50 km/s in the −xdirection), to estimate the reliability of the timing method.For each of the 1600 tetrahedra, the second step is performed 1000 times with the time delays adjusted with random errors, and the averaged relative error of the normal velocity component can be expressed as where V T and V M represent the calculated normal component of the velocity and the modeled one (50 km/s), respectively.
The relation between the relative error of the normal velocity and the tetrahedron geometry is shown in the left panel of Fig. 1.The sizes and the colors of the circles indicate the values of the relative errors: for those cases with the error less than 10%, the circles are very small, while the largest circles correspond to a relative error of greater than 100%.Note that the circles are sorted by size before being plotted, so there are no small circles hidden behind larger ones.
The uncertainty in the determination of the normal direction, on the other hand, is shown in the right panel of Fig. 1.The uncertainty is inferred by the 90-percent confidence interval (CI) of the angle between the calculated normal direction and the modeled one (along the −x-axis).The smallest circles correspond to the CIs of less than 10 • , which means that we are 90% confident that the angular error is less than 10 • , while the largest ones suggest the uncertainties of greater than 90 • .
It is clear that the errors are generally small except for the cases when the four spacecraft lie approximately along a straight line or within a plane.In these two extreme cases, the component equations in Eq. (1) degenerate to one and two component equations, respectively, leaving some unknowns undetermined.The result, therefore, is not surprising: for the cases of E close to 1, the separation vectors between each pair of satellites are in very similar directions, and therefore can hardly provide the information in the perpendicular directions; for the cases of high P, the separation vectors are generally coplanar, and the speed normal to the plane cannot be determined.
Note that here we use the discontinuity speed of 50 km/s which is ten times greater than the speed of the Cluster constellation (5 km/s).However, if the discontinuity moves much slower, say, 0.5 km/s (still in the −x-direction), the uncertainties of the timing method can be calculated to be quite different, as shown in Fig. 2.
The relative errors on the normal velocities shown in the left panel of Fig. 2, although not solely dependent on the values of E and P because of the differences in the orientations of these tetrahedra, are always greater than 20% even in the cases of regular tetrahedra.For E or P values greater than 0.8, the errors are generally greater than 100%, and even as large as 1000% for E or P close to 1. Similarly, the uncertainties in the normal directions (right panel of Fig. 2) are also much greater than those in the previous case.
The resulting much greater relative errors can be understood as an error amplification effect produced in the transformation process between two reference systems.If we are in the frame moving along with the Cluster constellation, as was suggested before, the relative position vectors R αβ in Eq. ( 2) could be replaced by the interspacecraft distance vectors D αβ .Therefore, if the previous two test runs are performed in the Cluster's frame, it is not necessary to adjust the D αβ values to obtain R αβ , and the latter run differs from the former one by having a discontinuity velocity of 5.5 km/s ten times smaller than that in the former run.Note that the D αβ values in the latter run are also ten times smaller, which is produced in the normalization process to have the averaged absolute value of the six time delays equal to 30 s.
The factors of ten, however, can be eliminated by substituting the V and D αβ values to Eq. ( 2), which suggests the relative errors ε 0 in the determination of the discontinuity velocity be the same for the two test runs if both performed in the Cluster's frame.As we know that the transformation to another frame of reference would not change the value of the absolute error, we have where V D0 , V D1 , and V C1 represent the discontinuity velocity in the frame moving with the Cluster constellation, the discontinuity velocity in the Earth's frame of reference, and the Cluster velocity in the Earth's frame, respectively.Therefore, if the discontinuity moves much slower than the Cluster constellation in the Earth's frame, the satisfaction of Eq. ( 4) results in significant amplification of the relative error ε 1 in the Earth's frame.In other words, since the velocity of the discontinuity should be cast in the Earth's frame for any physical usefulness, a small error in the velocity derived in the satellite frame can be translated into a very large relative error if the velocity of the discontinuity is small.(2004).The densities for an inbound and an outbound are shown with a factor-of-10 shift in the upper traces.The labeled speeds of the moving structures were derived from the timing method.
The error amplification effect, therefore, suggests that the uncertainties of the timing method may be large not only in the cases when the Cluster tetrahedron is greatly elongated or flattened, but also when the discontinuity moves much more slowly than the Cluster constellation.One example of the appearance of possible greater uncertainties occurs when the Cluster satellites move relatively faster with elongated configurations near their perigees and encounter slow moving plasmaspheric structures, and the applications of the timing method under these conditions should be carried out with special caution.

Example of applications
On 11 April 2002, between 04:15 and 06:30 UT, the Cluster satellites experienced a plasmasphere crossing and observed a plasmaspheric plume outside the plasmasphere (Darrouzet et al., 2004).Darrouzet et al. used the EFW (Gustafsson et al., 1997) and WHISPER (Décréau et al., 1997) measurements to obtain the electron densities and applied the timing method to determining the normal velocity of the plume.
However, it can be clearly seen in Fig. 3 that there are fluctuations superposed on the discontinuity, and the time differences between the observations of different satellites continuously evolve in time, especially during the Cluster inbound pass.Therefore, it might be relatively hard to determine the time delays accurately.
The small error in the determination of the time delays, on the other hand, can be significantly enlarged in the application of the timing method.As is shown in Fig. 3, the estimated velocity of the discontinuity is typically one order smaller than that of the Cluster constellation (4.6 km/s), and the averaged absolute value of the time delays between each spacecraft pair is around 30 s (see Fig. 5 of Darrouzet et al., 2004).Therefore, Fig. 2 can be used to estimate the error in the application of the timing method.During this time period, the Cluster tetrahedron elongation parameter E was 0.8, and the planarity parameter P was 0.2, which suggests that the uncertainties in the velocity determination could be as large as 100%, as the direction determination might be more reliable.
The uncertainty levels of ∼100%, therefore, suggest that one should avoid directly comparing two resulting velocities during different time intervals, although it is safe to claim that the plume moves very slowly.In other words, it is quite possible that some of the conclusions drawn by Darrouzet et al. (2004), e.g., the velocities measured along the inbound pass being larger than those along the outbound pass, are within the error bar of the data analysis if the obtained time delays have the uncertainties of ∼1 s or more.Therefore, it may be necessary to carefully reexamine the case, for example, by using the Eq.(1.7) of Sonnerup et al. (2008), to make sure that all of the time delays are accurately determined with the errors much smaller than 1 s despite the effect of the irregularities and the fluctuations shown in Fig. 3.
Besides estimating the time delay errors based on Eq. (1.7) of Sonnerup et al. (2008), another way to assess the reliability of the timing method is to calculate all of the four time delay residuals Fig. 4. Same format as in Fig. 2, but only those cases with all of the four residuals smaller than 1 s are selected in the simulation.
to higher reliabilities of the timing method.Here we perform the same simulation steps, again with the modeled discontinuity moving at 0.5 km/s in the −x-direction and the time delay standard deviation of 1 s, but select only those cases with all of the four residuals smaller than 1 s, to calculate the uncertainties of the timing method.The results are shown in Fig. 4, which may be directly compared with the uncertainties shown in Fig. 2.
It is shown that the uncertainties in Fig. 4 are reduced, especially those in the velocity determination.Given the E and P values of 0.8 and 0.2, respectively, the uncertainties in the velocity determination would be ∼30−40%, which are much smaller than the ∼100% ones suggested in Fig. 2. In other words, to ensure the appropriate usage of the timing method in this particular Darrouzet et al. (2004) event, one has to make sure that either the time delays are accurately determined (with uncertainties much smaller than 1 s), or all of the time delay residuals are adequately small (less than 1 s).Without these reexaminations, as we suggested before, some of the Darrouzet et al. (2004) conclusions, e.g., the plume velocities during the inbound pass being larger than those along the outbound pass, may be less convincing as they are likely to be within the error bar of the timing method.

Summary
In this study, we have performed simulations to estimate the uncertainty of the multi-spacecraft timing method, to investigate its dependence on the configuration of the spacecraft tetrahedron and the speed of the observed discontinuity in comparison with the constellation speed.By limiting the error sources of the timing method to the random errors produced in the determination of the time delays between each spacecraft pair, we have found that the uncertainty of the timing method could be very significant if the tetrahedron is largely elongated or flattened.Furthermore, if the discontinuity moves much more slowly than the constellation, a small error in the determination of the speed in the satellite frame of reference can be translated to a large relative error.The error estimation, therefore, provides a way for the users of the timing method to evaluate their result based on the estimated uncertainty levels in the determination of the time delays.For instance, in the study of slow-moving discontinuities, such as the plasmaspheric plumes (Darrouzet et al., 2004), even a minor error in the determination of the time delays may result in huge errors, a result that suggests that the timing method should be used with caution and the predictions from an analysis should be made with due attention to the uncertainties in the method.

Fig. 1 .Fig. 2 .
Fig.1.Influence of the tetrahedron configuration, parameter P and E, on the error of the timing method, denoted by the size and color of the circles, when the modeled discontinuity moves ten times faster than the Cluster constellation.Left and right panels correspond to the relative error in the determination of the normal velocity component and the 90-percent confidence interval of the angle between the normal direction and the modeled one, respectively.

Fig. 3 .
Fig.3.The plasmapause/plume crossings of the Cluster constellation, adapted fromDarrouzet et al. (2004).The densities for an inbound and an outbound are shown with a factor-of-10 shift in the upper traces.The labeled speeds of the moving structures were derived from the timing method.