This work aims to determine the spatial distribution and temporal evolution of soft X-ray emissivities (denoted as P) within the magnetosheath region. These emissivities are generated by charge exchange interactions between heavy solar wind ions and neutral hydrogen atoms from the terrestrial exosphere (see Eq. 1). The soft X-ray emissivity, also referred to as volumetric emission rate (VER), is expressed as: P =   α n H n s w ⟨ g ⟩ eV cm − 3 sec − 1 (2)

Here, α is the efficiency factor in units of [eV cm²], which includes solar wind charge exchange (SWCX) cross-sections, SWCX photon energies, and the ratio between the density of solar wind heavy ions to the density of solar wind protons (Whittaker and Sembay, 2016; Jung et al., 2022). Also, n _H indicates the density of exospheric H atoms in units of [atoms cm⁻³], n _sw is the solar wind proton density in units of [ions cm⁻³], and ⟨ g ⟩ = v s w 2 + v t h 2 [cm sec⁻¹] represents the relative velocity of neutrals and ions calculated with the plasma bulk speed v _sw and thermal speed vth under the assumption of negligible neutral velocity (Walsh et al., 2016; Sibeck et al., 2018).

An appropriate space-based photometric detector can provide routine observations of P such that these measurements can be used systematically to estimate their spatial distribution through inversion techniques. A single observation of P from a platform located at planetocentric distance r, viewing line-of-sight (LOS) look angle n ̂ , and acquired at time t can be expressed as: I r , n ̂ , t = 1 4 π ∫ 0 L max P r , t d l + I bkg n ̂ , t eV cm − 2 sec − 1 str − 1 (3) where I represents the incoming and directional photon flux, the term P(r, t) is the volumetric emissivity considered isotropic such that the factor 1/4π effectively extracts the photon flux along n ̂ , and the line-integral is evaluated from the origin of the LOS or spacecraft location (l = 0) to an appropriate boundary from where terrestrial soft X-ray emissions become negligible (l = L _max ). The term I _bkg denotes the astrophysical background along the LOS direction n ̂ that can be removed using sophisticated background models or alternative measurements, e.g., from Chandra X-ray Observatory (Weisskopf, 2003), ROSAT (Trümper, 1982), and XMM-Newton Jansen et al. (2001), among others.

The linearity between measured photon flux (I) and the volumetric emissivities (P), as presented in Eq. 3, enables the formulation of a discrete inverse problem (i.e., tomography) whose computational solution estimates the time-dependent 3-D soft X-ray emissivities at the magnetosheath from an ensemble of photon flux measurements of that region. In this study, we closely follow the steps presented by Cucho-Padin et al. (2022) in the context of exospheric tomography to formulate and solve this soft X-ray inverse problem.

First, we select the solution domain as the 3-D space that will be observed by the photometric sensor in a given period of time. Then, we discretize this solution domain into N non-overlapping 3-D voxels. The shape of the solution domain and voxels (e.g., cubic, spherical, cylindrical, or custom form) does not affect the mathematical formulation provided here; however, its selection should consider the desired spatial resolution of emissivities as well as the amount of computational resources used to solve the inverse problem. We illustrate the tomographic technique here using spherical shapes without loss of generality. Thus, we define the voxel size by its radial Δr, azimuthal Δϕ, and polar Δθ distances and assume that the soft X-ray emissivity is constant within the voxel volume.

Second, we define a [N × 1] column vector x containing all volumetric emissivities. We then create a [M × 1] column vector y of background-free measurements such that the mth element of y is y m = I m − I m bkg . Next, we calculate the intersection of each y measurement’s LOS with all voxels in the solution domain. This results in an ensemble of line sectors (dl) where each of them corresponds to a specific voxel. These dl values are used to form the [1 × N] row vector L[m], e.g., L[m] = [0, dl _n=2, …, 0, 0, dl _n=n , …, dl _n=N ], where a zero value indicates no intersection between measurement’s LOSs and the solution domain. This process is repeated with m LOSs to generate the [M × N] observation matrix L. For the sake of clarification, the L matrix can be generated a priori when (i) voxel sizes and (ii) the measurement’s LOS look angles ( n ̂ ) are known. The compact algebraic linear system that relates emissivities (x _n ) and measurements (y _m ) is given by y = L x + w , (4) where w denotes a [M × 1] measurement noise vector inherent to the optical acquisition. Eq. 4 is known as the forward emission model when the vector of emissivities, x, is provided. On the other hand, when values for photon flux measurements (y), their LOS’s look angles, and voxel sizes of the solution domain are provided, we have an inverse problem whose solution estimates values for emissivities. In the next subsection, we provide a technique to solve this inverse problem considering the time-dependent evolution of the magnetospheric region.

The rapid variability of IMF orientations may modify the global structure of the magnetosheath and magnetopause in timescales of tens of minutes. To investigate these dynamic processes, wide field-of-view soft X-ray technology is expected to have a high acquisition rate, e.g., one image every 5 min. Thus, the solution to the inverse problem defined in the previous section should be able to resolve volumetric emissivities from an ensemble of observations acquired during this short period of time.

To do so, we closely follow the statistical approach introduced by Norberg et al. (2023) in the context of ionospheric dynamic tomography. Thus, we first define a dynamic state-space framework that consists of two equations: y k = L k x k + w k (5) x k = F k − 1 x k − 1 + q k − 1 (6)

Here, Eq. 5 is the “measurement equation” and follows a similar structure as Eq. 4, the subscript k denotes the discrete time step, and the vector w _k is defined as a random vector with Gaussian distribution w k ∼ N ( 0 , R k ) (where N ( μ , Σ ) denotes a random vector that has a Gaussian probability distribution with mean vector μ and covariance matrix Σ). Furthermore, Eq. 6 is the “model evolution equation” where the [N × N] matrix F _k−1 models the temporal evolution of the states, typically derived from principle physics, and the [N × 1] vector q is the process noise with an assumed Gaussian distribution q k ∼ N ( 0 , Q k ) .

Based on the previous equations, we can define the probability distribution of the predicted state x ̂ at time k given a set of measurements y acquired at time k − 1 as p ( x k | y k − 1 ) = N ( x ̂ k , P ̂ k ) where x ̂ k = F k − 1 x ̄ k − 1 , (7) P ̂ k = F k − 1 P ̄ k − 1 F k − 1 T + Q k − 1 (8) and the superscripts (̄) and (̂) indicate the prior estimation and the future prediction of a given variable, respectively. The matrix P is known as the covariance matrix of states. The diagonal elements contain the variance of the soft X-ray emissivity in each voxel contained in vector x, and the off-diagonal elements indicate the covariance (or correlation) between a given pair of voxels.

Furthermore, if the probability distribution for the predicted state is (i) used as the prior distribution for the next time step k and (ii) updated with measurements also acquired at time k, we can obtain the posterior distribution defined as p ( x k | y k ) = N ( x ̄ k , P ̄ k ) where x ̄ k = P ̄ k L k T R k − 1 y k + P ̄ k − 1 x ̂ k P ̄ k = L k T R k − 1 L k + P ̂ k − 1 − 1 (9) Here, both x ̄ k and P ̄ k are presented in a space form and depend on their predicted versions. Also, the state-transition matrix F _k−1 is assumed to be the identity matrix I, which results in a purely random walk evolution model. For clarification, a random walk evolution model describes a simple process in which the next state only depends on the previous state and a fixed probabilistic variation which is provided in terms of variance. In our case, when F is equal to I, Eqs 7, 8 becomex _k = x _k−1 and P _k = P _k−1 + Q _k−1, respectively, which, in turn, shows that the variation from the previous state is given only by Q _k−1. Ideally, matrix F could incorporate magneto-hydrodynamics to predict the temporal evolution of the dayside magnetospheric region; however, the process of discretizing the MHD equations into a matrix would increase the complexity of the inverse problem. Previous studies demonstrated that a random walk evolution framework, along with continuous observations of the states, yields good results that are superior to static reconstructions (Zhang et al., 2005; Butala et al., 2010).

The iterative implementation of Eq. 9 provides time-dependent estimations of the state x ̄ k , and is known as Kalman filtering (KF). Typically, KF requires the use of the full covariance matrix P _k with N ² elements. The high memory allocation needed for this process fostered the development of the TomoScand Gaussian Markov Random Field (TS-GMRF) approach (Norberg et al., 2023; 2018). In TS-GMRF, Eqs 7, 8 are simplified using F _k−1 = I, and it is assumed that the predicted distribution (in Eq. 8) can be approximated with Q _k−1, i.e., P ̂ k ≈ Q k − 1 . Also, the TS-GMRF approach constructs the matrix Q ⁻¹, known as the “precision matrix,” using correlation lengths for the three dimensions (l _r , l _ϕ , l _θ ), a standard deviation function σ (r), and an ensemble of sparse differential matrices that impose smoothness on the estimated state vector x _k . Hence, when matrix Q ⁻¹ is provided, we can rewrite Eq. 9 as x ̄ k = L k T R k − 1 L k + Q k − 1 − 1 − 1 L k T R k − 1 y k + Q k − 1 − 1 x ̄ k − 1 (10) It is noteworthy that not only the precision matrix Q k − 1 , but also the observation matrix L _k is sparse due to the measurement geometry, such as the resulting system in Eq. 10 remains sparse and allows exploiting this property for efficient computation.

For the sake of clarification, the correlation length (l) is defined as the maximum distance between two points in the solution domain (e.g., x r 1 , x r 2 at locations r ₁ and r ₂) where the covariance of their physical values drops to 10% of the variance, i.e., l = argmin l Cov ( x r 1 , x r 2 ) − 10 % × Var ( x r 1 ) . The standard deviation function σ (r) used in constructing Q ⁻¹ is a [N × 1] vector that provides a standard deviation value to each voxel in the solution domain. Its mathematical definition is specific for the 3-D structure to be reconstructed and will be presented in Section 3. Also, the reader is referred to (Norberg et al., 2018; Cucho-Padin et al., 2022) for further details on implementing matrix Q ⁻¹.

To assess the confidence of our time-dependent, 3-D tomographic retrievals of soft X-ray emissivities, we use the residual error and the Structural Similarity index (SSIM). The residual error is defined by R e s . E r r o r = 100 % × ‖ y k − L k x ̄ k ‖ 2 ‖ y k ‖ 2 (11) and denotes the agreement between the 3-D estimated emissivities x ̄ k and the input radiance data y _k .

The SSIM index is defined as S S I M = 2 μ x g μ x ̄ + C 1 2 σ x g x ̄ + C 2 μ x g 2 + μ x ̄ 2 + C 1 σ x g 2 + σ x ̄ 2 + C 2 (12) where the terms μ x g and σ x g 2 denote the mean and variance of the ground truth vector x _g , the terms μ x ̄ and σ x ̄ 2 indicate the mean and variance of estimated emissivity vector x ̄ , and the term σ x g x ̄ represents the covariance between x _g and x ̄ . Also, the constant values C 1 = ( K 1 D ) 2 and C 2 = ( K 2 D ) 2 where D is the dynamic range of x, i. e., D = max(x _g ) - min(x _g ), and we use K ₁ = 0.01 and K ₂ = 0.03 as recommended in (Wang et al., 2023; Cucho-Padin et al., 2024). The index SSIM has values from 0 to 1, where SSIM = 1 indicates that the two input vectors are identical.

We follow the next steps to generate the ground truth vector x _g of soft X-ray emissivities. First, we use the magnetohydrodynamic (MHD) Open Geospace General Circulation Model (OpenGGCM) simulation reported in (Connor et al., 2021) to obtain the global spatial distributions of protons in the inner and outer magnetospheric region during controlled solar wind conditions. Using linear interpolation, we extract the OpenGGCM plasma data in the region of X _GSE , Y _GSE , Z _GSE ∼ [−20, 20] R _E with a spatial resolution of 0.1 RE. This spatial configuration defines a rectangular grid that contains N ^R = 64, 481, 201 voxels. In this proof-of-concept study, we aim to reconstruct the magnetosheath under varying solar wind magnetic conditions; thus, we conducted a 1-h MHD simulation with a 1-min resolution wherein the IMF suddenly varies from [B _x , B _y , B _z ] = [0, 0, 5] nT to [B _x , B _y , B _z ] = [0, 0, −5] nT at the 10th minute of the simulation. Values of solar wind density and velocity are assumed to be constants with values n _sw = 10 [cm⁻³], v _sw = 400 [km sec⁻¹].

Second, we identify and extract the proton population from the inner magnetospheric region simulated by the MHD model. Since this region does not contain sufficient heavy solar wind ions, it will not yield significant soft X-ray emission. For this task, we utilized the approach described by Samsonov et al. (2022), which determines the location of the inner magnetosphere through the following equations: p < p m s p + Δ p V x > V x s w × k v , (13) where p(msp) is the thermal pressure of the magnetospheric region, V _x (sw) is the solar wind along the x-axis, and the variables Δp and k _v are manually selected to adjust the comparison. For each 1-min MHD model, the thermal pressure of the magnetosphere is calculated, spanning values between 0.05 and 0.1 [nPa]. Also, V _x (sw) is set to −400 [km sec⁻¹], and we assume ΔP = 0.2 nPa and k _v = 0.15 as recommended in (Samsonov et al., 2022). If both conditions in Eq. 13 are matched, that location is considered part of the inner magnetosphere and should be removed from the analysis.

Third, we use a spherically symmetric neutral exosphere whose density distribution (n _H ) is given by (Connor and Carter, 2019): n H r = n 0 10 R E r 3 atoms cm − 3 (14) where n ₀ is the hydrogen (H) density at 10 R _E subsolar point and equal to 25 [atoms   cm⁻³], and the term r denotes the planetocentric distance in units of R _E . Also, we adopt an efficiency factor α = 1 × 10⁻¹⁵ [eV   cm²]. Finally, the global volumetric emissivity P can be directly calculated using Equation 2.

Figure 1 shows the temporal evolution of magnetosheath soft X-ray emissivities as a response to varying solar wind conditions, which will be denoted as x g k and will serve as the ground truth emissivities to be reconstructed. The top panel of Figure 1 shows the Bz component of the IMF. The bottom panel is subdivided into three rows displaying the ecliptic plane (XY _GSE ) and a 3-D visualization of the emissivities at the (A) 5^th, (B) the 20^th and (C) the 55^th minute after the start of the simulation.

Model-free, time-dependent, and 3-D tomographic reconstructions of soft X-ray emissivities in the magnetospheric region require simultaneous observations from distinct vantage points that can provide sufficient spatial information about the target. In this context, we evaluate tomographic reconstructions using a two-satellite configuration based on the results reported by Jorgensen et al. (2022), which indicate that using two soft X-ray imagers significantly reduces the estimation error obtained with a single instrument. Also, in a realistic scenario, using two satellites would reduce the cost of implementing additional space-based platforms and imaging sensors.

Our design of a spacecraft orbit closely follows that one proposed by Sibeck et al. (2023) for the Solar-Terrestrial Observer for the Response of the Magnetosphere (STORM) mission concept, whose objective is to continuously track and quantify the flow of solar wind energy through the magnetosphere using a set of optical instruments that measure energetic neutral atom (ENA) from the ring current, far ultraviolet (FUV) emission from the auroral region, and soft X-ray emission mainly from the dayside magnetosheath. Our study focused on the analysis of the soft X-ray measurements. The STORM mission concept comprises a single spacecraft in a circular orbit with a radius of 30 R _E and an inclination of 90° with respect to the ecliptic plane. The orbit’s orientation is nearly fixed in the inertial frame, allowing the Sun-Earth direction to rotate 360° relative to the orbit over 6 months. The spacecraft is three-axis controlled, with one body axis always to the nadir and the other two axes controlled to keep the sun on one side of the vehicle. Further, the soft X-ray imager onboard the STORM mission has a square field-of-view (FOV) of 23 [deg] and a constant boresight direction canted 18.5 [deg] towards the sun direction in the spacecraft body frame. Aside from periods of low orbital beta angle (Sun near the plane of the orbit), this design allows nearly continuous imaging of the subsolar point at 10 R _E (nominal location of the magnetopause).

The objective of this manuscript is to conduct a comparison study of tomographic reconstructions of the dayside magnetosheath from soft X-ray observations acquired by two satellites in STORM-like orbits. For this, we simulated seven (7) spacecraft with different phase angles but similar longitudinal locations that will be arranged in pairs for the tomographic process. Figure 2 shows the distributions of these satellites named spacecraft 0 (SC0) to spacecraft 6 (SC6) along the circular orbit (gray curve). The phase angles between satellites are fixed, and their values (in units of degrees) with respect to SC0 are SC1: 60 [deg], SC2: 75 [deg], SC3: 90 [deg], SC4: 105 [deg], SC5: 120 [deg], and SC6: 180 [deg].

In this work, we assert that tomographic reconstruction can be performed if and only if the two imaging instruments observe a given target from different vantage points simultaneously. Although this condition is ensured by our design of the imager’s boresight direction that allows nearly continuous observation of the subsolar point at 10 R _E , there are periods of time in which one or two imagers point sunward. To avoid contamination from direct sunlight, we define a Sun-Earth avoidance angle equal to 45 [deg]. In other words, if the angle created by the Sun-Earth line (X _GSE ) and the LOS direction ( n ̂ ) of a pixel in a soft X-ray imager is smaller than 45 [deg], we consider that tomography is not possible for this satellite’s 3-D position.

Based on the previous requirements, Figure 3 shows locations (in GSE spherical coordinates) where tomography is feasible using a pair of spacecraft for a period of 6 months. In each panel, the gray dots indicate the 1-h interval position of the SC0 (at r ≈ 30R _E ) during a complete translation around Earth. The blue dots show locations where SC0 and SCX (with X ∈ [1, 2, 3, 4, 5, 6]) simultaneously acquire soft X-ray radiance data free from sunlight contamination, i.e., tomography is possible. The number of locations to perform tomography significantly reduces as the phase angle between SC0 and SCX reaches 180 [deg], especially near the low-latitude region, as shown in the bottom right panel (SC0 and SC6). In order to conduct a comparison study of tomographic estimations using these six two-satellite configurations, we selected two vantage points for SC0 where tomography is feasible with all the remaining spacecraft. They are named P1 and P2 and are depicted as red stars in all six panels. Their specific positions in GSE Cartesian coordinates are P1 = [4.846, −14.143, −26.085] R _E and P2 = [3.941, 14.679, 25.827] R _E .

Figure 4 shows additional analysis for tomographic feasibility using these two-satellite configurations. Each panel shows a histogram of the number of hours when tomography is possible (using SC0 and SCX) with respect to the ecliptic longitude. The histogram dataset has been generated from data shown in Figure 3 that were accumulated along a 10-deg bin in the longitudinal dimension. The total number of hours when tomography can be performed is included on the top left of each panel. In a period of 6 months, tomography can be conducted between 1736 h (using SC0 and SC6) and 2,503 h (using SC0 and SC1). In addition, the number of hours when tomography is feasible remains nearly consistent (at ∼ 110 h) for all six cases when SC0 is located near the terminators (dawn and dusk regions or ecliptic longitudes around ±90 [deg]).

Furthermore, the total number of hours in a 6-month period is ∼4,320, which indicates that each spacecraft configuration can perform tomography between 40% and 57% of this period. Nevertheless, due to the geometry constraints, continuous acquisition along a longitudinal bin for ∼120 h and covering all latitudes is possible only near the terminators, as shown in Figure 3.

In this subsection, we define the experiments needed to evaluate the efficiency of the tomographic results using a two-satellite configuration. Also, we describe the generation of synthetic soft X-ray images and the details of implementing the dynamic tomographic approach.