Accuracy of multi-point boundary crossing time analysis

Recent multi-spacecraft studies of solar wind discontinuity crossings using the timing (boundary plane triangulation) method gave boundary parameter estimates that are significantly different from those of the well-established single-spacecraft minimum variance analysis (MVA) technique. A large survey of directional discontinuities in Cluster data turned out to be particularly inconsistent in the sense that multi-point timing analyses did not identify any rotational discontinuities (RDs) whereas the MVA results of the individual spacecraft suggested that RDs form the majority of events. To make multi-spacecraft studies of discontinuity crossings more conclusive, the present report addresses the accuracy of the timing approach to boundary parameter estimation. Our error analysis is based on the reciprocal vector formalism and takes into account uncertainties both in crossing times and in the spacecraft positions. A rigorous error estimation scheme is presented for the general case of correlated crossing time errors and arbitrary spacecraft configurations. Crossing time error covariances are determined through cross correlation analyses of the residuals. The principal influence of the spacecraft array geometry on the accuracy of the timing method is illustrated using error formulas for the simplified case of mutually uncorrelated and identical errors at different spacecraft. The full error analysis procedure is demonstrated for a solar wind discontinuity as observed by the Cluster FGM instrument.


Introduction
The analysis of discontinuities in space plasmas has received a lot of attention since the beginning of the space age.In the case of a planar discontinuity moving at constant velocity, its orientation can be estimated from magnetic field measurements using the minimum variance analysis (MVA) technique (Sonnerup and Cahill, 1967;Sonnerup and Scheible, 1998) based on the conservation law for magnetic flux.The MVA framework can also be applied to electric field measurements or plasma data if other conservation laws are used (e.g.Sonnerup et al., 2008), and it further allows to take into account physical or geometrical constraints.
Applications of the MVA technique to solar wind discontinuities were recently challenged in a comprehensive study based on data from ESA's Cluster satellites (Knetter et al., 2004;Knetter, 2005).Such multi-spacecraft missions offer an independent road to boundary parameter estimation through a crossing time analysis that effectively yields a boundary plane triangulation technique or, in brief, the socalled timing method.T. Knetter and colleagues found that the discontinuity normal vectors obtained with the timing approach differ from the MVA normals.Furthermore, the MVA normals at the individual spacecraft are often mutually inconsistent even though previously used quality criteria such as the intermediate-to-minimum eigenvalue ratio were met.The discrepancy became less pronounced when important quality thresholds like the required eigenvalue ratio and/or the change in magnetic field direction across the discontinuity were raised, and then the results turned out to be also more consistent with the timing normals.
The discrepancy of discontinuity normal vector estimates using the two principal methods has crucial implications on the physical interpretation of the measurements.In accordance with previous single-spacecraft studies of solar wind discontinuities (e.g.Tsurutani and Smith, 1979;Neugebauer et al., 1984;Lepping and Behannon, 1986;Söding et al., J. Vogt et al.: Accuracy of multi-point boundary crossing time analysis 2001), Knetter's MVA results gave significant normal magnetic field components B n for a large fraction of the solar wind directional discontinuities (DDs), and hence put them into the category of rotational discontinuities (RDs).The timing analysis, on the other hand, led to consistently small B n 's so that the DDs were either considered to be tangential discontinuities (TDs) or could not be clearly classified (Knetter et al., 2004;Knetter, 2005).The inconsistencies in the relative distributions of RDs and TDs based on different discontinuity analysis methods had been noted already earlier in a study of Horbury et al. (2001) using data from Geotail, IMP 8, and Wind.The state of the problem was summarized by Neugebauer (2006) who also reexamined a large number of solar wind discontinuities observed by the ISEE 3 spacecraft with inconclusive results, and discussed possible physical mechanisms.
The discrepancy in the distributions of boundary orientations obtained with different methods, and the resulting ambiguity in RD/TD classification may, for brevity, be termed discontinuity analysis inconsistency.In 2009, a team has formed at the International Space Science Institute in Bern to investigate the problem in detail.The present paper addresses one of the main team objectives, namely, the construction of a rigorous error analysis scheme for the timing method.Numerical experiments on this issue were carried out by Zhou et al. (2009).The purpose of our study is a fully analytical treatment of the problem.We start by introducing the reciprocal vector formalism in Sect. 2. In Sect.3, we show how the main variants of the multi-point timing method all allow to write the boundary slowness vector as a linear combination of crossing times and reciprocal vectors.The error analysis in Sect. 4 starts from the slowness vector formula to quantify the mean square errors in boundary orientation and speed in terms of the spacecraft configuration, the estimated parameters, and the uncertainties in crossing times and spacecraft positions.The crossing time uncertainty is studied further in Sect. 5.In Sect.6, the complete chain of boundary parameter estimation and error analysis is demonstrated using Cluster FGM measurements across a solar wind discontinuity.We conclude in Sect.7 by summarizing the important steps of the error analysis procedure.

Reciprocal vectors in multi-spacecraft analysis
To facilitate the use of vectors in the definition of dyads and for the purpose of matrix multiplication, we adopt the notation conventions described, e.g. in Paschmann and Daly (1998): vectors a,b,c,... are always understood as column vectors.They can be turned into row vectors by means of transposition denoted by the superscript T, e.g. a T .The hat symbol • indicates unit vectors, matrices are typeset in upright bold, and I denotes the identity matrix.
The positions of the spacecraft are given by r α (α = 1,...,S).Since we are mainly interested in the Cluster mis-sion with four spacecraft, we focus on S = 4 but consider the more general case where possible.Relative position vectors are written in the form r αβ = r β −r α .If the origin of our coordinate system coincides with the mean position (mesocenter) r * = (1/S) α r α of the spacecraft array, then obviously r * = 0 and the reference frame is called mesocentric.Spacecraft position vectors in a mesocentric frame are denoted as r * α , hence α r * α = 0.The so-called position tensor is defined through Throughout this paper, we assume that the number of spacecraft is at least four, and that they form a three-dimensional configuration that does not degenerate into a plane or a line.
Then the position tensor R * is a non-singular matrix (Vogt et al., 2008), and its inverse is the reciprocal tensor defined through where the (generalized) reciprocal vectors are given by Note the identity In the case S = 4, the vectors q α coincide with the reciprocal vectors of the spacecraft tetrahedron defined through (5) (Chanteur, 1998) where (α,β,γ ,λ) must be a cyclic permutation of (1,2,3,4).In this case the symbol K is used for the reciprocal tensor, i.e.
Useful algebraic identities for S = 4 are where δ αβ denotes the Kronecker delta symbol.
Reciprocal vectors can also be defined for three-spacecraft configurations, see Vogt et al. (2009).

Boundary parameter estimation from crossing data
The problem of computing boundary parameters from the crossing times and the positions of a multi-spacecraft array has been studied by a number of authors (e.g.Burlaga and Ness, 1969;Russell et al., 1983;Dunlop et al., 1988;Ann. Geophys., 29, 2239-2252, 2011 www.ann-geophys.net/29/2239/2011/Mottez and Chanteur, 1994;Schwartz, 1998;Dunlop and Woodward, 1998;Harvey, 1998;Soucek et al., 2004;Haaland et al., 2004;Vogt et al., 2008;Zhou et al., 2009).In its simplest and most popular form, the underlying model assumes a planar structure that varies only in the direction of the boundary unit normal vector n and propagates at a speed V along n relative to the spacecraft array.The model parameters n and V are conveniently combined into the boundary slowness vector The vector m can then be determined from a set of linear equations where the crossing times constitute the known data, and the spacecraft position vectors form the coefficient matrix that has to be inverted.There are essentially three variants of this procedure using (a) relative crossing data with respect to one reference spacecraft (e.g.Russell et al., 1983;Knetter et al., 2004), (b) all relative crossing data that are available (e.g.Soucek et al., 2004;Zhou et al., 2009), or (c) absolute crossing times and spacecraft positions (e.g.Haaland et al., 2004).To facilitate the comparison of these three options, we only consider four-spacecraft configurations and make use of algebraic identities for the tetrahedral reciprocal vectors k α (α = 1,2,3,4).Note that options (b) and (c) can be easily generalized to configurations with more than four spacecraft by means of the generalized reciprocal vectors q α .

Crossing data relative to one reference spacecraft
To uniquely determine the three-component slowness vector m, knowledge of three relative position vectors r ρα and the corresponding differences t ρα = t α − t ρ in crossing times are sufficient.Here the subscript ρ denotes the reference spacecraft, and α = ρ.The three conditions are thus V t ρα = r ρα • n which can be divided by the speed V to obtain Since any subset of three reciprocal vectors k β form a basis of three-dimensional space, the slowness vector can be expressed in the form m = β( =ρ) C ρβ k β .The symbol β( =ρ) indicates summation over all indices β except for β = ρ.To determine the three coefficients C ρα , we insert this ansatz into Eq.( 10) to obtain where the identity k β • r ρα = δ βα − δ βρ has been used.The slowness vector is thus given by The restriction α = ρ may be dropped by setting t ρρ = 0, and then m = α k α t ρα . (13)

Using all available relative crossing data
We can incorporate all available crossing data into the boundary analysis by assigning the role of reference spacecraft to each one of them in turn, then compute slowness vector estimates m (ρ) using Eq. ( 12), and, finally, average the results to obtain m = (1/4) ρ m (ρ) .The result can be easily rearranged to yield: where t ρρ = 0 as before.This form is completely equivalent to the least-squares result obtained by Harvey (1998) for a symmetrical treatment of relative crossing data.Using our notation, his Eq.(12.13) translates to and with t ρα = −t αρ the least-squares formula can be written in the form of Eq. ( 14).

Absolute crossing times
If a boundary crossing is unambiguously identified in the data of all spacecraft without reference to another, e.g. as the center time of a jump in one variable, or through correlation with a prescribed model profile, it is more appropriate to think in terms of absolute crossing times t α rather than their relative counterparts.We insert t ρα = t α − t ρ into Eq.( 14) and use the identity α k α = 0 to obtain the formula The result holds for arbitrary time offsets.
The following version of Eq. ( 17) was derived by Chanteur (1998) from spatial interpolation theory, as well as by Harvey (1998) and Vogt et al. (2008)  (1/4) α t α , was chosen as a reference point to center the time axis.It is straightforward to show that the linear combination of relative crossing times in Eq. ( 14) is equal to the respective absolute crossing time in the time frame centered at t * as (1/4) ρ t ρα = t * α .The case of an accelerated planar discontinuity was considered by Chanteur (1998), see Sect.14.5.2 of that publication.

Comments on implementation
In all variants of the boundary triangulation approach discussed above, the slowness vector m is given in terms of the crossing times and the reciprocal vectors.It is important to note that the latter have to be computed from the relative crossing position vectors, i.e. from r αβ = r β (t β ) − r α (t α ).
Here the spacecraft trajectories r α (t) and r β (t) have to be evaluated at the (absolute) crossing times t α and t β , respectively.This means that the set of relative crossing times alone is not sufficient to determine the solution uniquely but must be supplemented by at least one absolute time datum such as the crossing time of one reference spacecraft.
The discussion in the following Sect.4 shows that the relative crossing times which directly enter the slowness vector formulas should be known very precisely to yield accurate boundary parameter estimates.If δt denotes a reference value for the error in relative crossing times, then the additional absolute datum required for obtaining the crossing positions can tolerate an uncertainty δt that is somewhat larger than δt.The resulting positional inaccuracies are of the order n•u αβ δt where u αβ are the relative spacecraft velocities, corresponding to timing uncertainties n • u αβ δt /V .For the Cluster mission, n • u αβ /V is a very small quantity (of the order of 10 −3 or less).This means that as long as δt and δt are of the same order, we can disregard the contribution of the uncertainty in the additional absolute crossing time in the following error analysis.
If, however, instead of the actual crossing position vectors an instantaneous spacecraft configuration is used to compute the reciprocal vectors, another source of error comes into play that can no longer be neglected.Such a procedure yields additional timing inaccuracies of the order n • u αβ t αβ /V .Since the relative crossing times t αβ can be several orders of magnitude larger than their errors δt, the instantaneous configuration approximation may introduce significant inaccuracies and thus should be avoided.

Error analysis and array geometry
In this section we present the first part of the error analysis scheme for the timing approach to boundary parameter estimation.We give formulas for the errors in boundary orientation and speed, and assume that the primary uncertainties in crossing times and the positional inaccuracies are given in the form of error covariance matrices.The crossing time errors are quantified further in Sect. 5.

Analysis framework and general error formulas
The quality of boundary parameter estimation suffers from inaccuracies in crossing times and spacecraft positions.The problem was addressed, e.g. by Dunlop and Woodward (1998); Knetter (2005); Zhou et al. (2009).Analytical error formulas were derived by Gérard Chanteur for the absolute crossing times approach in various contexts (e.g.Chanteur, 1998Chanteur, , 2000Chanteur, , 2003;;Cornilleau-Wehrlin et al., 2003;Vogt et al., 2008), and a partial summary of Chanteur's results was also given in the PhD thesis of Knetter (2005).The analysis rests on Eq. ( 17) which is repeated here for convenience: Since Eq. ( 14) exhibits the same structure with t α being replaced by the expression (1/4) ρ t ρα , the line of reasoning can in principle be applied also to this case.
With only mild assumptions on error correlations, Chanteur arrived at the following formula for the unit normal covariance matrix and wrote the error in boundary speed V in the form see Eqs. (4.33) and (4.34) in Vogt et al. (2008).Here and δr δr T γ denotes the positional error covariance for spacecraft number γ .For the Cluster mission, the positional error covariance matrices are available through the Cluster Active Archive (Volpp and Sieg, 2010).
If ê is an arbitrary unit vector perpendicular to n, then êT δ nδ nT ê is a quadratic measure of the angular uncertainty of n in the plane spanned by the two vectors ê and n.If small compared to unity, this measure can be associated with the opening half angle (in ê-direction) of an elliptical cone of uncertainty for the boundary unit normal vector: where The formulas ( 19)-( 25) in this subsection provide the general framework for a rigorous error analysis of the timing method.They give the uncertainties in boundary orientation and speed in terms of the error covariance matrices of crossing times and spacecraft positions.For illustrating purposes, we proceed with order-of-magnitude estimates and simplifications that allow to highlight the influence of the spacecraft configuration and the boundary parameters on the overall accuracy.

Primary errors and relative importance
The error term C 1 is controlled by uncertainties in crossing time estimates, and C 2 originates from positional uncertainties (note that the errors in relative position vectors matter here).To assess the relative importance of the two error contributions, we assume crossing time uncertainties ∼ δt, relative positional errors ∼ δr, boundary speeds ∼ V , and inter-spacecraft separations ∼ L.
, and thus For the Cluster mission, δr is in the kilometer range or below (Volpp and Sieg, 2010).This has to be compared with the product V δt that is usually much larger in the geospace context.The Cluster FGM instrument with its high but finite time resolution of about 0.05 s (in normal mode) cannot be expected to yield discontinuity time uncertainties significantly smaller than δt ∼ 0.1 s in the presence of noise.Hence for boundary speeds of the order 100 km s −1 and above, C 2 is much smaller than C 1 , and the effects of positional errors can be safely neglected against those of timing uncertainties.
For instruments operating at spacecraft spin resolution, the timing error is expected to be in the range of the spin period, then this statement holds even for much smaller boundary speeds of order 10 km s −1 .If the timing uncertainties and positional covariances are the same for each spacecraft and mutually uncorrelated, the slowness error covariance matrix simplifies to where K is the reciprocal tensor.When in addition the positional error matrix is isotropic, i.e. then If positional inaccuracies can be ignored, the term in square brackets reduces to (δt) 2 .In the remainder of the present section this approximation is employed to study the influence of tetrahedron geometry on boundary parameter estimation accuracy.To recover the complete formulas with the effects of positional inaccuracies included, one may replace

Influence of the spacecraft array geometry
The error formulas for the simplified case are where ê is perpendicular to n, and δV /V is used as a shorthand notation for (δV ) 2 /V , i.e. the relative rms error in boundary speed.To study the influence of the spacecraft array geometry on boundary parameter accuracy, we write the reciprocal tensor in terms of its eigenvectors and the tetrahedron geometric parameters planarity P , elongation E, and the rms inter-spacecraft distance L. The expression for the resulting quadratic form c T Kc is given in Appendix A where also the parameters P , E, and L are defined.Note that P = 1 if all spacecraft are in one plane and E = 1 if they lie on a straight line (string-of-pearls configuration), and that spacecraft configurations close to an ideal tetrahedron correspond to small values of planarity and elongation: P ≈ 0 and E ≈ 0.
In the quadratic form c T Kc, we set c = n to compute δV , and c = ê for δθ.Hence the errors depend on the orientation of the spacecraft tetrahedron, on the length scale L, and on the shape parameters E and P .To assess the full range of possible errors, we assume that the boundary unit normal vector is aligned with the three eigenvectors ê(n) one by one.

Boundary unit normal aligned with the direction of elongation
The direction of elongation is given by the eigenvector ê(1) to the largest eigenvalue R (1) * of R * which corresponds to the smallest eigenvalue of K.If n = ê(1) , then ê = ê(2) yields the minimum angular uncertainty, and ê = ê(3) its maximum value: The error in boundary speed is given by This is the best case (highest accuracy) for the speed estimation but the worst case (lowest accuracy) for the computation of the boundary normal.

Boundary unit normal aligned with the direction of planarity
The direction of planarity is given by the eigenvector ê(3) to the smallest eigenvalue R (3) * of the position tensor, or the largest eigenvalue of the reciprocal tensor.If n = ê(3) , then ê = ê(1) yields the minimum angular uncertainty, and ê = ê(2) its maximum value: The error in boundary speed is given by This is the best case (highest accuracy) for the computation of the boundary normal and the worst case (lowest accuracy) for the speed estimation.

Boundary unit normal aligned with the eigenvector to the intermediate eigenvalue
The eigenvector ê(2) belongs to the intermediate eigenvalues both of the position tensor and the reciprocal tensor.If n = ê(2) , then ê = ê(1) yields the minimum angular uncertainty, and ê = ê(3) its maximum value: The error in boundary speed is given by This is the intermediate case for both boundary normal and speed estimation.
The results for the three cases presented above can be combined to yield representative errors for a particular spacecraft geometry.We average the squares of both the angular uncertainty and the error in speed to obtain where (41) The meaning of the term trace(K) = 4 α=1 |k α | 2 in the context of error amplification was recognized by Vogt and Paschmann (1998) in their study on the accuracy of spatial derivatives, and its importance was confirmed in the thorough analysis presented by Chanteur (2000) who further studied the dependence on L, P , and E. For further details the reader is referred to the original publications and to Vogt et al. (2008).Note that for planarity values close to one, the function A is well approximated as A (1−E) −1 (1−P ) −1 .

Reference error of the timing method
For geometrically ideal spacecraft configurations characterized by zero values of planarity and elongation, the directional inaccuracy δθ and the relative error in boundary speed δV /V are both given by For brevity, we refer to the term V δt/L as the reference error (of the timing method).
Figure 1 shows V δt/L as a function of inter-spacecraft length scale L and boundary speed V both in percent (blue solid contours, for δV /V ) and in degrees (red dashed contours, for δθ).The boundary speed values at the left y-axis are for δt = 0.1 s in which case the error formulas are For other values of δt, the numerical factors on the righthand side of the equation must be multiplied by δt/0.1 s.The numerical values of the boundary speeds for the case δt = 1 s have been added in Fig. 1 at the right y-axis for convenience.The smallest reference errors occur when L is large and both V and δt are small.In this sense magnetopause studies (boundary speed V ∼ 10 km s −1 ) using high-resolution Cluster FGM data (δt ∼ 0.1 s) from the year 2003 (L ∼ 5000 km) provide a best case scenario as directional inaccuracies could theoretically be as small as 0.01 deg, if in fact the magnetopause behaved as an ideal planar structure on the time scale of the transition (500 s) and on length scales close to one Earth radius.A worst case scenario for the reference error is the study of solar wind discontinuities (V ∼ 100 km s −1 ) using Cluster plasma measurements at spin resolution (δt ∼ 4 s) from the year 2002 (L ∼ 100 km) which gives a relative error in boundary speed of 400 %.High-resolution FGM data (δt ∼ 0.1 s) yield a value of 10 %.

Geometrical error amplification
Spacecraft array geometries that deviate from an ideal regular configuration are characterized by non-zero values of elongation and planarity.The errors δθ and δV /V are products of the reference error V δt/L and functions that depend Numerical experiments on the accuracy of the timing method were carried out by Zhou et al. (2009).They generated a reservoir of spacecraft tetrahedra with a homogeneous distribution in elongation and planarity, and then simulated crossings of planar discontinuities with timing errors that were identical at all four spacecraft and mutually uncorrelated.The resulting distributions of errors in boundary orientation and speed are shown as function of elongation and planarity in Figs. 1, 2, and 4 of Zhou et al. (2009).They compare nicely with the corresponding contour plots of the present study (Figs. 2 and 3), in particular with regard to the sharp increase in errors for values of elongation and planarity close to unity.We take this as a consistency check of our analytical error formulas that are easier to apply to actual data.

Crossing time errors
The timing method in boundary parameter estimation rests on the crossing times t α .Our error analysis scheme presented in the previous Sect.4 requires the crossing time error covariances δt α δt β as input parameters.The present section aims at quantifying these error covariances through a pattern www.ann-geophys.net/29/2239/2011/Ann.Geophys., 29, 2239-2252, 2011 matching approach: we study the similarity of a signal s with shifted versions of a pattern function p.For notational convenience, we write τ for the crossing time difference and refer to it also as the lag (time).An overbar ••• denotes the time averaging operation.The averaging window is assumed to be fixed with respect to the pattern function.Thus the window would have to move with the lag time τ if we defined the association measures in terms of the pair s(t) and p(t − τ ).Instead we choose to write the formulas using the pair s(t + τ ) and p(t) so that the averaging window is fixed and symmetric with respect to the time origin.

Association measures
A variety of association measures can be employed to quantify the similarity of a time series s and a pattern function p at a time shift τ .Here we choose the mean square deviation because it allows to derive analytical error formulas, and it is sensitive also to linear variations in the data.The latter statement is not true for Pearson's correlation coefficient where The coefficient γ is designed to measure linear correlation: if p(t) and s(t) are both linear functions, then γ (τ ) = 1 irrespective of the lag time τ .The correlation coefficient can be made sensitive to linear variations in the data by replacing s • and p • in the formula with its non-centered counterparts s and p.
Further association measures such as the mean absolute deviation |s(t + τ ) − p(t)| were tested by means of numerical experiments using synthetic data first, and then actual Cluster FGM measurements of solar wind discontinuities.In the noise-free limit the mean absolute deviation allows for an easier identification of crossing times than the mean square deviation.However, in the presence of noise or for actual measurements the differences turned out to be minor.
It may also be noted that mean deviation measures exhibit a misleading dependence on window width in graphical displays of association measures.This problem can be easily rectified by a multiplication with the number of data points in the averaging window to yield cumulative deviation measures (R. Wicks, private communication).The analytical error analysis given below is valid for both the mean square deviation as well as its cumulative counterpart.

Analytical error formulas
The following formula for the crossing time error covariances is derived in Appendix B: The underlying assumptions can be summarized as follows.
-Crossing times estimates are assumed to be based on the mean square deviations I α (τ ) = |s α (t + τ ) − p(t)| 2 between a pattern p and the signals s α at lag time τ .
-The pattern p(t) is expected to show a transition at the origin between approximately constant levels at both ends of the data window.Examples are an ideal step (Heaviside) function or a hyperbolic tangent profile.When timing is only relative, a windowed portion of a discontinuity crossing observed at a reference spacecraft ρ may also serve as a pattern function for the signal measured at another spacecraft α.
-The residuals h α are estimated through where τα, * is the lag time at the minimum of the mean square deviation I α .
-The residuals h α are assumed to be time-stationary random signals that are well characterized by their means and their correlation functions.Angular brackets ••• denote the ensemble averaging with respect to the residuals.
-An overbar ••• indicates time averaging, T w is the time interval used for averaging, and N is the number of data points in the time window.
-The time difference parameter in the sum runs from −T w to +T w , at least in principle.In boundary analysis practice, when pattern functions p(t) are characterized by constant levels left and right of the transition, the factor p (t)p (t + ) effectively cuts off the summation.
Note that although both and τ denote time differences, we prefer to use two different symbols as they appear in different contexts.
A second version of the error covariance formula can be obtained through normalization of the correlation functions.We define As explained in Appendix B, the factor (1 − | |/T w ) comes into play through the autocorrelation function p (t)p (t + ) and was thus included in the definition of G( ).The mean square (time) average of h α (t) is the minimum value of the mean square deviation I α : The crossing time error covariances can thus be written in the form: All information on the right-hand side of this equation can be constructed from the measurements.Implementation of the crossing time error covariance formulas is discussed below in Sect.7 where the complete error analysis scheme is summarized.
In some special cases the sum in Eq. ( 52) collapses to a single term.
-If the pattern function p(t) is an ideal step (Heaviside) function, the transition between the two states at both ends of the data window occurs within a sampling interval, then p (t) = 0 for t = 0 and thus G( ) = 0 for = 0.The error formula then simplifies to -If the residuals h α (t) and h β (t) can be represented as mutually uncorrelated white noise, then H αβ ( ) = 0 for α = β or = 0, so its only non-zero value is H αα (0) = 1, and we obtain The result is consistent with a formula derived by Alexander Khrabrov (private communication; see also Eq. 1.7 in Sonnerup et al., 2008) for this idealized case.
In nonlinear and turbulent space plasmas such as the solar wind, however, such correlations in the measurements cannot be disregarded.

Example
To demonstrate the error analysis scheme presented in this paper, Cluster FGM measurements of a solar wind discontinuity are considered.After computing the boundary parameter estimates, we proceed with the crossing time error covariance formula (52).Particular emphasis will be on the functions G and H αβ that quantify the effect of correlations in the set of residuals.Then the slowness vector covariances are computed and, finally, the errors of boundary speed and direction.

Fig. 4.
Crossing time error analysis using a tanh pattern function: FGM data (solid) and pattern functions (dashed, magenta) used for illustrating the crossing time error formulas.Shown are the B y components of the interplanetary magnetic field measured by the four Cluster spacecraft (S/C 1: black, S/C 2: red, S/C 3: green, S/C 4: blue) when they crossed a directional discontinuity shortly before and after 19:11:30 UTC.

Boundary parameter estimates
The data are shown in Fig. 4. We are looking at the magnetic field signature of a directional discontinuity in the solar wind crossed by all four Cluster spacecraft at around 19:11:30 UTC on 3 February 2003.Superposed are shifted versions of the hyperbolic tangent pattern function where the parameters B off , B amp , and T dis were obtained from visual comparison with the data.In principle, the model parameters could also be determined through a least-squares fit to a composite profile.The length of the pattern time window was chosen to be 10 s.The actual crossing times t α were determined from minima of the mean square deviations I α (τ ), see Eq. ( 45).Relative to the reference time 19:11:00 UTC, their numerical values are 35.345s, 29.324 s, 30.216 s, and 26.068 s for S/C 1 through S/C 4, respectively.The numbers are given with three digits as the sampling interval is T sam = 44.6 ms.
Then the reciprocal vectors k α were computed from the spacecraft positions at t α , and, finally, the boundary slowness vector m from Eq. ( 17): m Boundary speed and normal unit vector: The geometrical parameters of the spacecraft configuration are L = 3300 km, P = 0.19, and E = 0.27.Geometrical error amplification is small.The simplified error analysis presented in Sects.4.3-4.5 suggests that the relative error in boundary speed should be of the order of percent, and the directional error should be in the range of one degree.

Crossing time error covariances
To evaluate the crossing time error covariance formula (52), the functions G and H αβ have to be computed.For the hyperbolic tangent pattern function p(t) chosen here, the (normalized) autocorrelation function G( ) of the time derivative dp/dt drops to small values < 0.1 beyond ±50 samples away from the origin (| | > 2.2s), thus effectively limiting the summation in Eq. ( 52).
The residuals h α (t) in the upper panel of Fig. 5 were constructed according to Eq. ( 48).Note that by construction the residuals are centered around the actual crossing time.Welldeveloped structures in the diagram already suggest that the ideal white noise model does not apply.Products of correlation functions G( ) and H αβ ( ) are displayed in lower panel of Fig. 5 for α = 1.The results for α = 2,3,4 are qual- itatively very similar.The numerical values of the sum in Eq. ( 52) are the integrals under the curves.In our example they are in the range from 12 (α = 2,β = 3) to 31 (α = β = 1).
All variables on the right-hand side of the crossing time error formula ( 52) are now available and the crossing time error covariance matrix δt α δt β can be computed which in turn yields the slowness error covariances and all other uncertainties through the formulas in Sect. 4. The resulting matrix elements are given in Table 1.The diagonal elements (α = β) of the table can be understood as the square inaccuracies of the timing method at the respective spacecraft.
Taking the square root yields values δt 2 α around 0.1 s, i.e. in the range of 2-3 sampling intervals.
Note that also the off-diagonal elements (α = β) of the error covariance matrix δt α δt β are positive which reflects the fact that the residuals h α are positively correlated.In other words, the residuals share common substructures.When these common features are incorporated in the pattern function, it should in principle be better adapted to this particular data set.We constructed such an empirical pattern function by first averaging the four residuals h α (t) and then adding the resulting profile to the initial tanh pattern function.Using the empirical pattern function and a new set of residuals, the cross correlation (α = β) functions H αβ ( ) were found to be predominantly negative around the origin.Furthermore, the elements of the crossing time error covariance matrix δt α δt β turned out to be smaller in magnitude.Refining the pattern function may thus help to improve the accuracy of the timing method.A detailed study of the effects of empirical pattern functions on crossing time errors would be beyond the scope of the present paper and is left for future work.

Boundary parameter errors
As explained in Sect.4.1, the errors in boundary speed and normal unit vector are computed from the slowness error covariance matrix δmδm T which in turn is found from Eq. ( 21).Following the discussion in Sect.4.2, we disregard the contribution from the positional inaccuracies and consider the crossing time error covariances only.For the Cluster discontinuity crossing studied here, the numerical values of the slowness error covariances are given in Table 2.The resulting uncertainty in boundary speed is δV = 5 km s −1 , and the range of the directional error δθ is 0.6-0.9deg.

Summary
The principle variants of the timing approach to boundary analysis discussed in Sect. 3 require slightly different parameter estimation and error analysis strategies.In the present study we concentrated on absolute crossing times determined through minima of the mean square deviation I (τ ) of the data from a predefined pattern function such as a hyperbolic tangent profile.If the relative crossing time approach is employed, we recommend to construct an effective pattern function p(t) for the error analysis as follows: first apply time shifts to the signals so that the transitions all occur at the origin, and then average to obtain p(t).
The advantages of absolute crossing times over their relative counterparts are not only of technical nature.Comparing the data with a predefined pattern means that we are in explicit control of the features in the data that we wish to associate.In the relative crossing time method one compares segments of two time series around a transition (that has usually been identified by eyeballing) but the result can be distorted by substructures in the data that may move at different speeds than the boundary itself, and that may have been identified by some (pairs of) sensors but not by others.Furthermore, substructures that are moving in the plasma frame have different effects on the two main types of directional discontinuities which motivated our error analysis in the first place: TDs are stationary in the plasma frame whereas RDs propagate through the plasma.In the absolute crossing time method with a predefined pattern, such substructures become part of the residuals and are thus taken care of in the error analysis.Alternatively, they can be made explicit through an empirical pattern function as explained at the end of Sect.6.
To implement the multi-point crossing time method to boundary parameter estimation and the error analysis scheme presented in this paper, we recommend to proceed as follows.

Crossing times t α and boundary parameters m, n,V
Choose a pattern function p(t), construct the mean square deviations I α (τ ) of p(t) and the shifted signals s α (t +τ ), and identify the crossing times t α as the lag values at the minima of the I α 's.Take the spacecraft positions at t α to compute the reciprocal vectors k α .Obtain the boundary slowness vector m from Eq. ( 17), then compute V = 1/|m| and n = V m.

Boundary parameter errors δ nδ nT and (δV ) 2
Equation ( 19)-( 25) give the mean square errors of the boundary parameters n and V for general crossing time error covariances δt α δt β and spacecraft position covariance matrices δr δr T γ .To check if the latter make a significant contribution, carry out an order-of-magnitude assessment similar to the one in Sect.4.2.If the assessment is negative, positional inaccuracies can be disregarded and the error formulas simplify considerably.
In the second and third step, it is essential to construct the full crossing time error covariance matrix δt α δt β , and then to use the general formula for δ nδ nT .Correlations in the set of residuals are particularly important.If they are disregarded and an oversimplified white noise model is used to estimate δt α δt β , the crossing time errors may come out far too small.
••• denote the ensemble averaging with respect to the residual to be specified in more detail further below.
Let τ * denote the numerical value of the lag time at the minimum of the mean square deviation for the (hypothetical) "noise-free" case, and τ * the estimated lag time based on a "noisy" measurement.Note that in this context the term "noise" refers to all contributions to the signal other than the given pattern function.If noise was absent, then τ * = τ * .In the presence of noise, the mean square deviation I (τ ) is nonzero for all values of τ , and the estimated lag time τ * (i.e. the minimum of the empirical mean square deviation) differs from the true lag time τ * .
The mismatch of pattern and signal at time shift τ * defines the residual: (B2) To accomplish the error analysis, we wish to translate the mean square deviation I ( τ * ) into a function J (δt) where δt denotes the deviation of the estimated lag from its true value: We assume the residual and its derivatives to be sufficiently small compared to the derivatives of the pattern function so that only the dominant contributions need to be kept, namely, p h in the linear term, and p 2 in the quadratic term.
To arrive at an estimate for the error covariance matrix δt α δt β , we think of the residuals h α and h β as realizations of time-invariant and ergodic random processes that are well characterized by their means and correlation functions.Then ••• is the average with respect to the random functions h α and h β .Since the denominator does not depend on the residuals, it is a constant in the ensemble averaging procedure.
In principle, the time difference runs from −T w to T w .In practice, its scope is limited by the effective range of the correlation functions.

Fig. 1 .
Fig.1.Influence of array geometry on the accuracy of the timing method: reference error V δt/L.Blue and solid contour lines give the relative error in boundary speed δV /V in percent.Red and dashed contour lines give the directional uncertainty δθ in degrees.Control variables are the inter-spacecraft length scale L and the boundary speed both for crossing time inaccuracies δt = 0.1 s (annotation at the left y-axis) and δt = 1 s (annotation at the right y-axis).The spacing of the contour lines is logarithmic.

Fig. 2 .
Fig. 2.Influence of array geometry on the accuracy of the timing method: geometrical error amplification for the worst-case relative orientation of the boundary normal vector with respect to the spacecraft configuration.In this case the geometrical error amplification function is given by 1/[(1−E)(1−P )] where E is elongation and P is planarity.The contour lines in the plot are spaced logarithmically.

Fig. 3 .
Fig. 3. Influence of array geometry on the accuracy of the timing method: logarithmically spaced contours of the average geometrical error amplification function A(E,P )/ √ 3 in terms of the elongation E and the planarity P of the spacecraft tetrahedron.

Fig. 5 .
Fig. 5. Crossing time error analysis using a tanh pattern function.Top: residuals computed from the Cluster FGM data and tanh pattern functions in Fig. 4 as functions of centered array subscripts.The total time window is 10 s (sampling interval 44.6 ms).Bottom: products of correlation functions G( ) and H αβ ( ) that enter the sum in the crossing time error formula (52).Shown are the profiles for α = 1.The four colors correspond to spacecraft no.β as listed in the caption of Fig. 4.

Table 1 .
Crossing time error analysis using a tanh pattern function: elements of the crossing time error covariance matrix δt α δt β for the solar wind discontinuity observed by the FGM instruments onboard the Cluster spacecraft shortly before and after 19:11:30 UTC on 3 February 2003.The matrix elements are given in units of T 2 sam where T sam = 44.6 ms is the sampling interval.

Table 2 .
Crossing time error analysis using a tanh pattern function: elements of the slowness error covariance matrix δmδm T for the solar wind discontinuity observed by the FGM instruments onboard the Cluster spacecraft shortly before and after 19:11:30 UTC on 3 February 2003.The matrix elements are given in units of 10 −10 (s m −1 ) 2 .