Reconstruction of precipitating electrons and three-dimensional structure of a pulsating auroral patch from monochromatic auroral images obtained from multiple observation points

. In recent years, aurora observation networks using high-sensitivity cameras have been developed in the polar regions. These networks allow dimmer auroras such as pulsating auroras (PsAs) to be observed with a high signal-to-noise ratio. We reconstructed the horizontal distribution of precipitating electrons using computed tomography with monochromatic 20 PsA images obtained from three observation points. The three-dimensional distribution of the volume emission rate (VER) of the PsA was also reconstructed. The characteristic energy of the reconstructed precipitating electron flux ranged from 6 to 23 keV, and the peak altitude of the reconstructed VER ranged from 90 to 104 km. We evaluated the results using a model aurora and compared the model’s electron density with the observed one. The electron density was reconstructed correctly to some extent, even after a decrease in PsA intensity. These results suggest that the horizontal distribution of precipitating electrons 25 associated with PsAs can be effectively reconstructed from ground-based optical observations


Introduction
Aurora computed tomography (ACT) is a method for reconstructing the three-dimensional (3-D) volume emission rate (VER) of auroral emission based on monochromatic auroral images obtained from multiple observation points (e.g., Aso et al., 1990).
The horizontal distribution of precipitating electron flux can be simultaneously obtained using ACT without rocket or satellite observations (Tanaka et al., 2011).Previous studies have applied ACT to bright and well-shaped discrete auroras, such as the quiet arc during the substorm growth phase and multiple auroral arcs (Aso et al., 1990(Aso et al., , 1993(Aso et al., , 1998;;Frey et al., 1996;Nygrén et al., 1997;Tanaka et al., 2011).However, ACT has not been applied to pulsating auroras (PsAs).
A PsA is a type of diffuse aurora that appears as irregular patches showing quasi-periodic on-off switching of its intensity with a periodicity of ~2-20 s (Yamamoto, 1988).The intensity is somewhat dimmer than that of a typical discrete aurora (some hundreds of R up to tens of kR at 557.7 nm; a few hundred R to ~10 kR at 427.8 nm) (McEwen et al., 1981).It has been difficult to apply ACT to PsAs because the signal-to-noise ratio (SNR) of PsA images is lower than those of discrete aurora images.However, remote operation of many high-sensitivity cameras via the internet and an archive system capable of storing a massive amount of aurora data make it possible to observe PsAs with a high SNR.
The Magnetometers Ionospheric Radars All-sky Cameras Large Experiment (MIRACLE) network consists of nine all-sky cameras (ASCs) located in the Fennoscandian region.Two of the ASCs with intensified charge-coupled devices (ICCDs) were replaced with cameras possessing the newer technology of electron-multiplying CCDs (EMCCDs) in 2007 (Sangalli et al., 2011).Ogawa et al. (2020) developed a low-cost multi-wavelength imaging system for aurora and airglow studies and installed Watec monochromatic imagers (WMIs) at several locations in the north and south polar regions.A WMI consists of a highly sensitive CCD camera made by Watec Co., Ltd (Japan).These cameras are suitable for studying very faint auroral structures such as PsAs.In this study, we attempted to use these high-sensitivity cameras and ACT methods to reconstruct the 3-D VER of a PsA patch and the horizontal distribution of precipitating electrons for the first time.the European incoherent scatter (EISCAT) radar operates is also shown.We selected 427.8 nm auroral images in which PsA patches were detected at the EISCAT radar observation point.The reconstructed results were compared with the electron density observed using the EISCAT radar in Sect.3.4.Figure 1a shows 427.8 nm auroral images obtained by the three ASCs from 00:53:30 UT to 00:53:42 UT.The temporal resolution of each ASCs was 2 s.A median filter of 3 × 3 pixels was applied to auroral images to improve the SNR.We also composited auroral images obtained from four WMI CCD cameras of the same type at SKB.The auroral images at SKB were calibrated to derive the absolute emission intensity using the optical facilities at the National Institute of Polar Research, Japan (Ogawa et al., 2020b), and have a time ambiguity of ~1-2 s.We determined the time when the auroral image was obtained by aligning the temporal changes in the PsA patch as shown in auroral images from ABK, KIL, and SKB.The ACT method used in this study is based on the method proposed by Tanaka et al. (2011).We adopted an oblique coordinate system with the origin (O) at coordinates of (69.4°N, 19.2°E).The X-axis was anti-parallel to the horizontal component of the geomagnetic field, the Y-axis was eastward, and the Z-axis was anti-parallel to the geomagnetic field and perpendicular to the Y-axis (see Fig. 2 in Tanaka et al. (2011)).The simulation region ranged from −75 to 75 km, from −100 to 100 km, and from 80 to 180 km for the X-, Y-, and Z-axes, respectively.We set the energy (E) range to extend from 300 eV to 100 keV.The energy axis contains the information on the auroral emission altitude.This is because 65 the higher the electron energy, the lower the stopping height of the precipitating electrons.This reconstruction region was divided linearly into nx × ny × nz voxels along the X-, Y-, and Z-axes and logarithmically into nE bins in the E direction.We set the parameters (nx, ny, nz, nE) to (75, 100, 50, 50), corresponding to a spatial mesh size of 2 × 2 × 2 km.These parameters were selected so that each voxel has at least one line-of-sight crossing from the pixels in the auroral images.

MIRACLE
The differential flux of precipitating electrons was reconstructed by maximizing the posterior probability (| &), where f is a vector of the differential flux of precipitating electrons and  & is a vector of gray levels at pixels in the auroral images obtained with ASCs.According to Bayes' theorem, the posterior probability (| &) is given by (Tanaka et al., 2011) where () is a vector of grey levels obtained by line-integrating the VER in the line-of-sight direction from each pixel (Eq. (8) in Tanaka et al. (2011)).The VER was derived from model  using the aurora emission model (Eq.(3) in Tanaka et al. (2011)).Σ $! is the inverse covariance matrix,  is the variance of ∇ " , and the second-order derivative of  is taken with respect to x, y, and E. The second term in Eq. ( 1) represents the smoothness of  in space and energy directions.We determined Σ $! as the standard deviation calculated from each auroral image.To determine the standard deviation at each pixel in an auroral image, we calculated the mean value and standard deviation in the central 5 × 5-pixel regions for 120 images.Next, we derived a regression line between the mean value and the standard deviation.Finally, we converted the grey level at each pixel to the standard deviation using the equation of the derived regression line.To maximize the posterior probability, it is necessary to minimize the function where ,  ) , and cj are the so-called hyperparameters, which are constants corresponding to the weighting factors for the spatial () and energy ( ) ) derivative terms and the correction factors for the relative sensitivity between cameras (cj), respectively.
The subscript j signifies the three observation points (ABK, KIL, and SKB).The parameter  ./0was fixed at 1.The summation was conducted for the first term in Eq. ( 2) because cj and  * $! were different for the three ASCs.
We carried out the change of variables  = exp() to take advantage of the non-negative constraint on the differential flux f (i.e.,  ≥ 0).We then minimized the function (; ,  ) ,  10/ ,  /23 ) by implementing the Gauss-Newton algorithm with the initial value  (5) = log/ (5) 0, where The hyperparameters were determined using the fivefold cross-validation method (Stone, 1974).First, elements of the vector ) with a trial-and-error method.In addition, the number of iterations for the Gauss-Newton algorithm was also simultaneously determined to be 200 to minimize  ̅ .
The PsA patches shown in Fig. 1a are embedded in the background diffuse auroral emission.We found that a horizontally uniform diffuse aurora causes ambiguity in the reconstruction result because the altitude of the uniform auroral structure cannot be determined from the single-wavelength images.Thus, we subtracted the background emission from the images prior to ACT reconstruction.We created the background emission image by assuming the same value for all voxels.The background VER was taken to be 75 cm −3 s −1 , corresponding to the spatially averaged observed background emission intensity.In this analysis, we assumed that the diffuse aurora showed a uniform VER in all voxels; however, the VER of the diffuse aurora depends generally on the altitude.We note that the altitude dependence of the VER did not affect the analysis result.This is because the horizontal distribution of the background emission intensity (i.e., diffuse aurora) in the auroral image was dependent only on the zenith angle  (∝ cos ).This distribution is not dependent on the altitude distribution of the VER, if the VER is horizontally uniform.

Reconstruction of a pulsating aurora patch model
We reconstructed a PsA patch model from pseudo auroral images to evaluate the analytical error of ACT before reconstructing the PsA patch from the observed auroral images.To create the pseudo auroral images, we prepared the horizontal distributions of the total energy, Q0, and the characteristic energy, Ec.We then derived the 3-D VER, L, as shown in Fig. 2a.The total energy was assumed to have a Gaussian shape in horizontal directions with a maximum value of 6 mW m −2 .The energy distribution was considered to be a Maxwellian distribution with a uniform characteristic energy of 15 keV.Pseudo auroral images were obtained from L by solving the forward problem (Fig. 2b).We added random noise from a normal distribution with a mean value of 0 and the standard deviation determined from observed auroral images.
Figure 2c shows Q0, Ec, and L reconstructed from the pseudo auroral images.The values of Q0 were calculated as When we assume the energy distribution to be a Maxwellian distribution, the characteristic energy can be written as . We calculated the errors between the model and the result for Q0, Ec, and L (Fig. 2d).The median values of the errors were −5% for Q0, −21% for Ec, and −11% for L. The northwestern part of Q0 was overestimated by at most 23%, the edge part was underestimated by at most 29%, and the central part was underestimated by ~8%.The central part of Ec was reconstructed correctly.In comparison, the edge part (especially the northwestern part) was underestimated by at most 56%.The underestimation of Ec was caused by the overestimation of the emission altitude (Fig. 2d).patch is vertically thin and horizontally wide.In addition, the SNR at the edge part is lower than at the central part since we assumed a Gaussian shape for the horizontal distribution of Q0.These factors would tend to reduce the accuracy at the edge part.In this particular event, the optical observation points of ABK, SKB, and KIL were located at the southeastern direction from the PsA patch at the EISCAT radar observation point (TRO) (see Figure 1).Therefore, an accurate reconstruction of the 130 VER peak altitude was more difficult at the northwestern part than at other parts.This was because of the insufficient information on the northwestern part of the PsA patch owing to the bias of the optical observation point distribution.
Notably, the reconstructed results using the hyperparameters determined by the cross-validation method revealed unexpected fine structures.To avoid this phenomenon, we set the lower limit of  by a different method, namely by minimizing the sum of the squares of the residuals between the model and the reconstructed result of Q0 and Ec.The lower limit on  makes it challenging to reconstruct actual fine-scale structures in the patches.

Precipitating electrons
Figures 3a and 3b show Q0 and Ec as reconstructed from the observed auroral images (Fig. 1a).The maximum value of Q0 was ~6 mW m −2 .The reconstructed Ec ranged from 6 to 23 keV.These energies are consistent with observation results from sounding rockets and low-altitude satellites (e.g., McEwen et al., 1981;Miyoshi et al., 2015).The horizontal distribution of Ec was neither uniform nor stable in the patch during the pulsation.Thus, the ACT method is helpful for investigating PsAassociated temporal variations in the horizontal distribution of precipitating electrons without rocket or satellite observations.In particular, the southwestern part of Ec was enhanced at 00:53:38 UT.It should be noted that the edge and northwestern parts of Ec are expected to be underestimated owing to analytical error, as shown in Fig. 2d.These temporal variations indicate changes in the cyclotron resonance energy of whistler-mode chorus waves during the pulsation.The chorus waves scatter electrons into a loss cone near the magnetic equator.The cyclotron resonance energy of chorus waves depends on the background magnetic field, electron density, and wave frequency (e.g., Kennel and Petschek, 1966).Therefore, the observed temporal variations indicate changes in the magnetospheric source region's background magnetic or plasma environment during the pulsation.
The characteristic energy can also be enhanced by the field aligned current (FAC).The absence of FAC in the PsA patch has been observed because the PsA patch have no shear motion and the energy spectra of precipitating electrons observed by rockets and satellites did not show a field-aligned acceleration, which are typically observed in discrete aurora (Davis, 1978).On the other hand, several studies have reported the FAC associated with PsA patches (Fujii et al., 1985;Gillies et al., 2015;Hosokawa and Ogawa, 2010).If the upward and downward FACs flow at the edge of the PsA path, the potential drop due to the upward FAC can accelerate precipitating electrons and enhance their characteristic energy (Sato et al., 2004;Shepherd and Fälthammar, 1980).However, we cannot suggest the existence of the FAC at the PsA patch's edge for this evet because the total energy was not enhanced at the edge (Figure 3).If the FAC exists, both total and characteristic energy should be increased.
Multievent analysis is needed to examine the FAC in PsA patches.The 3-D current structure in the PsA patch is beyond the scope of this study.Its reconstruction using ACT and EISCAT_3D radar is planned in the future.

Volume emission rate
The 3-D distributions of the VERs derived from the reconstructed electron flux was obtained by solving the forward problem described in Section 2 (see Figure 4a).Cross-sections in the horizontal plane at an altitude of 94 km are shown in Fig. 4b.The peak altitude ranges from 90 to 104 km (Fig. 4c).The full width at half maximum is almost uniform with a median value of ~20 km (Fig. 4e).The errors of the peak altitude and the altitude width in Figures 4d and 4f were derived using the model and reconstructed VER shown in Figure 2. The high peak altitude at the northwestern part is expected to be overestimated by at most 8% owing to analytical error, as shown in Fig. 2d.For the altitude width, the most part is expected to be overestimated by ~2% (Fig. 4f).The reconstructed peak altitude and width are consistent with those determined in studies using stereoscopic observations or an incoherent scatter radar (Brown et al., 1976;Jones et al., 2009;Kataoka et al., 2016).Stenbaek-Nielsen and Hallinan (1979) reported the existence of thin (<1 km vertical extent) PsA patches based on stereoscopic observations, but our results do not support their results.
The peak altitude of the PsA patch was also estimated by a different method (Appendix A).We projected the observed auroral images at altitudes ranging from 80 km to 120 km with an interval of 2 km (Movie A1, http://doi.org/10.5446/57558).The emission altitude was determined to be the altitude at which the sum of the squares of the residuals between the two projected images reached a minimum value (Fig. A1).The estimated peak altitude range was 92 to 106 km from 00:53:30 UT to 00:53:40 UT (Fig. A2).These altitudes closely match those determined by ACT.

Electron density
The altitude profiles of VER at the EISCAT radar observation point shown in Fig. 4g were converted to the ionospheric electron density and compared with the actual data observed by the EISCAT radar.The continuity equation for the electron density can be written as where ne [m −3 ] is the electron density, L [m −3 s −1 ] is the VER, k is a positive constant for converting VER to the ionization rate (see Appendix B), and  9GG [m 3 s −1 ] is the effective recombination rate.We derived the electron density from the VER by solving Eq. ( 4) with the Runge-Kutta method.The initial value was derived from Eq. (4) under steady-state conditions (i.e., ∂ 9 / = 0) using reconstructed L at 00:53:36 UT.The VERs were interpolated linearly in time to use the Runge-Kutta method.The altitude profile of  9GG has been investigated by several studies using rocket-and ground-based measurements.185 Vickrey et al. (1982) summarized many of these results and proposed the following best fit parameterization:  GH8 = 2.5 × 10 $!" exp(−/51.2) [m 3 s -1 ], (5) where z [km] is the altitude.Semeter & Kamalabadi (2005) used the effective recombination coefficients  IJ $ and  J ! $ for NO + and O2 + , respectively (Walls and Dunn, 1974), as upper and lower bounds on  9GG : Here Tn [K] is the neutral temperature.The red lines in Fig. 5 show the derived electron densities using these three recombination coefficients.We note that these values are underestimated compared to the electron densities observed by the EISCAT radar (black lines in Fig. 5).This underestimation probably comes from the background emission subtraction from the auroral images prior to ACT and from ambiguity in the effective recombination coefficients.The electron densities reconstructed from auroral images without background emission subtraction are shown as blue lines for reference in Fig. 5.
The reconstruction results from the images, including background emission, approached the electron density profile observed with the EISCAT radar.We note that the electron density was reconstructed correctly to some extent even after the auroral emission intensity decreased at 00:53:40 UT.This correct reconstruction considered the time change in the continuity equation.
The electron density would seem to have rapidly decreased after 00:53:40 UT, if the time change was not considered.This result suggested that the time change should be considered ( 9 / ≠ 0 in Eq. ( 4)) when using the continuity equation to derive the electron densities associated with PsAs.It should be noted that the electron density is underestimated at higher altitudes (>~140 km) even if the background emission was included.This underestimation can be improved by reconstructing low-energy electron flux from auroral images of various wavelengths (e.g., 844.6 nm).

Conclusions
We applied the ACT method to a PsA patch for the first time and reconstructed the horizontal distribution of precipitating electrons from 427.8-nm auroral images obtained from three observation points.We improved the previously proposed ACT method by adding the following processes: the subtraction of the background diffuse aurora from the auroral images prior to ACT, the estimation of the relative sensitivity between ASCs, and the determination of the hyperparameters of the regularization term.The characteristic energies of the reconstructed electron fluxes (6 to 23 keV) and the peak altitudes of the reconstructed VERs (90 to 104 km) were consistent with those found in previous studies.We determined that the horizontal distribution of the characteristic energy was neither uniform nor stable in the patch during the pulsation, further underlining the shortcomings of rocket and satellite observations for investigating PsAs.These spatial and temporal variations imply changes of the electron density and magnetic field in the magnetosphere.ACT error was evaluated using an auroral patch model.The characteristic energy of electron flux was correctly reconstructed at the center part of the patch but underestimated at the patch edge by at most 56%.The reconstructed electron flux will be improved in future work by incorporating auroral images of various wavelengths.
Although we reconstructed the differential flux of precipitating electrons from auroral images using ACT, Tanaka et al. (2011) extended ACT to a method called generalized-ACT (G-ACT).G-ACT uses multi-instrument data, such as ionospheric electron density from incoherent scatter radar, cosmic noise absorption from imaging riometers, and the auroral images.They demonstrated that the incorporation of the ionospheric electron density from the EISCAT radar improved the accuracy of the reconstructed electron flux.Furthermore, 3-D ionospheric observation by EISCAT_3D (http://www.eiscat3d.se) is scheduled to begin in 2023.In the future, we will improve the reconstructed electron flux by conducting G-ACT using electron density data from the EISCAT or EISCAT_3D radar.

Figure 2 :
Figure 2: (a) The horizontal distribution of the prepared total energy Q0 and the characteristic energy Ec of precipitating electrons and the vertical cross-section of the volume emission rate L along the dashed lines in the left and middle panels.We derived L from the prepared Q0 and Ec values using the aurora emission model.Q0 and Ec are not shown for Q0 values less than 1 mW m −2 .(b) Pseudo auroral images obtained by line-integrating the VER.(c) The horizontal distribution of Q0 and Ec and the vertical cross-section of L reconstructed by aurora computed tomography from the pseudo auroral images.(d) The errors of Q0, Ec, and L, calculated as (Error) = [(Result) -(Model)] / (Model) × 100.

Figure 3 :
Figure 3: (a) Total energy, Q0, and (b) characteristic energy, Ec, of the precipitating electron flux reconstructed from the observed auroral images.Results of Q0 and Ec where Q0 is less than 1 mW m -2 are not shown.

Figure 4 :
Figure 4: (a) Reconstructed 3-D distribution of volume emission rates (VERs) L. VERs less than 1 cm −3 s −1 are not shown.(b) Cross-sections in the horizontal plane at an altitude of 94 km.VERs are not shown for Q0 values less than 1 mW m −2 .(c) Peak altitudes of the reconstructed L and (d) their errors determined using the model aurora.(e) Altitude widths of the reconstructed L and (f) their errors determined using the model aurora.(g) Altitude profiles of L at the European incoherent scatter radar observation point as indicated by black plus marks in Figs.4b-4f.

Figure 5 :
Figure 5: Electron density altitude profiles, ne, converted from the reconstrucred volume emission rates with the subtraction of the background emission (BGE) (red lines), those without the BGE subtraction (blue lines), and those observed by the European incoherent scatter radar (black lines).Details of effective recombination coefficients   ,    " , and   " are explained in the text.The measurement uncertainties are represented by error bars.