Errors on derived results

1. Precise current density determination for the study of the magnetopause current layer with multiple spacecraft J. De Keyser, F. Darrouzet, E. Gamby Belgian Institute for Space Aeronomy, Brussels, Belgium johan.dekeyser@aeronomie.be • Errors on derived results • The result of LSGC, applied to a scalar field f, is f and the associated error covariance matrix C(δfi, δfj). This matrix is usually not diagonal: Because of the distribution of points and the coupled way in which the gradient components in space and time are computed, the errors on the gradient components are not statistically independent. When applying LSGC to vector fields, as in the curlometer, there is – in addition – a coupling between the errors on the components of the gradients of Bx, By, and Bz, in particular because of the div B = 0 condition. It is important to take these non-vanishing cross-correlations into account while establishing the error estimates on the current components: • µ02 ‹ δjx2› • = ‹ ( δ∂yBz - δ∂zBy )2 › • = ‹ (δ∂yBz)2 › + ‹ (δ∂zBy)2 › - 2 ‹ δ∂yBz•δ∂zBy› • µ02 ‹ δjy2› • = ‹ ( δ∂zBx - δ∂xBz )2 › • = ‹ (δ∂zBx)2 › + ‹ (δ∂xBz)2 › - 2 ‹ δ∂zBx•δ∂xBz› • µ02 ‹ δjz2› • = ‹ ( δ∂xBy - δ∂yBx )2 › • = ‹ (δ∂xBy)2 › + ‹ (δ∂yBx)2 › - 2 ‹ δ∂xBy•δ∂yBx› • where the colored terms are cross-correlations between errors on different components. Of course, there are also cross-correlations between the errors on the current components: • µ02 ‹ δjx • δjy› • = ‹ ( δ∂yBz - δ∂zBy ) ( δ∂zBx - δ∂xBz )› • = ‹ δ∂yBz• δ∂zBx› + ‹ δ∂zBy• δ∂xBz› • - ‹ δ∂yBz• δ∂xBz› - ‹ δ∂zBy• δ∂zBx› • µ02 ‹ δjx • δjz› • = ‹ ( δ∂yBz - δ∂zBy ) ( δ∂xBy - δ∂yBx )› • = ‹ δ∂yBz• δ∂xBy› + ‹ δ∂zBy• δ∂yBx› • - ‹ δ∂yBz• δ∂yBx› - ‹ δ∂zBy• δ∂xBy› • µ02 ‹ δjy • δjz› • = ‹ ( δ∂zBx - δ∂xBz ) ( δ∂xBy - δ∂yBx )› • = ‹ δ∂zBx•δ∂xBy› + ‹ δ∂xBz•δ∂yBx› • - ‹ δ∂zBx•δ∂xBy› - ‹ δ∂xBz•δ∂xBy› • For the divergence, which one may a priori impose to vanish, one has • ‹ (δ(•B))2› • = ‹ ( δ∂xBx + δ∂yBy +δ∂zBz ) 2› • = ‹ (δ∂xBx)2 › + ‹ (δ∂yBy)2› + ‹ (δ∂zBz)2› • + 2 [ ‹δ∂xBx•δ∂yBy› + ‹δ∂xBx•δ∂zBz› + ‹δ∂yBy•δ∂zBz› ] • Implementation • We are currently preparing a free and portable stand-alone implementation of LSGC, in the form of a Matlab library. • This library consists of computing routines (apply to scalar and vector fields) producing f and C(δfi, δfj) • lsgc_avds, lsgc_avd, lsgc_avs, lsgc_ads, lsgc_av, lsgc_ad, lsgc_as, lsgc • computes the gradients in a set of points, with automatic determination of velocity (v), homogeneity directions (d) and/or scales (s), when the other parameters are given • lsgc_inst_ads, lsgc_inst_ad, lsgc_inst_as, lsgc_inst • computes instantaneous gradients (uses spatial homogeneity scales for estimating total error and weighting the system) • lsgc_cgc • classical gradient computation (no weighting); for more than 4 simultaneous data, a least-squares solution is returned • and of routines for retrieval of the results • lsgc_f • returns the field + error margins • lsgc_gradx, lsgc_gradt • returns the spatial/temporal gradient + error margins • lsgc_curl, lsgc_div • returns the curl or divergence of a vector field + error margins • lsgc_vframe, lsgc_lc, lsgc_sign • returns the frame velocity, homogeneity scales, curvature signs with which the end result was obtained • lsgc_neq, lsgc_svd, … • number of equations used, singular values or the problem … • Conclusion • A robust curlometer can be obtained with least-squares gradient computation, with the divergence-free constraint built into the method. It provides error estimates on the current density vector: This is essential to assess the significance of the current densities, particularly for narrow structures as the magnetopause, or for weak currents. • References • Darrouzet, F., J. De Keyser, P. M. E. Décréau, J. F. Lemaire, and M. W. Dunlop. ‘Spatial gradients in the plasmasphere from Cluster’. Geophys. Res. Lett. 33, L08105, 2006. • De Keyser, J., F. Darrouzet, M. W. Dunlop, and P. M. E. Décréau. ‘Least-squares gradient calculation from multi-point observations of scalar and vector fields: Methodology and applications with Cluster in the plasmasphere’. Ann. Geophys. 25, 971–987, 2007. • De Keyser, J., ‘Least-squares multi-spacecraft gradient calculation with automatic error estimation’, Ann. Geophys., 26, 3295–3316, 2008. • Abstract • In this contribution we discuss the abilities of a least-squares curlometer (an application of least-squares gradient computation techniques, LSGC) to provide precise current density vectors, together with error estimates, something that is of considerable importance for the study of the magnetopause. The success of such a curlometer depends on the reliable determination of the motion of the current layer, as that helps to ensure that the optimal set of data is used for the computation of the currents. Another aspect is the use of constraints to enforce a zero magnetic field divergence. Additional geometric constraints may be imposed. The technique can be applied without any a priori limitation on the number of spacecraft or their configuration. • Classical 4-point gradient computation • Classical gradient computation (CGC) needs 4 simultaneous measurements to find the spatial gradient. • - Homogeneity: all sampling positions must be within the gradient region • Large relative errors are inevitable with numerical differentiation, so you • need accurate and well-calibrated data • Ill-conditioning of the problem if sampling positions are nearly co-planar • Least-squares gradient computation • We compute the gradient of a field f in x0 from the measurements fi at xi (red dots) by a number of spacecraft. In the neighborhood of x0, f is approximated by • where f0 and are the value and the space-time gradient at x0, with Δx = x - x0. This neighborhood is the homogeneity domain, described by orthogonal unit vectors U = […uj…] and homogeneity scales lj. The error δfi on each fi comes from the measurement error δfi,meas and the approximation error (quadratic upper bound for a linear approximation) • δfi,approx ≈ δfi,meas║Δxi║² • where the norm is based on the homogeneity properties: • ║Δxi║² = Δxi┬UL-2U┬Δxi • We now compute the gradient by solving the over-determined system in a weighted least-squares sense. The weights are the inverse of • δf²i = δf²i,meas ( 1 + ║Δxi║4[ max( 1, exp(║Δxi║²) ]² ) • the total error on each measurement. The second term estimates the error of the local Taylor approximation and depends on the homogeneity properties. Data are irrelevant for measurements far from x0 because of the growing approximation error. In practice, we include all data in a neighborhood (the light green ellipse) by setting a limit on the total error. Depending on the spacecraft configuration and the sampling rate, the number of points can be large, allowing a more precise evaluation of the gradient. • Such least-squares gradient computation (LSGC) produces error estimates on the gradient that include the effects of both the measurement and the approximation errors. In a recent development of the method, the homogeneity properties are automatically computed. As a result, the error estimates are fairly reliable indicators of the actual total error on the result. • A least-squares curlometer • The description given above for the gradient of a scalar field f can easily be extended to the case of a vector field, such as the magnetic field B : One can simply perform LSGC on each of the individual components. Ignoring time variations, the local current density vector is given by • j = curl B / m0 • thus creating a curlometer technique. • In the context of LSGC it is fairly easy to impose strict conditions (linear constraints) on the solution. An obvious constraint is to require explicitly that div B = 0. In addition, one can impose linear geometric constraints of the form • and specify that the gradient must be perpendicular to direction d (note that several such constraints can be imposed simultaneously). Imposing constraints tends to lead to a better numerical conditioning and/or avoiding singular cases. By adding a priori knowledge about the solution, one guarantees that the solution physically makes sense, so that the actual quality of the solution improves. • Determining the homogeneity parameters • The detection of magnetospheric currents is often quite difficult, since currents typically occur in sheets, like the magnetopause, in which case one has to deal with small-scale structures that are rarely sampled very well, or as distributed current systems, like the ring current, in which case the current density is rather low. It is therefore important to be able to reliably estimate the error bars on the obtained gradients and current density vectors, in order to be able to assess the physical significance of the gradient and current density vectors that one obtains. • The error bars depend on the measurement errors, but also on the spacecraft configuration relative to the scale sizes of the magnetic field structures, i.e., on the homogeneity parameters. In practice, you must be able to compute the homogeneity parameters automatically. • Determination of scales • Assume that the directions uj are given. We then try to determine the length scales lj automatically by looking at the residuals of the overdetermined system and • requiring the ratio A = ² / (N-M) → 1 • or by analyzing the pattern of residuals to improve the lengths with a heuristic so that Aj = ljnew/lj → 1 • In both cases, we use BFGS multi-dimensional optimization to minimize the target function • F = A + A-1 or F = Σj ( Aj + Aj-1 ) • This target function usually has several minima; we have devised an algorithm to find an initial guess close to the global minimum. Various strategies can help to reduce BFGS computer time consumption, which is especially welcome in the multi-dimensional case. • There is often not enough information in the sampled data set to really determine all lj. We use simple heuristics to deal with this. • Determination of homogeneity directions • While there is no general rule giving a priori knowledge about the spatial homogeneity directions, the following is a good choice, especially for the magnetopause: • direction 1 : the local magnetic field direction B: this is appropriate when the structures are roughly aligned with the magnetic field. • direction 2 : the gradient direction f itself (i.e. the second direction is a vector in the plane formed by B and f, perpendicular to B); at thin current sheets such as the magnetopause, the gradient changes rapidly in the direction of the gradient itself. • direction 3 : completes the orthonormal set • Since this set of homogeneity directions depends on the gradient itself, it has to be computed iteratively. • Determination of frame motion • The motion of the frame in which one best determines the gradient is usually the frame in which the structures do not seem to change much (equivalent to the choice of the homogeneity direction in time). That choice amounts to find the phase speed vframe (speed of the structures) and to compute the gradient in a frame moving at that speed. • in a fixed frame in a moving frame • Depending on the spacecraft configuration and frame motion, the number of data points actually inside the homogeneity domain can be somewhat larger. If the structures convect rather than change in time - as for individual magnetopause crossings - the approximation error in time can be considered small, leading to sharper error estimates. The phase speed is obtained from • dfi/dt = ∂fi/∂t + vframe.fi = 0 • Scalar case: gives the speed along the gradient direction. • Vector case: there are 3 equations, fixing the vector vframe. • As the phase speed depends on the gradient itself, it has to be computed iteratively. • Note that there are some similarities with M. Hamrin’s GALS technique, but we include a full adaptivity of the homogeneity length scales and a determination of the error on the gradient.