Remote Sensing Estimates of CME Density in the Ecliptic Using the STEREO Heliospheric Imagers

We present a method to estimate electron densities in the ecliptic plane using the STEREO Heliospheric Imagers (HIs). The nature of Thomson scattering of photospheric light by solar wind electrons is such that visible‐light observations by HI provide a means to infer information about plasma density from afar. This is achieved using discrete tomography, whereby line‐of‐sight integrals from HI are used to estimate electron density in the heliosphere over a predefined grid. The technique is applied to the Earth‐impacting coronal mass ejection (CME) launched on 12 December 2008. The two vantage points afforded by STEREO are insufficient to reproduce the density structure of the CME in detail; however, the technique is successful in locating the presence of a density enhancement associated with the CME. When applied to consecutive images, we are able to use the technique as a means to track the CME through interplanetary space. From these observations we make estimates of the CME radial velocity and density profiles, the results of which are consistent with measurements made in situ by Wind at L1. This method is presented as a new approach to determine CME propagation from HI observations, one that avoids many of the assumptions of CME morphology and dynamics that are often applied when tracking CMEs in such data. We also expect that this technique will prove useful, and a test of the long‐standing heliospheric reconstruction technique employed by using single images over time (Jackson et al., 2006, https://doi.org/10.1029/2004JA010942) when views from many wide‐angle imagers become available.


Introduction
Coronal mass ejections (CMEs) are a means by which the Sun releases large quantities of its energy into interplanetary space. Such occurrences are frequent (Gopalswamy, 2010, show respective rates of <0.5 and >6 per day at solar minimum and maximum in Solar Cycle 23), and, should they hit Earth, can be the cause of severe geomagnetic disturbances. Knowledge of the state of the Earth-directed solar wind, particularly CMEs, is therefore essential to anticipating such occurrences. While solar wind monitors, such as ACE and Wind, orbiting the stable Sun-Earth L1 point provide us with accurate in situ measurements of solar wind parameters, their location means that they are only able to give about 1-hr advanced warning. Spacecraft positioned closer to the Sun, by making use of an electronic solar sail, have been proposed (Janhunen, 2010); however, these would only improve forecast times by about twofold. An alternative approach is to attempt to infer solar wind and CME parameters from afar using remote sensing instrumentation. This is often achieved using coronal observations, for example, the Large Angle Spectroscopic Coronagraph (LASCO, Brueckner et al., 1995) on board SOHO or the STEREO Coronagraphs (COR, Howard et al., 2008). These may then used to drive models that simulate plasma transport in the heliosphere (e.g., ENLIL, Odstrcil et al., 2005, or EUHFORIA, Pomoell & Poedts, 2018, which are in turn used to predict the future state of the solar wind at interplanetary distances. The development of wide-angle imagers, such as the Solar Mass Ejection Imager (SMEI, Eyles et al., 2003) on board Coriolis and the STEREO Heliospheric Imagers (HIs, Eyles et al., 2009) bridge this observational gap between the Sun and ∼1 AU and allow CMEs to be tracked continuously in this range. Each STEREO spacecraft hosts two HI cameras, HI-1 and HI-2, with respective fields of view (FOVs) of 20 • and 70 • and image cadences of 40 and 120 min. The former is off pointed from the Sun by 14 • and the latter by 53.7 • , and both are centered in the ecliptic. This configuration provides ecliptic elongation ranges of 4-24 • and 18.7-87.7 • , which is far beyond the range observed by coronagraphs. Due to difficulties in subtracting the background F-corona and residual stray light from HI images, we do not make use of HI-2 data in this study. Each HI camera uses a 2,048 × 2,048 pixel CCD array; however, these are binned 2 × 2 before transmission to Earth to produce 1,024 × 1,024 pixel science images. In this paper we are concerned with interpreting the science images and as such we present our analysis in terms of image pixels, which refer to the latter.
There are two principle means by which the three-dimensional positions of CMEs are located via visible-light remote sensing: stereoscopy and tomography. The former is a means of triangulation, whereby an observed front is tracked through the instrument FOV, which represents the leading edge, or leading surface, of the CME or observed feature. Viewpoints from multiple spacecraft then allow the position of the feature to be inferred. The latter, and potentially more informative approach, is the tomographic inversion of visible-light images, which combines projections from multiple, spatially separated vantage points of the optically thin plasma to constrain its three-dimensional density structure (e.g., Davila, 1994;Jackson & Froehling, 1995). From the perspective of space weather forecasting, speed, density, and magnetic field are the most important properties for estimating the geo-effectiveness of the solar wind and CMEs and so accurate density measurements from remote sensing would provide significant improvements to our forecasting capabilities. The triangulation approach may in fact be practiced with a single observing spacecraft using assumptions about CME morphology and propagation (e.g., Davies et al., 2012;Kahler & Webb, 2007;Lugaz et al., 2009;Rouillard et al., 2008;Rollett et al., 2016). Due to the nature of Thomson scattering (Tappin & Howard, 2009a), CMEs observed in coronagraphs are often assumed to propagate in the sky plane. In wide-angle imagers it is instead often assumed that CMEs have reached a constant speed, whereby their observed angular motion may be used to infer their speed and direction (e.g., Davies et al., 2009;Sheeley et al., 1999). Gopalswamy et al. (2009) acknowledge that the sky plane assumption produces inaccurate speed measurements and show also that some CMEs are still accelerating well into distance range covered by HI-1, thus suggesting the constant speed assumption may also be inaccurate in some cases. Stereoscopic approaches combine two or more viewpoints to track CMEs (e.g., Byrne et al., 2010;Davies et al., 2013) to remove the need for some of these assumptions and better constrain CME position. Both the stereoscopic and tomographic approaches allow the evolution of CME kinematic properties to be tracked by observing successive sets of images. The stereoscopic approach may be practiced with just two viewpoints to allow triangulation of the CME leading edge; however, its weakness is that it relies upon assumptions of CME morphology and therefore can reveal no information of the CME structure. The strength of the tomographic approach is that it need not rely on this assumption and, instead, produces a plasma density distribution over the observed region from which the presence, or absence, of a CME may be determined. However, the method requires a large number of viewpoints (a minimum of four according to Davila, 1994) to be capable of resolving detailed and small-scale density structures.
The tomographic technique is well established and has been performed since the first coronagraphs were launched (e.g., Kastner et al., 1977, with OSO 7 andAltschuler, 1979, with Skylab). These techniques were first applied to background coronal structures, which were assumed to be quasi-static, and so multiple observation angles are provided by solar rotation. This technique is known as Solar rotational tomography (SRT) and is the favored approach for determining coronal density structures when Thomson scattering polarization brightness is available. Further SRT studies were performed by Zidowitz et al. (1996Zidowitz et al. ( , 1997Zidowitz et al. ( , 1999 using both the Mauna Loa MK3 Coronameter and C1, the inner LASCO camera on board SOHO, which have FOVs extending to 2.44 and 3 R ⊙ , respectively. Frazin (2000) and Frazin and Janzen (2002) applied the technique to LASCO C2, which extends the observational range to 6 R ⊙ . Moran and Davila (2004)  ratio of polarized to unpolarized light observed from the single viewpoint of LASCO to achieve the polarimetric reconstruction of the three-dimensional structure of two CMEs. The launch of STEREO in 2006 allowed the technique to be performed from multiple vantage points throughout the heliosphere. Kramar et al. (2009) applied SRT to COR1 observations to reconstruct coronal densities between 1.5 and 4 R ⊙ over two 14-day periods. de Patoul et al. (2013) use a combination of stereoscopic and SRT methods to reconstruct density structures in the region between 1.01 and 1.39 R ⊙ from the three points of view of SOHO/EIT and STEREO/EUVI. While SRT has been successfully applied to quasi-static structures, it is not applicable to CMEs, which evolve on much smaller timescales than the solar rotation frequency. Antunes et al. (2009) combine both inverse modeling and forward modeling using the graduated cylindrical shell model (Thernisien et al., 2006 to constrain the structure of a CME from the two points of view of COR. Frazin et al. (2009) demonstrated that combined C2, C3, and COR1 data, with the addition of forward modeling, could be used to reconstruct the density structure of a CME, and Frazin (2012) goes further to include MHD simulations in the reconstruction.
The first tomographic approach to attempt CME reconstruction, and the first making use of multiple vantage points, was by Jackson and Froehling (1995). They combine observations from the Zodiacal Light Photometers on board the HELIOS 2 spacecraft and the Solwind coronagraph to reconstruct a CME that occurred on 7 May 1979. A velocity profile is fitted to the observations as a means to bridge the gap between the FOVs of each instrument. The result was a density distribution with a large three-dimensional extent comprising two lobes of dense material. This solution was found to be consistent with the observed brightness but did not resemble the expected shell or loop-like CME structure. A method was developed by Asai et al. (1998), Jackson et al. (1997Jackson et al. ( , 1998, and Kojima et al. (1998), which uses the SRT technique and radio interplanetary scintillation (IPS) as a means to successfully constrain a density model of the heliosphere. The Jackson et al. (1998) model also employs physical principles in order to help provide a solution, assuming radial solar wind outflow and enforcing conservation of mass and mass flux. Density and velocity are projected outward from a solar wind inner boundary, or source surface, such that they can be compared to IPS line-of-sight (LOS) measurements, and observational gaps are filled in by means of a Gaussian filter. This model was further developed to permit time-dependent three-dimensional (heliospheric latitude, longitude, and time) tomography (Jackson & Hick, 2000). Further development allowed the use of Helios photometer data as an input and was able to reproduce the three-dimensional density distribution associated with both corotating and fast-moving structures (Jackson & Hick, 2002, which are found to be in agreement with in situ observations of the same structures from five different spacecraft by Burlaga et al. (1980). The time-dependent IPS model was later developed in order to incorporate Thomson scattering measurements from SMEI . Radial and tangential magnetic fields from the Current Sheet Source Surface magnetic field model (Zhao & Hoeksema, 1995) have also been integrated since 2005. The Current Sheet Source Surface model provides the magnetic field components on the source surface, which is propagated outward by the flow lines of velocity from the tomography model, based on the assumption of flux freezing (Dunn et al., 2005). The IPS model continues to be used to make real time predictions of global solar wind parameters using near-real-time data from the Solar-Terrestrial Environment Laboratory, Japan. Near-Sun results from this modeling effort have also been used to drive the inner boundary of MHD models such as ENLIL to provide real-time forecasts Yu et al., 2015).
In the majority of multispacecraft studies in STEREO era, the tomographic technique is typically restricted to coronagraph or other near-Sun observations. This is because the contribution from dust-scattered light (the F-corona) is far less dominant near the Sun when compared to the light scattered by electrons (the K-corona). However, at elongation angles covered by HI or SMEI, resolving any CMEs becomes an effort in accurate background subtraction of the F-corona. The typical K-corona, including CMEs, at the edge of the HI-1 FOV (24 • elongation in the ecliptic) may have an observable radiance of a few 10 −15 B ⊙ , while the F-corona is typically 10 2 times greater. In the FOV of HI-2, this ratio is closer to 10 3 . Therefore, techniques used to simply track CMEs require a precision of only ∼ 10 −2 when subtracting the F-corona. Photometric analysis, such as tomography, requires a far higher accuracy in order to achieve accurate results. DeForest et al. (2011) are able to calibrate HI-2 data to 10 −17 B ⊙ using a combination of techniques, and this processing has been used to study small-scale density structures (e.g., Howard et al., 2013Howard et al., , 2012Howard & DeForest, 2014. The state-of-the-art of these techniques is summarized in DeForest et al. (2016); however, the techniques are only currently applied to HI data from STEREO-A, resulting from difficulties with the STEREO-B instruments. Unfortunately, we require data from HI-1 on both spacecraft for this study,

10.1029/2019JA027175
and as such, we must rely on more primitive image processing techniques. However, as we show here, a solution is possible using both HI-1A and HI-1B with relatively simple corrections to the data. Various studies have used polarized light measurements from coronagraphs to determine CME structure (e.g., Crifo et al., 1983;Moran et al., 2010;Poland & Munro, 1976), in addition to polarimetric reconstruction discussed above. DeForest et al. (2016) have argued the case for it a requirement for the next generation of HIs, and it will soon be realized by the upcoming Polarimeter to Unify the Corona and Heliosphere (PUNCH) mission. Section 2 presents a summary of Thomson scattering in the context of HI observations and a mathematical formulation of this as an inverse problem that is to be solved for electron densities. Section 3 then presents a description of the image processing techniques we apply to STEREO/HI data and how we translate these data into electron density estimates.

Thomson Scattering in the Heliosphere
A thorough summary of Thomson scattering in the heliosphere is presented by Minnaert (1930), in which the effects of polarization and solar limb darkening are also detailed. These are revisited by Tappin and Howard (2009a) who present a clearer pictorial description of the scattering process with emphasis on how the observed intensity is dependent on the Sun-electron-observer scattering angle, . The authors focus particularly on how the theory is applied to observations at large elongation angles, for the purpose of interpreting HI data. Accounting for integration over the photosphere, over the point spread function of each pixel and along the pixel LOS, Tappin and Howard (2009a) arrive at an equation that is equivalent to where e is the differential cross section for perpendicular scattering, R ⊙ is the solar radius, u is the solar limb darkening coefficient, and N e is electron density. The equation applies to an observer at a distance r 0 from the Sun and an LOS elongation of , where is the electron-Sun-observer angle. Equation (1) differs from that of Tappin and Howard (2009a) in that the van de Hulst coefficients have been replaced by their approximate values when scattering occurs at a sufficient distance from the Sun.

The Inverse Problem
The key to tomographic inversion of HI images lies in formulating, and solving, an inverse equation of the form The radiance values, in each image pixel, are contained in the array y, and x is an array of densities over a predefined grid, each element, n, of which represents the density in a single grid cell. The matrix H is simply a physical operator that relates the discrete values of x and y. Each HI-1 image consists of 1, 024×1, 024 pixels. The work in this paper is concerned only with densities in the ecliptic plane, and therefore, our observational data are the row of pixels, in each instrument, whose pointing direction is closest to the ecliptic. The array y therefore possesses M = 2048 elements, each representing the image pixel nearest to the ecliptic plane; 1,024 pixels from each spacecraft. The number of elements, N, of x is dependent on our choice of grid resolution and the size of the region where the HI FOVs overlap, which varies throughout the STEREO mission. H is then an M × N matrix, each element of which has a value of or 0, if LOS m does not intersect grid cell n. Tappin and Howard (2009b) present an algorithm for determining the elements of H over a spherical, heliocentric grid, which is as follows: 1. Starting at the observer, compute the distance to the next grid cell boundary in each dimension. 2. The distance through the grid cell is the shortest of these. 3. Compute the integrand at the midpoint of the path through the cell. 4. Add the contribution to the total. 5. Subtract the distance used from the distance in the other two dimensions and recompute the distance to the next grid cell in the remaining dimension. 6. Repeat Step 2 until z ≥ 2(1 + sin( 2 )).
Here we choose a spherical, heliocentric grid with dimensions 1 160 AU × 1.5 • × 1.5 • in radial distance, azimuthal angle, and polar angle. Although the polar angle is not required for our two-dimensional ecliptic reconstruction, we include it in the calculations as a sanity check. For a given pair of HI-1A and HI-1B images, the position of each spacecraft and the pointing direction of each image pixel is known. It is therefore possible to determine the matrix, H, from equation (3). Combined with the observations, y, equation (2) presents an inverse problem, which may be solved for x yielding the electron density distribution in the ecliptic plane.
Due to the nature of the problem, the matrix, H, is found to contain a large number of zeros. It therefore cannot be dealt with analytically, and, instead, an iterative algorithm, the Conjugate Gradient (CG) method, is employed in order to determine x. The CG method in fact refers to a group of methods, each of which may be applicable depending on the nature of the matrix in the inverse equation. These methods are iterative algorithms for solving large systems of linear equations, such as equation 2. Because the method is iterative, it may be stopped after a finite number of steps, after a sufficient level of accuracy has been reached. Several approaches to solving linear problems using iterative methods, including the CG algorithm and its variants, are covered by Barrett et al. (1994). The BiConjugate Gradient Stabilised (BiCGSTAB) algorithm (van der Vorst, 1992) is used here because it may be applied to nonsymmetric matrices, such as H, and is known to converge quicker than other methods. By inputting an initial guess, x (0) , of the density, the algorithm will determine an updated value of x (i) after each iteration, which serves to reduce the residual, y − Hx (i−1) .

In Situ Measurements
The aforementioned tomographic technique is applied to the 12 December 2008 CME, studied comprehensively by such authors as Byrne et al. (2010), Davis et al. (2009), andLiu et al. (2010). The reasons for choosing this CME are discussed later. This near Earth-directed CME impacted the spacecraft at L1, with which we can measure the density in situ. Because the HI-1 FOVs cover less than 0.5 AU on the Sun-Earth line (see Figure 3) at the time when this CME occurred, we have no direct means directly compare the HI and in situ observations to assess the fidelity of our density reconstruction. Therefore, we make assumptions about how the CME density and velocity evolve beyond this distance, as it travels toward Earth, such that we may compare our result to the in situ signature measured by Wind. Richardson et al. (2005) take in situ measurements from Helios 1, 2, ACE, Wind, and Ulysses and determine the average evolution of CME parameters as a function of heliocentric distance. They show density to follow a power law equaling N e (r) = 6.2r −2.3 , where the density decrease with increasing r is greater than the ambient solar wind due to CME expansion. CME velocities are shown to be approximately constant at around 450 km s −1 , which is consistent with previous studies. While a great deal of CME acceleration takes place in the low corona, CMEs have often established an approximately constant speed by the time they enter the HI-1 FOV. Gopalswamy et al. (2009), for example, show that 87% of CMEs observed by the LASCO (Brueckner et al., 1995) had stopped accelerating within the 32 R ⊙ FOV. We therefore make the assumptions that the radial velocity of our Earth-impacting CME is constant and that its density follows a power law. This allows us to extrapolate our HI-1 observations out to L1 and compare our results with those measured in situ by the Wind spacecraft that occupies this position. Although part of the motivation for the inversion of HI data was to avoid these kind of assumptions, we wish to emphasize that they are only used as a means to validate the results against in situ data and are not required to perform the inversion itself.

Method
Before tomography may be performed on HI data, there are two issues that must be addressed. First, the observed radiances in HI include not only Thomson scattered photospheric light but also light scatted from the dust of the F-corona and the contribution from background stars, planets, other solar system bodies and any residual stray light. These must therefore be separated out to provide a data array y containing only the contribution from free electrons. The second issue is that we do not expect equation (2) to possess a unique solution, and therefore, some regularization must be applied to arrive at a physically consistent result. The first of these issues amounts to corrections applied to the data, y, in equation (2) and the second to corrections applied to the equation itself.

Image Processing and Background Subtraction
The K-coronal brightness is typically a few percent of the total observed brightness in HI at elongation angles of a few degrees, and this reduces when moving away from the Sun ). This is the reason The left-hand images show L1 data, the middle images show L2 data where the F-corona and stray light have been removed, and the right-hand plots have had a robust mean applied to remove stars. The right-hand plots therefore represent an estimate of the radiance due to the K-corona alone. The radiance of the L1 data has been scaled by 0.01 relative to the other two, and all six plots are shown on a log scale due to the large range in observed radiances. that tomographic reconstruction of coronal density is more easily applied to coronagraph observations close the Sun. Using the methods described in this section, however, it is possible to separate the excess density associated with transient CMEs from the quasi-static F-corona and background K-corona, and from other sources in order perform tomography using the HI-1 instruments on STEREO.
The Level-0, or L0, data from each HI camera undergo several stages of processing after being received by ground stations in order to arrive at a value for the K-corona. These are described in detail by Eyles et al. (2009) and  and are summarized here. The results of the main stages of processing are shown in Figure 1. First, the data are corrected for shutterless readout ) and for a large-scale flat-field (Bewsher et al., 2010;Tappin et al., 2015). The flat-field correction accounts for uneven responses of pixels across the CCD array by applying a gain map, derived from analyzing the starfield . Moreover, a pointing calibration is applied to the data, by comparing the locations of stars in the HI image with their positions predicted using a star catalogue (Brown et al., 2009). These processes are standard calibration techniques performed in order to produce the L1, data from the L0 data. An additional stage of processing is to remove the contribution from the F-corona and residual stray light to produce L2 data (Figure 1), which is also a standard calibration performed by the PI team (STFC/RAL Space). These methods are based on the assumption that the F-corona does not vary on time scales as small as the HI image cadences. A sequence of images is taken over a period of m days. In nominal mode, with no missing images, this consists of 36 images per day for HI-1. The F-corona is then approximated, for each individual pixel, as the average of the lowest nth percentile of all values in that image pixel. As part of the official HI-1 pipeline, RAL Space provides data with the mean of the lowest quartile over 1 and 11 days, that is, n = 25 and m = 1 or 11. The 1-day lower quartile mean has been subtracted to produce the L2 data in Figure 1. These background values are then subtracted from an individual image to reveal structures in the K-corona. Finally, the starfield must also be removed. This is achieved, here, by simply replacing each pixel by the Figure 2. The mean ecliptic radiance observed by STEREO HI-1 over a 2-week (28 November to 11 December 2008) period with varying levels of processing applied. The top panels show the L1 data in blue with one standard deviation indicated by dashed lines. The black line is the theoretical observed radiance assuming a density profile of 5r −2 electrons cc −1 . The bottom panel is plotted with a smaller range on the axis and shows the same theoretical K-corona in black with three methods of background subtraction applied: The mean of the lowest quartile over 1 day (red), the mean of the lowest quartile over 11 days (green), and the minimum value over 1 day (purple). In each case the dashed lines indicate one standard deviation from the mean. robust mean of its neighborhood; from an array of p × p pixels centered on each image pixel the median and median absolute deviation are calculated. Any pixel values exceeding q deviations from the median are excluded and the mean value of the remaining pixels is returned. This procedure is used to exclude outliers that are presumed to result from the presence of background stars. In this paper, we use values of p = 21 and q = 2 to produce the final (right hand) pair of images in Figure 1, which we will refer to colloquially as L3 data.
The bottom panels in Figure 2 show three examples of L3 data that have had different background subtraction methods applied to them in the second stage of image processing. In each case these L3 data have been averaged in each pixel over a two week period (28 November to 11 December 2008) to show the mean radiance and variability in each pixel as a function of elongation. In each case, the background subtraction of methods are indicated in the figure: The mean of the lowest quartile over 1 day is subtracted to produce the red data, the mean of the lowest quartile over 11 days is subtracted to produce the green data, and the minimum value over 1 day is subtracted for the purple data. These subtraction methods are all based on the assumption that transient solar wind structures, such as CMEs, are superimposed on a static background comprised of the F-corona. However, in practice, these methods also remove some of the quasi-static K-corona and so the L3 values, plotted in the bottom panels Figure 2, should be considered as excess density observed above the background level. Subtracting a background over a shorter period of 1 day is performed to remove the streamer belt contribution to the signal and reveal fast-moving transient structures such as CMEs. The background calculated over a longer period of 11 days is intended to retain more of the streamer belt structure. Figure 2 shows all the L3 methods to be far lower than the theoretical background K-corona (calculated using 5r −2 electrons cc −1 Cox, 2000) at low elongation angles. The 1-day minimum and 1-day lower quartile mean subtraction methods produce data that are an order of magnitude fainter than the expected K-coronal background. The 11-day lower quartile mean subtraction method removes less of the K-corona at small elongation angles than the other two methods but is still much lower than the expected value. This demonstrates the problem that in subtracting the F-corona, we are also subtracting a large BARNES 7 of 19 10.1029/2019JA027175 contribution from the K-corona as well. In the case of the 1-day lower quartile mean subtraction the relative intensity observed by HI-1A to HI-1B is most consistent with the theoretical K-corona. The reason that the data have higher values at a given elongation for HI-1A than HI-1B is because the spacecraft are at respective solar distances of 0.97 and 1.04 AU. Equations (1) and (3) both show that, for a given LOS elongation, observed brightness is inversely proportional to the Sun-spacecraft distance, r 0 . For the 1-day minimum subtraction method the HI-1A observations are much greater compared to the HI-1B observations and therefore this method produces inconsistent results between the two spacecraft. For the 11-day lower quartile subtraction we see that the data from HI-1B are brighter than from HI-1A at a given elongation. This is likely to be due to inconsistencies in the background subtraction over the longer 11-day period. However, in this paper we are concerned with the fast-moving structure revealed by the 1-day subtraction and as such the data produced from the 11-day subtraction are not used in our analysis. Because of the inconsistencies between data from each spacecraft using the 1-day minimum method, we instead opt for the 1-day lower quartile subtraction method in order to carry out our following analysis. The top panel of Figure 2 shows that the variability of the L1 observations becomes increasingly large at greater elongation angles, and so the assumption that this background is static breaks down for pixels in the outer edge of the HI-1 FOV. As such, we restrict our further analysis to pixels with an elongation within 22 • of the Sun.
After the launch of STEREO, it was found that the pointing of the HI-1 instrument on STEREO-B made many small movements. This was originally thought to be due to the impact of dust particles, due to the fact that the cameras on STEREO-B are facing its direction of motion (Davis et al., 2012). More recently,  have shown that this effect is in fact the result of vibrations caused when the attitude-control reaction wheels on STEREO-B are spinning faster than 350 revolutions per second. This is thought to be because the HI-1B instrument became loose during launch, and the same effect is not seen in HI-1A. The scale of these movements is usually very small and translates to ∼ 1 3 of an image pixel. It is, however, enough to severely affect the standard method of background subtraction and so the motion must be accounted for. The pointing of the HI cameras is computed using the known starfield, and if the offset in successive images changes by more than one pixel, then this motion is accounted for before performing the background subtraction.
The presence of planets, comets, and other solar system bodies in the FOV during some periods of HI data contributes to observed radiances via reflected sunlight. There are further potential contributions that occur due to dust impact signatures (Davis et al., 2012) and internal reflections (Halain et al., 2011). In addition there is an approximately biannual passage of the galactic plane through the FOV, which is considerably denser than the typical background starfield and causes the robust mean correction to fail. These contributions make it very difficult to produce L3 data and so render many periods of HI data unusable for the purpose of tomography.

Regularization
While the previous section dealt with corrections that must be made to the data y, here we deal with corrections to the inverse equation, y = Hx, as a whole in order to ensure a physically realistic result. Specifically, we require a spatially smooth density distribution. In order to reduce the contribution of observed noise to the resulting density reconstruction, Frazin (2000) introduces a condition that ensures that the density variation is smooth. This is achieved by rewriting the inverse problem in terms of a new matrix, A, that satisfies the equation where b and A are given by The array b is created by extending the array y, with a number of additional elements equal to 0. The matrix, H, ensures the estimate,x, is consistent with the data, and the matrices, D 2 r and D , suppress large gradients in the resulting estimate. D 2 r and D are each a discrete approximation of the respective gradient operators, 2 ∕ r 2 and ∕ , and the regularization parameter, , may be tuned to control their relative contribution to the solution. The elements of the regularization matrices, D 2 r and D , are determined using a finite difference scheme. The value of used here is chosen such that the weighting of the regularization matrix is approximately 5% of the observation matrix. This is found to ensure smoothness in the resulting density, without preventing the solving algorithm from converging toward a value that satisfies the data.
Due to the large range in observed radiances (approximately 2 orders of magnitude in HI-1), equation (2) is multiplied by a weighting factor, w, to increase the contribution from LOSs that lie at large distances from the Sun. Kramar et al. (2009) Both element y m and row H m in equation (2) are multiplied by w m , which still delivers the solution, x, corresponding to the original system, y = Hx. Davila (1994) show with very basic simulated data that the optimum separation when performing tomography with two spacecraft is around 90 • . This seems evident if we consider that the technique is intended to reduce LOS ambiguities. Here, it is the HI-1 LOSs that we require to intersect at 90 • for optimal reconstruction and this occurred at around the end of 2008. However, this coincides with the depths of solar minimum and therefore few CMEs are observed. By the time we began to see a significant number of CMEs in HI (from the beginning of 2011), the spacecraft separation was close to 180 • and therefore offered little use in reducing LOS ambiguities in the density estimates. Because we are interested in observations in the region common to the field FOV of both spacecraft, this means that the CMEs to which we may apply the technique are restricted to Earth-directed events. The difficulty in background subtraction when the galactic plane is present in the FOV of either HI-1 camera also further reduces the periods to which we are able to test the technique. For these reasons, we have selected just a single event, the well-studied 12 December 2008 CME, to provide a proof of concept for the method.

Data Selection Criteria
A final consideration is the contribution to our observed radiance that comes from the regions that are not common to the FOV of both spacecraft. This problem is illustrated by Figure 3, where panel (a) shows the respective FOVs of HI-1A and HI-1B in the ecliptic plane. The solid line originating from STEREO-A is the LOS of the HI-1A boresight (14 • elongation). Panel (b) shows the contribution to the total radiance in this pixel, as a function of LOS distance if we assume a typical K-coronal density profile of 5r −2 cc −1 , which is taken from Cox (2000). The area bounded by vertical dashed lines is the region of overlap with the HI-1B FOV. However, as was discussed earlier, our background subtraction method removes the stationary K-corona and so we are observing excess density above this level. We therefore must assume that this excess density is completely contained within the region common to the FOV of both spacecraft, which is reasonable because we are observing an Earth-directed CME. Additionally, the positions of the STEREO spacecraft at this time also mean that this region is close to the Thomson surface, where observed scattering is maximized. While it is not ideal to include this assumption, it is unavoidable when using this technique. Figure 3a illustrates the reason that we choose to avoid the inclusion of LASCO data in our inversion. The largest of the LASCO coronagraphs, C3, has an elongation coverage of approximately 1-8 • , and there is therefore an overlap in the fields of view of the C3 instrument and both HI-1 cameras. The CME studied in this paper enters the C3 FOV at approximately 13:00 UT on 12 December 2008 and is observed as a faint halo that fades after several frames. Therefore, while there is some overlap in observations between HI-1 and C3, we do not therefore believe that observations from C3 would contribute any useful information to solving densities in the region covered by the HIs.

Tomographic Inversion
The CME we study was launched at around 03:00UT on 12 December 2008 (Byrne et al., 2010); it passed over Wind approximately 4 days later and was observed by both HI-1A and HI-1B (e.g., Davis et al., 2009;Liu et al., 2010) and also the HI-2 cameras (e.g., DeForest et al., 2011). One hundred eight pairs of HI-1 images from both spacecraft are selected covering the 72-hr period beginning at 00:09 UT on 12 December 2008. The data are processed to arrive at the corrected L3 K-corona, which uses a 1-day lower quartile mean background subtraction, centered on the image of interest, and a 21×21 robust mean excluding image pixels that exceed two median absolute deviations. The elements of the observation matrix, H, are computed for each LOS over a grid with a resolution of 1∕160 AU radially and 1.5 • in longitude, and an initial guess at the density, x (0) , is taken to be n e = 0.5r −2 , which is one tenth of the assumed background K-corona. The solving algorithm is then performed on equation (2) for 100 iterations for each image pair. This number is chosen because, due to the regularization terms, the algorithm typically ceases to converge beyond this point.
The solution from the pair of HI-1 images at 05:29 UT on 13 December is shown in Figure 4c. It should be emphasized that here, and in subsequent figures, the densities are multiplied by r 2 in order to account for the steep decrease in density with increasing distance from the Sun. For illustrative purposes the solution is also shown with various levels of regularization present. Panel (a) contains no regularization, while panel (b) shows the solution when the solving algorithm is required to enforce spatial smoothness by means of the gradient terms in equation (5). Finally, panel (c) also includes the use of the weighting term in equation (6). It can be seen that the regularization helps us to arrive at a more physically consistent result. The negative densities are permitted because we are observing the excess density above the K-coronal background. The gradient and weighting terms serve to provide a smooth density and remove a lot of the LOS smearing that is seen in (a). The result in (c), however, is still rather amorphous and does not conform to the conventional three-part CME model of a dense shock/leading edge preceding a flux rope cavity of low electron density followed by a dense core. Instead, the density is concentrated in single, large droplet. This is consistent with the findings of (Davila, 1994), who used rather primitive simulated data to study the effect of spacecraft number and separation on solar tomography. They found that a reconstruction based on data from two spacecraft, separated by 90 • , lost all detail regarding "fine-scale structure" but were still capable of locating the "centroid of each major emission region". This is seen too in the results presented here, which means that we learn little of CME morphology from a two-spacecraft reconstruction. However, the major emission region in this case is a CME and so tomography provides a means to track its position in successive images. Figure 5 shows six examples of the density solution derived from tomographic inversion of HI-1 images, beginning at 11:29 UT on 12 December 2008 and at 6-hr intervals thereafter. This covers the period before, and during, the time that the CME is present in the common FOV. Each solution is achieved using the same regularization and weighting constraints and so Figure 5d is therefore identical to panel (c) in the previous figure. There are some general trends that are observed in each solution. When the CME is observed by HI-1 Figure 5. Reconstructed electron density distributions from six pairs of HI-1 images, each separated in time by 6 hr. The time of observation is shown above each panel. Panel (a) shows the reconstructed density distribution prior to the CME entering the common HI-1 FOV. In panel (b) the CME is seen by both HI-1A and HI-1B, but we are unable to identify it in our density reconstruction. In the following panels (c-f), the CME is distinguishable, to some extent, and we are able to track its passage through the heliosphere out to beyond 0.3 AU.
(panels b-f), we see a significant amount of LOS smearing. This is particularly apparent on the left-hand side (toward STEREO-A) of panels (b) and (c). We attribute this to inconsistent background subtraction between data from the two observing spacecraft. This is good evidence that more thorough background separation methods are required in order to perform this technique more successfully. A further feature common to all panels is an inconsistently low or high density observed in the lower left side of the approximately trapezoidal region. This results from the use of an explicit forward difference scheme to calculate the longitudinal gradient operator, D , whereby the cells on the anticlockwise "edge" of the solution region have no gradient constraints because they are adjacent to "empty" cells.
The CME we are observing first appeared in the respective HI-1A and HI-1B FOVs at 15:29 and 13:29 UT on 12 December. Hence, panel (a) shows a measurement of the ambient solar wind density at 11:29 UT, preceding the first observations of the CME. We generally see a low density throughout the observed region, with higher density in the edges that correspond to the outer limit of the HI-1A and HI-1B FOVs. Because our L3 data are a measure of excess density, we do not expect them to produce a good reconstruction in the absence of a CME and this is apparent in panel (a). Panel (b) shows the tomographic inversion of HI-1 images at 17:29 UT on 12 December, shortly after the CME entered both FOVs. The density appears less smooth than in panel (a); however, we see nothing that resembles a CME, but instead a large density enhancement smeared along the LOS of the inner edge of HI-1A. Panel (c) shows the tomographic inversion of the HI-1 data at 23:29 UT on 12 December, the pair of images shown in Figure 1. Although Figure 1 clearly shows a CME, this is not reproduced by the tomography, which again contains a large amount of LOS smearing. Again, we attribute this to poor background separation in the region of the HI-1A FOV near the Sun. By 05:29 UT on 13 December, in panel (d), we can clearly distinguish a distinct density enhancement that is concentrated in the central region of the plot, with a peak at 0.2 AU along the Sun-Earth line. This dense central structure is associated with the CME observed in HI-1 images; however, again, we see some LOS smearing. There is also a secondary, smaller and less clear density structure that is present at approximately x =0.05 AU, =0.22 AU, which is potentially associated with the unresolved CME leading edge. This is not entirely distinguishable from LOS smearing toward STEREO-A, however. Six hours later, at 11:29UT in panel (e), the central density structure has moved further through the HI-1 FOVs to 0.25 AU, still on the Sun-Earth line. The secondary feature is still visible and has also moved with the central structure. By 17:29 UT, in panel (f), these structures are still visible but are beginning to become indistinguishable from the background solar wind. The central density structure is now at 0.3 AU on the Sun-Earth line and the secondary structure has reached the edge of the observable region of space. These 6-hourly intervals are but a few examples of our density reconstructions; the HI-1 image cadence is 40 min and so a 6-hr period contains nine pairs of images. For each pair of HI-1 images between 00:49 and 22:09UT on 13 December, the CME is distinguishable as a significant density enhancement above the background level. This corresponds to 22 hr and 40 min of observations, or 34 pairs of images, for which we can record the location and value of the peak in density as a means to track the CME. Figure 6 shows the HI-1A and HI-1B observations and their reconstructions presented in the form of time-elongation profiles, often referred to as J-maps (e.g., Davies et al., 2009), corresponding to the ecliptic plane. This serves to illustrate the fidelity between the HI-1 L3 data and the values computed from each density solution. The latter is achieved by applying equation (2) to each density solution (i.e., the computed radiance, y c , is equal to y c = Hx, wherex is the density solution). The fractional residual between the observed and computed radiances is also plotted. The observed radiance J-maps (top) show a bright front, corresponding to the CME, entering the near-Sun edge of the HI-1A and HI-1B FOVs at 15:29 and 13:29 UT on 12 December, respectively. These propagate outward for a period of just over 24 hr, until they become too faint to distinguish at the outer edges of the HI-1 FOVs. In both the HI-1A and HI-1B data, we observe a bright band that falls under the black diamonds. This is followed by a darker region behind the CME. These 10.1029/2019JA027175 features are thought to correspond to the dense buildup of material in front of the CME and the rarefied region behind it in the magnetic cavity.
A significant difference between the observed and reconstructed radiances is that the latter appear much smoother. This results from the use of the regularization terms, which prevent the solution from converging exactly toward the data. However, as shown in Figure 4, this regularization is necessary to arrive at a physically meaningful solution. Despite this effect, the root-mean-square value of the residual over all image pixels is 19% in HI-1A and 23% in HI-1B. However, this is generally lower during the time in which the CME is observed, because our method is tailored toward observing excess density. We also expect there to be some smoothing present in the observed radiances values due to CME motion. A typical CME traveling at 400 km s −1 will travel through approximately nine HI-1 pixels during the 20-min exposure time in a single image. However, this is a problem that is difficult to avoid when using HIs because the long exposure time is required when observing at such great distances from the Sun. There are two vertical structures in the observed HI-1B data at around 01:29 and 08:09 UT on 12 December, and a third at around 14:09 UT on 13 December, that are due to erroneous background subtraction in HI-1B images. This prevents the solving algorithm from converging because there is no physically meaningful density solution that satisfies these data. As a result, we observe large residual values at these times in the reconstruction for both spacecraft. However, two of these occur before the CME enters the FOV of either HI-1A or HI-1B. In the computed radiance plots, we see similar features to those seen in the observed radiance plots. The bright feature associated with this CME can be seen and is at the correct elongation. However, it is smoothed out somewhat due to the regularization. During the period between 16:09 and 19:29 UT on 12 December we see the residual in HI-1B increases noticeably; the mean value is 29% during this period. This is shortly after the CME enters the FOV common to HI-1 instruments but before we are able to reconstruct it as a clear density enhancement.
As noted previously, the CME is visible as a distinct density enhancement for a period of 22 hr (34 image pairs) from 00:49 UT on 13 December to 22:09 UT the same day. We assume the location of the peak density of this region to be an approximate measure of the CME apex location. This is of course a gross oversimplification. The black and red diamonds overplotted in Figure 6 show the elongation along which our CME apex lies for each of the 28 image pairs. In both the HI-1A and HI-1B plots we see that this corresponds to the bright region following the dense CME leading edge. In the previous Figure 5 we saw that the HI-1 reconstruction was unable to reproduce any of the expected CME structure. However, we might expect the densest part of the CME appear ahead of the CME flux rope, that is, in the vicinity of the diamonds in Figure 6. In the middle row of Figure 6, we see that the absolute residual values are greater in the region of the CME apex. This is because we are forcing the solution to have a smooth density gradient when in fact we would expect a sharp change in density between the dense CME front and both the ambient solar wind ahead and the rarefied flux rope behind. Figure 7 is a time-height profile of electron density between the Sun and Earth out to approximately 0.4 AU, produced using all 108 of the tomographic reconstructions performed on the HI data. It is constructed by taking the density values in the row of grid cells that lie on the Sun-Earth line for each time step and stacking them horizontally in sequence. The path of the CME through this figure is again represented by diamonds and we see a very flat gradient that conforms to our expectation that the CME possesses a near constant speed.

In Situ Observations at Earth
While the consistency of the reconstruction with the data in Figure 6 is shown to be good, the consistency of the density estimates are more difficult to establish. We know from the lack of structure in the density reconstruction that the solution is not representative of the structure of the CME. However, we have quantified the density and inferred the motion of the CME from successive pairs of HI-1 images. Additionally, we know this CME to be Earth impacting and therefore have a means to validate our estimates of density and speed against data recorded in situ by Wind at L1. The CME is not observable in HI-1 images beyond approximate 0.35 AU, and so we must extrapolate our observations to 1 AU. The radial evolution of CME speed and density were addressed earlier in this paper and as such we can, in this case at least, make the assumptions that the CME is moving at a constant speed and that its density obeys a power law. Therefore, we make a linear fit to our CME position, as a function of time, from which the gradient reveals the CME speed. We likewise fit a power law to the density, as a function of distance, to arrive at an estimate of the density profile. In Figure 8, the extrapolated velocity and density profiles are plotted. Panel (a) shows the heliocentric distance of the CME density peak as a function of time, during the 22-hr period over which it is tracked. The overplotted red line shows a linear fit to CME position as a function of time. In fact, if we apply a second order polynomial we arrive at a value of a = −2.5 m s −2 for acceleration. It is clear, however, that the motion is indeed very linear, which is consistent with established CME behavior in this region of the heliosphere. As such, we instead use the linear fit for simplicity, which gives a constant speed of v = 375 ± 5 m s −1 . This speed is quite consistent with the speed of the sheath observed at L1, 355 ± 9 km s −1 (Möstl et al., 2014). Kinematic studies of this CME using HI reveal speeds of 411 ± 23 km s −1 (HI-A) and 417 ± 15 km s −1 (HI-B) ), while Liu et al. (2010 combine observations from both spacecraft to arrive at a speed of 363 ± 43 km s −1 . Davies et al. (2012) apply single-spacecraft tracking techniques using HI-1 and HI-2. They apply a range of CME half-widths between 0 • and 90 • to arrive at a range of CME speeds of 378 to 387 km s −1 for STEREO-A and 404 to 413 km s −1 for STEREO-B. Davies et al. (2013) apply stereoscopic tracking using the HI cameras on both spacecraft and determine speeds that are consistent with their single-spacecraft Figure 8. (a) CME heliocentric distance as a function of time. The red lines represents a linear fit, which corresponds to a constant CME speed of 375 ± 5 km s −1 . (b) CME peak density as a function of heliocentric distance. The red lines overplotted represents a power law, which gives a density of 14.7 ± 0.8r −1.90±0.04 . measurements. Figure 8b shows the radial evolution of the CME density peak, during the 22 hr of CME measurements. We observe an electron density profile of (14.7 ± 0.8)r (−1.90±0.04) cc −1 , which is reasonable compared to the average CME value of 6.2r −2.3 cc −1 derived by Richardson et al. (2005). From our profile we would, of course, expect to observe a density very close to 14.7 ± 1 cc −1 at L1, which is somewhat lower than the sheath density measured in situ of 16 ± 4 cc −1 . This in situ value, from Möstl et al. (2014), is derived by taking the mean density within the sheath, while the peak value seen by Wind was above 25 cc −1 . Our value represents the peak the density excess of the CME, on top of the K-coronal background, as determined by HI-1 tomography. If we add a typical background of 5r −2 electrons cc −1 (Cox, 2000) to our value then it is reasonably consistent with that measured in situ. Finally, Möstl et al. (2014) quote the arrival time of the shock front in situ at L1 to be 06:36 UT on 16 December. Extrapolating the position of the CME derived from HI-1 to L1, we find an arrival time of 22:53 UT on 16 December, which is over 16 hours later than observed. However. This is explained by the failure of the tomographic inversion to reproduce any detailed structure of the CME, specifically the shock. This can be understood more clearly by comparison with the Wind in situ data. Figure 9 shows the solar wind parameters of the CME on 16-17 December 2008 measured in situ at L1 by Wind. The faster, dense CME sheath is indicated by the first two vertical lines, while the flux rope is the magnetized and rarefied region of plasma situated between the second and third vertical lines. The velocity, and its uncertainty, derived from the HI-1 tomography is overplotted in black, with uncertainties in light blue (third panel). In the final panel of Figure 9 we see the in situ plasma density with a HI-derived estimate over plotted. This is achieved by taking the HI-derived values of density along the Sun-Earth line, as with Figure 7, and time shifting them to L1 according to our estimated speed of 375 km s −1 . Each density value 10.1029/2019JA027175 is also scaled by r −2 to account for the radial density profile, and we only include values from times where the CME is present, that is, the columns in Figure 7 where the white diamonds are overplotted (00:49 to 22:09 UT on 13 December). This provides estimated values of density at L1 as a function of time. These values are averaged hourly to produce the mean and standard deviation. Because this represents excess density, above the background K-corona, we add a baseline of 5 electrons cc −1 in order to account for the subtracted K-corona, which is the background that we have assumed throughout our analysis based on Cox (2000). These values are plotted in the bottom panel of Figure 9.
Due to the limited amount of density information that we infer from HI-1 tomography, our speed estimate represents a single value corresponding to the CME and we learn nothing of the background solar wind. Therefore, we see no variation of velocity with time, as is observed by Wind. However, the dense region of plasma in the CME sheath is exactly what we are observing in HI and our speed estimate is quite consistent with the sheath velocity seen in Wind data. The lack of detailed structure in our derived density distribution is again apparent when compared to the complexity observed in situ, which is attributed to the lack of information available from the two vantage points of STEREO. Our density is much smoother as a function of time, while conversely there are very sudden changes observed by Wind, particularly in the region of the CME sheath. The smoothness is in fact enforced on the density solution by the gradient terms that we include in the solving algorithm. The reason for including these terms is to force the algorithm to reduce LOS smearing in the solution, caused by the limited points of view. However, using HI we estimate a peak density at Earth approximately 22:00 UT on December 16, which is very close to the peak density observed by Wind. This explains the discrepancy in arrival time compared with the value from Möstl et al. (2014), where the authors were in fact measuring the arrival time of the initial high-density front, and we are identifying the combined total density within the CME. The HI-1 density also shows a rarefied region of plasma in the region in which the flux rope is measured by Wind, although all of the density structure ahead of the sheath and behind the flux rope is missed. Of course, the main reason for attempting to extrapolate these densities out to 1 AU was as a means to test the fidelity of the tomographic reconstruction. However, it is of note that an attempt to make this kind of measurement with HI, and to compare it to in situ data near Earth, has not been undertaken previously.

Conclusions
We have demonstrated a method for modeling a CME in the heliosphere using tomographic inversion of visible-light images from the HI-1 instruments on board each STEREO spacecraft. The method is capable of locating the presence of a density enhancement associated with the CME launched on 12 December 2008; however, it is unable to reproduce any small-scale structure. While it is known from previous studies (e.g., Davila, 1994;Jackson & Froehling, 1995) that this method is ineffective using classical inversion techniques from just two viewpoints, ours is the first attempt to apply the method to HI data. This is not to say that inversion methods are not possible using just two viewpoints, however. The models developed by Jackson et al. (1998) and Kojima et al. (1998) were successful and allowed density reconstructions to be achieved using both HELIOS (Jackson & Hick, 2002 and SMEI observations (Jackson et al., , 2006 with the inclusion of some assumptions of the physics of the solar wind. Indeed, in most fields outside of medical imaging, it is commonplace to include some modeling and physical assumptions to account for the lack of available viewpoints (e.g., Ma et al., 2005 in ionospheric mapping and Li et al., 2012 in seismic imaging). However, our method presents an alternative to the typical approach of simply tracking CME features in HI observations, whereby the observed leading edge from one or more spacecraft is used to infer the propagation of these features. Our approach negates the need for some of the assumptions that these methods often rely on, such as fixed CME morphology. Additionally, it has the potential to gain some inference of CME internal structure in the form of density distribution. When compared to values observed by Wind at L1, we find that our velocity estimate for this CME is indeed very good as we track it through our fields of view. Our density estimate is also reasonably consistent, albeit smoothed, because we have an insufficient number of vantage points with which to reconstruct the detailed structure of the CME. Our reconstructed density distribution is smeared out along the two spacecraft LOS directions and results in much lower values than those which are observed in situ. This is made worse by the fact that we must enforce smoothing on the solving algorithm to limit this LOS smearing. The most significant improvement to this method would therefore be to have more observing spacecraft that are well separated spatially, which would allow the solution to be better constrained and the regularization to be relaxed. Observations of this CME from LASCO, situated at L1, are available and would provide a third vantage point to potentially improve the solution. While we have since lost contact with STEREO-B, the recent launch of Parker Solar Probe, and forthcoming launch of Solar orbiter, will both provide wide-angle imagers (WISPR and SoloHI, respectively). Additionally, proposed future heliospheric imaging missions such as at L1, L5, and the planned PUNCH mission in Earth orbit mean that we have the potential to perform this technique with much greater success in the near future.
A real difficulty in observing CMEs with HI is the background subtraction discussed throughout this paper. This is less of an issue for the feature-tracking methods, which require only the angular position of the observed front, as a function of time, in order to track CME motion. Photometric techniques, such as tomographic inversion, require a precise measurement of K-coronal brightness from every observing spacecraft in order to return a valid density solution. More sophisticated image processing has been applied in other studies to HI-1A data (e.g., Howard et al., 2012Howard et al., , 2013Howard & DeForest, 2014, which are able to clearly reveal small-scale structure in the K-corona. The potential to develop these techniques on future wide-angle imagers means that we may soon have multiple spacecraft capable of observing the K-corona to high degree of precision. In this study we have only used HI-1 data, because our background subtraction in HI-2 is too inaccurate. These improvements toward image processing techniques therefore also have the potential to extend the reconstruction of CMEs in the heliosphere to greater distances still. Effects such as the presence of the galactic plane, onboard image processing problems and the presence of planets, comets, and other bright bodies make many periods of HI-1 data unusable for tomography. Additionally, the optimal spacecraft separation during the time period studied in this paper does not exist for the majority of the rest of the STEREO mission. For these reasons we have only presented one CME as a single case study to demonstrate the tomographic inversion technique. It is not clear how effective this technique would be for faster CMEs, or whether it would be possible at all for non-Earth-directed CMEs. While the analysis presented here produces two-dimensional density estimates, the theory and method are equally applicable in three dimensions. Performing such three-dimensional analysis would have the potential to provide more information on the CME structure; however, this would be unlikely to improve the inability to reproduce small-scale structure. Performing the analysis in three dimensions, however, greatly increases the number of image pixels in the array and the number of grid cells in the array x in equation (2) and therefore increases the size of the matrix H significantly. Further to this, a third gradient term is required in equation (5) to ensure smoothness in the polar direction, which increases the size of the equation further still and makes the computational demands of solving the equation considerably higher. A solution is still certainly possible if we greatly reduce the resolution of both the grid and of the image. However, for the sake of this study, because we were concerned with the Earth-impacting CME signature, we instead opted to perform the analysis in the ecliptic only.
A further means by which this method would be greatly improved is by polarized brightness measurements light. The form of the Thomson scattering equations is such that the contribution to total observed brightness from the polarized and unpolarized components of scattered light have a different dependence on where along the LOS an electron lies. The consequence of this is that polarized brightness measurements provide more information with which to constrain the spatial distribution of electron density. This kind of polarimetric reconstruction is routinely applied to coronagraph observations. While none of the past and current HIs (and neither of SoloHI or WISPR) possess this technology, the planned PUNCH mission will do. This is an important new technology that further improves the potential to perform successful tomographic inversion with wide-angle imagers.