Magnetic Field Reconstruction for a Realistic Multi-Point, Multi-Scale Spacecraft Observatory

Future in situ space plasma investigations will likely involve spatially distributed observatories comprised of multiple spacecraft, beyond the four and five spacecraft configurations currently in operation. Inferring the magnetic field structure across the observatory, and not simply at the observation points, is a necessary step towards characterizing fundamental plasma processes using these unique multi-point, multi-scale data sets. We propose improvements upon the classic first-order reconstruction method, as well as a second-order method, utilizing magnetometer measurements from a realistic nine-spacecraft observatory. The improved first-order method, which averages over select ensembles of four spacecraft, reconstructs the magnetic field associated with simple current sheets and numerical simulations of turbulence accurately over larger volumes compared to second-order methods or first-order methods using a single regular tetrahedron. Using this averaging method on data sets with fewer than nine measurement points, the volume of accurate reconstruction compared to a known magnetic vector field improves approximately linearly with the number of measurement points.


INTRODUCTION
Plasmas, which are ubiquitous throughout the universe, are readily available for study in the natural laboratory of space. Electromagnetic fields play a fundamental role in the transport, heating, and acceleration of charged particles that compose plasmas. In order to characterize fundamental processes governing heliospheric plasmas, the space plasma community has utilized in-situ spacecraft measurements of electromagnetic fields and charged particles. These in-situ measurements include the characterization of the vector magnetic field B at a spacecraft via magnetometers; see §2.4 of Verscharen et al. (2019).
Knowledge of B from a single magnetometer is limited; single-point measurements can not construct the full three-dimensional structure characteristic of processes such as magnetic reconnection and plasma turbulence. To avoid this shortcoming, ESA's CLUSTER (Escoubet et al. 2001), NASA's THEMIS (Angelopoulos 2008) and MMS (Burch et al. 2016) missions have employed four-and five-spacecraft configurations, where each spacecraft is equipped with an instrument suite that includes a magnetometer. These missions study the boundaries of the Earth's magnetosphere, including how magnetic reconnection transfers magnetic energy into kinetic energy of plasma particles.
Analysis techniques have been created for multispacecraft missions, such as CLUSTER, which search for specific types of plasma waves (Constantinescu et al. (2006)) and which analyze current sheet structure (Narita et al. (2013)) for a four spacecraft configuration. Knowledge of the direction of wave propagation allows us to use multi-spacecraft filtering (Pinçon & Motschmann (1998)) to determine the general polarisation properties of any multi-point measurement of a wave field in space plasmas. Measurements from exactly four spacecraft (e.g. a tetrahedron of spacecraft) can be used to estimate current density via the Curlometer technique (Robert et al. (1998)).
The Cluster and MMS missions have also utilized the Curlometer technique to interpolate the value of the magnetic field over a region near the tetrahedron's barycenter, regardless of the field's geometry. However, these interpolations are limited to measuring fluctuations on a scale on the same order as that of their interspacecraft distances (e.g. Robert et al. 1998;Forsyth et al. 2011). To study multi-scale processes, such as plasma turbulence, with structures on characteristic length scales that cover many orders of magnitude, we must employ measurements from more than four spacecraft. Therefore, we develop a method which extends the magnetic field reconstruction technique Curlometer to configurations of more than four spacecraft.
Many such multi-spacecraft missions have been proposed, e.g. Cross-scale (Schwartz et al. 2009), AME (Dai et al. 2020) and HelioSwarm ), but in order to optimize such missions, it is urgent to robustly quantify the impact of particular spacecraft configurations on multi-point analysis methods, capturing the effects of the physical scales spanned by the spacecraft in the observatory and the geometry of the polyhedra that can be drawn from the constituent spacecraft. Such quantification will help demonstrate that a proposed mission will be able to usefully analyze a large number of magnetometer measurements made in the pristine solar wind, magnetosphere, and magnetosheath. It will also assist in the optimization of spacecraft configurations and quantification of errors derived from multi-point, multi-scale measurements. In this paper, we focus on the fidelity of the reproduction of the magnetic field using a sparsely sampled set of measurements whose spatial configuration is based upon realistic configurations of the proposed nine-spacecraft HelioSwarm observatory, described for instance by Plice et al. (2020).
The reconstruction method is described in §2, the results are applied to two magnetic field models, including a numerical simulation of turbulence, in §3, with a concluding discussion in §4.

Geometrical Definitions
Given N spacecraft, we identify C(N, k) polyhedra with k vertices. As spatial divergence analysis methods, (e.g. Paschmann & Daly 1998, 2008Dunlop et al. 1988) require at least four vertices to resolve three-dimensional structure, we only consider polyhedra with at least four vertices, known as tetrahedra. For N = 9, there are 126 (i.e. 9 choose 4) tetrahedra, 126 polyhedra with 5 vertices, 84 with 6 vertices, 36 with 7 vertices, 9 with 8 vertices, and 1 with 9 vertices, for a total of 382 polyhedra with at least 4 vertices.
Each polyhedron is characterized in terms of its size and shape. Because measurements from all d spacecraft are weighted equally, we define the barycenter of the q th polyhedron with set D of d vertices drawn from the N ≥ d spacecraft positions x i as Given the barycenter, we then define the volumetric tensor of the q th polyhedra with set D of d vertices as Here x ij represents the j th ∈ {x, y, z} component of the position vector for the i th spacecraft. The eigenvectors of the tensor R q,d represent the three semi-axes of the polyhedra and are associated with the eigenvalues and c q,d = R q,d 3 (minor axis), where a ≥ b ≥ c (a more detailed analysis of the eigenvalues can be found in Chapt 12 of Paschmann & Daly 1998).
To provide a useful geometric interpretation of these shapes, we define a characteristic size L, as well as an elongation E and a planarity P (see chapter 16.3 of Paschmann & Daly 1998) 1 :

First-Order Method
In a first-order Taylor series expansion, we use the values of the magnetic field, B, measured at four spacecraft positions x i to estimate the value of B (and its corresponding directional derivatives) at any other point in space, ξ (Fu et al. 2015(Fu et al. , 2020. The Taylor expansion is: ∀i ∈ {1, 2, 3, 4}, m ∈ {x, y, z}.
In this equationB i m is the measured m th component of B at the i th spacecraft, B m is the computed m th component of B at ξ, ∂ k B m is the computed derivative of the m th component of B with respect to the k th direction at ξ, and r i k is the relative position of spacecraft i with respect to ξ. In other words, if x ik is the k th component of spacecraft i's location, then r i k := x ik −ξ k . This is a system of 12 equations with 12 unknowns, where the 12 equations represent the x, y, and z components of B for each of the four spacecraft. The 12 unknowns are the x, y, and z components of B at ξ and the nine terms in the Jacobian of B at ξ.
This system can be reformatted into linear (Ax = b) form and solved with a common linear system solver. This 12-dimensional linear system (shown in full detail in Appendix A, Equation A8) comprises the first-order reconstruction method.
This magnetic field reconstruction method is related to the Curlometer method (Dunlop et al. 1988;Robert et al. 1998), which utilizes Ampère's law to calculate the current density J as the curl of B. The Curlometer solves the same set of equations, but uses the partial derivatives to estimate the current density at the center of each tetrahedron. The Curlometer method has been widely applied to four-spacecraft magnetic field measurements made for instance by Cluster and MMS, (c.f. Chapter 16.2 of Paschmann & Daly 1998). Future missions, such as the proposed HelioSwarm Observatory , will have more than four spacecraft. Therefore, for every reconstructed point, ξ, we can apply this reconstruction method for each of the C(N, 4) tetrahedra and average the reconstructed values, yielding a statistically larger base of estimates and improving the accuracy of the reconstruction.

Second-Order Method
Because the proposed HelioSwarm Observatory has nine spacecraft, we can use measurements of B from all nine spatial points simultaneously to apply a secondorder reconstruction method. This method, also based on a Taylor series expansion, is more accurate for values located near the center of the expansion (i.e. near the barycenter of the nine-spacecraft constellation) than a single implementation of the first-order method. Following the work of Torbert et al. (2020), we write: These terms are the same as in the first-order method, with the addition of ∂ j ∂ k B m , the second derivative of the m th component of B with respect to the k th and j th directions, at ξ. This is a system of 31 equations with 30 unknowns. 27 of these equations are associated with the x, y, and z components of B from the nine spacecraft. There are four additional constraints, imposed by the magnetic field having zero divergence, as well as the divergence of the magnetic field having zero gradient. The 30 unknowns are the x, y, and z components B at ξ, the nine terms in the Jacobian of B at ξ, and the 18 terms in the Hessian of B at ξ (excluding the 9 redundant terms).
This system can be reformatted into linear (Ax = b) form where A is a 31 × 30 matrix (shown in full detail in Appendix B, equation B16). This system is over-determined, therefore in general, an exact solution does not exist. However, we can find an approximate solution via the method of ordinary least squares. This method finds the solution to the problem Ax = b which minimizes the two-norm of the error, i.e.
This second-order reconstruction method is referred to as M 2 throughout this paper.

Quantifying Error
We define the error at any point in space, ξ, as: where B true (ξ) is the magnetic field vector at point ξ and B calc (ξ) is the reconstructed magnetic field vector at point ξ.
Given that we can determine the value of this error at all points in a simulation or for a given analytic field, we also define (θ) as the proportion of the volume that is reconstructed with less than θ% error. For a sufficiently dense grid of N x ×N y ×N z uniformly-spaced points, (θ) can be estimated as To define the physical volume in which a given magnetic field reconstruction is accurate, (θ) is translated into a dimensional quantity by multiplying it by the total volume covered by the N x × N y × N z grid.

Error Minimization Techniques
As the first-order method ( §2.2.1) only requires a single tetrahedron of spacecraft to estimate B, in this paper we will test four selection methods for using a subset of the 126 tetrahedra to improve the reconstruction. These methods combine the statistically large set of tetrahedra with our knowledge of the spacecraft positions relative to ξ and the geometry of all 126 tetrahedra.
For the first method, M 1.1 , at each point in space we reconstruct the magnetic field using all 126 tetrahedra to produce 126 estimates for B(ξ). We then average over these B(ξ) values component-wise to estimate B(ξ).
For the second method, M 1.2 , we perform the same averaging as method one, but only include tetrahedra whose barycenter are within a characteristic distance of ξ. i.e. for each reconstructed point ξ, only include tetrahedra j in the average which satisfy where L j is the characteristic size and (r 0 ) j is the barycenter of the j th tetrahedron. For the third method, M 1.3 , we perform the same selection as method two, but with the added restriction that the shape of tetrahedron j must be quasi-regular. In terms of the geometric quantities of the spacecraft configuration (defined in equation 3), this translates to elongation E and planarity P being sufficiently small. Because E and P are symmetric with respect to orientation, we will define a composite geometric parameter χ j Small χ j implies that both the elongation and planarity of tetrahedron j are small. For method M 1.3 , we restrict our averaging to only include tetrahedra where For the fourth method, M 1.4 , we perform the same selection of tetrahedra as method three, but require the tetrahedra included in the averaging to be more regular. For method M 1.4 , our shape and position requirement is The value of 0.6 was selected because page 408 of Paschmann & Daly (1998) shows it be a threshold value for elongation and planarity which separates the well performing 'pseudo-sphere type' and 'potato type' spacecraft configurations from the poorer performing 'knife blade type', 'cigar type', and 'pancake type' configurations.
The first-order methods M 1.1 , M 1.2 , M 1.3 , M 1.4 will be compared to the second-order method M 2 as well as the first-order method applied to a single regular (i.e. χ = 0) tetrahedron of spacecraft. This single regular tetrahedron will have the same characteristic scale as the nine-spacecraft configuration it is compared to.

Models
To validate and quantify the errors of our reconstruction, we implement our reconstruction methods on two magnetic field models, a simple current sheet and a numerical simulation of turbulence.

Simple Current Sheet
For our first model, we define a magnetic field where B is analytically defined at all spatial points. This field, which represents a simple current sheet, can be described in cylindrical coordinates as The variable σ represents the current sheet characteristic width and J 0 represents the magnitude of the current at its center.

Turbulence Simulation
Physically realistic fields, such as those generated by turbulence in the solar wind, are significantly more complex than the simple current sheet model of equation 13. We therefore test our reconstruction techniques on magnetic fields drawn from numerical simulations of turbulence. In particular, we utilize the magnetic fields from a fully developed turbulence simulation performed with the five moment, multi-fluid solver within the Gkeyll simulation framework (Hakim et al. 2006;Wang et al. 2015Wang et al. , 2020. This turbulence simulation is designed to represent plasma behavior in the pristine solar wind at 1AU. We use the five moment (n s , u s , p s ), two fluid (s = p, e) plasma model to evolve a proton-electron plasma. We note that the five moment, two fluid model formally reduces to Hall MHD in the limit m e → 0 and 0 → 0 (Srinivasan & Shumlak 2011), where 0 is the vacuum permittivity. We use a reduced (proton to electron) mass ratio of m p /m e = 100, a temperature ratio of T p /T e = 1, Alfvén velocity of v A /c = B/ µ 0 n p m p c 2 = 0.02, plasma beta (ratio of plasma thermal pressure to magnetic pressure) of β p = 2µ 0 n p T p /B 2 = 1, and adiabatic index γ = 5/3. We employ an elongated domain L x = L y = 0.2L z = 100πρ p with resolution n x = n y = n z = 448. Lengths are normalized to the proton gyroradius ρ p = v tp /Ω p , the ratio of the proton thermal speed v tp = 2T p /m p and the proton cyclotron frequency Ω p = q p B/m p . We choose a uniform background density and magnetic field, B 0 = B 0ẑ , and initialize the simulation with the three dimensional extension of the Orszag-Tang vortex (Orszag & Tang 1979)   where z ± = δu ± δB/ √ µ 0 ρ 0 are the Elsasser variables (Elsasser 1950), k x,y = 2π/L x,y , and k z = 2π/L z . The initial amplitude, z 0 = 0.2, is chosen to satisfy the critical balance condition, The simulation is run for one Alfvén crossing time, t A = 1500/Ω p , at which point the turbulence has fully developed and reached a steady state. In Fig. 1, we plot the trace magnetic energy spectrum as a function of k ⊥ = k 2 x + k 2 y , with a k −5/3 ⊥ dashed line plotted for reference. The steep roll-over in the spectrum at k ⊥ ρ p 1 is due to numerical diffusion from the finite volume scheme employed by Gkeyll.
To compare the simulation to the selected spacecraft configurations with separations in physical units, we note that the proton gyroradius can be written as With the constants in the turbulence simulation of m p /m e = 100, β p = 1, ω pe = c de = 5.64 × 10 4 √ n, we set n e = 0.2829 cm −3 , so that ρ p = 100 km.
We extract from this simulation a 3-dimensional grid of values representing the plasma's physical parameters at different points in space. From this grid, we use trilinear interpolation to estimate the value of B at any point in the simulation volume.

Spacecraft Configurations
To illustrate our reconstruction methods for realistic spacecraft configurations, we study these methods using Elongation ( three different nine-spacecraft configurations. The spacecraft configurations are selected from the phase A design reference mission (DRM) of the proposed HelioSwarm Observatory concept, corresponding to hours 94, 144, and 205 of the science phase. These hours are selected because they represent a selection of spacecraft tetrahedra that have significantly different distributions of their elongation, planarity, and length. In Table 1 we note the geometric characteristics of the overall nine-vertex polyhedra for each of the three configurations.
We also calculate the size, elongation, and planarity of all 126 tetrahedron in each configuration and display them in Fig. 2, noting the minimum and maximum values of these three parameters for each configuration in Table 1. The wavelengths associated with the overall, minimum, and maximum scales, kρ p = 2πρ p /L, are overlaid on Fig. 1, using a fiducial value of ρ p = 100 km.

APPLICATION OF RECONSTRUCTION
To find the expected error at all points in space near a particular spacecraft configuration, we take a Monte Carlo approach and place the barycenter of each nine-spacecraft configuration into a known magnetic field at random locations. We then reconstruct the magnetic field on a grid of points centered at the barycenter of the nine-spacecraft configuration using the first-and second-order reconstruction methods. The location of each point in the reconstructed grid is constant with respect to the spacecraft configuration. Therefore, we find the average of the errors, θ, at all reconstructed grid points for all elements of the Monte Carlo ensemble, allowing the calculation of the expected value of error at each point on the grid.
Additionally, we compare the divergence found on a grid of points sampled from the baseline current sheet and turbulence simulation magnetic fields with that of the same points sampled from the fields reconstructed using our first-order reconstruction methods. This comparison yields divergence values of similar magnitude in the baseline and reconstructed fields, which indicates that our reconstruction methods do not introduce nonphysical values of divergence.

Current Sheet
We present an example magnetic field reconstruction of the simple current sheet model ( §2.3.1) in Figure 3. Here, we use the first-order method M 1.3 to reconstruct the magnetic field in the z = 0 plane for the simple current sheet, Eqn13 with σ = 2000 km using the hour 94 spacecraft configuration. There is little difference between the reconstructed and original fields near the center of the spacecraft configuration, and the difference in vectors increases with distance from the spacecraft configuration's center.
We perform 200 Monte Carlo iterations of reconstruction using each method, observing that 200 was sufficient to point-wise converge in error. The characteristic width of the current sheet is chosen as a uniform random variable σ ∼ U [500, 5000] km, while the barycenter of the nine-spacecraft configuration is selected as a 3D uniform random variable r 0 ∼ U [−1000, 1000] 3 km. We reconstruct a 30 × 30 × 30 grid of points ξ that extends 100 km past the furthest spacecraft in all directions.
The errors computed for the reconstruction of the simple current sheet are displayed in Figures 4, 5, and 6 for the hour 94, 144, and 205 configurations respectively. These figures illustrate the ensemble-averaged errors along a 2D plane orthogonal to the current intersecting a given nine-spacecraft configuration's barycenter. The first four panels correspond to the four first-order reconstruction methods, M 1,1 , M 1,2 , M 1,3 and M 1,4 , the fifth panel corresponds to the second-order method M 2 , and the final panel corresponds to the reconstruction obtained from the standard first-order method applied to a single regular tetrahedron, with E = P = 0. This single tetrahedron has the same characteristic size L as the overall nine-spacecraft configuration, calculated as twice the major axis of the volumetric tensor, Eqn 2, evaluated using all nine points. With four spacecraft, we cannot reconstruct the magnetic field with the secondorder method, nor can we select subsets of tetrahedra with advantageous geometric characteristics, so only the first-order reconstruction method from a single tetrahedron is used. We see that near the barycenter of each of the nine-spacecraft configurations (located at the origin of Figures 4, 5, and 6) the magnetic field can be reconstructed to within 1% accuracy. By comparing method M 1.1 with methods M 1.2 and M 1.3 in these figures, we also conclude that leveraging knowledge of the tetrahedral shapes and positions expands the region of high-accuracy reconstruction. Unfortunately, overly restrictive conditions limit the number of tetrahedra available to average over, limiting the size of the reconstructed region. In fact, the bottom left panel of Figure 6 is empty because none of the 126 tetrahedra in the hour 205 configuration satisfy the geometric requirement that χ j ≤ 0.6 demanded by M 1.4 . Additionally, the second-order reconstruction method M 2 is accurate for only a small volume when  Figures 4-10, we see that the behavior of methods M 1.1 , M 1.2 , M 1.3 , M 1.4 , and M 2 is distinct to that of the reconstruction using a single regular tetrahedron. The single regular tetrahedron only accurately reconstructs the magnetic field of the current sheet near each of the four spacecraft. Due to the angular symmetry in the current sheet and the fact that none of the four spacecraft are positioned on the z = 0 plane, the area of most accurate reconstruction appears to be a ring on the bottom right panel of Figures 4, 5, and 6.

Turbulence Simulation
We present an example magnetic field reconstruction of the turbulence simulation ( §2.3.2) in Figure 7. Here, we use the first-order method M 1.3 to reconstruct the magnetic field in the z = 0 plane in the turbulence simulation. Mirroring the behavior described in Figure  3, there is little difference between the reconstructed and original fields near the center of the spacecraft configuration. We perform 50 Monte Carlo iterations of reconstruction using each method, observing that 50 was more than enough to point-wise converge in error. The barycenter is chosen as a uniform random variable so that all spacecraft remained in the 31415 × 31415 × 157079 km simulation cube. We then construct a 30 × 30 × 30 grid of points ξ. Each dimension of this grid is selected so that the overall size of the grid extends 100 km past the furthest spacecraft in all directions.

Simulation Reconstruction
The errors computed from the turbulence simulation reconstruction are displayed in Figures 8, 9, and 10 for the configurations at hours 94, 144, and 205 respectively. The panels shown are organized in the same order as the simple current sheet reconstruction. In Table 2 we show the volume (in units of 10 6 km 3 ) of the magnetic field that can be reconstructed with errors less than 1%, 5%, and 10%. This is done for all three of the investigated spacecraft configurations, and using all five of the ninespacecraft reconstruction methods, M 1.1 , M 1.2 , M 1.3 , M 1.4 , and M 2 . In the bottom half of this table, we compare the volume reconstructed using a single regular tetrahedron to that of our five reconstruction methods.
Near the barycenter of each of the nine-spacecraft configurations, located at the origin of Figures 8, 9, and 10, the magnetic field can be reconstructed to within 1% accuracy. The second-order method, M 2 , can only reconstruct the magnetic field to within 10% accuracy in a small region near the barycenter of the configuration, while the first-order methods can reconstruct the magnetic field within 10% over a much greater area. This is the case because the second-order Taylor series expansion diverges quadratically with distance away from the barycenter of the spacecraft configuration, while the first-order Taylor series only diverges linearly with distance. Since our goal is to maximize the volume of accurate reconstruction, the first-order methods are superior. However, the secondorder method may be more accurate at reconstructing the values of the magnetic field very close to the barycenter of a spacecraft configuration.
The largest disparity compared to the current sheet simulations occurs for the single regular tetrahedron case, on the bottom right panels of each figure. The magnetic field is again only reconstructed accurately near each of the spacecraft, but because the turbulence simulation lacks angular symmetry, these regions manifest as spheres centered around each spacecraft. These four spheres appear to be the same size on the bottom right panel of each figure because the spacecraft are equidistant from the z = 0 plane. Shown in Table 2, the single regular tetrahedron reconstructs the largest volume with less than 10% error, however the first-order methods reconstruct larger volumes with smaller errors. To maximize the volume reconstructed with less than 1% error, it appears it is best to use the first-order method M 1.4 , detailed in   Table 2. Volumes (in units of 10 6 km 3 ) with reconstructed magnetic field error less than 1%, 5%, or 10% for the three configurations using the four first-order methods and the second-order method discussed in §2.2. These volumes are compared to the equivalent regions reconstructed from a single regular tetrahedron with the same characteristic size as the overall ninespacecraft configuration. §2.2.4 (if a sufficient number of quasi-regular tetrahedra can be formed from the nine spacecraft configuration).

Sensitivity to Number of Spacecraft
We analyze how the volume reconstructed with less than 5% error varies as a function of the number of spacecraft. This analysis was completed using the Monte Carlo sampling of the turbulent simulation as described in §3.
For N ∈ {4, 5, 6, 7, 8, 9} spacecraft, we reconstructed the value of the magnetic field at all 30×30×30 points ξ using all C(N, 4) tetrahedra. We then use the M 1.3 firstorder reconstruction method of §2.2.4 to reconstruct B at all points ξ. The errors everywhere are computed using equation 7, and the volume where the error is less than 5% is computed using equation 8 multiplied by the total reconstructed volume. We visualize the errors of this method for the hour 94 configuration in Figure 11. In this example, we find that as the number of spacecraft is increased, the area which is reconstructed with a high accuracy also increases. As this result depends on which particular subset of spacecraft are chosen for a given N, we next investigate whether this increase is holds for an arbitrary selection of spacecraft.
We start by choosing 4 out of the 9 spacecraft of the hour 94 configuration. These spacecraft measurements are used to estimate the value of B everywhere via the first-order reconstruction method M 1.3 . We find the volume over which we can reconstruct B with an error less than 5%. This process is repeated for all 126 possible choices of 4 spacecraft. We repeat all of these volume calculations, initializing the spacecraft configuration at 50 different locations within the simulated turbulent B field. Finally, we take the mean of all 126 × 50 volume values and plot them in Figure  12. In these averages, we omit the instances where no tetrahedra pass the selection criteria of method M 1.3 . We repeat this process for N = 5, 6, 7, 8, and 9 spacecraft from the hour 94 configuration, as well as for the hour 144 and 205 configurations.
As shown in Figure 12, we see that increasing the number of spacecraft measurements available increases the volume of the magnetic field reconstructed with less than 5% error. The variance of this reconstructed volume is smallest for the hour 94 configuration, which contains the most tetrahedra which are quasi-regular (χ j ≤ 1). However, it is not the case that the hour 94 configuration has the highest average volume which is reconstructed with less than 5% error.
We also track the instances where zero of the available tetrahedra in the set of N spacecraft meet the shape threshold of χ j ≤ 1 for the M 1.3 method. The  Table 3. From the nine-spacecraft configurations of hours 94, 144, and 205, we select a subset of N ∈ {4, 5, 6, 7, 8, 9} spacecraft. We determine the probability that this N spacecraft configuration does not contain a tetrahedron which passes the threshold shape requirements of first-order reconstruction method M1.3.
percentage of arrangements where this occurs is shown in Table 3 as a function of spacecraft configuration (hour) and number of spacecraft, N . We see from this table that for the analyzed configurations, there must be at least seven spacecraft measurements to guarantee that at least one tetrahedron passes the previously stated shape criteria.

DISCUSSION
We have demonstrated that our reconstruction methods are an effective way to leverage magnetometer measurements from a configuration consisting of more than four spacecraft. We have defined a shape metric, χ, for a tetrahedron of spacecraft which can be used as a threshold criterion. Estimates of magnetic field derived from tetrahedron which do meet the threshold value of χ will be discarded, as they are misshapen and therefore more likely to produce erroneous estimates. Finally, we have shown that increasing the number of spacecraft in a configuration will increase the volume over which the magnetic field can be accurately reconstructed, as well as increase the likelihood that some tetrahedra of spacecraft in the configuration are well shaped.
In Table 2 we demonstrated that our second-order reconstruction method M 2 does not reconstruct the magnetic field with high accuracy over a large volume. However, we have shown that methods M 1.3 and M 1.4 , which average over a subset of the many available tetrahedra formed by nine spacecraft, improves the field reconstruction. This work indicates that the subset of tetrahedra which should be averaged over needs to consider each tetrahedron's spacial proximity to the reconstructed point as well as its geometric properties. By comparing results from spacecraft configurations with different tetrahedral geometric configurations, we find that designing spacecraft trajectories which maximize the number of tetrahedra that are quasiregular (i.e. χ ≤ 1) is essential to improving the accuracy of the reconstructed magnetic field.
This work can help optimize future multi-spacecraft missions, such as HelioSwarm.
The selection of  Hour 94 Hour 144 Hour 205 Figure 12. Mean values of volume which were reconstructed with less than 5% error for the three nine-spacecraft configurations analyzed.
The dashes above/below the markers represent one standard deviation away from the mean volume for each configuration.
tetrahedra which are included in the calculation of B can be tuned to maximize the volume over which the field is reconstructed accurately, or it can be tuned to recreate B as accurately as possible over a small volume. The first-order methods discussed here can be applied to reconstruct any vector field which is sparsely sampled by in-situ measurements, as no assumptions are made about the physical properties of the field.
The first-order reconstruction method applied to a single tetrahedron reconstructs the magnetic field perfectly at each spacecraft location. However, using any of our proposed composite first-order reconstruction methods, which average over many of these reconstructions, negates this behavior. In future work, we plan to construct a weight function which, when introduced into the tetrahedral averaging, returns this desired limiting behavior. Additional future work could include characterizing methods of predicting the surface insideof-which we have less than a prescribed error value for an arbitrary configuration of spacecraft.
The authors would like to thank the HelioSwarm Science and Flight Dynamics teams for discussions and comments during the execution of this project, in particular Laura Plice and Jonathan Niehof. This material is based upon High Performance Computing (HPC) resources supported by the University of Arizona TRIF, UITS, and Research, Innovation, and Impact (RII) By writing out all terms in the sum, we see that this iŝ We now write this in vector notation asB If we repeat this process for all spacecraft i ∈ {1, 2, 3, 4} (only using component m), we find that we can combine the four vector equations into one matrix equation Now we need to generalize this to include all values of m. We first note that the matrix from equation A4 is independent of m. This independence means that if we replace m in equation A4 with x, y, or z, this matrix will stay the same. By doing this procedure, we see the three systems are: We now combine the three linear equations of A5, A6, and A7 into a single linear system of Ax = b form. This single system's independent variable x (and dependent variable b) will be a concatenation of the independent (and dependent) variables of systems A5, A6, and A7. The matrix of this combined linear system is a block matrix where the matrix from equation A4 is a block on the diagonal of larger (12 × 12) matrix. These three matrix blocks on the diagonal ensure that the matrix from equation A4 is being applied to the x, y, and z components independently. This system written out is:  x r 1 y r 1 z 0 0 0 0 0 0 0 0 1 r 2 x r 2 y r 2 z 0 0 0 0 0 0 0 0 1 r 3 x r 3 y r 3 z 0 0 0 0 0 0 0 0 1 r 4 x r 4 y r 4 z 0 0 0 0 0 0 0 0 0 0 0 0 1 r 1 x r 1 y r 1 z 0 0 0 0 0 0 0 0 1 r 2 x r 2 y r 2 z 0 0 0 0 0 0 0 0 1 r 3 x r 3 y r 3 z 0 0 0 0 0 0 0 0 1 r 4 x r 4 y r 4 z 0 0 0 0 0 0 0 0 0 0 0 0 1 r 1 x r 1 y r 1 z 0 0 0 0 0 0 0 0 1 r 2 x r 2 y r 2 z 0 0 0 0 0 0 0 0 1 r 3 x r 3 y r 3 z 0 0 0 0 0 0 0 0 1 r 4 x r 4 y r 4 We now solve this linear system and find the values of B x , ∂ x B x , ∂ y B x , ... to apply our first-order reconstruction method.

B. SECOND-ORDER SYSTEM
Making a few adjustments, the second order system can be derived in exactly the same manner as the first. We see that for spacecraft i and component m, equation 5 iŝ By writing out all of the terms in the sums, we seê We now use the fact that second derivatives are symmetric to see that there are three redundant terms in this equation. This is because Therefore equation B10 can be simplified tô We write equation B11 in vector notation We repeat this calculation for all nine spacecraft i and combine the vector equations in the following matrix equation: 2r 1 x 2r 1 y 2r 1 z r 1 x r 1 x 2r 1 x r 1 y 2r 1 x r 1 z r 1 y r 1 y 2r 1 y r 1 z r 1 z r 1 z 2 2r 2 x 2r 2 y 2r 2 z r 2 x r 2 x 2r 2 x r 2 y 2r 2 x r 2 z r 2 y r 2 y 2r 2 y r 2 z r 2 z r 2 z 2 2r 3 x 2r 3 y 2r 3 z r 3 x r 3 x 2r 3 x r 3 y 2r 3 x r 3 z r 3 y r 3 y 2r 3 y r 3 z r 3 z r 3 z 2 2r 4 x 2r 4 y 2r 4 z r 4 x r 4 x 2r 4 x r 4 y 2r 4 x r 4 z r 4 y r 4 y 2r 4 y r 4 z r 4 z r 4 z 2 2r 5 x 2r 5 y 2r 5 z r 5 x r 5 x 2r 5 x r 5 y 2r 5 x r 5 z r 5 y r 5 y 2r 5 y r 5 z r 5 z r 5 z 2 2r 6 x 2r 6 y 2r 6 z r 6 x r 6 x 2r 6 x r 6 y 2r 6 x r 6 z r 6 y r 6 y 2r 6 y r 6 z r 6 z r 6 z 2 2r 7 x 2r 7 y 2r 7 z r 7 x r 7 x 2r 7 x r 7 y 2r 7 x r 7 z r 7 y r 7 y 2r 7 y r 7 z r 7 z r 7 z 2 2r 8 x 2r 8 y 2r 8 z r 8 x r 8 x 2r 8 x r 8 y 2r 8 x r 8 z r 8 y r 8 y 2r 8 y r 8 z r 8 z r 8 z 2 2r 9 x 2r 9 y 2r 9 z r 9 x r 9 x 2r 9 x r 9 y 2r 9 x r 9 z r 9 y r 9 y 2r 9 y r 9 z r 9 z r 9 We repeat the matrix in equation B13 three times as blocks on the diagonal of a larger matrix, one for each component m ∈ {x, y, z}. However, because the matrix from equation B13 is 9 × 10, this results in a matrix which is 27 × 30. A linear system with a matrix of dimension 27 × 30 is under-determined, and therefore, cannot be solved. We remedy this by using the known physical properties of the magnetic field B. We know that magnetic fields have no monopoles, therefore the divergence of B and its gradient is zero, i.e. ∇(∇ · B) = 0 ∇ · B = 0.
Applying this fact, we must satisfy the following four equations: Because these constraints are in terms of the same independent variables as the overall linear system (ie B and its partial derivatives), we can include these four constraints directly into our linear system. This is done by concatenating four rows onto the bottom of our 27 × 30 matrix, and adding four zeros onto the bottom of our dependent variable vector. These four rows are: (B15) These four rows are constructed so that 1's are in positions corresponding to the partial derivatives in constraints B14, and there are 0's everywhere else.
In conclusion, by solving the linear system