Application of generalized aurora computed tomography to the EISCAT_3D project

EISCAT_3D is a project to build a multi-site phased-array incoherent scatter radar system in northern Fenno-Scandinavia. We demonstrate via numerical simulation how useful monochromatic images taken by a multi-point imager network are for auroral research in the EISCAT_3D project. We apply the generalized aurora computed tomography (G-ACT) method to modelled observational data from real instruments, such as the Auroral Large Imaging System (ALIS) and the EISCAT_3D radar. G-ACT is a method for reconstructing the three-dimensional (3D) distribution of auroral emissions and ionospheric electron density (corresponding to the horizontal two-dimensional (2D) distribution of energy spectra of precipitating electrons) from multi-instrument data. It is assumed that the EISCAT_3D radar scans an area of 0.8° in geographic latitude and 3° in longitude at an altitude of 130 km with 10 × 10 beams from the radar core site at Skibotn (69.35° N, 20.37° E). Two neighboring discrete arcs are assumed to appear in the observation region of the EISCAT_3D radar. The reconstruction results from G-ACT are compared with those from the normal ACT as well as the ionospheric electron density from the radar. It is found that G-ACT can interpolate the ionospheric electron density at a much higher spatial resolution than that observed by the EISCAT_3D radar. Furthermore, the multiple arcs reconstructed by G-ACT are more precise than those by ACT. In particular, underestimation of the ionospheric electron density and precipitating electrons' energy fluxes inside the arcs is significantly improved by G-ACT including the EISCAT_3D data. Even when the ACT reconstruction is difficult due to the unsuitable locations of the imager sites relative to the discrete arcs and/or a small number of available images, G-ACT allows us to obtain better reconstruction results.

1 Introduction

EISCAT_3D is a multi-point phased array incoherent scattering radar system under construction in northern Fenno-Scandinavia as of November 2023 and is expected to be operational in winter 2023. The EISCAT_3D radar will be able to measure the three-dimensional (3D) distribution of ionospheric parameters, such as the electron density, electron temperature, ion temperature, and ion Doppler velocity, at a resolution that is more than 10 times higher than that of the existing EISCAT radar. Thus, it is expected to provide new insights into various science topics pertaining to auroral physics, ionospheric physics, magnetosphere–ionosphere coupling, and so on (McCrea et al., 2015; Wannberg et al., 2010).

The height distribution of the ionospheric electron density in the auroral region, which is related to the energy distribution of auroral precipitating electrons, is essential for clarifying if the precipitating electrons experienced acceleration and to determine their place of origin. In addition, we can estimate the ionospheric conductivity from the electron density by using empirical models (e.g., the Mass Spectrometer and Incoherent Scatter (MSIS) atmosphere model and the International Geomagnetic Reference Field (IGRF) model; Hedin, 1991; Alken et al., 2021). It is well known that the spatial distribution of ionospheric conductivity plays an essential role in the magnetosphere–ionosphere coupling process (e.g., Ellis and Southwood, 1983; Glaßmeier, 1984; Itonaga and Kitamura, 1988). It may be possible to deduce a 3D current system from the ionospheric conductivity distribution by using magnetic field data from a ground-based magnetometer array and/or ionospheric electric field data from radars (Kamide et al., 1981; Vanhamäki and Amm, 2007).

On the other hand, it is useful to utilize optical imaging observations to study auroral dynamics. Radars generally have a high range resolution; however, scanning a particular area with multiple beams is time-consuming. In contrast, an optical imager has a high angular resolution. It can measure angular distributions at a higher temporal resolution than a radar, even though it can detect only the integrated luminosity along the line of sight. In other words, the radar and optical imager are complementary. Therefore, it is essential to use image data with radar data effectively. Figure 1 shows a schematic illustration of the relationship between the radar and imager observations.

Figure 1 Schematic illustration of the relationship between the radar and optical imager observations.

There are some ground-based imager networks in northern Fenno-Scandinavia, e.g., the Aurora Large Imaging System (ALIS) (Brändström, 2003), the all-sky cameras in the Magnetometers – Ionospheric Radars – All-sky Cameras Large Experiment (MIRACLE) (Syrjäsuo, 2001), and Watec Monochromatic Imager (WMI) (Ogawa et al., 2020). In particular, the ALIS was designed to obtain the 3D distribution of the optical emissions in the mesosphere, thermosphere, and ionosphere, and was recently developed into ALIS_4D (https://alis4d.irf.se/, last access: 5 May 2024). By applying the auroral computed tomography (ACT) technique to the monochromatic images taken at some ALIS stations, it is possible to retrieve the 3D distribution of auroras that have a horizontal scale of several tens to hundreds of kilometers (Aso et al., 1998; Gustavsson, 1998; Gustavsson et al., 2001; Simon Wedlund et al., 2013). Conversely, the inverse problem of ACT is ill-posed and ill-conditioned, because the optical image data correspond to the emission intensity integrated along the line of sight and only a few images are usually available. Thus, assumptions need to be made to solve the inverse problem, which often makes it somewhat difficult to interpret the results of the data analysis from a physical point of view.

Aso et al. (2008) and Tanaka et al. (2011) have extended ACT to the generalized ACT (G-ACT). This method can reconstruct the spatial and energy distributions of auroral precipitating electrons by using multi-instrument data, such as the ionospheric electron density from an incoherent scatter radar, the cosmic noise absorption (CNA) from an imaging riometer, and optical monochromatic images. Tanaka et al. (2011) demonstrated via numerical simulation that the reconstruction from only auroral images can be improved by G-ACT by using the height profile of the electron density from the EISCAT radar.

In the current study, we investigated how effective the combination of the EISCAT_3D radar and the monochromatic imager network is for auroral research by conducting a simulation. We apply the G-ACT method to the modelled observational data, i.e., electron density from the EISCAT_3D radar and the multiple monochromatic images from the ALIS, and compare the reconstruction results with those obtained by the normal ACT and the radar's electron density data. We selected the ALIS (not ALIS_4D) as the monochromatic imager network for this simulation study, because the monochromatic images from the ALIS can be used for both the normal ACT analysis and the G-ACT analysis. It is possible to compare the auroral 3D distributions reconstructed by the two methods in the same region because one of the ALIS stations is located in Skibotn (69.35° N, 20.37° E), Norway, which is the core site of the EISCAT_3D radar, and the field of view (FOV) of the ALIS imagers covers the radar observation region.

2 Forward analysis

2.1 Observatories and instruments

Figure 2 shows the locations of the stations used in this simulation study and each instrument's fields of view (FOVs) at an altitude of 130 km. Blue quadrangles and red crosses correspond to the FOVs of the ALIS imagers and the beam position of the EISCAT_3D radar, respectively. The core site of the EISCAT_3D radar is located at Skibotn, Norway.

It was assumed that the EISCAT_3D radar scanned an area of geographic latitude from 68.6 to 69.4° N and longitude from 18.767 to 21.767° E at an altitude of 130 km with 10 × 10 beams, which corresponds to a spatial resolution of 0.08° (about 8.9 km) in latitude and 0.3° (about 12 km) in longitude. It was also assumed that the electron density was detected at altitudes between 90 and 170 km. The dot-dash line indicates the region where the reconstruction results were evaluated, which is the same as the dot-dash line in Fig. 3a.

Figure 2 Locations of the stations used in this study and the fields of view (FOVs) of the ALIS imagers and EISCAT_3D radar at an altitude of 130 km. It was assumed that the viewing direction of the imagers was set to the south of Skibotn and the radar scanned an area of geographic latitude from 68.6 to 69.4° N and longitude from 18.767 to 21.767° E with 10 × 10 beams.

Each ALIS station has a sensitive high-resolution (1024 × 1024 pixels) unintensified monochromatic CCD imager with a six-position filter wheel for narrowband interference filters (427.8, 557.7, 630.0, and 844.6 nm) (Brändström, 2003). The FOV of each imager is about 50 to 90°. It was assumed that the viewing direction was set to the south of Skibotn and the filter was fixed to the N${}_{\mathrm{2}}^{+}$ first negative band (427.8 nm) for all stations. We postulated that the image size was reduced to 256 × 256 pixels after 4 × 4 pixel binning.

Figure 3 (a) Horizontal distribution of the total energy flux (Q₀) of the incident auroral electrons. The top and right correspond to the geomagnetic northward and geomagnetic eastward directions, respectively. The thick line and dot-dash line indicate the region used for the inverse analysis and the region where the reconstruction results were evaluated, respectively. (b) Energy distribution of the incident electrons at the peak location of the discrete arcs.

2.2 Distribution of incident auroral electrons

Figure 3a indicates the horizontal distribution of the total energy flux (Q₀) of the auroral precipitating electrons that was assumed for the forward analysis. It was also assumed that two neighboring discrete arcs appeared over the southern sky of Skibotn. In this simulation, an oblique coordinate system was adopted with an origin at Skibotn, with the x axis pointing in the geomagnetic southward direction, the y axis pointing in the eastward direction, and the z axis anti-parallel to the geomagnetic field (cf., Fig. 2 of Tanaka et al., 2011). The inclination and declination angles of the geomagnetic field were 78 and 6°, respectively. The calculation ranges were −40 to 100 km, −150 to 150 km, and 90 to 190 km for the x, y, and z directions, respectively. The spatial mesh sizes (Δ_x, Δ_y, Δ_z) were 1, 2, and 2 km for the x, y, and z directions, respectively. The discrete arcs were assumed to have a sinusoidal shape in the y direction and a Gaussian shape in the x direction. The distance between the two arcs was 20 km. The area surrounded by a thick line is the reconstruction region used for the inverse analysis in Sect. 3. The dot-dash line indicates the region where the reconstruction results are evaluated in Sect. 4.

Figure 3b shows the energy distribution of the differential number flux of the precipitating electrons at the peak location of the arcs. The energy spectrum is represented by the sum of the Gaussian distribution ($f_{\mathrm{g}}\left(E\right)=A_{\mathrm{g}}\exp \left\{-{\left(E-E_{\mathrm{0}}\right)}^{\mathrm{2}}/W^{\mathrm{2}}\right\}$) and two power-law distributions for the low-energy tail ($f_{\mathrm{pl}}\left(E\right)=A_{\mathrm{pl}}{\left(E/E_{\mathrm{0}}\right)}^{-a}$, E≤E₀) and high-energy tail ($f_{\mathrm{ph}}\left(E\right)=A_{\mathrm{ph}}{\left(E/E_{\mathrm{0}}\right)}^{-b}$, E≥E₀), which has been introduced by Strickland et al. (1993) as a typical spectrum of discrete auroras. E₀, W, a, and b were set to 5 keV, 0.15×E₀, 1.0, and 3.0, respectively, for the entire simulation region. The energy range for the calculation was between 0.3 and 20 keV and divided logarithmically into 50 intervals.

Since the main purpose of this paper is to compare the results from the two analysis methods, ACT and G-ACT, we assumed a rough but typical auroral shape, size, and energy distribution. Actual auroras have a variety of shapes, including multiple arcs, structured shapes, patchy shapes, and very thin arcs less than 100 m thick. However, it is difficult to examine such a large number of auroral types in this paper due to a lack of space. If one wants to evaluate the accuracy of the tomographic analysis results for real auroras, simulations should be performed using modelled auroras that resemble them (e.g., Fukizawa et al., 2022).

2.3 Calculation of modelled data set

The formulation of our method is based on that used by Janhunen (2001). The forward problem was solved by using the distribution of the incident electrons described in Sect. 2.2. The height profile of volume emission rate ${\mathit{L}}_{x_{\mathrm{1}},y_{\mathrm{1}}}\left(z\right)$ along the field line at a certain horizontal coordinate (x₁, y₁) is calculated by

$$\begin{array}{}\text{(1)}& {\mathit{L}}_{x_{\mathrm{1}},y_{\mathrm{1}}}\left(z\right)={\mathit{m}}_{\mathrm{1}}{\mathit{f}}_{x_{\mathrm{1}},y_{\mathrm{1}}}\left(E\right),\end{array}$$

where ${\mathit{f}}_{x_{\mathrm{1}},y_{\mathrm{1}}}\left(E\right)$ is the energy distribution of differential number flux of incident electrons at the top of the ionosphere (x₁, y₁, z_max), and m₁ is a matrix operator for calculating ${\mathit{L}}_{x_{\mathrm{1}},y_{\mathrm{1}}}\left(z\right)$ from ${\mathit{f}}_{x_{\mathrm{1}},y_{\mathrm{1}}}\left(E\right)$. We adopted Rees' model (Rees, 1989) to obtain the energy deposition rate to the atmosphere from the differential flux and the method proposed by Sergienko and Ivanov (1993) to calculate the 427.8 nm volume emission rate from the energy deposition rate. The elements of m₁ are a function of the atmospheric parameters, which were calculated by using the MSIS-90 atmosphere model (Hedin, 1991). The derivation of m₁ is described in detail in the Appendix of Tanaka et al. (2011).

Assuming that m₁ is independent of x and y, Eq. (1) can be expanded in the x and y directions as follows:

$$\begin{array}{}\text{(2)}& \left(\begin{array}{c}{\mathit{L}}_{\mathrm{1},\mathrm{1}}\\ {\mathit{L}}_{\mathrm{1},\mathrm{2}}\\ \mathrm{⋮}\\ {\mathit{L}}_{n_x,n_y}\end{array}\right)=\left(\begin{array}{cccc}{\mathit{m}}_{\mathrm{1}}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{0}& {\mathit{m}}_{\mathrm{1}}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{\ddots }& \mathrm{0}\\ \mathrm{0}& \mathrm{\dots }& \mathrm{0}& {\mathit{m}}_{\mathrm{1}}\end{array}\right)\left(\begin{array}{c}{\mathit{f}}_{\mathrm{1},\mathrm{1}}\\ {\mathit{f}}_{\mathrm{1},\mathrm{2}}\\ \mathrm{⋮}\\ {\mathit{f}}_{n_x,n_y}\end{array}\right),\end{array}$$

$$\begin{array}{}\text{(3)}& \mathit{L}={\mathit{M}}_{\mathrm{1}}\mathit{f}.\end{array}$$

In Eq. (3), f is a function of x, y, and E and has a length of $n=n_E\times n_x\times n_y$, and L is a function of x, y, and z and has a length of $m=n_z\times n_x\times n_y$. M₁ is a large sparse matrix whose size is m×n.

In a similar manner to L, a square of ionospheric electron density $\mathit{D}(x,y,z)$ generated by the incident electrons is given by

$$\begin{array}{}\text{(4)}& \mathit{D}={\mathit{M}}_{\mathrm{2}}\mathit{f},\end{array}$$

where

$$\begin{array}{}\text{(5)}& {\mathit{M}}_{\mathrm{2}}=\left(\begin{array}{cccc}{\mathit{m}}_{\mathrm{2}}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{0}& {\mathit{m}}_{\mathrm{2}}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{\ddots }& \mathrm{0}\\ \mathrm{0}& \mathrm{\dots }& \mathrm{0}& {\mathit{m}}_{\mathrm{2}}\end{array}\right).\end{array}$$

M₂ has the same size as M₁. In Eq. (5), the electron density measured by the EISCAT_3D radar was assumed to be caused by the auroral precipitation only. For the derivation of m₂, we assumed that the ionospheric electron density is quasi-steady-state and the ionization loss in the E-layer is dominated by the recombination process. Again, refer to the Appendix of Tanaka et al. (2011) for more details.

A gray level g_i at a pixel i in the auroral image is approximated by a linear integration along a line of sight, as shown below:

$$\begin{array}{}\text{(6)}& g_i={\displaystyle \frac{c_g\left(\mathit{\theta },\mathit{\phi }\right)}{\mathrm{4}\mathit{\pi }}}\int L\left(r,\mathit{\theta },\mathit{\phi }\right)\mathrm{d}r,\end{array}$$

where r, θ, and φ are polar coordinates whose origin is located at the center of the camera lens, and c_g(θ,φ) is a sensitivity and vignetting factor. Equation (6) can be represented by

$$\begin{array}{}\text{(7)}& \mathit{g}={\mathit{P}}_{\mathrm{1}}\mathit{L}={\mathit{P}}_{\mathrm{1}}{\mathit{M}}_{\mathrm{1}}\mathit{f},\end{array}$$

where g is a gray-level vector which has l_g elements, and P₁ is a l_g× m matrix used to calculate g by integrating L along the line of sight.

The square of electron density observed by the EISCAT_3D radar d is expressed by the matrix expression:

$$\begin{array}{}\text{(8)}& \mathit{d}={\mathit{P}}_{\mathrm{2}}\mathit{D}={\mathit{P}}_{\mathrm{2}}{\mathit{M}}_{\mathrm{2}}\mathit{f},\end{array}$$

where P₂ is a l_d× m matrix that extracts data in the voxels corresponding to the radar observation locations from $\mathit{D}(x,y,z)$.

We added the noise to d and g and finally obtained modelled data, $\stackrel{\mathrm{̃}}{\mathit{d}}$ and $\stackrel{\mathrm{̃}}{\mathit{g}}$. Gaussian noise with a standard deviation of 5 % of the electron density was added to the electron density data. The offset of 300 R was added to the gray level data and then Gaussian noise with a standard deviation of $\sqrt{\mathit{g}+\mathrm{300}}R$ was added. Figure 4a and b show the modelled ionospheric electron density that should be obtained by the EISCAT_3D radar and the modelled auroral images at five ALIS stations.

Figure 4 (a) Modelled ionospheric electron density data obtained by the EISCAT_3D radar. (b) Modelled auroral images taken at five ALIS stations. Top and right of the images correspond to the northward and westward directions, respectively.

3 Inverse analysis

Our inverse analysis method is based on the Bayesian model. According to Bayes' theorem, the probability that model f is true after data $\stackrel{\mathrm{̃}}{\mathit{b}}$ were observed, i.e., the posterior probability $P(\mathit{f}\mathrm{|}\stackrel{\mathrm{̃}}{\mathit{b}}$), is expressed by

$$\begin{array}{}\text{(9)}& P\left(\mathit{f}\mathrm{|}\stackrel{\mathrm{̃}}{\mathit{b}}\right)={\displaystyle \frac{P\left(\stackrel{\mathrm{̃}}{\mathit{b}}\mathrm{|}\mathit{f}\right)P\left(\mathit{f}\right)}{P\left(\stackrel{\mathrm{̃}}{\mathit{b}}\right)}}\propto P\left(\stackrel{\mathrm{̃}}{\mathit{b}}\mathrm{|}\mathit{f}\right)P\left(\mathit{f}\right),\end{array}$$

where $P(\stackrel{\mathrm{̃}}{\mathit{b}}\mathrm{|}\mathit{f}$) is the likelihood, which is the probability of observing data $\stackrel{\mathrm{̃}}{\mathit{b}}$ given model f; P(f) is the prior probability of model f; and $P(\stackrel{\mathrm{̃}}{\mathit{b}}$) is the marginal probability of $\stackrel{\mathrm{̃}}{\mathit{b}}$. In this study, P(f) and $P(\stackrel{\mathrm{̃}}{\mathit{b}}\mathrm{|}\mathit{f}$) are given by

$$\begin{array}{}\text{(10)}& {\displaystyle }P\left(\mathit{f}\right)\propto \exp \left(-{\displaystyle \frac{{∥{\mathrm{\nabla }}^{\mathrm{2}}\mathit{f}∥}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}\right),\text{(11)}& {\displaystyle }\begin{array}{rl}& P\left(\stackrel{\mathrm{̃}}{\mathit{b}}\mathrm{|}\mathit{f}\right)\propto \exp \mathit{\{}-{\sum }_j{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\stackrel{\mathrm{̃}}{\mathit{b}}}_j-{\mathit{b}}_j\left(\mathit{f}\right)\right)}^T\\ & {\mathbf{\Sigma }}_j^{-\mathrm{1}}\left({\stackrel{\mathrm{̃}}{\mathit{b}}}_j-{\mathit{b}}_j\left(\mathit{f}\right)\right)\mathit{\}},\end{array}\end{array}$$

where σ² is the variance of ∇²f, where ${\mathbf{\Sigma }}_j^{-\mathrm{1}}$ is the inverse covariance matrix, and where j means the kind of data. It was assumed that the modelled data are independent from each other, so ${\mathbf{\Sigma }}_j^{-\mathrm{1}}$ has zero off-diagonal elements and the inverse of the variance of ${\stackrel{\mathrm{̃}}{\mathit{b}}}_j$ in the diagonal elements. ${\stackrel{\mathrm{̃}}{\mathit{b}}}_j$ corresponds to the modelled data $\stackrel{\mathrm{̃}}{\mathit{g}}$ and $\stackrel{\mathrm{̃}}{\mathit{d}}$ for j=1 and 2, respectively, and they include the noise. b_j(f) corresponds to g and d in Eqs. (7) and (8). Equation (10) indicates the smoothness constraint on f with respect to x, y, and E. In Eq. (11), it was assumed that the modelled data ${\stackrel{\mathrm{̃}}{\mathit{b}}}_j$ has Gaussian errors. By substituting Eqs. (10) and (11) into Eq. (9), $P(\mathit{f}\mathrm{|}\stackrel{\mathrm{̃}}{\mathit{b}}$) is given by

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& P\left(\mathit{f}\mathrm{|}\stackrel{\mathrm{̃}}{\mathit{b}}\right)\propto \exp [-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}\mathit{\{}\sum _j{w}_j^{\mathrm{2}}{\left({\stackrel{\mathrm{̃}}{\mathit{b}}}_j-{\mathit{b}}_j\left(\mathit{f}\right)\right)}^T\\ & {\mathbf{\Sigma }}_j^{-\mathrm{1}}\left({\stackrel{\mathrm{̃}}{\mathit{b}}}_j-{\mathit{b}}_j\left(\mathit{f}\right)\right)+{∥{\mathrm{\nabla }}^{\mathrm{2}}\mathit{f}∥}^{\mathrm{2}}\mathit{\left\}}\right],\end{array}\end{array}$$

where w_j is a hyper-parameter, which is a constant corresponding to the weighting factor for each instrument data.

Maximization of the posterior probability is equivalent to minimization of the function inside the curly brackets of Eq. (12), which is given by

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& \mathit{\varphi }(\mathit{f};w_j)=\sum _j{w}_j^{\mathrm{2}}{\left({\stackrel{\mathrm{̃}}{\mathit{b}}}_j-{\mathit{b}}_j\left(\mathit{f}\right)\right)}^T\\ & {\mathbf{\Sigma }}_j^{-\mathrm{1}}\left({\stackrel{\mathrm{̃}}{\mathit{b}}}_j-{\mathit{b}}_j\left(\mathit{f}\right)\right)+{∥{\mathrm{\nabla }}^{\mathrm{2}}\mathit{f}∥}^{\mathrm{2}}.\end{array}\end{array}$$

Here, we define $\mathit{r}(\mathit{f};w_{\mathrm{1}},w_{\mathrm{2}})$ as follows;

$$\begin{array}{}\text{(14)}& \mathit{r}\left(\mathit{f};w_{\mathrm{1}},w_{\mathrm{2}}\right)=\left(\begin{array}{c}{w}_{\mathrm{1}}{\mathbf{\Sigma }}_{\mathrm{1}}^{-\frac{\mathrm{1}}{\mathrm{2}}}\left(\stackrel{\mathrm{̃}}{\mathit{g}}-\mathit{g}\left(\mathit{f}\right)\right)\\ w_{\mathrm{2}}{\mathbf{\Sigma }}_{\mathrm{2}}^{-\frac{\mathrm{1}}{\mathrm{2}}}\left(\stackrel{\mathrm{̃}}{\mathit{d}}-\mathit{d}\left(\mathit{f}\right)\right)\\ {\mathrm{\nabla }}^{\mathrm{2}}\mathit{f}\end{array}\right).\end{array}$$

Then, Eq. (13) is given by

$$\begin{array}{}\text{(15)}& \mathit{\varphi }(\mathit{f};w_{\mathrm{1}},w_{\mathrm{2}})={∥\mathit{r}(\mathit{f};w_{\mathrm{1}},w_{\mathrm{2}})∥}^{\mathrm{2}}.\end{array}$$

Here, we change the variables by letting f=exp (x) to take advantage of the non-negative constraint of f (i.e., f≥0). Then, the minimization of $\mathit{\varphi }(\mathit{x};w_{\mathrm{1}},w_{\mathrm{2}})$ becomes a non-linear least squares problem with respect to x, so we solved it by the Gauss–Newton algorithm.

In the Gauss–Newton method, the parameter x proceeds by iteration, ${\mathit{x}}^{(k+\mathrm{1})}={\mathit{x}}^{\left(k\right)}+\mathrm{\Delta }{\mathit{x}}^{\left(k\right)}$, where the increment Δx^(k) at the kth step is a solution of the following equation:

$$\begin{array}{}\text{(16)}& \left({\mathit{J}}^T\left({\mathit{x}}^{\left(k\right)}\right)\mathit{J}\left({\mathit{x}}^{\left(k\right)}\right)\right)\mathrm{\Delta }{\mathit{x}}^{\left(k\right)}=-{\mathit{J}}^T\left({\mathit{x}}^{\left(k\right)}\right)\mathit{r}\left({\mathit{x}}^{\left(k\right)}\right),\end{array}$$

where J(x) is the Jacobian matrix of r(x) with respect to x. Since Eq. (16) is a normal equation with a large sparse matrix, we solved it by the conjugate gradient (CG) method. The initial values, x⁽⁰⁾, were obtained in advance from only gray level data, $\stackrel{\mathrm{̃}}{\mathit{g}}$, by solving the minimization of ϕ(f) with respect to f. We solved the linear least squares problem by the simultaneous iterative reconstruction technique (SIRT) method (Aso et al., 1998) with f⁽⁰⁾=10⁷ [m⁻² s⁻¹ eV⁻¹] and used the solution f^* for the initial value of the Gauss–Newton algorithm (i.e., ${\mathit{x}}^{\left(\mathrm{0}\right)}=\log \left({\mathit{f}}^{*}\right))$. The hyper-parameters (w₁, w₂) were determined by using the 5-fold cross-validation (Stone, 1974).

The flow of the inverse analysis is summarized as follows.

1. Calculate the initial value of x, x⁽⁰⁾; x⁽⁰⁾ is given by ${\mathit{x}}^{\left(\mathrm{0}\right)}=\log \left({\mathit{f}}^{*}\right)$, where f^* is the solution to minimize ϕ(f), which is solved by the SIRT method. Only auroral images are used for this step (i.e., $\mathit{\varphi }\left(\mathit{f}\right)={∥{\mathbf{\Sigma }}_{\mathrm{1}}^{-\mathrm{1}/\mathrm{2}}\left(\stackrel{\mathrm{̃}}{\mathit{g}}-\mathit{g}\left(\mathit{f}\right)\right)∥}^{\mathrm{2}}$), and the initial value of f is set to 10⁷ [m⁻² s⁻¹ eV⁻¹].

2. Determine the hyper-parameters (w₁, w₂) so as to minimize $\mathit{\varphi }(\mathit{x};w_{\mathrm{1}},w_{\mathrm{2}})$ ($={∥\mathit{r}(\mathit{x};w_{\mathrm{1}},w_{\mathrm{2}})∥}^{\mathrm{2}}$) by the 5-fold cross-validation method. In this step, the values of w₁ and w₂ are selected from pre-created lists, and the same algorithm as shown in step 3 is used to solve the minimization of $\mathit{\varphi }(\mathit{x};w_{\mathrm{1}},w_{\mathrm{2}})$.

3. Solve $min\left[\mathit{\varphi }(\mathit{x};w_{\mathrm{1}},w_{\mathrm{2}})\right]$ with respect to x using w₁ and w₂ determined in the step 2 by the Gauss–Newton algorithm. In the Gauss–Newton algorithm, x proceeds by iteration, ${\mathit{x}}^{(k+\mathrm{1})}={\mathit{x}}^{\left(k\right)}+\mathrm{\Delta }{\mathit{x}}^{\left(k\right)}$, where Δx^(k) is obtained by solving the normal Eq. (16) by the conjugate gradient method. The reconstructed differential flux is obtained by substituting the solution into f=exp (x).

4 Results from the inverse analysis

In this chapter, we show some results reconstructed by the normal ACT and the G-ACT methods. To quantify the performance of these methods, we calculate the mean absolute error (MAE) and the mean absolute percentage error (MAPE). The MAE and MAPE are defined by

$$\begin{array}{}\text{(17)}& \mathrm{MAE}={\displaystyle \frac{\mathrm{1}}{N}}\sum _i^N\left|{\widehat{\mathit{\xi }}}_i-{\mathit{\xi }}_i\right|\end{array}$$

and

$$\begin{array}{}\text{(18)}& \mathrm{MAPE}={\displaystyle \frac{\mathrm{1}}{N}}\sum _i^N\left|{\displaystyle \frac{{\widehat{\mathit{\xi }}}_i-{\mathit{\xi }}_i}{{\mathit{\xi }}_i}}\right|\times \mathrm{100},\end{array}$$

respectively, where ${\widehat{\mathit{\xi }}}_i$ is the reconstruction and ξ_i the true (input) value. The MAE and MAPE are used for different types of data.

The inverse analysis was performed for the reconstruction region shown in Fig. 3a (−40 km < x < 100 km, −70 km < y < 70 km) by using the same spatial and energy grids as those for the forward analysis. Figure 5a shows the precipitating electrons' total energy flux (Q₀). Figure 5b indicates Q₀ as reconstructed only from five auroral images by the ACT method. In this paper, the ACT method uses only ALIS images and solves the minimization of ϕ(x;w₁). The results are displayed for the region of −20 km < x < 80 km and −50 km < y < 50 km. It appears that Q₀ was reconstructed well; however, there are two points to be noted: one is an underestimation of Q₀ at the peak location of each discrete arc and the other is an overestimation between the two arcs. The energy flux at the center of the reconstructed arcs is slightly smaller than the input flux. On the other hand, the energy flux between the two arcs is greater than the input flux, particularly at y<0. For example, Q₀ at (x, y) = (45, −20 km) is 1.47 mW m⁻² for the input flux and 7.30 mW m⁻² for the reconstructed one by ACT.

Figure 5 (a) Horizontal distribution of the incident auroral electrons' total energy flux (Q₀). Points A ([x,y] = [35, −22 km]), B ([x,y] = [23, 16 km]), C ([x,y] = [52, −8 km]), and D ([x,y] = [38, 16 km]) indicate the locations where the reconstructed height profiles of the electron density and energy spectra of the incident electrons are shown in this paper. (b) Q₀ reconstructed from five ALIS images using the ACT method. (c) Q₀ reconstructed by the G-ACT method using five ALIS images and the electron density from 10 × 10 beams of the EISCAT_3D radar. (d) Q₀ reconstructed by the G-ACT method using five ALIS images and the electron density from 21 × 21 beams of the EISCAT_3D radar. (e) Q₀ reconstructed only from the electron density from 21 × 21 beams of the EISCAT_3D radar.

The MAE was used for evaluating the reconstruction of the total energy flux, because the total energy flux includes values close to zero and thus the MAPE was unsuitable for it. The MAE was calculated by using all data in the evaluation area (−20 km < x < 80 km, −50 km < y < 50 km). The MAE values are shown in each panel of Fig. 5. The MAE for the ACT reconstruction was 2.11 mW m⁻².

Figure 5c shows Q₀ as reconstructed by G-ACT using both the ALIS images and the electron density data from the 10 × 10 beams of the EISCAT_3D radar. In this panel, the underestimation of Q₀ at the center of each arc and the overestimation of Q₀ between the two arcs brought by the normal ACT were significantly improved (MAE = 1.87 mW m⁻²). To more clearly show the impact of the electron density data on the improvement, we tested the case that the radar scanned the same area with 21 × 21 beams. Figure 5d presents Q₀ reconstructed by G-ACT using the electron density from the 21 × 21 beams. It is evident that Q₀ reconstructed by G-ACT is more accurate than that reconstructed by ACT (MAE = 1.68 mW m⁻²). The Q₀ value at (x, y) = (45, −20 km) between the two arcs was improved to 4.36 mW m⁻² (2.29 mW m⁻²) by the G-ACT method with the electron density data from the 10 × 10 beams (21 × 21 beams) of the EISCAT_3D radar.

Figure 5e shows Q₀ derived from only the electron density data from the 21 × 21 beams of the EISCAT_3D radar. A larger spatial grid size (${\mathrm{\Delta }}_x={\mathrm{\Delta }}_y=\mathrm{5}$ km, Δ_z=3 km) was used for this inverse analysis. Since the spatial distribution of the electron density data from the radar was much sparser than the grid size for the ACT and G-ACT cases (i.e., Δ_x=1 km, ${\mathrm{\Delta }}_y={\mathrm{\Delta }}_z=\mathrm{2}$ km), a larger grid size was required to collect enough electron density data to solve the inverse problem, even for the 21 × 21 beam scan. The two discrete arcs were roughly reconstructed; however, the horizontal resolution was too low to resolve the fine-scale structure of the arcs.

Figure 6 shows height profiles of the electron density along the field lines at the locations (A, B, C, and D) shown in Fig. 5a. These locations were selected to emphasize the difference in the reconstruction between the normal ACT and G-ACT methods. The black line represents the true profile of the electron density, which was derived by the forward analysis using the incident electrons described in Sect. 2.2. The red squares show the modelled electron density data from the 10 × 10 beams of the EISCAT_3D radar. Several radar data exist along the field line at these locations. It is difficult to estimate the energy distribution of the precipitating electrons as well as the height distribution of the electron density from such few data; this is why the large grid size was used for the inversion from the EISCAT_3D radar data (Fig. 5e). The green crosses and blue circles correspond to the electron density reconstructed by the normal ACT and G-ACT methods, respectively. The spatial distribution of the electron density data obtained from only the optical images by ACT is much denser than those from the EISCAT_3D radar. However, the electron density is smaller than the true values, especially above the height of the peak density, which is consistent with the underestimation of Q₀ (Fig. 5b). The underestimation of the electron density was significantly modified by G-ACT using the EISCAT_3D data. What is most important here is that the electron density can be interpolated at a much higher spatial resolution than that expected from only the EISCAT_3D radar.

The MAPE was used for evaluating the reconstructed electron density because it has a wide scale from 10¹⁰ to 10¹² m⁻³, and the MAE was unsuitable for it. The MAPE values are shown in each panel of Fig. 6. The MAPE values for the electron density reconstructed by G-ACT are smaller than those by ACT at all the locations.

Figure 6 Height profile of the ionospheric electron density at the points (a) A, (b) B, (c) C, and (d) D, which are shown in Fig. 5a. The black curve and solid red squares represent the electron density derived from the incident auroral electrons and modelled observational data from the EISCAT_3D radar (10 × 10 beams), respectively. Green crosses and blue circles show the electron density reconstructed by the ACT and G-ACT methods, respectively.

Figure 7 shows the energy distribution of the differential number flux of precipitating electrons at A, B, C, and D, as reconstructed by the ACT and G-ACT methods. The modelled electron density data from the 10 × 10 beams was used for the G-ACT analysis. This figure indicates that the differential flux reconstructed by the normal ACT tends to be underestimated in the energy range lower than the peak energy (E₀). The underestimation of the differential flux was modified by G-ACT, particularly at the energy corresponding to the altitude where the electron density was obtained from the radar. In the assumed situation, the reconstruction results from G-ACT tends to be better also at the other points (except for A, B, C, and D) than those from ACT. Again, the MAPE was used for evaluating the reconstructed differential flux. Although the MAPE values for the differential flux are greater compared with those for the electron density profiles, the MAPE values for the G-ACT reconstruction are smaller than those for the ACT reconstruction at all the locations.

Figure 7 Energy distribution of the differential number flux of the incident electrons at points (a) A, (b) B, (c) C, and (d) D, which are shown in Fig. 5a. The black curve shows the energy distribution of the incident electrons. Green crosses and blue circles show the differential number flux reconstructed using the ACT and G-ACT methods, respectively.

Figure 8a and b show Q₀ as reconstructed by the ACT and G-ACT methods using data from three ALIS stations (Kiruna, Silkkimuotka, and Tjautjas). All three stations are located to the south of the discrete arcs. Under this condition, it was difficult for ACT to reconstruct the neighboring multiple arcs precisely from the images, because they overlap and cannot be distinguished from each other (MAE = 4.46 mW m⁻²). However, it was demonstrated that the G-ACT method using the electron density from the radar is capable of reconstructing the Q₀ of multiple arcs. The underestimation of Q₀ for both of the two arcs was greatly improved (MAE = 3.19 mW m⁻²), although it still remained because the radar beams were somewhat sparse. Again, we tested the 21 × 21 beam scan case on a trial basis. The reconstructed Q₀ was better improved for both arcs (Fig. 8c; MAE = 1.86 mW m⁻²). Figure 8d and e show the reconstructions made by ACT and G-ACT when using data from two ALIS stations (Kiruna and Silkkimuotka). In this case, ACT was not able to separate the two discrete arcs, and the northern arc disappeared (MAE = 5.87 mW m⁻²). The northern arc was partially reconstructed by the G-ACT method; however, the reconstruction was still difficult in the 10 × 10 beam scan case (MAE = 3.97 mW m⁻²). Even in such a case, if a sufficient number of electron density data were available, the G-ACT method was able to reconstruct the Q₀ of the two arcs very well (Fig. 8f; MAE = 1.80 mW m⁻²).

Figure 8 Upper panels show the total energy flux (Q₀) of the incident electrons reconstructed by using three ALIS images (Kiruna, Silkkimuotka, and Tjautjas). (a) Q₀ reconstructed by the normal ACT, (b) by G-ACT using the electron density data from 10 × 10 beams, and (c) by G-ACT using the electron density data from 21 × 21 beams. Lower panels show Q₀ reconstructed by using two ALIS images (Kiruna and Silkkimuotka). Panels (d), (e), and (f) were obtained by the same method as panels (a), (b), and (c), respectively.

Figure 9 shows the height profiles of the electron density at A and B. The upper and lower panels show the reconstruction results obtained by using three and two ALIS stations, respectively. It can be confirmed that the electron density data obtained from the EISCAT_3D radar were effectively used to improve the reconstructed result by the normal ACT, particularly around the altitude where the radar data exist. Although the electron density reconstructed by ACT with a few images is much lower than the true value, G-ACT enables the underestimation to be corrected. It was confirmed that the MAPE values for the reconstruction results by G-ACT were much smaller than those by ACT at all the cases.

Figure 9 Height profile of the ionospheric electron density. The format of this figure is similar to that of Fig. 6. Panels (a) and (b) indicate the electron density reconstructed by using three ALIS images (Kiruna, Silkkimuotka, and Tjautjas) at A and B, respectively. Panels (c) and (d) show the electron density reconstructed by using two ALIS images (Kiruna and Silkkimuotka) at A and B.

5 Discussion

The EISCAT_3D radar can observe the ionospheric parameters at a much higher spatiotemporal resolution than the existing EISCAT radar. However, if one is interested in the auroral phenomena that have a horizontal scale larger than several tens of kilometers (such as growth-phase arcs, multiple arcs, spirals, westward-traveling surges, and omega bands), the spatial distribution of the ionospheric electron density data obtained by the beam scan of the EISCAT_3D radar may be too sparse to study the fine-scale structures inside the auroras. It is evident that the horizontal spatial resolution is too low to capture both the entire structure and fine-scale structure of the aurora (e.g., Fig. 5e).

The G-ACT method that combines the electron density data with the optical images may enable us to interpolate the electron density data at a much higher spatial resolution than that observed by the EISCAT_3D radar. In particular, this method is effective for the reconstruction of the 3D fine-scale structure of an aurora over a wide horizontal area at high temporal resolution. For instance, this method can provide the fine horizontal structure of the height-integrated ionospheric conductivity of mesoscale (10–1000 km) auroral phenomena at short sampling intervals. This indicates that it is possible to estimate the 3D current system of such auroral phenomena by using the magnetic field measured by a ground-based magnetometer network or the ionospheric electric field from radars (e.g., Vanhamäki and Amm, 2007).

Since the auroral images usually include observational noise, it is often difficult to reconstruct the auroral 3D distribution precisely by using the ACT method. As for the multiple arcs assumed in this study, the total energy flux of the precipitating electrons reconstructed by ACT was underestimated inside the discrete arcs. This is because the two neighboring arcs overlapped when viewed from several imagers and were difficult to perfectly separate. In Figure 5b, the reconstructed electron flux between the arcs was greater than the modelled flux; instead, the flux inside the arcs decreased. When the multiple arcs overlapped from all imagers, it was quite difficult for ACT to distinguish them from each other (Fig. 8). Of course, the reconstruction result depends on the condition such as the relative position of the aurora and the imagers, the noise level, and the shape of the aurora, and different conditions cause the reconstructed electron flux to be overestimated. We demonstrated that G-ACT can significantly reduce the reconstruction errors caused by ACT.

Here, we discuss the timescale of auroral phenomena to which the G-ACT method is applicable. We estimated the integration time required for observing the ionospheric electron density with the EISCAT_3D radar by Eqs. (63), (65), and (66) of Virtanen (2011). It was assumed that the range resolution is 2 km, the full width at half maximum (FWHM) of the beam is 1.4°, the observation frequency is 233 MHz, and the transmitter power is 3.5 MW, which corresponds to the power in the first stage of the EISCAT_3D radar. Then, the integration time needed to achieve the standard deviation of 5 % is less than 0.05 s per beam in the altitude range of 90–200 km when the electron density is greater than 1.0 × 10¹¹ m⁻³. Thus, the 10 × 10 beam scan assumed in this study takes about 5 s.

The interpulse period (IPP) between the pulses of the EISCAT_3D radar is about 5 ms for the E-region ionosphere observation. Thus, if a 16-bit pulse code is used, the minimum temporal resolution becomes 0.08 s per beam, resulting in 8 s for the 10 × 10 beam scan. In practice, the temporal resolution depends on the pulse code and background electron density; therefore, the 10 × 10 beam scan of the electron density in the E-region ionosphere may be made in less than 5 s.

The temporal resolution of the optical imager depends on the performance of the imager, the wavelength of the filter, the auroral emission intensity, etc. Since the monochromatic images are required for the G-ACT analysis, the temporal resolution of high-sensitivity imagers (e.g., electron-multiplying CCD (EMCCD) imagers) with the monochromatic filters is usually a few seconds or less and can be higher than that of the 10 × 10 beam scan of the radar. For example, Fukizawa et al. (2022) reconstructed the 3D distribution of pulsating aurora every 2 s by ACT using the 427.8 nm auroral images. Furthermore, 10 Hz sampling monochromatic all-sky imagers observing 427.8 nm auroral emission have been operative in Tromsø, Norway (Hosokawa et al., 2023).

In addition, the steady state of the electron density was assumed in this study, as given by

$$\begin{array}{}\text{(19)}& {\displaystyle \frac{\partial N_{\mathrm{e}}}{\partial t}}=q-\mathit{\alpha }{N}_{\mathrm{e}}^{\mathrm{2}}\approx \mathrm{0}\end{array}$$

in the E-region ionosphere. Here, q is the ion production rate, N_e is the electron density, and α is the effective recombination coefficient. The steady-state condition is satisfied when the incident electron precipitation does not change over timescales longer than the ion recombination time constant, $\mathit{\tau }=\mathrm{1}/\mathit{\alpha }{N}_{\mathrm{e}}$ (e.g., Semeter and Kamalabadi, 2005). It is well known that α has a large uncertainty (Penman et al., 1979). By using α used by Semeter and Kamalabadi (2005), τ is between 16 and 50 s in the altitude of 90–190 km when the electron density is 1.0 × 10¹¹ m⁻³, and τ decreases as the electron density increases. Thus, the reconstruction results by G-ACT using the current model are valid if the auroral arcs are stable for a longer time than τ. However, it is straightforward to add the time derivative term of the electron density to our model (i.e., $\stackrel{\mathrm{̃}}{\mathit{d}}={\mathit{N}}_e^{\mathrm{2}}+\left(\mathrm{1}/\mathit{\alpha }\right)\left(\partial {\mathit{N}}_e/\partial t\right))$, because this term can be estimated by the continuous observation of the electron density. Such a modified model is available if $\partial {\mathit{N}}_e/\partial t$ is stable during the data acquisition interval. We will examine the modification of the model in the near future.

The mesoscale auroras that we mentioned here have the following typical drift speeds; 70–170 m s⁻¹ for the equatorward drift of the growth-phase arcs (Karlsson et al., 2020), 1–2 km s⁻¹ for the westward-traveling surges (Kamide and Baumjohann, 1993), and 200–800 m s⁻¹ for the eastward drift of the omega bands (Vokhmyanin et al., 2021) at the ionospheric altitude. Pulsating auroras are also mesoscale diffuse aurora, which switch on and off with a quasiperiodic oscillation period of 2–20 s (Lessard, 2012). To study auroral phenomena with relatively fast temporal variations, the number and direction of the beams need to be adjusted. In such situations, simulation studies as shown in this study may be useful in planning observations.

6 Conclusions

We demonstrated via numerical simulation that the combination of optical imagers and the EISCAT_3D radar is very powerful for the study of aurora physics, since they have a complementary relationship with each other. G-ACT, which was used to reconstruct the 3D distribution of auroras (corresponding to the horizontal 2D distribution of the electron energy spectra) from multi-instrument data, can be applied as a technique to take advantage of the optical image data effectively. It has the capability to interpolate the electron density observed by the EISCAT_3D radar at a higher spatial resolution, in particular for mesoscale (10–1000 km) auroral phenomena. G-ACT enables auroral phenomena to be better reconstructed than when the normal ACT is used. Even if ACT cannot reconstruct the auroral distribution precisely, G-ACT may allow us to reduce the reconstruction error. Therefore, it is important to construct multi-point monochromatic imager networks that cover the observation region of the EISCAT_3D radar in the near future.