Empirical reconstruction and long-duration tracking of the magnetospheric boundary in single- and multi-spacecraft contexts

The magnetospheric boundary is always moving, making it difficult to establish its structure. This paper presents a novel method for tracking the motion of the boundary, based on in-situ observations of the plasma velocity and of one or more additional observables. This method allows the moving boundary to be followed for extended periods of time (up to several hours) and aptly deals with limitations on the time resolution of the data, with measurement errors, and with occasional data gaps; it can exploit data from any number of spacecraft and any type of instrument. At the same time the method is an empirical reconstruction technique that determines the one-dimensional spatial structure of the boundary. The method is illustrated with single- and multi-spacecraft applications using data from Ampte/Irm and Cluster.


Introduction
The magnetopause (MP) and the boundary layer (BL) are constantly moving.A major cause of this motion is the changing solar wind dynamic pressure: The magnetospheric boundary moves so as to re-establish the equilibrium between total solar wind and magnetospheric pressure.Total pressure changes of only a few percent cause the MP/BL to move inward or outward over 1000 km on time scales ranging from seconds to hours; extreme compressions/decompressions of the magnetosphere correspond to inward/outward displacements of several R E .Boundary motion can also be induced by a plasma instability, such as the Kelvin-Helmholtz instability (see Sckopke et al., 1981;Fitzenreiter and Ogilvie, 1995;Otto and Fairfield, 2000; Correspondence to: J. De Keyser (Johan.DeKeyser@bira-iasb.oma.be)Owen et al., 2004c;Hasegawa et al., 2004a, and references therein) in which the free energy implied in the shear flow along the MP/BL drives a surface wave.This mechanism is expected to occur at the magnetospheric flanks where the flow shear is high, when the magnetosphere is embedded in a fast solar wind stream, and for the northern interplanetary magnetic field when the stabilizing effect of the magnetic field line curvature pressure is minimal.Typical surface wave frequencies are 0.001-0.05Hz.
Knowing the time-varying position of the boundary is essential for converting a time series of in-situ observations into spatial information, an operation known as "empirical reconstruction".A well-known example is the classical single-spacecraft determination of magnetopause thickness (Berchem and Russell, 1982): crossing duration times the average plasma speed in the boundary normal direction gives magnetopause thickness.This example highlights the ingredients needed for time-to-space conversion: 1) It is assumed that the structure is not dynamically changing, but only convected.2) An appropriate reference frame has to be established, where boundary motion is measured in the boundary normal direction.3) The boundary is taken to be a planar tangential discontinuity, through which there is no plasma flow, so that the plasma velocity approximates the boundary velocity.4) The time scale can then be translated into a spatial scale by using boundary velocity information.The present paper deals with a more sophisticated empirical reconstruction method.It should be noted that magnetic field-based reconstruction methods also exist (e.g., Hau and Sonnerup, 1999;Hu and Sonnerup, 2003;Hasegawa et al., 2004b); we will briefly compare both approaches in the Conclusions section.Paschmann et al. (1990), and later Phan and Paschmann (1996), proposed to integrate the plasma velocity over time to obtain the time-varying position of the boundary.Straightforward application of this idea has limited success because the errors grow rapidly as the analysis interval becomes longer.De Keyser et al. (2004) have applied this idea in the multi-spacecraft case: Combining plasma velocity data from all spacecraft can reduce the errors.It has also been shown that problems due to inaccurate boundary velocity information can be remedied by using an optimization approach, in which a model boundary motion profile is matched to the observed boundary motion proxy while simultaneously fitting other parameters to prescribed spatial profiles (De Keyser et al., 2002).The present paper introduces a new empirical reconstruction technique that is based on the optimization approach.Section 2 explains the rationale behind the method.Section 3 describes an efficient algorithm for solving the optimization problem.In Section 4 the technique is applied to single-and multi-spacecraft observations of the MP/BL by Ampte/Irm and Cluster.Section 5 summarizes the main features of the method and suggests a number of physical problems to which the method could be applied.

Observations
Consider a time interval [t start , t stop ].
Let the position x k sc of spacecraft k=1, . .., K be known.Some of these spacecraft (at least one of them) measure the local plasma velocity v k (t k,v i ) with uncertainty δv k (t k,v i ), at times t k,v i ∈[t start , t stop ] for i=1, . .., N k,v .The times of measurement on the different spacecraft are, in general, not synchronized.Some or all of these spacecraft provide magnetic field observations B k (t k,B i ), i=1, . . ., N k,B .The magnetic field instruments usually sample the medium at different times than the plasma instruments.
In-situ observations are assumed to be available for physical quantities l=1, . .., L, which we will call "guiding variables".Any observable can be used as a guiding variable, although it is most useful to consider variables that distinctly change across the MP/BL, since the values of such variables provide clear information about the position of the spacecraft relative to the boundary.All guiding variable observations f k,l (t k,l i ) with error estimates δf k,l (t k,l i ), i=1, . . ., N k,l do not have to be available on all spacecraft.There are no restrictions on the sampling times, as different instruments typically operate at their specific sampling rate.Observations from different spacecraft have to be intercalibrated first; an absolute calibration is not required, but remains desirable for a proper interpretation of the physical quantity.

Reference frame
The standard boundary normal coordinate system for analyzing individual MP/BL crossings (Sonnerup and Cahill, 1967) is valid for short time intervals only, as it postulates a fixed boundary orientation.For a long-duration analysis that includes multiple and/or partial boundary crossings, we resort to a generalization (De Keyser et al., 2002): If the boundary structure is magneticly field-aligned (a tangential discontinuity), the local normal can be obtained as . Applying variance analysis to the set of normals n, possibly adding the constraint n z =0 (assuming the structure to be at most twodimensional, so that the surface normal remains confined to the xy plane), and ordering the principal axes of the covariance matrix from highest to lowest variance, results in an intermediate frame x y z.We consider here the case where the minimum variance is well separated from the two other covariance matrix eigenvalues, indicating that n remains confined to the plane x y and that the MP/BL is at most twodimensional.Frame x y z is then rotated around z over an angle θ to obtain a frame xyz, where x points along the "average" outward normal.This direction (angle θ ) is identified by trial and error.It is the direction along which one can best measure the inward-outward motion of the boundary; for instance, at the magnetospheric flanks, a useful first estimate of the "average" outward normal direction can be obtained by requiring the tailward magnetosheath flow to be perpendicular to it, so that only the relatively modest and variable flow component associated with the inward-outward motion of the boundary remains in the "average" outward normal direction.The method described below is not very sensitive to the precise orientation, unlike earlier methods.The xyz frame is the desired generalization.

Boundary position and motion
The position x mpbl (t) of the MP/BL can be obtained by integrating its motion v mpbl (t) along x in time: where x mpbl (t c ) defines a reference position at time t c , the center of the time interval.It is implicitly assumed that x mpbl (t) is a single-valued function, which might not be true in the nonlinear stage of the evolution of surface waves (Otto and Fairfield, 2000) or in the presence of isolated magnetosheath plasma entities inside the magnetosphere (Lemaire and Roth, 1991).If the plasma near the MP/BL moves collectively back and forth as a planar structure, without plasma flow across the MP/BL, the locally measured plasma velocity v x (or v ⊥x ) can be used as a substitute for the boundary velocity (Paschmann et al., 1990;De Keyser et al., 2004).Alternatively, one may use the electric drift or the velocity of a single plasma component.Integrating the oscillatory and sign-varying v mpbl , however, leads to an ever-increasing relative error on the result as the integration interval becomes larger.This is aggravated by the limited precision and time resolution of plasma velocity measurements; a time resolution of a few seconds at most appears to be required in practice.Significant data gaps pose an insurmountable problem.
In order to address these difficulties, the method described below regards v x as a proxy for v mpbl rather than a substitute.It allows both to differ somewhat, depending on the error estimate δv x .The guiding variables provide the additional information needed to determine v mpbl .

Optimization problem
Tracking MP/BL motion can be formulated as an optimization problem: Determine boundary motion v mpbl (t) and find a spatial profile f l (x) for each guiding variable so that v mpbl matches the proxy v k x and so that each f l fits the corresponding f k,l data.Below, we cast this optimization problem into a mathematical expression.
Each guiding variable profile f l (x) is defined over a prespecified interval [x min , x max ].The same interval is used for all guiding variables.A profile is represented by a linear spline function with equidistant nodes x i , i=1, . . ., µ (spatial resolution x=(x max −x min )/(µ−1)): with φ(s)=1− |s| for |s| <1 and zero elsewhere.The values f l lef t and f l right at x min and x max , respectively, are prespecified.Strictly positive guiding variables (e.g.densities and temperatures) should be replaced by their logarithms to avoid the need for positivity constraints on their spline representation.The time profile v mpbl (t) is represented by an equidistant linear spline over [t start , t stop ] with nodes t i , i=1, . .., γ (time resolution t=(t stop −t start )/(γ − 1)): The boundary position x mpbl (t) is computed from v mpbl (t) by Eq. (1).The optimization domain consists of the Lµ spline coefficients f l i and the γ coefficients v mpbli .While a fairly low number of spatial points µ is often sufficient (tens of points), the number of time points γ can be large, depending on the length of the time interval (typically an hour or longer) and the desired time resolution (from minutes down to a few seconds).
The deviation between guiding variable data and their spatial fit can be expressed in a least-squares sense by The fit f l is evaluated at x k (t k,l i )=x k sc (t k,l i )−x mpbl (t k,l i ), the position of spacecraft k relative to the boundary.In principle, the guiding variable observations must be taken in the moving frame.For quantities that are not invariant under a (Galilean) transformation, such as velocities or electric fields, the above formulation should be extended to include the frame transformation (which depends on the unknown v mpbl ).In practice, however, the errors introduced by ignoring oscillatory frame motion can often be regarded as a kind of measurement noise.The match between the boundary velocity and the proxy is quantified, again in a least-squares sense, by The smoothness of the spatial profiles f l (x) and of the time profile v mpbl (t) is measured by the mean square of their second order derivatives The overall target function that must be minimized can then be expressed as where the τ l f , τ l n , τ v , and τ s are normalization constants (defined in Appendix A) whose purpose is to make sure that all contributing terms have a fixed order of magnitude that is essentially independent of the number of spacecraft, the number of guiding variables, and the number of data samples.
The normalized weights w l in Eq. ( 2) can be used to bias the optimization toward certain guiding variables ( L l=1 w l =1).The guiding variable profile smoothness parameter 0<λ<1 and the boundary velocity profile smoothness parameter 0<σ <1 are dimensionless constants that determine the importance of the second derivatives of the guiding variable profiles and of the boundary velocity profile, respectively, in the overall target function.The proxy confidence parameter, the constant 0<α<1, fixes the trade-off between matching the boundary velocity proxy and fitting the guiding variable observations.If the proxy is thought to be a highly reliable indicator of boundary motion, its value should be taken close to 1, so that the guiding variable terms F f and F n become negligible.If the proxy is only a crude indicator, a smaller α will ensure that all terms in the target function contribute to defining the optimal solution.
The G l n and G s terms are regularization terms that guarantee the uniqueness of the solution.Occasionally, there may be no observations in a subinterval [x i−1 , x i+1 ], so that the values f l i , l=1, . .., L drop out of the expressions for G l f ; hence, their values would be undetermined if it were not for the G l n terms.A similar problem arises for the v mpbli when there are no observations in [t i−1 , t i+1 ].This happens whenever t is smaller than the time resolution of the data, for instance, in data gaps.Here, the G s term solves the problem.

Solving the nonlinear problem
Minimizing target function Eq. ( 2) is difficult for three reasons.First, it is nonlinear as the boundary velocity coefficients, through the boundary position, appear in the argument of the nonlinear profiles f l (x).Second, the dimension of the search space is large, especially if the time interval is long and if the time resolution is high.Third, evaluating the target function is expensive, particularly if there are several spacecraft and guiding variables.Nevertheless, an effective optimization strategy can be constructed by exploiting the particular properties of the optimization problem.
The problem has at least one minimum (there is a lower bound F ≥0) and it obviously cannot have local maxima.Finding all the minima is equivalent to solving ∂F /∂v mpbli = 0, i=1, . . ., γ , ∂F /∂f l j = 0, j =1, . . ., µ; l=1, . . ., L. (3) This nonlinear system may have multiple solutions, from which the global minimum has to be selected.The following properties help to find that global minimum:

Property 1 Fitting guiding variables only
If the boundary motion is given (α=0), the problem consists of L independent minimizations of that is, the f l (x) are least-squares smoothing spline fits of the guiding variable observations, corresponding to the second set of equations in Eq. ( 3).This set consists of L linear systems in the f l j , so that each fit is the unique global minimum of the F l .These linear systems and their solution are discussed in Appendix B.

Property 2 Fitting boundary velocity only
For α=1 the target function becomes that is, v mpbl (t) is a spline fit of the proxy observations, corresponding to the first set of equations in Eq. (3).This system is linear in the v mpbli coefficients, so that there exists a unique solution, which necessarily is the global minimum.If α<1 but not too small, when the proxy is believed to be reliable, a spline fit or a simple interpolation of the proxy should provide a good initial solution for the v mpbli .

Property 3 Smoothness of guiding variable profiles
Strong smoothing of the guiding variable profiles (λ→1) implies monotonous changes from f l lef t to f l right across the spatial domain.The profiles can then locally be approximated by linear functions so that each x mpbli depends linearly on (some of the) v mpblj through integration (1).System ( 3) is then linear in the vicinity of the solution: A unique extremum exists there, which must be the global minimum.Enforcing a strong degree of smoothness, therefore, prevents the optimization from being trapped in local minima.Strong smoothing, however, obscures the spatial structure.

Exploring the solution space
Minimization algorithms search for solutions with progressively lower target function values.While Property 3.1 ensures that there are no local minima in the vicinity of the solution when λ→1, such local minima may exist for the lower values of λ that are used in practice.The optimization process can avoid such local minima by starting with a good initial guess and by exploring alternatives not only close to the current best solution, but more distant ones, too.

Producing good starting solutions
Good starting solutions can be generated by considering p=1, . . ., P subproblems defined on nested subintervals symmetric around the center time t c .Each subproblem is simpler to solve because the number of unknowns is lower (γ is smaller) and because target function evaluation is cheaper (less data during a shorter time interval).Solving the p=1 problem begins by choosing v 1 mpbl as an interpolant or a fit of the proxy (Property 3.1) and computing the corresponding spatial profiles (Property 3.1).There is not much error accumulation in integration (1) for a short time interval, so that this initial guess is quite good.Problem p=1 is also small: Finding the global minimum does not require much work.It is always possible to add interpolated proxy data in mpbl , so that it can serve as an initial guess for the p-th problem.Both extensions are rather short so that integration errors remain acceptable; the initial guess is a good one.By repeated application of this idea, the time interval can be extended without being hindered by error accumulation for as long as the physical assumptions warrant.We use subintervals whose length increases exponentially with p, so that a constant fraction of new information is added when going to the next problem.
Alternative approaches include first solving the problem with a large t (which is cheap) and progressively refining it, or beginning with λ=1 (for which Property 3.1 guarantees a unique solution) and gradually lowering it.

Exploring variant solutions
An optimization procedure should explore variant solutions that are not only in the immediate vicinity of the current best solution, in order to avoid getting stuck in a local optimum.Since the current best solution is most likely to be close to the global optimum (and this probability increases as the optimization proceeds), however, most of the effort should be invested in examining the local neighborhood.A good balance between local and nonlocal searching is offered by the multi-level technique outlined below.
The solution space is explored by repeatedly generating a variant ṽmpbl (t)=v mpbl (t)+ v mpbl (t) of the current best solution.Then ṽmpbl (t) is integrated to obtain xmpbl (t).The modified spatial profiles are found as outlined by Property 3.1 and they are interpolated at the spacecraft positions relative to xmpbl at the times of guiding variable observations.The spatial profiles and ṽmpbl (t) are differentiated twice.Target function ( 2) is then evaluated and compared to the original value to either accept or reject the variant.Evaluating the quality of a variant is clearly computationally expensive.
We examine variants for which v mpbl <ζ .The optimization procedure starts with ζ t=x max −x min , so that the search is indeed nonlocal.Subsequently, ζ is halved when no more improvements can be made with the current value.This is repeated until ζ reaches a prescribed precision ζ * in the v mpbli search space.We consider changes v mpbl that are nonzero in a subinterval [t i−1 , t i+m+1 ].If also then x mpbl (t) =0 in that subinterval alone: Only a few terms in the target function are affected, so that evaluating the effect of the change is much cheaper.We consider changes on a binary hierarchy of time scales m t that range from the order of the tracking duration, for m=2 b (b= log 2 γ is the number of binary levels) down to the time resolution, for m=1.This is yet another aspect of nonlocal exploration.
The simplest change is the one depicted in Fig. 1, where v mpbl has two nonzero values, with opposite sign, at t i and t i+m .The net effect of such a change is to displace the points in [t i−1 , t i+m+1 ]; the displacement is ξ =ζ t, except at the edges of that subinterval.One should, however, respect the integration constant in Eq. ( 1): x mpbl (t c ) = 0.A constant must therefore be added to x mpbl when t c ∈[t i−1 , t i+m+1 ].Each change v mpbl affects the F v terms only in two subintervals of length 2 t around t i and t i+m (which merge when m≤2).The effect on F s is cheap to evaluate as d 2 v mpbl /dt 2 is nonzero only at t i−1 , t i , t i+1 , t i+m−1 , t i+m , and t i+m+1 , the changes being 1, −2, 1, −1, 2, and −1 times ζ /( t) 2 , respectively.The changes are applied by scanning over overlapping intervals: Starting from [t c−m , t c ] or [t c , t c+m ], the interval is repeatedly shifted left or right, respectively, over m/2 t (or over t if m=1).
There are two options for exploring the f l j dimensions.The standard procedure is to recompute the spatial profiles for each ṽmpbl by solving the L least-squares problems (Property 3.1).This is the most rigorous way to proceed, but also the most computationally intensive.When it is known a priori that the profiles are not affected much by the change, e.g. when the boundary displacement is smaller than the spatial resolution, ξ < x, it is sufficient to update the profiles once, after each full scan of possible changes.This is a much cheaper alternative.Each v mpbl then implies a change in F l f (with unmodified f l j ) that includes only terms in [t i−1 , t i+m+1 ] as xmpbl =0 there; as the profile is not updated, F n does not change at all.The update of the profiles at the end of a scan leads to a further decrease of F as F l f and F l n are minimized.Evaluating a few terms rather than the whole target function speeds up the computation.Note also that the target function change F becomes much smaller than F as the optimization proceeds, so that a significant number of digits is lost when comparing solutions.By comparing just a few terms, however, this loss of precision is avoided.θ =26 • (see Sect. 2.2), giving x=[+0.594−0.803+0.050]y=[+0.640+0.509+0.575],z=[−0.487−0.309+0.817] in GSE coordinates, so that x is the average outward normal, y points roughly sunward, and z is the invariant direction and pointing northward.
Since the three-dimensional plasma electrostatic analyzer (Paschmann et al., 1985) and the fluxgate magnetometer (Lühr et al., 1985) observations are given with 4.3-s resolution, a reconstruction with t=5 s is computed.The guiding variables are the ion density n i and the magnetic field component B z , with equal weights w l =0.5.A relative error of 10% on n i and an absolute error of 0.1 nT on B z were specified.The proxy for boundary motion used here is v ⊥x .Considering this proxy to be reliable, it is given a high relative importance by setting the proxy confidence parameter α=0.8.The spatial profiles are defined by µ=60 points and x=500 km, with a little spatial smoothing λ=0.0001.A minor amount of boundary velocity profile smoothing σ =0.0001 is also imposed on v mpbl (t) to prevent irregular variations where data points are missing.The optimization problem is solved in 10 steps by starting with a small time interval at the center and extending it by 55% in each step to arrive at γ =541 points for the final problem, following the procedure outlined in Section 3.2.1.The subproblems are solved to a relative precision on the v mpbli of 10 −3 ; the final solution is computed with 10 −5 precision.Once v mpbl is computed, spatial 1-D model profiles can be fitted to the observations for any physical quantity, representing MP/BL spatial structure.Given these model profiles and v mpbl (t), the time variations at each spacecraft can be "predicted".
The time profiles (Fig. 2, top panels) show the observations as discrete data points and the model predictions as continuous curves.The model matches the observations very well, not only for the guiding variables n i and B z , but also for the ion temperature T i , the full magnetic field vector B, and the plasma velocity v. Single-valued spatial profiles for these quantities do indeed exist (Fig. 2, bottom left).Because of the rotation of the magnetic field through the magnetopause, there is a distinct change of B y coincident with each B z reversal.To the extent that boundary motion is due to surface waves with rather short wavelengths, the associated curvature produces B x and B y variations that are not accounted for by the model.For the same reason, the fits for v x and v y are worse than for v z .The fifth panel of Fig. 2 compares v mpbl to the v ⊥x data.In spite of observational errors, boundary curvature effects, and compressibility of the plasma, v mpbl (t) closely follows the proxy, except when the boundary moves back-and-forth rapidly.Boundary motion (seen at 5-s resolution) can be extremely fast, occasionally in excess of 250 km/s, and variable, implying accelerations of up to 10 km/s 2 .Boundary position and spacecraft trajectory are shown in the sixth panel; the boundary moves over ∼2 R E during the time interval considered.
The fairly small scatter of the data around the spatial profiles (Fig. 2, bottom left) confirms that the structure of the boundary is essentially one-dimensional.Although only n i and B z are used as guiding variables, the other physical pa-rameters follow similarly well-defined profiles (less so for B x , B y , v x and v y , as explained above).The high magnetic shear magnetopause is identified by the reversal in B z .Magnetic field strength appears to be enhanced just earthward of the magnetopause.The boundary layer is evident in the n i , T i , v y , and v z profiles.These spatial profiles represent the average structure over the whole time interval; precise instantaneous magnetopause thicknesses should be obtained by analyzing individual crossings with specific techniques (e.g.Dunlop et al., 2002;Haaland et al., 2004a,b).
The convergence plots (Fig. 2, bottom right) the progress of the optimization procedure.The optimization is restarted for each of the 10 subproblems.Requiring a relative precision (v mpbl )≤10 −5 , corresponding to ∼0.005 km/s, might seem absurdly small, but because integrating such a small v mpbl over ∼1 h implies a ∼20 km change on x mpbl , a succession of a number of such changes could modify x mpbl by hundreds of kilometers (in a worst-case scenario).The relative deviation (F −F * )=(F −F * )/F * of the target function value F for the current solution from the (unknown) optimum F * was estimated by comparing the solutions to the last one found.The bottom panel shows how the target function and its contributing terms evolve.The final solution corresponds to a ratio F f :F v :F n :F s of about 8:4:3:1, indicating that fitting the guiding variables and the proxy dominate the optimization process; the v mpbl smoothness constraint is of minor importance.
The result obtained here is definitely better than the reconstruction by De Keyser et al. (2002).Considering the plasma density, for instance, shows that the relative heights of the high density plateaus given by the older method (De Keyser et al., 2002, , Fig. 1) are not fully consistent with what is observed (e.g. the model density around 08:27 is too low, while that around 08:32 UT is too high), while there is complete consistency with the present method (Fig. 2, top).This is not surprising: The older method prescribed the shape of the spatial profiles; the profiles also had to satisfy 1-D MHD equilibrium conditions.The present method relaxes these constraints, thus enlarging the solution space so that a better result can be found.As a bonus, the new method determines the shape of the spatial profiles rather than prescribing it.
We have also computed a reconstruction (not shown) for the same event with a time resolution t=30 s; all other parameters were the same.It is found that the time profiles match the data in a more crude but still consistent way, as features on a time scale < t can no longer be represented.At this time resolution, the v mpbl profile is barely able to follow the v ⊥x proxy.There is also a larger scatter of the data points around the spatial model curves.While it can be more difficult to find a good initial guess, the optimization process requires less computations per iteration (γ is smaller), so that the solution is found more quickly for larger t.The optimization can be further accelerated by first averaging the guiding variable observations over a time scale on the order of t.Empirical reconstruction for Ampte/Irm MP/BL pass on 6 December 1984, 08:00-08:45 UT (x is the average outward normal, y is roughly sunward, and z is the invariant direction) computed with a time resolution t=5 s, using v ⊥x as the boundary motion proxy and ion density n i and magnetic field B z (with equal weights w l =0.5) as guiding variables, and with proxy confidence parameter α=0.8, with µ=60 points in the spatial profiles, with spatial resolution x=500 km, with guiding variable profile smoothness parameter λ=0.0001, and with boundary velocity profile smoothness parameter σ =0.0001.Top panels: Time profiles with observations (markers) and predictions from the reconstruction (curves).Bottom left: Spatial profiles.Bottom right: Convergence history; the problem was solved in 10 steps, with a subproblem precision of 10 −3 and a final precision of 10 −5 .

Extending the time interval
The as guiding variables with equal weights w l =0.5.The proxy confidence parameter is α=0.8, µ=80 points with a spatial resolution x=500 km are used to represent the guiding variable spatial profiles.To avoid local minima, more spatial smoothing is needed than before (guiding variable profile smoothness parameter λ=0.001).The boundary velocity profile smoothness parameter is σ =0.0001.Because of the longer time interval, 15 subproblems are used (a 45% extension of the time interval per step), each solved up to 3 digits in precision, while the final solution is computed up to 5 significant digits for v mpbl .The overall agreement in Fig. 3 between model and observations over the entire time interval is remarkable.The scatter of data points around the spatial profiles cannot be less than in Fig. 2, as this reconstruction is not an improvement but an extension of that earlier calculation.Note the abnormally low plasma densities observed around 06:45 UT that contribute to that scatter; these are probably invalid data points; the method is robust enough to tolerate the presence of a few such points.When there is no valuable information in the proxy and when the guiding variables are more or less constant (e.g. in the magnetosheath, away from the MP/BL), the optimization does not constrain v mpbl : The boundary motion found for the last 10-15 min of the time interval is not very reliable.
A huge amount of data is involved in this calculation: The search space is large (γ =2161) and there are ∼2500 data points per guiding variable.The optimization requires a lot of computational effort.The result, however, is rewarding: a common underlying structure for three hours of multi-instrument observations.The MP/BL appeared to be undulating with a small amplitude-over-wavelength ratio, and magnetosheath conditions seemed to be rather constant.1984, 17:30-19:15 UT, computed with a time resolution t=5 s, using v ⊥x as the boundary motion proxy and ion density n i (weight 0.4), ion temperature T i (weight 0.4), and magnetic field B z (weight 0.2) as guiding variables, and with proxy confidence parameter α=0.8, µ=80 points in the spatial profiles, a spatial resolution x=500 km, guiding variable profile smoothness parameter λ=0.001, and boundary velocity profile smoothness parameter σ =0.0001.Top panels: Time profiles with observations (markers) and the predictions from the reconstruction (curves); data gaps are highlighted.Bottom left: Spatial profiles.Bottom right: Convergence history; the problem was solved in 5 steps, with a subproblem precision of 10 −3 and a final precision of 10 −5 .

Dealing with data gaps
In order to demonstrate the robustness of the method, it is applied here to the inbound Ampte/Irm MP/BL pass on 29 December 1984, 17:30-19:15 UT, shown in Fig. 4: Several data gaps occur during this time interval, as marked out in the time profiles in the figure; the magnetic field data gaps are slightly smaller than the plasma data gaps.This pass is also interesting because it encompasses the historically first empirical reconstruction by Paschmann et al. (1990, Fig. 4), for 18:04:33-18:18:15 UT, which covers only part of the boundary layer.
The procedure of Section 2.2, with θ =21 • , leads to a reference frame x=[+0.729−0.643 −0.233], y=[+0.629+0.496 +0.599], z=[−0.269−0.583 +0.766] in GSE coordinates.The reconstruction is computed with high time resolution ( t=5 s), as there is a strong time variability in the data, especially in the interval 08:20-08:40 UT; this variability might be related to unsteady magnetosheath conditions (some interplanetary magnetic field variability was observed prior to the time interval considered here).We use v ⊥x as a proxy for the boundary motion.The guiding variables are n i (weight 0.4), T i (weight 0.4), and B z (weight 0.2).The magnetic shear is about 92 • for this event, so that the magnetic field changes are less pronounced than in the previous examples.It is nevertheless worthwhile to include B z as a guiding variable because it partly fills in the plasma data gaps.The proxy confidence parameter is α=0.8.The spatial profiles are defined by µ=80 points, a spatial resolution x=500 km, and a guiding variable profile smoothness parameter λ=0.001.The boundary velocity profile smoothness parameter is σ =0.0001.The problem was solved in 5 steps, each computed up to 3 digits precise, while the final solution has 5 significant digits.
The time profiles (Fig. 4, top panels) show that there is a reasonably good fit.Because of the smoothness imposed on v mpbl (t), a smoothly varying boundary motion is obtained in the data gaps.The few B z data points in the middle of the largest plasma data gap do constrain the boundary position there, in spite of the absence of v ⊥x , n i , and T i data.The model predictions are less reliable in regions where there are observations for only some guiding variables (e.g.no plasma data but only B z data) and meaningless where there are no observations at all (common plasma and field data gaps).In spite of the data gaps, the method is able to connect both sides of a data gap because of the assumption of a common underlying spatial structure (Fig. 4, bottom left): While the exact behavior of v mpbl in the data gaps strongly depends on the somewhat arbitrary choice of the boundary velocity profile smoothness parameter σ , its time-integral (the change in position of the boundary across the gap) is well-determined.Note, again how, the data-driven F f and F v terms dominate the optimization result (Fig. 4, bottom right).

Multi-spacecraft application
The Cluster four-spacecraft mission conducts detailed observations of the MP/BL.The first, the third, and the fifth plot in Fig. 5 show the electron density obtained by the PEACE instruments (see, e.g., Owen et al., 2001) on the four spacecraft, the magnetic field strength from the four FGM magnetometers (Balogh et al., 2001), and the plasma velocity magnitude observed by the CIS/HIA instrument (Rème et al., 2001), only operating on C1 and C3, during an outbound pass on 6 June 2001, 01:10-03:30 UT.While C1, C2, and C4 are close together and see similar profiles, C3 is about 3000 km earthward, systematically recording lower densities and higher field strength.
The procedure outlined in Section 2.2, now using magnetic field data from the four spacecraft, with θ =41 • , produces a reference frame x=[+0.616−0.503 −0.606], y=[+0.640−0.130 −0.758], z=[+0.460+0.854 +0.242] in GSE coordinates.An empirical reconstruction is computed using v ⊥x obtained from CIS/HIA on C1 and C3 as a proxy for boundary motion.The electron density n e and the magnetic field B z are used as guiding variables (equal weights w l =0.5): Density or magnetic field variations are observed by at least one of the spacecraft throughout most of the interval.The time resolution is t=20 s.The proxy confidence parameter is α=0.5.The guiding variable pro-files are defined by µ=60 points, a resolution x=500 km, and a smoothness parameter λ=0.001.The boundary velocity time profiles are defined by γ =421 points and a smoothness parameter σ =0.0001.The problem is solved in 10 steps (problem size increases 51% per step) with an intermediate precision of 3 digits and a final precision of 5 digits.Note that there are about 3200 guiding variable observations to be accounted for in the target function, the large number being due to the 4ss time resolution of the data and the fact that there are four spacecraft.The "predicted" time profiles are shown in the second, the fourth, and the sixth panel in Fig. 5.There is a nice overall agreement, but some differences are clearly visible.Interestingly, the method "predicts" the plasma velocity profiles at C2 and C4, for which no actual observations are available.There are considerable differences between the v ⊥x measurements from C1 and C3, where v mpbl , being a kind of 20-s average, follows their overall behaviour.The figure also shows the position of the four spacecraft relative to the boundary position x mpbl .
The spatial profiles of n e , B, and v (Fig. 5, bottom) show a significant scatter of the data points rather than being truly single-valued relations.Note, however, that this scatter remains within each spacecraft data set, while the curves from all the spacecraft are consistent with each other.It can therefore be concluded that the observed time variations and the differences between the time profiles from different spacecraft can be explained to a large extent by the back-and-forth motion of the boundary and by the relative positions of the four spacecraft, but there are differences that may be due to boundary curvature or to time dependent effects (for instance, the higher-than-average density observed at times of rapid boundary motion).

Conclusions
This paper describes a robust method for tracking the position of the magnetospheric boundary for an extended time period, given observations of a proxy for boundary motion and measurements of one or more observables that can be used as guiding variables.The method is based on an optimization approach that simultaneously reproduces the position of the boundary as a function of time and the onedimensional spatial profiles of the guiding variables.This boundary tracking method therefore inherently is also a onedimensional empirical reconstruction method.The method determines the boundary motion by correcting the proxy observations so that the corresponding guiding variable data show minimal scatter around their one-dimensional spatial profiles.Whatever the precision of the proxy and the guiding variable observations, one should aim at a reasonable balance between the goals of fitting guiding variable observations and proxy observations.
The optimization problem defined here is, in its original form, hardly tractable.By exploiting the properties of the problem, however, an effective optimization algorithm has been developed that includes a procedure to find good :10-03:30 UT, computed with a time resolution t=20 s, using v ⊥x from C1 and C3 as the boundary motion proxy, with electron density n e and magnetic field B z (equal weights w l =0.5)) as guiding variables, and with proxy confidence parameter α=0.5, with µ=60 points in the spatial profiles, with a spatial resolution x=500 km, with guiding variable profile smoothness parameter λ=0.001 and boundary velocity profile smoothness parameter σ =0.0001; the optimization used 5 steps, with a subproblem precision of 10 −3 and a final precision of 10 −5 .Top panels: Time profiles of electron density observations (PEACE, C1-C4); n e predictions from the reconstruction (C1-C4); magnetic field strength observations (FGM, C1-C4); |B| predictions (C1-C4); plasma velocity magnitude observations (CIS/HIA, C1 and C3); |v| predictions (C1-C4); boundary motion proxy v ⊥x (C1 and C3) and v mpbl ; and boundary position x mpbl and spacecraft trajectories (C1-C4).Bottom panels: Spatial profiles of electron density (C1-C4), magnetic field (C1-C4), and plasma velocity (C1 and C3).
starting' solutions, a particular choice for the boundary velocity alternatives that are explored, an efficient strategy for updating the guiding variable profiles, and a cheap way to evaluate target function changes.
There are two extreme situations to which this method can be applied.If the time interval is short and the proxy is accurate enough to obtain nearly single-valued spatial profiles by integrating the proxy directly, using Eq. ( 1), the method modestly improves the quality of the spatial profiles and of the boundary velocity as a function of time, whatever the spatial structure implied by the profiles (not necessarily monotonous).If, however, the time period is extended or if the proxy is rather inaccurate, a solution can always be found with smooth guiding variable profiles (Property 3.1) but any substructure in the boundary layer will be blurred as the spatial profiles are forced to be essentially monotonous.Practical applications are situated in between those two extremes.
We have shown with Ampte/Irm examples that the method compares favourably to older methods, that the time interval can be expanded to include a full MP/BL pass, that the method is robust and that it can tolerate data gaps.The Cluster four-spacecraft application illustrates how the computed spatial profiles and the relative position of the spacecraft in the average boundary normal direction can account for the differences in the time profiles measured by the individual spacecraft.Measuring the structure of the boundary simultaneously at different points increases the volume of data that constrains the boundary motion and structure, and adds confidence to a proper interpretation of the model by allowing one to verify the hypotheses of time-stationarity and onedimensionality.
The proposed method has a number of limitations that are related to the assumptions: -The tangential discontinuity assumption: a) This assumption is made when we try to establish an appropriate boundary-aligned reference frame.This is not essential for the method; one could have used other criteria to find such a frame.A small nonzero B n would not affect the orientation of the frame significantly.Moreover, the method is not very sensitive to the precise frame orientation.b) This assumption was also used to justify taking the locally measured v n as a proxy for boundary motion.Clearly, a small nonzero flow through the boundary would not dramatically alter the quality of the proxy.And the method can easily deal with a significant uncertainty margin on the proxy, for instance, by appropriately choosing the weight factors in the target function.
-The assumption of planarity: a) This assumption is made to justify taking v n as a boundary motion proxy but, as said above, the method is fairly tolerant for inaccurate proxies.b) This assumption is used to justify that the structure can be represented by a spatial profile which depends on one variable only.c) It is also used to argue that different spacecraft, separated in distance along the boundary, would observe the same structure at the same distance from the nominal position of the boundary.As we have illustrated in the examples, no strict planarity is required, but it is sufficient that the amplitude-to-wavelength ratio is small.
-The structure is not dynamically changing: If there would be significant dynamical changes, it makes no sense to search for a common long-term structure.The time interval can then be split, and the analysis can be done separately for subintervals during which no dynamic changes take place.Although it is a serious limitation to assume that there are no major dynamic changes, it is much less restrictive than requiring that the magnetospheric boundary is not accelerating or that its acceleration would be constant, as in most, if not all, other magnetopause analysis methods.
While these assumptions have to be made to justify the method, the examples indicate that the method actually is rather tolerant to modest violations of these assumptions.
In any case, the method can be used as a first step in the analysis to give an impression of the overall structure and of the changing boundary position.If one finds that the method is not able to properly represent the structure with a one-dimensional spatial profile, one can then turn to other methods that introduce additional dimensions (time, a second and/or a third spatial dimension).
We have already mentioned the existence of magnetic field-based reconstruction methods (e.g., Hau and Sonnerup, 1999;Hu and Sonnerup, 2003;Hasegawa et al., 2004b).The most evolved of these methods use the measured magnetic field and plasma pressure as boundary conditions to solve the two-dimensional magnetohydrostatic Grad-Shafranov equation, thus providing a two-dimensional picture of the magnetic field configuration during a magnetopause crossing.This requires the existence of a de Hoffman-Teller frame; roughly speaking, the structure needs to be non-accelerating, and that limits the method to rather short time intervals, up to a few minutes at most.As a consequence, this method can exploit multi-spacecraft data only for relatively short spacecraft separation distances.The empirical reconstruction method discussed in the paper has a different goal: It primarily aims at determining the motion of the boundary, while the one-dimensional reconstruction is a by-product.Rather than excluding time periods of boundary acceleration, the intent is to study long time intervals that include acceleration, and therefore multiple partial and/or full crossings of the magnetospheric boundary.The method differs from Grad-Shafranov reconstruction on the following points: -It does not rely on the total plasma pressure, which may be difficult to determine when not all plasma components are measured.
While Grad-Shafranov reconstruction is aimed at a detailed study of individual magnetopause crossings, empirical reconstruction studies the long-term motion of the boundary and how that results in multiple boundary crossings.Both types of methods complement each other.In summary, the method proposed here seems to be, as far as we know, the only method able to track the boundary position for a very long time interval.Applying the method with care, and remaining aware of its limitations, it can give valuable information about the moving magnetospheric boundary.It therefore seems to be an interesting addition to our arsenal of single-and multi-spacecraft analysis techniques for in-situ observations of the magnetospheric boundary.It could prove useful in the following physical applications: -Being able to assess that, at least in a number of cases, much (if not all) of the observed time variability can be explained in terms of the motion of a boundary with a fixed structure, suggests that no other time-dependent physical processes are at work that significantly alter the structure of the boundary.This is an important physical result.
-Knowing the position of the magnetospheric boundary is an essential aspect of understanding the solar wind/magnetosphere interaction.Correlating boundary position with upstream solar wind conditions could add to our insight into how solar wind pressure drives boundary motion, and whether there may be other reasons for boundary motion (e.g.surface instabilities).Much of our current knowledge is based on correlating individual magnetopause crossings with solar wind data, while the long-duration reconstructed boundary position obviously holds much more information.
-The method gives information about the thickness of and possible density structure inside the LLBL (see also De Keyser et al., 2004).This is important for understanding mass transport across the magnetopause.
-In studying MP/BL passes on the dayside and on the magnetospheric flanks, we find that the magnetospheric boundary is often undulating rather gently, with the undulations propagating tailward as the back-and-forth motion is combined with the magnetosheath flow.To the extent that this propagation speed can be considered constant, the time axis in the plots of x mpbl (t) can be rescaled to a distance along the MP/BL, so that the curves represent the spatially undulating shape of the boundary.The method can thus be used to study the shape of the boundary layer, including surface wave amplitudes and wavelengths (if the amplitude over wavelength ratio is sufficiently small).

Appendix A Normalization of the target function
In order to control the optimization problem, the target function is based on three dimensionless constants: the proxy confidence parameter α, the guiding variable profile smoothness parameter λ, and the boundary velocity profile smoothness parameter σ .As these constants vary between 0 and 1, the relative importance of the terms in the target function change.One can only assign a meaning to the particular values of these constants if the target function terms F f , F n , F v , and F s are normalized.To this end, the normalization constants τ l f , τ l n , τ v , and τ s were introduced in expression (2) for the target function: with f k,l max = max i f k,l (t k,l i ); f k,l min , v k x,max , and v k x,min are similarly defined.
These constants make the order of magnitude of the terms independent of the number of spacecraft K, of the number of guiding variables L, of the number of sampling points, and of the overall measurement precision.Each constant τ l f measures the sum of squares of the range of a guiding variable relative to the measurement error.Normalizing G l f by this constant then produces a term that describes the deviation of observations from the guiding variable spatial profile relative to the range of that variable.Note that the measurement errors play a role in determining the relative importance of individual guiding variable observations, but the normalization removes the overall order of magnitude of measurement precision from the target function.Constant τ v similarly normalizes the proxy fitting term.The constants τ l n and τ s define reference values for the average squared second order derivatives, so that G l n /τ l n and G s /τ s are of order unity.The second order derivatives in the expressions for G l n , G s , and τ s are approximated by finite differences.The expression for the τ l n is an explicit estimate for the average squared second order derivative of a smooth transition from the minimum to the maximum guiding variable value across the spatial domain.

Appendix B Computing guiding variable profiles
The optimization problem requires that one can find, for any given boundary velocity profile v mpbl (t), the best fits to the guiding variable observations.With target function (2), this reduces to the well-known linear least-squares problem of determining a smooth spline fit.The construction of such a fit requires some computational effort and has to be repeated often in the overall optimization process; we therefore opted for using piecewise linear splines, for which such a fit can be computed relatively fast.In the present context, the linear µ × µ least-squares problem has the following form: Let the sets S k,l i group all observations for which the position of spacecraft k relative to the boundary x k,l j =x k (t k,l j ) is inside the interval [x i , x i+1 ].One can then express the matrix elements by Matrix A l f and vector b l f correspond to the F l f term in the target function.Matrix A n stems from the smoothness term F l n .Matrix A b and vector b l b force the fit to have the prescribed values at either end of the spatial interval.Note that system (B1) is band-diagonal; it can be solved efficiently.It is well-conditioned for λ→1, but may be singular if λ=0 (see also Property 3.1).

Fig. 1 .
Fig. 1.The optimization considers changes v mpbl (t) whose timeintegral is zero, the net effect of which is to displace the boundary in t i−1 , t i+m+1 .(a) Change v mpbl and the corresponding x mpbl .(b) Change over an interval that includes reference time t c ; adding an offset preserves x mpbl (t c )=0.
Fig. 2.Empirical reconstruction for Ampte/Irm MP/BL pass on 6 December 1984, 08:00-08:45 UT (x is the average outward normal, y is roughly sunward, and z is the invariant direction) computed with a time resolution t=5 s, using v ⊥x as the boundary motion proxy and ion density n i and magnetic field B z (with equal weights w l =0.5) as guiding variables, and with proxy confidence parameter α=0.8, with µ=60 points in the spatial profiles, with spatial resolution x=500 km, with guiding variable profile smoothness parameter λ=0.0001, and with boundary velocity profile smoothness parameter σ =0.0001.Top panels: Time profiles with observations (markers) and predictions from the reconstruction (curves).Bottom left: Spatial profiles.Bottom right: Convergence history; the problem was solved in 10 steps, with a subproblem precision of 10 −3 and a final precision of 10 −5 .

Fig. 3 .
Fig.3.Reconstruction for the full Ampte/Irm MP/BL pass on 6 December 1984, 06:30-09:30 UT, with a time resolution t=5 s, using v ⊥x as the boundary motion proxy and n i and B z (equal weights w l =0.5) as guiding variables, and with proxy confidence parameter α=0.8, with µ=80 points in the spatial profiles, with a spatial resolution x=500 km, with spatial profile smoothness parameter λ=0.001, and with boundary velocity profile smoothness parameter σ =0.0001.Top panels: Time profiles with observations (markers) and predictions from the reconstruction (curves).Bottom left: Spatial profiles.Bottom right: Convergence history; the problem was solved in 15 steps, with a subproblem precision of 10 −3 and a final precision of 10 −5 .

Fig. 4 .
Fig. 4.Reconstruction for the full Ampte/Irm MP/BL pass on 29December 1984, 17:30-19:15 UT, computed with a time resolution t=5 s, using v ⊥x as the boundary motion proxy and ion density n i (weight 0.4), ion temperature T i (weight 0.4), and magnetic field B z (weight 0.2) as guiding variables, and with proxy confidence parameter α=0.8, µ=80 points in the spatial profiles, a spatial resolution x=500 km, guiding variable profile smoothness parameter λ=0.001, and boundary velocity profile smoothness parameter σ =0.0001.Top panels: Time profiles with observations (markers) and the predictions from the reconstruction (curves); data gaps are highlighted.Bottom left: Spatial profiles.Bottom right: Convergence history; the problem was solved in 5 steps, with a subproblem precision of 10 −3 and a final precision of 10 −5 .

Fig. 5 .
Fig.5.Empirical reconstruction for an outbound MP/BL pass by Cluster on 6 June 2001, 01:10-03:30 UT, computed with a time resolution t=20 s, using v ⊥x from C1 and C3 as the boundary motion proxy, with electron density n e and magnetic field B z (equal weights w l =0.5)) as guiding variables, and with proxy confidence parameter α=0.5, with µ=60 points in the spatial profiles, with a spatial resolution x=500 km, with guiding variable profile smoothness parameter λ=0.001 and boundary velocity profile smoothness parameter σ =0.0001; the optimization used 5 steps, with a subproblem precision of 10 −3 and a final precision of 10 −5 .Top panels: Time profiles of electron density observations (PEACE, C1-C4); n e predictions from the reconstruction (C1-C4); magnetic field strength observations (FGM, C1-C4);

----
It computes accelerated boundary motion; -It covers a long time interval (hours); It considers multiple crossings of the magnetospheric boundary; It can use a variety of data (not necessarily magnetic field and plasma pressure); It can exploit multi-spacecraft data even for moderately large spacecraft separation;