Observing electric field and neutral wind with EISCAT 3D

Measurements of height-dependent electric field (E) and neutral wind (u) are important governing parameters of the Earth's upper atmosphere, which can be used to study, for example, how auroral currents close or how energy flows between the ionized and neutral constituents. The new EISCAT 3D (E3D) incoherent scatter radar will be able to measure a three-dimensional ion velocity vector (v) at each measurement point, which will allow less stringent prior assumptions about E and u to be made when estimating them from radar measurements. This study investigates the feasibility of estimating the three-dimensional electric field and neutral wind vectors along a magnetic field-aligned profile from E3D measurements, using the ion momentum equation and Maxwell's equations. The uncertainty of ion drift measurements is estimated for a time and height resolution of 5 s and 2 km. With the most favourable ionospheric conditions, the ion wind at E region peak can be measured with an accuracy of less than 1 m/s. In the worst case, during a geomagnetically quiet night, the uncertainty increases by a factor of around 10. The uncertainty of neutral wind and electric field estimates is found to be strongly dependent on the prior constraints imposed on them. In the lower E region, neutral wind estimates have a lower standard deviation than 10 m/s in the most favourable conditions. In such conditions, also the F region electric field can be estimated with uncertainty of about 1 mV/m. Simulated measurements of v are used to demonstrate the ability to resolve the field-aligned profile of E and u. However, they can only be determined well at the heights where they dominate the ion drift, that is above 125 km for E and below 115 km for u. At the other heights, the results are strongly dependent on the prior assumptions of smoothness.

1 Introduction

One of the main parameters that incoherent scatter radars (ISRs) measure is ion drift velocity v. This can be related primarily to electric field E and neutral wind u, making it possible to use ISR measurements for estimating these parameters, which are of interest when, for example, studying the electrodynamics of the aurorae borealis (e.g. Takahashi et al., 2019), determining how auroral currents close within the ionosphere, and studying how energy is transferred between the ionosphere and the neutral atmosphere (e.g. Aikio and Selkälä, 2009; Kosch et al., 2011; Cai et al., 2016).

The method for simultaneously estimating E and u for the auroral ionosphere using an incoherent scatter radar was first described by Brekke et al. (1973). Since then, this technique has been used and improved (see Nygrén et al., 2011, and references therein). The velocity of both the ion and neutral wind is related to other ionospheric parameters, such as ion-neutral collision frequency and electric field through the momentum equation of the ions. Some parameters can be measured, while others need to be taken from models. The terms with pressure gradients and gravitation are commonly neglected. The electric field can be deduced from measurements higher up where ion-neutral collisions are negligible and used further down by assuming that electric field along the magnetic field line is constant (Brekke et al., 1994; Heinselman and Nicolls, 2008).

Currently ion drifts are measured with a monostatic radar by pointing the transmit beam in three or more different directions and measuring the projection of the ion-velocity vector onto these directions. If the ion velocity then is assumed constant or slowly changing in the horizontal direction for all of the pointing directions, an estimate of the ion drift can be made (Heinselman and Nicolls, 2008; Nicolls et al., 2014b). If the observations are made with slowly moving dish-based radars, making the measurements in different pointing directions might take several minutes (Williams et al., 1984). With modern phased array systems, the scanning time can be reduced (Heinselman and Nicolls, 2008).

Multistatic radars can contribute to the measurement challenge by simultaneously measuring a common scattering volume from multiple different directions. Provided there are at least three linearly independent ion velocity measurements, the full vector can be determined. With dish antennas, the ion velocity is obtained only from one intersection volume at a time (Williams et al., 1984; Risbeth and Williams, 1985). In order to obtain measurements along a line or over a volume of space, the antennas need to be steered, which again takes time. In 2008, only the EISCAT UHF system could offer multistatic measurements (Heinselman and Nicolls, 2008). Later, the receiver antennas in Kiruna and Sodankylä were converted to receive the signals from EISCAT VHF instead (Kero, 2014; Mann et al., 2016). Earlier, there were also other ISRs that were multistatic (Williams et al., 1984).

One of the capabilities of the upcoming EISCAT3D (E3D) incoherent scatter radar is that it can simultaneously observe scatter from ionospheric plasma from at least three different geographically separated receiver sites, each using multiple simultaneous receiver beams that intersect the transmit beam at multiple heights; see Fig. 1. This is made possible by the use of phased array technology (e.g. Wirth, 2001), which allows for fast beam steering and allows the receivers to form multiple simultaneous beams. A measurement of the ion velocity vector along the radar transmit beam is then possible without any other assumptions than homogeneity of ion velocity within the common scattering volumes where the transmit and receive beams intersect (McCrea et al., 2015; Virtanen et al., 2014).

Figure 1 Locations of the initial three E3D sites in Skibotn, Karesuvanto, and Kaiseniemi (a) and geometry for measuring ion velocity vectors (b). All velocity measurements shown are carried out with one transmit beam. For this to work, the receivers are able to form all the needed receiver beams simultaneously.

In this article, we study the ion velocity measurement capability of EISCAT3D. Since ion velocity depends on both electric field, through the Lorentz force, and neutral wind, through collisions, the ion velocity contains information on both electric field and neutral wind. At each range along the transmitted beam, we therefore have six unknowns, three electric field and three neutral wind components, but only three observed ion velocity components. This makes the problem slightly underdetermined. The traditional approach to handle this problem has been to assume a constant electric field along the magnetic field line. Here, we relax this assumption and investigate the physical nature of the problem, showing that physics-based constraints give us additional equations similar to first-order Tikhonov regularization. The ion velocity estimations and the corresponding uncertainty calculations are described in Sect. 2. The inverse problem of determining electric field and neutral wind from the ion velocities is described in Sect. 3. In Sect. 4, model calculations are used to illustrate the resolution and accuracy that are possible.

2 Ion wind

ISR measurements mainly provide four ionospheric parameters: electron density n_e, electron temperature T_e, ion temperature T_i, and the ion velocity component along the Bragg scattering vector w_p. The connection between w_p and the Bragg scattering vector k_p for the transmit–receive pair p and the ion velocity v is

$$\begin{array}{}\text{(1)}& w_p={\mathit{k}}_p\cdot \mathit{v}/\left|{\mathit{k}}_p\right|+{\mathit{\epsilon }}_p,\end{array}$$

where ε_p is a random variable that models the velocity measurement errors.

The basis for finding the velocity of the uncharged wind and electric field is through the ion velocity. The measurements w_p and the unknown velocity v can be set up as a linear inverse problem (Heinselman and Nicolls, 2008; Nygrén et al., 2011):

$$\begin{array}{}\text{(2)}& \mathit{w}=\mathbf{K}\mathit{v}+\mathit{\epsilon },\end{array}$$

where ${\mathit{w}}^{\top }=[w_{\mathrm{1}},\mathrm{\dots }{w}_P]$ is a vector containing independent measurements of ion-line Doppler shift, ${\mathbf{K}}^{\top }=\left[{\mathit{k}}_{\mathrm{1}}\mathrm{\dots }{\mathit{k}}_P\right]$ is the theory matrix, and ε is the noise vector. We assume that the noise is independent and identically normal distributed with zero mean and variance of ${\mathit{\sigma }}_w^{\mathrm{2}}$, which means that we assume that the line-of-sight (LOS) ion velocity measured with different receivers has the same uncertainty.

With the first stage of E3D, there will be P=3 measurements of the ion velocity at every range, as shown in Fig. 1. Since the theory matrix then is quadratic, it will be possible to find v with low uncertainty and without restrictions if the measurements are sufficiently independent linearly (see Aster et al., 2013; Risbeth and Williams, 1985). The Bragg scattering vectors can be calculated from the preliminary positions of E3D as mentioned by Kero et al. (2019). We assume a target in the direction of the magnetic field at ionospheric range, extending outwards from the Skibotn transmitter site.

As it later will become an advantage to have the ion velocity in magnetic field coordinates, we have to transform the scattering vector matrix K. The transformation matrix from geographic to local magnetic coordinates is

$$\begin{array}{}\text{(3)}& {\mathbf{R}}_{\mathrm{geo}\to \mathrm{gmag}}=\left[\begin{array}{lll}\cos \mathit{\delta }& -\sin \mathit{\delta }& \mathrm{0}\\ \sin I\sin \mathit{\delta }& \cos \mathit{\delta }\sin I& \cos I\\ -\cos I\sin \mathit{\delta }& -\cos I\cos \mathit{\delta }& \sin I\end{array}\right],\end{array}$$

where δ is the declination, and I is the dip angle of the magnetic field (Heinselman and Nicolls, 2008). We note that the transform is a rotation, meaning that it is orthonormal, and ${\mathbf{R}}_{\mathrm{geo}\to \mathrm{gmag}}{\mathbf{R}}_{\mathrm{geo}\to \mathrm{gmag}}^{\top }$ is equal to the identity matrix. The LOS velocities in geographic coordinates are then related to the ion velocity in local magnetic coordinates as follows:

$$\begin{array}{}\text{(4)}& \mathit{w}={\mathbf{KR}}_{\mathrm{geo}\to \mathrm{gmag}}^{\top }\mathit{v}+\mathit{\epsilon }.\end{array}$$

We find the ion velocity by solving the inverse problem by using the linear least-squares method:

$$\begin{array}{}\text{(5)}& \widehat{\mathit{v}}={\mathbf{R}}_{\mathrm{geo}\to \mathrm{gmag}}{\left({\mathbf{K}}^{\top }\mathbf{K}\right)}^{-\mathrm{1}}{\mathbf{K}}^{\top }\mathit{w}.\end{array}$$

The uncertainty of the ion velocity estimate is quantified using the following covariance matrix:

$$\begin{array}{}\text{(6)}& {\mathbf{\Sigma }}_v={\mathbf{R}}_{\mathrm{geo}\to \mathrm{gmag}}{\left({\mathbf{K}}^{\top }\mathbf{K}\right)}^{-\mathrm{1}}{\mathbf{R}}_{\mathrm{geo}\to \mathrm{gmag}}^{\top }{\mathit{\sigma }}_w^{\mathrm{2}}.\end{array}$$

We will use this uncertainty later when estimating electric field and neutral wind.

2.1 Uncertainty of ion wind velocity

Vallinkoski (1989) describes a method for finding the uncertainty of ISR parameter estimates. Our procedure is similar. Like the other ionospheric parameters, the ion velocity component along the Bragg scattering vector w_p is estimated from the autocorrelation function (ACF) that is measured by the radar. The ACF $\mathit{\rho }={\left[\mathit{\rho }\left({\mathit{\tau }}_{\mathrm{0}}\right),\mathit{\rho }\left({\mathit{\tau }}_{\mathrm{1}}\right),\mathrm{\dots },\mathit{\rho }\left({\mathit{\tau }}_{T-\mathrm{1}}\right)\right]}^{\top }$ is described by the theory for incoherent scatter (Kudeki and Milla, 2011). Here, τ is the time lag. The theory provides a non-linear relationship between the parameters θ and the ACF ρ.

$$\begin{array}{}\text{(7)}& \mathit{\rho }=\mathit{f}\left(\mathit{\theta }\right)+\mathit{\epsilon }\end{array}$$

The parameters θ include ionospheric plasma parameters and parameters specific to the radar experiment. In this relationship, there are also measurement errors, which are modelled with a random variable ε.

To simplify the uncertainty calculations, we linearize the relationship between ACF and the parameters measured with the ISR. The first-order Taylor polynomial for the ACF around a parameter estimate θ^′ is

$$\begin{array}{}\text{(8)}& \mathit{\rho }=\mathit{f}\left({\mathit{\theta }}^{\mathbf{\prime }}\right)+\mathbf{J}(\mathit{\theta }-{\mathit{\theta }}^{\mathbf{\prime }})+\mathit{\epsilon },\end{array}$$

where

$$\begin{array}{}\text{(9)}& \mathbf{J}={\left[\begin{array}{lll}\frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{0}}}{\mathrm{d}{\mathit{\theta }}_{\mathrm{1}}}& \mathrm{\cdots }& \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{0}}}{\mathrm{d}{\mathit{\theta }}_n}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \frac{\mathrm{d}{\mathit{\rho }}_{T-\mathrm{1}}}{\mathrm{d}{\mathit{\theta }}_{\mathrm{1}}}& \mathrm{\cdots }& \frac{\mathrm{d}{\mathit{\rho }}_{T-\mathrm{1}}}{\mathrm{d}{\mathit{\theta }}_n}\end{array}\right]}_{\mathit{\theta }={\mathit{\theta }}^{\mathbf{\prime }}}\end{array}$$

is the Jacobian of f evaluated at $\mathit{\theta }={\mathit{\theta }}^{\mathbf{\prime }}$. We move the constant parameters over to the left side and get an inverse problem with the solution

$$\begin{array}{}\text{(10)}& \widehat{\mathit{\theta }}={\left({\mathbf{J}}^{\mathrm{H}}{\mathbf{\Sigma }}_{\mathit{\rho }}^{-\mathrm{1}}\mathbf{J}\right)}^{-\mathrm{1}}{\mathbf{J}}^{\mathrm{H}}{\mathbf{\Sigma }}_{\mathit{\rho }}^{-\mathrm{1}}\left(\mathit{\rho }-\mathit{f}\left({\mathit{\theta }}^{\mathbf{\prime }}\right)+\mathbf{J}{\mathit{\theta }}^{\mathbf{\prime }}\right),\end{array}$$

where the superscript H denotes Hermitian transpose. The linearized covariance matrix quantifying the uncertainty of the estimate is

$$\begin{array}{}\text{(11)}& {\mathbf{\Sigma }}_{\widehat{\mathit{\theta }}}={\left({\mathbf{J}}^{\mathrm{H}}{\mathbf{\Sigma }}_{\mathit{\rho }}^{-\mathrm{1}}\mathbf{J}\right)}^{-\mathrm{1}}.\end{array}$$

The uncertainty is dependent on how well the ACF is measured. This depends on the signal-to-noise ratio. The signal strength is dependent on the ionospheric plasma parameters as mentioned above, the radar equation, and the experiment design with pulse length and coding, etc. The noise level is determined by the system noise temperature, which depends on the implementation of the receiver electronics and the sky noise temperature at the radar frequency.

To determine the ACF, we calculate the ISR spectrum as described by Kudeki and Milla (2011) and take its inverse Fourier transform. We then multiply it with the signal strength, which we take from the radar equation

$$\begin{array}{}\text{(12)}& P_{\mathrm{S}}={\displaystyle \frac{P_{\mathrm{t}}{G}_{\mathrm{t}}{G}_{\mathrm{r}}{\mathit{\lambda }}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathit{\chi }}{(\mathrm{4}\mathit{\pi }{)}^{\mathrm{3}}{R}_{\mathrm{t}}^{\mathrm{2}}{R}_{\mathrm{r}}^{\mathrm{2}}}}\cdot \mathcal{V}{n}_{\mathrm{e}}\cdot {\displaystyle \frac{\left(\mathrm{4}\mathit{\pi }{r}_{\mathrm{e}}^{\mathrm{2}}\right)}{\mathrm{1}+T_{\mathrm{e}}/T_{\mathrm{i}}}},\end{array}$$

where P_t is transmit power, G_t is transmit gain, G_r is receive gain, λ is the radar wavelength, R_t and R_r are the distance between the target and transmitter and receiver, χ is the polarization angle, and r_e is the classical electron radius. The scattering volume 𝒱 is approximated as a spherical sector

$$\begin{array}{}\text{(13)}& \mathcal{V}={\displaystyle \frac{\mathrm{2}\mathit{\pi }\mathrm{\Delta }r}{\mathrm{3}}}\left(\mathrm{1}-\cos {\displaystyle \frac{\mathit{\phi }}{\mathrm{2}}}\right)\left({\displaystyle \frac{\mathrm{\Delta }{r}^{\mathrm{2}}}{\mathrm{4}}}+\mathrm{3}{R}_{\mathrm{t}}^{\mathrm{2}}\right),\end{array}$$

where Δr is the resolution in range direction, and φ is the one-way half-power beamwidth of the radar. For bistatic cases when receiving in Karesuvanto or Kaiseniemi, we assume that the receiver sees the whole scattering volume such that we do not need to include possible losses because the radar beams do not overlap completely. The noise power P_N is obtained using the Nyquist–Johnson noise model as follows:

$$\begin{array}{}\text{(14)}& P_{\mathrm{N}}=k_{\mathrm{B}}{T}_{\mathrm{sys}}{f}_{\mathrm{B}},\end{array}$$

where k_B is the Boltzmann constant, T_sys is the system noise temperature, and f_B is the bandwidth of the signal, which is given by $f_{\mathrm{B}}=c/\left(\mathrm{2}\mathrm{\Delta }r\right)$, where c is the speed of light.

The calculations require that measurements of the different lags of the ACF are uncorrelated. This means that the covariance matrix of the measurement errors Σ_ρ is diagonal. We can make this assumption if the different lags are measured using a coded long pulse with a low signal-to-noise ratio (Lehtinen and Häggström, 1987). A solution for cases with signal-to-noise ratio over 1 can be to shorten the baud length and so increase the resolution in range direction (Lehtinen and Damtie, 2013). This results in a weaker signal from every range but provides more independent measurements that can be averaged to obtain the desired range resolution.

Figure 2 Ionospheric parameter profiles we used to calculate the ion velocity errors.

We can use this outline to calculate the uncertainty in w_p for several representative cases. For the radar parameters of E3D, we use frequency f=233 MHz, one-way half-power beamwidth 2^∘, both transmit and receive gain equal to 38 dB, transmit power 5 MW, and a noise temperature of 200 K. These are, to the best of our knowledge, the performance parameters of the latest revision of the EISCAT 3D design, which may still change before the final implementation. We use a scattering angle of 90^∘, even if it is not absolutely correct when receiving in Karesuvanto and Kaiseniemi.

In order to investigate the performance of the radar in different geophysical conditions, we have studied three different cases: (1) daytime; (2) night-time without auroral precipitation, as modelled by IRI (quiet night-time); and (3) night-time with auroral precipitation. Each of these cases have different ionospheric plasma parameter profiles consisting of n_e, T_e, T_i, and m_i. The key parameter that influences observability is n_e, as the signal-to-noise ratio is to first order proportional to this parameter. The temperatures and the ion mass also influence the uncertainty somewhat but much weaker than the electron density. For the representative cases, we used the plasma parameters for 20 February 2014 at three times: 14:00, 23:00, and 21:20 UTC. The profiles are calculated by the IRI-2016 model (Bilitza et al., 2017), except for the aurora case at 21:20, where we used data from EISCAT UHF for electron density and the temperatures. We integrated the EISCAT data over 10 min in order to obtain plasma parameters with smooth profiles. For calculating the magnetic field, we use the international geomagnetic reference field (see Thébault et al., 2015). The ionospheric parameter profiles for the three representative cases are shown in Fig. 2. In this article, these profiles are only used for calculating the uncertainty in the ion velocity measurements.

For the analysis, we assumed an experiment where the baud length is 15 µs, the pulse consists of 51 bauds, and there is an interpulse period of 5 ms and an integration time of 5 s. We use an analysis range resolution of 2250 m, corresponding to the baud length. This range resolution we see as a compromise for an experiment investigating both the E and F region. For simplicity, we kept the range resolution constant.

Figure 3 1σ uncertainty in line-of-sight velocity (a), and the hence following three components of the ion velocity vector in magnetic field coordinates: perpendicular east (b), perpendicular north (c), and field-aligned direction (d).

The uncertainty of the LOS ion velocities is shown in Fig. 3. According to the figure, the uncertainty at daytime and auroral night-time is considerably lower than at night-time without aurora. While the uncertainty varies from about 5 m/s at 100 km to 20 m/s at 140 km altitude in the non-aurora night case, for daytime and auroral night-time conditions, the uncertainty is smaller than 3 m/s. In general, the uncertainty is smaller where the signal-to-noise ratio is high. This occurs primarily at E region heights, where the electron density is comparatively high. At F region heights, the electron density is also high. However, this is about twice as far as the E region, and the backscattered signal is therefore weaker.

It is worth noting that the test case is close to a solar maximum, which means that the electron density is comparatively high. At solar minimum, the electron density in the ionosphere is in general about a factor of 2 lower (e.g. Brekke, 2013), and the uncertainty in ion velocity will be higher. One can compensate for this by integrating the LOS ion velocity over a larger number of range gates, leading to a reduced range resolution. For example, Nygrén et al. (2011) used 10 km range resolution at E region heights in an experiment with EISCAT UHF.

Using Eq. (6), we obtain the uncertainties of the ion velocity components, which are plotted in Fig. 3b–d. The uncertainty in magnetic field-aligned component is very similar to the LOS uncertainties. This is expected because all lines of sight do not differ much from the magnetic field line direction. Therefore, the uncertainty of the ion velocity components perpendicular to the magnetic field line is a factor of 3–5 times higher. At the highest altitudes, the scattering vectors are even more similar, which leads to an increased uncertainty.

3 Neutral wind and electric field

The velocity of ion and neutral wind are coupled through collisions as described by the ion momentum equation. This can be found by taking the first moment of the Boltzmann equation (e.g Inan and Gołkowski, 2011). We assume that we can treat the ions as a single fluid. The momentum equation is

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{n}_{\mathrm{i}}{m}_{\mathrm{i}}\left[{\displaystyle \frac{\mathrm{d}\mathit{v}}{\mathrm{d}t}}+\left(\mathit{v}\cdot \mathrm{\nabla }\right)\mathit{v}\right]& =-\mathrm{\nabla }{\mathbf{P}}_{\mathrm{i}}+n_{\mathrm{i}}{m}_{\mathrm{i}}\mathit{g}+\\ & q_{\mathrm{i}}{n}_{\mathrm{i}}\left(\mathit{E}+\mathit{v}\times \mathit{B}\right)-\\ & \sum _k{n}_{\mathrm{i}}{m}_{\mathrm{i}}{\mathit{\nu }}_{\mathrm{i}k}\left(\mathit{v}-{\mathit{v}}_k\right),\end{array}\end{array}$$

where n_i is the number density of ions, m_i is the ion mass, P_i is the ion pressure tensor, g is the gravitational acceleration, q_i is the ion charge, E is the electric field, B is the background magnetic field, ν_ik is the momentum transfer collision frequency between ions and particle species k, and v_k is the velocity of particle species k. We assume that spatial variations of the ion velocity are small such that we can neglect the term (v⋅∇)v. Further we assume that the pressure is isotropic, so we can write the pressure tensor as a scalar p_i. Only collisions between ions and neutrals are of importance to change the ion velocity (Brekke, 2013); other collision terms can be neglected. If the ions obey the ideal gas law, the ion pressure p_i can be written as p_i=n_ik_BT_i. Additionally, we neglect local temperature variations such that ∇T_i=0. Finally, as in previous work, we also neglect the contribution from pressure gradients and gravity. With all these assumptions, Eq. (15) can be rewritten as

$$\begin{array}{}\text{(16)}& n_{\mathrm{i}}{m}_{\mathrm{i}}{\displaystyle \frac{\mathrm{d}\mathit{v}}{\mathrm{d}t}}=q_{\mathrm{i}}{n}_{\mathrm{i}}\left(\mathit{E}+\mathit{v}\times \mathit{B}\right)-n_{\mathrm{i}}{m}_{\mathrm{i}}{\mathit{\nu }}_{\mathrm{in}}\left(\mathit{v}-\mathit{u}\right),\end{array}$$

where u is the neutral wind velocity.

For steady-state conditions $\left(\frac{\mathrm{d}\mathit{v}}{\mathrm{d}t}=\mathrm{0}\right)$, the ion velocity in the magnetic field coordinate system becomes (see Brekke, 2013; Heinselman and Nicolls, 2008)

$$\begin{array}{}\text{(17a)}& {\displaystyle }{v}_x=u_x+{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+{\mathit{\kappa }}_{\mathrm{i}}^{\mathrm{2}}}}\left[{\displaystyle \frac{{\mathit{\kappa }}_i}{B}}{E}_x-{\mathit{\kappa }}_{\mathrm{i}}\left(u_y+{\displaystyle \frac{{\mathit{\kappa }}_i}{B}}{E}_y\right)-{\mathit{\kappa }}_{\mathrm{i}}^{\mathrm{2}}\left(u_x\right)\right]\text{(17b)}& {\displaystyle }{v}_y=u_y+{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+{\mathit{\kappa }}_{\mathrm{i}}^{\mathrm{2}}}}\left[{\displaystyle \frac{{\mathit{\kappa }}_i}{B}}{E}_x+{\mathit{\kappa }}_{\mathrm{i}}\left(u_x+{\displaystyle \frac{{\mathit{\kappa }}_i}{B}}{E}_x\right)-{\mathit{\kappa }}_{\mathrm{i}}^{\mathrm{2}}\left(u_y\right)\right]\text{(17c)}& {\displaystyle }{v}_z=u_z+{\displaystyle \frac{{\mathit{\kappa }}_i}{B}}{E}_z.\end{array}$$

Here, κ_i is the ion mobility

$$\begin{array}{}\text{(18)}& {\mathit{\kappa }}_i={\displaystyle \frac{q_{i}B}{m_i{\mathit{\nu }}_{in}}},\end{array}$$

where the subscript z denotes the direction along (antiparallel to) the magnetic field, x horizontally towards east, and y perpendicular to the other two directions, giving a right-handed system. Since in this article we are only considering the ion mobility, we will drop the subscript i from now on and just write κ instead.

The component equations can be combined into a compact matrix equation (see Heinselman and Nicolls, 2008).

$$\begin{array}{}\text{(19)}& \mathit{v}={\displaystyle \frac{\mathit{\kappa }}{B}}\mathbf{C}\mathit{E}+\mathbf{C}\mathit{u},\end{array}$$

where C is the matrix

$$\begin{array}{}\text{(20)}& \mathbf{C}=\left[\begin{array}{ccc}\frac{\mathrm{1}}{\mathrm{1}+{\mathit{\kappa }}^{\mathrm{2}}}& \frac{\mathit{\kappa }}{\mathrm{1}+{\mathit{\kappa }}^{\mathrm{2}}}& \mathrm{0}\\ \frac{-\mathit{\kappa }}{\mathrm{1}+{\mathit{\kappa }}^{\mathrm{2}}}& \frac{\mathrm{1}}{\mathrm{1}+{\mathit{\kappa }}^{\mathrm{2}}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right].\end{array}$$

When estimating the neutral wind and electric field at a certain altitude, Eq. (19) has to be solved. This is an underdetermined inverse problem with six unknowns, which are all components of both the electric field and the neutral wind velocity. For measurements, we only have the three components of the ion velocity. To resolve this, some a priori assumptions or constraints are required.

The original solution of Brekke et al. (1973) was to use the fact that κ≫1 at F region altitudes; therefore the ion drift is determined only by the electric field. Then this is assumed to be constant along the magnetic field line. However, the electric field may not be constant in reality (e.g. Sangalli et al., 2009). Such an assumption then affects the neutral wind estimates.

It is possible to assume that the neutral wind and electric field vary smoothly in the whole range of interest and use the full profile of all ion wind measurements to obtain estimates of the neutral wind and electric field. We will outline a procedure to specify a smoothness constraint based on Maxwell's equations in order to give a physically feasible solution.

We start by discretizing the problem as follows: we have a set of ion wind velocity vectors v_i…v_H which are measurements of Eq. (19) integrated over a height range defined by the weighting functions d_i(h)

$$\begin{array}{}\text{(21)}& {\mathit{v}}_i=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\left[{\displaystyle \frac{\mathit{\kappa }}{B}}\mathbf{C}\mathit{E}+\mathbf{C}\mathit{u}\right]d_i\left(h\right)\mathrm{d}h+{\mathit{\epsilon }}_i,\end{array}$$

where ε_i is the noise in measurement i and assumed to be normally distributed with zero mean and covariance described by Eq. (6) (in magnetic field coordinates). We assume that the unknowns can be described by a set of basis functions

$$\begin{array}{}\text{(22)}& \mathit{E}\left(h\right)=\sum _{j=\mathrm{1}}^{N_E}\mathit{\kappa }\left(h\right){\mathit{\beta }}_j{\mathit{b}}_j\left(h\right)/B\left(h\right)\end{array}$$

and

$$\begin{array}{}\text{(23)}& \mathit{u}\left(h\right)=\sum _{j=N_E+\mathrm{1}}^{N_E+N_u}{\mathit{\beta }}_j{\mathit{b}}_j\left(h\right).\end{array}$$

This allows Eq. (21) to be written as

$$\begin{array}{}\text{(24)}& {\mathit{v}}_i=\sum _{j=\mathrm{1}}^N{\mathit{a}}_{ij}\left(h\right){\mathit{\beta }}_j+{\mathit{\epsilon }}_i,\end{array}$$

where

$$\begin{array}{}\text{(25)}& {\mathit{a}}_{ij}=\left\{\begin{array}{lr}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\frac{\mathit{\kappa }\left(h\right)\mathbf{C}\left(h\right){\mathit{b}}_j\left(h\right)d_i\left(h\right)}{B\left(h\right)}\mathrm{d}h& \mathrm{for}\mathrm{1}\le i\le N_E\\ \underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\mathbf{C}\left(h\right){\mathit{b}}_{\mathit{j}}\left(h\right)d_i\left(h\right)\mathrm{d}h& \mathrm{for}{N}_E<i\le N\end{array}\right.,\end{array}$$

which can be calculated before solving the problem and therefore can be regarded as constants.

We assume that the weighting functions d_i(h) for the ion wind measurements are boxcars with centre at a certain height and extending exactly halfway to the centre of the nearest box in both directions. At the ends, the measurement height boxes are symmetric around their centre. The basis functions for the unknowns b_j(h) are also boxcars. Other basis functions could also be used. We further assume that κ, B, and the rotation matrix C are constant throughout our measurement height boxes d(h).

Equation (24) in matrix form then can be written as follows:

$$\begin{array}{}\text{(26)}& \mathit{V}=\mathbf{A}\mathit{x}+\mathit{\xi },\end{array}$$

where ${\mathit{V}}^{\top }=[{\mathit{v}}_{\mathrm{1}}^{\top },\mathrm{\dots },{\mathit{v}}_H^{\top }]$, $\mathbf{A}=\left[\begin{array}{lll}{\mathit{a}}_{\mathrm{1},\mathrm{1}}& \mathrm{\cdots }& {\mathit{a}}_{\mathrm{1},N}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ {\mathit{a}}_{H,\mathrm{1}}& \mathrm{\cdots }& {\mathit{a}}_{H,N}\end{array}\right]$, ${\mathit{x}}^{\top }=[{\mathit{\beta }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\beta }}_N]$, and ${\mathit{\xi }}^{\top }=[{\mathit{\epsilon }}_{\mathrm{1}}^{\top },\mathrm{\dots },{\mathit{\epsilon }}_H^{\top }]$.

In order to regularize the problem, we use Gauss' and Faraday's laws. Faraday's law for a static magnetic field, $\mathbf{\nabla }\times \mathit{E}=\mathbf{0}$, gives us three equations for the gradient of the electric field:

$$\begin{array}{}\text{(27a)}& {\displaystyle \frac{\mathrm{d}{E}_y}{\mathrm{d}z}}-{\displaystyle \frac{\mathrm{d}{E}_z}{\mathrm{d}y}}=\mathrm{0}\text{(27b)}& {\displaystyle \frac{\mathrm{d}{E}_x}{\mathrm{d}z}}-{\displaystyle \frac{\mathrm{d}{E}_z}{\mathrm{d}x}}=\mathrm{0}\text{(27c)}& {\displaystyle \frac{\mathrm{d}{E}_y}{\mathrm{d}x}}-{\displaystyle \frac{\mathrm{d}{E}_x}{\mathrm{d}y}}=\mathrm{0}.\end{array}$$

Gauss' law for a charge-neutral plasma, $\mathbf{\nabla }\cdot \mathit{E}=\mathrm{0}$, can be written as

$$\begin{array}{}\text{(28)}& {\displaystyle \frac{\mathrm{d}{E}_x}{\mathrm{d}x}}+{\displaystyle \frac{\mathrm{d}{E}_y}{\mathrm{d}y}}+{\displaystyle \frac{\mathrm{d}{E}_z}{\mathrm{d}z}}=\mathrm{0}.\end{array}$$

Equations (27a), (27b), and (28) are added to the theory matrix A to regularize the electric field. The derivatives $\frac{\mathrm{d}{E}_x}{\mathrm{d}z}$, $\frac{\mathrm{d}{E}_y}{\mathrm{d}z}$, and $\frac{\mathrm{d}{E}_z}{\mathrm{d}z}$ are approximated with finite differences, with dz equal to the range step.

The horizontal gradients $\left(\frac{\mathrm{d}{E}_x}{\mathrm{d}x}+\frac{\mathrm{d}{E}_y}{\mathrm{d}y},\frac{\mathrm{d}{E}_z}{\mathrm{d}x},\frac{\mathrm{d}{E}_z}{\mathrm{d}y}\right)$ are not specified by our measurements. We therefore treat them as Gaussian random variables ${\mathit{\xi }}_{j,(x,y,z)}$ with zero mean and some variance ${\mathit{\alpha }}_{j,(x,y,z)}^{-\mathrm{2}}$. Eqs. (27a), (27b), and (28) then result in

$$\begin{array}{}\text{(29)}& {\displaystyle \frac{E_{j,(x,y,z)}-E_{j+\mathrm{1},(x,y,z)}}{\mathrm{\Delta }{h}_E}}={\mathit{\xi }}_{j,(x,y,z)}.\end{array}$$

These equations are added to the theory matrix. This implies that we assume these three derivatives of the electric field to be smaller than $\mathrm{2}/{\mathit{\alpha }}_{j,(x,y,z)}$ 95 % of the time. For a box size of Δh_E, this means that $(E_{j,(x,y,z)}-E_{j+\mathrm{1},(x,y,z)})\sim \mathcal{N}(\mathrm{0},\mathrm{\Delta }{h}_E^{\mathrm{2}}{\mathit{\alpha }}_{j,(x,y,z)}^{-\mathrm{2}})$, which is similar to first-order Tikhonov regularization but with a regularization constant ${\mathit{\alpha }}_{j,(x,y,z)}$ that varies with both height and electric field component (see, for example, Roininen et al., 2011). It is worth pointing out that the constraints are obtained from Maxwell's equations and therefore have a physical interpretation.

Constraints, such as Eq. (29), will favour smoother solutions that are closer to being constant-valued (Aster et al., 2013). Throughout this paper, we will loosely use “flatness” to describe how close a function is to a constant value, as the magnitude of the left-hand side of Eq. (29) is minimized when the function is constant.

For the neutral wind, we also use first-order Tikhonov regularization as described for the electric field above. In addition, we use zeroth-order Tikhonov regularization to restrict the neutral wind to smaller magnitudes. This corresponds to the following statistical assumptions

$$\begin{array}{}\text{(30)}& \begin{array}{ll}{u}_{j,(x,y,z)}-u_{j+\mathrm{1},(x,y,z)}& ={\mathit{\zeta }}_{\mathrm{1},j,(x,y,z)}\\ u_{j,(x,y,z)}& ={\mathit{\zeta }}_{\mathrm{0},j,(x,y,z)},\end{array}\end{array}$$

where ${\mathit{\zeta }}_{\mathrm{1},j,(x,y,z)}\sim \mathcal{N}(\mathrm{0},\mathrm{\Delta }{h}_u^{\mathrm{2}}{\mathit{\gamma }}_{\mathrm{1},j,(x,y,z)}^{-\mathrm{2}})$ and ${\mathit{\zeta }}_{\mathrm{0},j,(x,y,z)}\sim \mathcal{N}(\mathrm{0},{\mathit{\gamma }}_{\mathrm{0},j,(x,y,z)}^{-\mathrm{2}})$. The first row regularizes to the flatness of the profile, and the second constrains the magnitude.

This procedure can be interpreted as adding equations for the derivatives of the unknowns to the theory matrix, where these equal to zero with some uncertainty variance justified by physics. This gives us a problem with smooth, well-behaved solutions provided that the constraints are strong enough.

The regularized linear least-squares solution of the inverse problem is then

$$\begin{array}{}\text{(31)}& \widehat{\mathit{x}}={\left({\mathbf{A}}_R^{\top }{\mathbf{\Sigma }}_m^{-\mathrm{1}}{\mathbf{A}}_R\right)}^{-\mathrm{1}}{\mathbf{A}}_R^{\top }{\mathbf{\Sigma }}_m^{-\mathrm{1}}\mathit{m},\end{array}$$

where m is the extended measurement vector ${\mathit{m}}^{\top }=\left[{\mathit{V}}^{\top }{\mathbf{0}}^{\top }\right]$, and A_R is the theory matrix A extended with the constraints (29) and (30). We will discuss the measurement error covariance matrix Σ_m in the next subsection.

Figure 4 1σ uncertainty of estimates of electric field (a, c, e) and neutral wind (b, d, f). The left column (a, b) shows the perpendicular east components, the middle column shows the perpendicular north components (c, d), and the right column shows the field-aligned components (e, f). The solid lines shows the results for the daytime profile, the dashed lines are for the night profile, and the dashed–dotted lines show results for the night-time profile with aurora. The colours show different regularization parameters. The cyan lines use the numbers derived from Sangalli et al. (2009), the yellow line shows results where the variation in electric field is one-tenth, and the blue lines show a thousandth of these.

3.1 Uncertainty calculations

The measurement uncertainty of the ion wind vector estimate at a range i is quantified by the covariance matrix in Eq. (6). When we combine measurements from different heights to a single vector, the covariance matrix becomes a block matrix with all individual covariances ${\mathbf{\Sigma }}_{v_i}$ along the diagonal,

$$\begin{array}{}\text{(32)}& {\mathbf{\Sigma }}_V=\left[\begin{array}{lll}{\mathbf{\Sigma }}_{v_{\mathrm{1}}}& \mathrm{\cdots }& \mathbf{O}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathbf{O}& \mathrm{\cdots }& {\mathbf{\Sigma }}_{v_H}\end{array}\right],\end{array}$$

where O is the zero matrix. This assumes that measurements from different heights do not correlate.

When expanding the theory matrix to include the regularizations, we also have to expand the covariance matrix. The inverse problem is regularized with a set of values ${\mathit{\alpha }}_{j,(x,y,z)}$, and ${\mathit{\gamma }}_{(\mathrm{0},\mathrm{1}),j,(x,y,z)}$ which control the smoothness of electric field and neutral wind as a function of height. The values we use for the regularization also form the uncertainty of the added measurements. They are, however, not assumed to be co-varying, and therefore these only add diagonal terms:

$$\begin{array}{}\text{(33)}& {\mathbf{\Sigma }}_L=\mathrm{diag}\left\{{\mathit{\alpha }}_{\mathrm{1},x}^{-\mathrm{2}},{\mathit{\alpha }}_{\mathrm{1},y}^{-\mathrm{2}},{\mathit{\alpha }}_{\mathrm{1},z}^{-\mathrm{2}},\mathrm{\cdots },{\mathit{\alpha }}_{N,x}^{-\mathrm{2}},{\mathit{\alpha }}_{N,y}^{-\mathrm{2}},{\mathit{\alpha }}_{N,z}^{-\mathrm{2}}\right\}.\end{array}$$

The covariance matrix of the regularized measurements then becomes

$$\begin{array}{}\text{(34)}& {\mathbf{\Sigma }}_m=\left[\begin{array}{ll}{\mathbf{\Sigma }}_V& \mathbf{O}\\ \mathbf{O}& {\mathbf{\Sigma }}_L\end{array}\right].\end{array}$$

As the inverse problem then should be solvable using Eq. (31), the uncertainty of the solution is given by

$$\begin{array}{}\text{(35)}& {\mathbf{\Sigma }}_{\widehat{x}}={\left({\mathbf{A}}_R^{\top }{\mathbf{\Sigma }}_m^{-\mathrm{1}}{\mathbf{A}}_R\right)}^{-\mathrm{1}}.\end{array}$$

This can be considered as the a posteriori estimation error covariance for the electric field and neutral wind.

3.2 Regularization parameters

Before calculating the electric field and neutral wind estimate uncertainties by inserting values into the equation, assumptions must be made on how strongly the problem should be regularized. With ISR, the variation in electric field and the neutral wind have typically been measured in their own height ranges, neutral wind up to around 140 km, and electric field above that. Knowledge on the variation at the other heights is sparse, and it is therefore not obvious what good choice for the regularization constants α for the electric field or γ for the neutral wind would be.

Simultaneous observations of electric field and neutral wind have been made with sounding rockets. However, there are not sufficiently many of such measurements to fully characterize the statistics of the altitude variation of electric fields and neutral winds. Altitude profiles of electric field and neutral wind can still be used for estimating typical magnitudes of their gradients in order to find suitable values for the regularization parameters α and γ. Here, we will use measurements from the Joule II rocket campaign, where altitude profiles from 85 to 210 km of electric field and neutral winds below 130 km were derived (Sangalli et al., 2009). Since the rocket did not travel exactly along the magnetic field line, the variation of the electric fields along the trajectory is larger than along the magnetic field. Therefore, the variance of the electric field gradient will be overestimated, leading to a softer regularization.

At higher altitudes, the electric field is expected to be constant along the magnetic field because of the high field-aligned conductivity. We therefore use two estimates of the variation of the electric field, one for high and one for low altitudes. We assume that the variance is the same for the three components. Based on the Sangalli et al. (2009) measurements, we estimate that the largest electric field variation is 20 mV/m over a 2.5 km range at about 90 km and 5 mV/m over the same range at 190 km altitude. We set the regularization parameters to match these variations. This means that we assume that the largest variations in electric field measured by the rocket experiment are relatively rare (occur 5 % of the time). Our regularization of the field-aligned gradient is then ${\mathit{\alpha }}_{j,(x,y,z)}^{-\mathrm{1}}$= 1 µV/m² at 190 km and 4 µV/m² at 90 km altitude. In between these, we interpolate the variation linearly. We choose our measurement region to be similar, between 80 and 200 km height, and can also extrapolate the variation linearly. This we will call the “measurement-based” regularization. Additionally, we have calculated the uncertainty for two cases where we constrain the electric field more strongly towards flatness. This can be seen as more similar to the commonly used assumption that the electric field is constant. We do this by dividing the regularization for the E field by 10 and 1000.

We assume a 1σ variation of the neutral wind gradient of 20 m/s/km for all heights. In addition, we add an assumption that the neutral wind estimates follow a normal distribution with zero mean and standard deviation of 200 m/s, which corresponds to using 0.005 s/m as the zeroth-order Tikhonov regularization parameter.

3.3 Ion-neutral collision frequency

Use of the correct ion-neutral collision frequencies is crucial for calculating the ion mobilities correctly. Therefore accurate collision frequencies are necessary for estimating the electric field and neutral wind. The ion-neutral collision frequency can be calculated theoretically (see Schunk and Nagy, 2009) or measured with ISR (Nicolls et al., 2014a; Davies et al., 1997). Both methods will result in uncertainty of the collision frequency on the magnitude of 50 % but somewhat lower for the ISR measurements. In this study, we will ignore this uncertainty. Any uncertainty in the collision frequency will add to the error budget.

In this study, we have calculated the collision frequencies using

$$\begin{array}{}\text{(36)}& {\mathit{\nu }}_{in}={\displaystyle \frac{{\sum }_{i,j}{s}_{i,j}{n}_i{n}_j}{{\sum }_i{n}_i}},\end{array}$$

where s_i,j is the collision frequency coefficient (CFC) between ion species i and neutral species j. We use the CFCs from Schunk and Nagy (2009) for the most usual neutral and ion species N₂, O₂, and O and NO⁺, O⁺, and ${\mathrm{O}}_{\mathrm{2}}^{+}$. Where the collision is resonant, we simply assume a reduced temperature of 400 K to calculate the CFC. The particle densities are calculated by the MSIS atmospheric model (see Picone et al., 2002).

3.4 Electric field and neutral wind uncertainty

We can now investigate the expected performance of E3D for estimating electric fields and neutral winds as a function of height. The variances of the estimates are the diagonal of the a posteriori covariance matrix, Eq. (34). As the performance depends on ionospheric conditions, we study the same three ionospheric conditions as for the ion velocity uncertainty (see Fig. 2). The performance also depends on the a priori smoothness constraints. Figure 4 shows these 1σ uncertainties for the different ionospheric conditions and regularization constraints. The ionospheric conditions are indicated with line style, and the different smoothness assumptions are indicated with colour. Cyan is the measurement-based regularization α, which is defined in Sect. 3.2. We also use two increasingly stronger regularization constraints for the electric field; the yellow line uses 10α, and the blue uses 1000α. This means that yellow and blue lines are assuming the horizontal gradients of the electric field to be a factor of 10 or 1000 smaller in magnitude than the cyan line.

The uncertainties of the perpendicular electric field (Fig. 4a and c) can be divided into two regions: above and below approximately 125 km. Above 125 km, the electric field uncertainty is primarily defined by measurement uncertainty. Below this height, it is primarily constrained by the regularization as ion velocity is less dependent on electric field due to ion demagnetization. For the parallel E field (Fig. 4e), the ionospheric conditions play a smaller role.

At low altitudes the uncertainty is above 10 mV/m for the measurement-based regularization. For the higher altitudes, its size depends on the ionospheric conditions but is around 1 mV/m in the perpendicular directions and approximately a factor of 3 lower in the field-aligned direction. With stronger constraints towards flatness, the uncertainty decreases, but one has to remember that this comes at the cost of blurring out smaller scale variations. For the lowest range, the yellow and cyan lines indicate estimates of the electric field with too large an uncertainty to be useful. This means that we can not measure electric field with a useful accuracy below 125 km unless we can make assumptions of horizontal gradients being less than approximately 4 nV/m² (blue line).

The uncertainties of the neutral wind components are shown in Fig. 4b, d, and f. The neutral wind is best estimated to an accuracy of approximately 10 m/s between 90 and 125 km. Below 90 km, the electron density is typically lower, which increases LOS ion velocity measurement errors. Above 125 km, the ion-neutral collision frequency decreases rapidly, which makes the ion drift increasingly independent of the neutral wind. At highest altitudes, the uncertainty is merely constrained by our assumptions on neutral wind amplitude (200 m/s). Best estimates are obtained at around 100 km altitude.

The usable range of neutral wind measurements depends strongly on the prior assumption on the smoothness of the electric field. The strongest regularization, corresponding to the smallest horizontal electric field gradient assumption, indicated with the blue line, leads to a neutral wind uncertainty of less than 30 m/s up to 150 km. However, this altitude is greatly reduced with less strict prior assumptions on the electric field gradient; see yellow and cyan lines.

Typical values for the perpendicular electric fields are on the order of tens of millivolts per metre (e.g. Nygrén et al., 2011; Sangalli et al., 2009). Sometimes they are smaller, as reported by Nygrén et al. (2012), but can also be an order of magnitude larger during active auroral conditions, as seen by Dahlgren et al. (2011). Such electric fields are larger than the measurement uncertainties of the possible E3D estimates at high altitudes, even without restrictive regularization. Even for our worst case, quiet night-time conditions, the E3D ion velocity uncertainties are smaller than typical ion drifts. At lower altitudes, the electric fields grow.

Typically the horizontal components of the neutral wind are on the order of tens of metres per second. However, it is not uncommon with stronger neutral winds (Heinselman and Nicolls, 2008; Nygrén et al., 2011, 2012; Brekke, 2013). At collision-dominated altitudes below 115 km, the uncertainties are smaller than these typical values, so accurate neutral wind measurements should be achievable. The vertical component of the neutral wind is smaller but can in extreme cases reach 100 m/s (Nygrén et al., 2011, 2012; Brekke, 2013). Only in the lowest part of the E region is the expected uncertainty smaller than typical vertical winds.

It is important to remember that the results and their uncertainties presume that all assumptions of the flatness of the electric field or neutral wind profile are true. If our assumptions on the magnitude of the electric field gradients or neutral wind gradients are too small, the uncertainties presented are overly optimistic.

We can compare our results with earlier measurements of Dahlgren et al. (2011), which used the tristatic EISCAT UHF to measure the electric field at 220 km altitude under similar conditions as we used for our aurora case. The experiment setup was similar, except for the radar itself. If we look at the time period between 19:28 and 19:36, the horizontal electric field components had a magnitude of up to 250 mV/m but mostly around 30 mV/m. Typical standard deviations are tens of millivolts per metre. With our model and E3D, such electric fields should be measurable with a factor of 10 improvement of uncertainty down to approximately 125 km.

The earlier mentioned Joule II rocket experiment was accompanied with ion velocity measurements at PFISR, which were used to estimate the neutral wind at the same heights (Heinselman and Nicolls, 2008). Also here, the ionospheric conditions look most like our aurora example, but the radar pointed in seven directions to find the different ion wind components. The neutral wind profiles were integrated over 15 min, and their uncertainties are similar to those in Fig. 4b and d at the highest altitudes but somewhat higher further down.

Figure 5 Example simulation of electric field and neutral wind estimates. The estimates were calculated from simulations of ion velocity based on measurements of Sangalli et al. (2009). The layout of the figure is as in Fig. 4. The colours show different regularization parameters. The cyan line uses the numbers derived from Sangalli et al. (2009), the yellow line shows results where the variation in electric field is one tenth, and the blue line shows a thousandth of these. The black line shows the values which were used to simulate the ion velocity. We note that the axes on the plots are different.

4 Simulated measurement

In order to demonstrate what a electric field and neutral wind estimate profile could look like, we have simulated a E3D measurement and analysed it. We based the simulations on the Joule II rocket measurement presented in Sangalli et al. (2009). During the downleg flight, the rocket measured neutral wind at altitudes 90–130 km by tracing chemical releases. The electric field was measured already from 210 km altitude. Since Sangalli et al. (2009) did not include field-aligned components, we used a synthetic profile. We used the electric field and neutral wind profiles to simulate E3D ion velocity measurement with noise added from Eq. (26). These simulated ion velocity measurements were used to estimate electric field and neutral wind. By comparing these with the original data set, we can visualize how good the E3D estimates are. The results are shown in Fig. 5. We use the same regularization schemes as for Fig. 4.

The results confirm that the electric field is estimated well above 125 km, as predicted by the uncertainty estimates in Sect. 3.4. Below 125 km the electric field is not estimated well. For all regularization schemes, the behaviour at lower altitudes is similar. The electric field is estimated to be a constant value corresponding approximately to the value at 125 km altitude, indicating that the regularization contributes with all information of the electric field where ions are demagnetized.

The neutral wind is in general better estimated below 120 km altitude, where it is the largest influencer of the ion wind. Above approximately 125 km, the neutral wind is not well measured in any of the cases. This is not surprising as the neutral wind has little effect on the ion velocity at higher altitudes.

When using the strongest flatness constraints on the electric field, this causes the estimates of the neutral wind to be more fluctuating than the original values (see, for example, Fig. 5b). We believe that the reason is that the model tries to fit the unknowns to the ion wind measurements but “knows” a priori that the electric field is constant, so all the variation in ion velocity must be explained by the neutral wind instead of the electric field. If the constraints on the electric field are relaxed, the estimate of all unknowns is closer to the original values.

At the heights where the regularization plays a smaller role, the deviation from the original values seems similar to the predicted uncertainties shown in Fig. 4.

5 Discussion

Earlier ISR studies on neutral wind have assumed that the electric field is exactly constant along the magnetic field – mainly due to the lack of three-dimensional ion vector velocity measurements along the whole radar transmit beam. The technique presented in this paper allows us to relax this assumption with a scheme that arises from Maxwell's equations and assumption of horizontal smoothness of electric field. A special case of our regularization scheme is the case where the electric field is approximately constant as a function of height along a magnetic field line. This corresponds to a very strong smoothness assumption on horizontal gradients of electric field (see dark blue line in Figs. 4 and 5). The technique presented in this study can be thus seen as a generalization of the commonly used technique for estimating electric field and neutral wind.

In addition to studies of electric fields and neutral winds separately, a use case of the technique presented here is investigation of Joule heating. For Joule heating, both electric fields and neutral winds are necessary (e.g. Aikio et al., 2014, and references therein). With improved estimates of electric field and neutral wind, the Joule heating can be calculated with higher accuracy.

Our results indicate that it will be possible to observe an altitude profile of electric field and neutral velocity using E3D. However, it is only possible to reconstruct either the electric field or the neutral wind at any given altitude region. This is ultimately due to the fact that above an altitude of approximately 125 km, the ion drift is to a large extent determined by electric field and nearly unaffected by neutral velocity. Similarly, below 125 km, the ion drift is primarily determined by neutral wind.

For future measurements, one important question to solve is what regularization parameters should be used. If the constraints are too weak, the problem is underdetermined and the solution noisy. The classical approach is to assume that the electric field is constant along the magnetic field line. However, the electric field may not always be the constant. Then, as the example shows, the neutral wind estimates have to compensate for variations in the ion drift measurements due to under-resolved electric field variability. Relaxing the assumption of a constant electric field will in these cases improve the results.

Adjusting the regularization constants must be done with caution, since the problem easily becomes underdetermined. Therefore it is important to justify the choices of regularization. For the electric field, we used regularization parameters, which are estimated from in situ rocket measurements. However, the optimal values of the regularization parameters for general use are still to be found. For the values we used, only constraining the electric field was not enough, and we also constrained the gradient of the neutral wind in the same way. Here too, the exact values can be discussed. Forcing the neutral wind velocity gradients to be too small causes the estimates of the neutral wind to fit worse to the ion wind at collision-rich heights. This then increases the noise in the electric field here. If the variation is allowed to be too large, the problem is not solvable. In order to allow for higher variations in the neutral wind but also to use all information we have about it, we added a size constraint of 200 m/s. As can be seen in the uncertainty plots, this restricts the size of the neutral wind components to become smaller.

In future work, the model can, for example, be improved in one of the following ways. If somehow measurements of neutral wind or electric field exist, these can be added additionally as constraints. Such measurements could, for example, be the movement of meteor smoke, polar mesospheric summer echoes, or other measurements of events in the ionosphere that imply size or direction of neutral wind or electric field. An independent measurement of mean neutral wind can often be obtained up to about 100 km using meteor radars (e.g. Stober et al., 2018).

In this work, we used the same resolution in time for both electric field and neutral wind. The large mass in the neutral atmosphere causes the neutral wind to vary more slowly than the electric field. Nygrén et al. (2011) took advantage of this to use different time resolutions for the different parameters. In the future, it would be an advantage to include this for our model as well.

The technique discussed in this study can be extended further. With the help of phased array technology, E3D will allow fast beam scanning to be used to measure how ion vector velocity and electron density vary within a volume of space. This type of measurement may potentially result in improved estimates of electric field and neutral wind, as more physics-based regularization can be added. We can use Gauss' and Faraday's laws without the need to treat the horizontal gradients as unknown random variables, as they will be determined by the measurements. We can also introduce constraints that are not possible for a one-dimensional profile. It will be possible to apply Ampere's law to enforce current continuity. We can also apply the Navier–Stokes equations to enforce that the neutral wind is approximately consistent with anelastic flow. Estimating electric field and neutral wind within a volume is a topic of future work.