Reconstruction of the magnetopause and low-latitude boundary layer topology using Cluster multi-point measurements

Multi-spacecraft data from Cluster allow for a more detailed magnetopause and boundary layer structure determination than ever before. Reconstruction methods, in which the time variability observed during a pass is interpreted as being due to boundary motion and/or spatial structure, are particularly well suited for this task. Such methods rely on the availability of plasma and field data and adopt the tangential discontinuity hypothesis to determine the motion, acceleration, boundary structure, boundary curvature and surface wave speed over an extended time interval. In this paper one- and two-dimensional reconstruction methods are applied to multi-spacecraft data for the first time.


Introduction
The solar wind−magnetosphere interaction leads to the formation of a magnetospheric boundary that consists of the magnetopause (MP) and often, but not always, a boundary layer (BL). This position of this boundary is known to change with time. If the boundary moves gradually, it may be considered locally to be a moving planar slab. Rapid boundary motion, together with the tailward flow of magnetosheath plasma near the boundary, leads to surface waves. The ratio of the amplitude of the motion relative to its wavelength determines whether a slab (one-dimensional, 1-D) or a surface Correspondence to: J. De Keyser (johan.dekeyser@bira-iasb.oma.be) wave (two-or three-dimensional, 2-D or 3-D) geometry is appropriate.
This paper deliberately does not address the question of the origin of boundary motion, but focuses on methods for determining that motion from in situ spacecraft observations. It is then possible to deconvolve the data to obtain the spatial structure of the boundary. Such methods (e.g. De Keyser et al., 2002;De Keyser and Roth, 2003) are known as empirical reconstruction techniques. They should be distinguished from magnetic field-based reconstructions (Walthour et al., 1993;Hau and Sonnerup, 1999;Hu and Sonnerup, 2003). They are fairly recent; we apply them here in their multispacecraft form to Cluster data for the first time. This is especially interesting because the four spacecraft provide a rather dense coverage of the MP/BL structure. Empirical reconstructions put data from different instruments in the same topological context. This provides us with a remarkably complete picture of boundary structure.

Reconstruction techniques
Spacecraft observe a lot of variability when they cross the magnetospheric boundary, as exemplified by the fact that many passes show multiple MP crossings. A major goal of the Cluster mission is to disentangle space and time variations. In the reconstructions discussed here, it is assumed that the observed structure has no intrinsic time variability: The observed variability is due to spatial structure that convects across the spacecraft. The result of the reconstruction process is a picture of this spatial structure. The essence of reconstruction consists of locating observations in a reference frame appropriate for that structure. A first step is to orient the frame so that it is more or less aligned with the structure. A second step is to consider motion of the reference frame relative to the structure.
Traditionally, an MP/BL-aligned frame is given by the "boundary normal coordinates" based on the minimum variance analysis (MVA) of some vector quantity, most often the magnetic field (Sonnerup and Scheible, 1998). The idea behind this type of frame is that locally (in both space and time) an individual MP crossing can be considered planar. Reconstruction, however, typically uses data from an extended period of observations, corresponding to relatively large spatial structures. A generalization of the traditional boundary normal coordinates is therefore needed. One way to obtain such a frame is the following (De Keyser et al., 2002). Given the time series of magnetic field observations B(t k ), and assuming that the structure is fieldaligned as in a tangential discontinuity, the local normal to the structure is n(t k )=B(t k )×δB(t k / ||B(t k ×δB(t k ||, where δB(t k )=(B(t k+1 )−B(t k−1 ))/(t k+1 −t k−1 ) is the central difference of the magnetic field time series. MVA of n(t k ), possibly adding the constraint <n z >=0, and sorting the variance directions x , y , z in descending order of variance magnitude, often shows that the variance along z is significantly smaller than in the other directions, implying that the structure can be regarded as 2-D: The structure can then be projected onto the x y plane. As will be shown in the examples, it is usually possible to find an appropriate rotation around z that transforms the x y axes into xy, where x is the average outward boundary normal and y is along the average tangential direction (while z is always tangential to the boundary). The xyz frame is an appropriate generalized boundary normal coordinate frame.
The following choices concerning reference frame motion are now possible: -For 1-D structure: Let the xyz frame comove with the boundary along x. A proxy for the boundary position is where v x is the in situ plasma velocity in the normal direction, and where ξ (t o ) defines a reference position at an arbitrary time t o ; this proxy is based on the idea that the plasma near the boundary comoves with it (Paschmann et al., 1990). This integration is not easy to carry out in practice because of the accumulation of errors; some difficulties have been discussed by De Keyser et al. (2002). Note that there are several ways to verify boundary motion: One can check whether the electric drift in stationary uniform field regions (from the Cluster EFW or EDI instruments) matches the measured plasma velocity (from the Cluster CIS/HIA ion spectrometer), or one can compare with boundary velocities obtained from the relative timing of sharp boundary crossings by the four spacecraft.
-For a 2-D structure: Let the reference frame move along y with the surface wave speed v wave . Rather than interpreting time variability in terms of motion along x, it is then ascribed to spatial structure being convected past the spacecraft in the y direction.
-For a 2-D structure with global motion: A more realistic reconstruction can be obtained by combining the two previous techniques: Let the frame move in y with v wave and in x by following only the low frequency global motion found from ξ (t) (De Keyser and Roth, 2003). The rationale behind this is that a simultaneous global compression or expansion of the magnetosphere can accompany surface waves. This approach is particularly useful for determining the structure of periodic surface waves, where one can easily assess how the boundary moves from period to period.
In practice, a 1-D reconstruction is first attempted. Using the 1-D comoving frame defined above, one can check whether the structure can be represented with profiles that depend only on x, that is, if a single-valued profile is obtained when plotting all observations as a function of x. If not, at least one additional dimension is needed to accurately represent the structure. Assuming stationarity and invariance along z, we therefore introduce variability along y by locating the observations in a 2-D comoving frame as discussed above.

An example of 1-D reconstruction
As an example of 1-D reconstruction we consider a dawn side inbound pass by the Cluster spacecraft on 23 April 2001, between 14:00 and 17:00 UT ( Fig. 1), at about 07:30 local time. The spacecraft separations were of the order of 1000 km. While the magnetopause is crossed around 14:26 UT, several transient variations are seen both in the magnetic field recorded by the fluxgate magnetometer (FGM) (Balogh et al., 1997) and in the plasma density measured by the ion spectrometer (CIS/HIA) (Rème et al., 1997). The nature of these transients is not immediately clear. They could be the signature of detached entities of magnetosheathlike plasma inside the magnetosphere, of structures traveling along the MP/BL, such as surface waves or flux transfer events, or they could simply be successive partial MP/BL encounters due to in-and outward motion of the overall boundary.
The most pronounced transient occurs between 16:18 and 16:40 UT; the ion density observations n i (t) are shown in Fig. 2a. MVA of the local normals n(t k ) for the four spacecraft, with the constraint <n z >=0, results in the x y z frame. Consider a series of xyz frames obtained by rotation in the x y plane. For each rotation angle, the boundary position ξ (t) is evaluated by integrating the CIS/HIA plasma velocity (from spacecraft 1 and 3 on which CIS/HIA is operational) and the spatial ion density profile n i (x sc (t)−ξ (t)) is plotted. For a particular rotation angle this curve appears to be more or less single-valued (Fig. 2b): The transitions from magnetosphere-like into magnetosheath-like plasma and back again overlap, and the curves for both spacecraft overlap as well. The structure can therefore be considered one-dimensional. The corresponding reference frame is x= (  estimates the surface wave speed to be half of the 150 km/s tailward magnetosheath flow. The amplitude over wavelength ratio therefore is about 1:15, sufficiently small to justify the planar approximation. The single-valued spatial profile reflects the 1-D structure of the boundary. An empirical model profile is obtained by fitting a curve to these data, thus smoothing out small-scale spatial or temporal variability and reducing the effects of nonsystematic instrument errors. Given this empirical model, the time profiles that would re-sult from the known boundary motion can be predicted (not shown); the good agreement with the observed time profiles indicates that the difference between both can be fully explained in terms of the spacecraft separation along x only, confirming again the 1-D geometry of the MP/BL. Spatial profiles can be obtained in the same way for any instrument on board any of the spacecraft.

An example of 2-D reconstruction
A major coronal mass ejection on 22-23 November 2001 had its main shock hitting the Earth on 24 November at 06:00 UT, with solar wind pressure up to 80 nP and interplanetary magnetic field (IMF) up to 65 nT, and leading to a magnetic storm with D st = -250. One day later the geomagnetic activity had subsided to quiet conditions. At that time, on 25 November 2001, 00:00-04:00 UT, the Cluster configuration (separations ∼2500 km) made an outbound pass through the dusk flank MP/BL around 19:00 local time; the Earth was then immersed in a tenous, cold, low β, fast solar wind, with a velocity of 780 km/s. The IMF was strong and steady around 20 nT. Observations show rather regular, quasi-periodic, large-amplitude variations in all observed plasma and field parameters that persisted for hours (Fig. 3). Spectral analysis shows a period of about 4 min, especially early in the pass, while larger amplitude variations seem to recur every 10 to 12 minutes.
We employ the 2-D reconstruction method that accounts for global motion, first proposed by De Keyser and Roth (2003) but now applied with Cluster multi-spacecraft observations. As in the 1-D example, the ion velocity measured by CIS/HIA is used to obtain the motion of the MP/BL. Although those data are available only on spacecraft 1 and 3, they allow us to find ξ (t) during the time interval of interest and to use it as the basis for a reconstruction with data from all four satellites. The CIS/HIA data start from about 01:55 and show about 10 large-scale recurrent variations before the final MP passage around 03:11 UT, corresponding to an average wave period of about 11 min, spanning a time interval of about an hour and a half. Consider the two wave periods (N=2) in the time interval 01:53:00 to 02:15:00 UT (wave period T is 11 min.). Reconstruction with a small number of periods is less likely to be disturbed by some non-periodic phenomenon (or phenomena with a different period); the wave can more reliably be regarded as a stationary, convecting structure. Reconstruction with a large number of periods leads to a better coverage of the wave structure, but as the observation time interval becomes longer, the procedure to obtain the boundary position ξ (t) by integrating the normal velocity will fail at some point due to error accumulation.
Constrained MVA of the local surface normal is again used to establish an appropriate reference frame. We subsequently rotate this frame around z so as to point x along the average outward normal. A good first estimate for this rotation is found by trying to obtain single-valued curves as in the 1-D example (note that no clear single-valued 1-D profile is found, which is why the 2-D reconstruction method is applied). A rotation over −56 • works best (x=(0.279, 0.926, 0.253), y=(−0.825, 0.366, −0.430), z=(−0.491, −0.089, 0.867) in GSE).
A low-frequency filter has to be applied to ξ (t) to separate slower motion (to be treated as a global expansion or compression) from motion at frequencies at or higher than the wave frequency (interpreted as being due to spatial structure). As outlined by De Keyser and Roth (2003), such a filter is obtained by defining the low frequency motion ξ lf (t) to be the linear interpolant of the set ξ (t 0 +mT), m=0, ..., N. Indeed, all those points correspond to the same phase of the wave, so that this interpolant follows the overall boundary motion from period to period. The xyz frame therefore moves with speed dξ lf (t)/dt along x, and with a constant wave speed v wave along y. Every measurement made by any of the spacecraft can then be located in this reference frame. Because of the periodicity, all the y coordinates can be taken modulo L, where L=v wave T is the wavelength, to superpose the data from each spacecraft obtained during each pass through a wave cycle. Table 1 lists the typical magnetospheric and magnetosheath parameters in the comoving frame that will be used in the reconstruction. Note the very fast magnetosheath flow along the magnetospheric boundary, and hence the large flow shear across it (650 km/s and more). Figure 4 shows CIS/HIA ion density (only Cluster 1 and 3), PEACE electron density (Johnstone et al., 1997), FGM magnetic field strength, WHISPER electric wave intensity in the 2-80 kHz band (Décréau et al., 1997), STAFF perpendicular magnetic wave intensity in the 0.3-10 Hz band (Cornilleau-Wehrlin et al., 1997), EFW spacecraft potential (Gustafsson et al., 1997), and the CIS/HIA plasma velocity along x and y in the nonmoving frame. All data were used at spin resolution, about 4 s. (The 213 ms resolution WHISPER data were resampled after applying a 0.25 Hz low band-pass filter.) The periodicity is most evident in the ion density, magnetic field strength, and wave intensity profiles. Some data, like the STAFF data, vary over more than 4 orders of magnitude. These large variations occur as the spacecraft passes alternately from the magnetosphere through the boundary layer into the magnetosheath and back. Such crossings through the MP/BL are not always complete: In the first period the spacecraft never make it to the magnetosheath proper, as the magnetosheath ion density is not reached. One can verify that the MP (manifest as a modest change in B y and B z ; not shown) is not seen during the first period, but it is crossed briefly in the second period.
The wave speed v wave must lie between the magnetospheric and the magnetosheath y velocity, that is, between 0 and 600 km/s. The reconstruction method offers a way to determine its value more precisely, as will be discussed later on; for the moment v wave =200 km/s will be adopted. This implies a wavelength of about 21 R E . Locating the data points from all spacecraft in the xy plane leads to an irreg-ularly distributed set of points (De Keyser and Roth, 2003). These data are supplemented with the magnetospheric and magnetosheath values from Table 1, imposed at a distance of 4000 km away from the inner-and outermost data points, respectively. A Delaunay triangulation of this set of points is then created (this is a particular triangulation that covers the xy plane with triangles that are as "small" as possible, where "small" has to be taken in a specific geometrical sense). Linear interpolation per triangle is then used to obtain the data on a fine mesh. Finally, those data are averaged onto a coarser mesh with physically relevant spatial resolutions, in this case 2000 km in x and 4000 km in y. We work with two different triangulations (Fig. 5): one for four spacecraft quantities and one for the two spacecraft CIS data. The four spacecraft triangulation in particular gives a pretty good coverage of the wave structure. Figure 6 plots the reconstructions in the xy plane. These reconstructions must be interpreted as a synoptic map of the observations, rather than as a spatial structure per se: Given the long wave period and the high wave speed, the wavelength turns out to be so long that boundary curvature cannot be neglected; for shorter wavelengths, the reconstruction could be considered as a representation of the actual wave geometry. Figure 6a gives the reconstructed ion density from CIS/HIA (2 spacecraft). Figure 6b shows the PEACE electron density. As these data are available from the four satellites, the reconstruction is correspondingly more detailed. The intensity of electric waves in the 2-80 kHz band from WHISPER and of perpendicular magnetic waves in the 0.3-10 Hz range from STAFF is given in Fig. 6c and 6d, respectively. The EFW spacecraft potential, which is correlated with the ambient plasma density, is reconstructed in Fig 6e. Figure 6f shows the magnetic field strength reconstruction; Fig. 6g shows the magnetic field strength overplotted with the xy projections of the field vectors. At this location at the dusk side boundary the magnetic field strength in the magnetosheath (38 nT) is actually higher than that in the magnetosphere (21 nT). The surface wave indentation is clearly visible; note that some magnetic field compression takes place inside the wave indentation as the field goes up to 43 nT there. Note also the magnetic field depression inside the wave structure that corresponds to the MP current layer; it can be identified in Fig. 6i as the place where the magnetic field changes orientation abruptly. Figure 6h plots the reconstructed CIS/HIA xy velocity vectors superposed on its color-coded v z component. (The v xy vectors are transformed to the comoving frame: v x has been corrected for low frequency motion and v y has been corrected by the wave speed). This plot allows to determine the actual wave speed by requiring that the flow should be tangential to the boundary in the comoving frame. Indeed, in this plot the flow vectors are more or less parallel to the isolines for v z . Note in particular how this condition is fulfilled for the flow vortex in the indentation of the magnetosheath into the magnetosphere. Moreover, the location of the flow reversal in y coincides with the flow shear in the z direction. This suggests that the choice v wave =200 km/s is reasonable. Qualitatively similar reconstructions are found for wave speeds from 150 to 300 km/s; it is therefore hard to determine the phase speed more accurately. For shorter wavelengths and/or for a larger separation between the spacecraft in the y direction the above criteria would become more strict. The wave speed can also be found from deHoffmann-Teller analysis for subintervals (Khrabrov and Sonnerup, 1998); no deHoffman-Teller frame would exist for the extended time interval if there is an important overall boundary motion. Note also how the reconstructions, despite the completely different nature of the various kinds of data, all reflect the same overall wave geometry. As in the 1-D case, the empirical models of Fig. 6 can be used to reproduce the time profiles that would be seen by the spacecraft. Comparison of these model profiles with the actual observations turns out to be quite satisfactory (not shown).
We conducted an independent a posteriori analysis of the structure by means of a local analysis technique. This technique determines the normal direction (in the xy plane) and the velocity and the acceleration along that normal, of planar features observed by all spacecraft. It is not always easy to identify corresponding features because often smaller scale structures seem to be superposed on the overall structure. Also, the assumption of planarity on the spacecraft separation scale is not always satisfied. About 20 features were identified in the magnetic field strength and the spacecraft potential data ; those where the technique reported a significant acceleration (possibly due to non-planarity) were discarded. The normal directions agree with the overall wave picture. The boundary normal velocity can be converted into a tangential speed (if the features move only along y): Values are found of 30 km/s in the magnetosphere, around 150 km/s near the boundary layer inner edge, and up to 250 km/s and more on the magnetosheath side. This fits in the surface wave picture, but suggests that v wave =150 km/s might have been a better choice.

Conclusion
We have demonstrated the feasibility and utility of reconstruction methods that use multi-spacecraft data for the first time. Such methods rely on the possibility to track the motion of the boundary by measuring the plasma velocity in the vicinity of the boundary; because of error accumulation while integrating the normal velocity, the applicability of these methods is limited. Also, these methods assume that the structure of the boundary is intrinsically stationary. The 1-D reconstruction method described here is fairly general. It obviously fails if the structure cannot be considered to be planar. Improvements are possible when one uses modelbased reconstruction methods, such as the one described by De Keyser et al. (2002). The 2-D reconstruction method separates low frequency motion from the periodic motion, in the assumption that there is no motion of the boundary at comparable or higher frequencies (unless "noise" that is filtered out by spatial averaging); the validity of this assumption may be questionable as well.
Reconstruction techniques with multi-spacecraft observations are particularly interesting because they combine a vast amount of information into a global picture of the boundary. Including the data from several instruments provides an additional bonus in terms of completeness of the picture. This may be a starting point to improve our understanding of the physical processes at work.