New method for solving inductive electric fields in the non-uniformly conducting ionosphere

We present a new calculation method for solving inductive electric fields in the ionosphere. The time series of the potential part of the ionospheric electric field, together with the Hall and Pedersen conductances serves as the input to this method. The output is the time series of the induced rotational part of the ionospheric electric field. The calculation method works in the time-domain and can be used with non-uniform, time-dependent conductances. In addition, no particular symmetry requirements are imposed on the input potential electric field. The presented method makes use of special non-local vector basis functions called the Cartesian Elementary Current Systems (CECS). This vector basis offers a convenient way of representing curl-free and divergence-free parts of 2-dimensional vector fields and makes it possible to solve the induction problem using simple linear algebra. The new calculation method is validated by comparing it with previously published results for Alfv én wave reflection from a uniformly conducting ionosphere.


Introduction
In this paper we present a new method for calculating inductive electric fields in the ionosphere. It is well established that on large scales the ionospheric electric field is well approximated by a potential field (e.g. Untiedt and Baumjohann, 1993). This is understandable, since the temporal variations of large-scale current systems are generally quite slow, in the time scales of several minutes, so inductive effects should be small. However, studies of Alfvén wave reflection, made, e.g. by Yoshikawa and Itonaga (1996); Buchert and Budnik (1997); Buchert (1998); Yoshikawa and Itonaga (2000); Correspondence to: H. Vanhamäki (heikki.vanhamaki@fmi.fi) Lysak and Song (2001) and Sciffer et al. (2004), have indicated that in some situations inductive phenomena could well play a significant role in the reflection process, and thus modify the nature of the ionosphere-magnetosphere coupling.
Recently Vanhamäki et al. (2005) showed by approximate calculations that inductive electric fields associated with some very dynamic ionospheric phenomena, including Westward Travelling Surge (WTS), -bands and Giant Pulsations, are locally very significant. They calculated the inductive fields caused by self-induction in the ionosphere (primary process) and also by the ground-induced currents flowing in the conducting ground (secondary process). They concluded that the inductive electric fields are indeed small at large scales, but may locally be almost as large as the associated potential fields. These local "hot-spots" tended to occur in those areas where the field-aligned currents (FAC) were largest, so in these areas the inductive processes could well contribute to the ionosphere-magnetosphere coupling. Vanhamäki et al. (2005) also concluded that at ionospheric altitudes the secondary contribution from ground induction is always very small and smoothly distributed, and in practice negligible when compared to the larger and more concentrated primary contibution from ionospheric self-induction. However, the calculation method used by Vanhamäki et al. (2005) was rather approximate, giving only order of magnitude estimates. The induced electric fields in the ionosphere were calculated as vacuum fields, i.e. the currents driven by the induced fields themselves were neglected. This approximation probably gives induced electric fields that are too large, as the effect of the neglected current should tend to decrease the induced fields according to Lenz's law. These preliminary results indicate that the inductive processes in the ionosphere are significant enough to merit a closer study. As a first step in this direction we present in this paper an exact and self-consistent method for calculating the ionospheric inductive electric fields.
The new method that we introduce is not based on the reflection of incident Alfvén waves at the ionospheric boundary, as the previous models by Yoshikawa and Itonaga (1996); Buchert (1998); Lysak and Song (2001) and Sciffer et al. (2004). Instead, we divide the ionospheric electric field into potential and rotational parts, assume the potential part to be known and calculate the induced rotational part. One may think that the potential and rotational electric fields are created by incident and reflected Alfvén, and fast magnetosonic waves, repectively, but this is not necessary in our calculation method. This kind of approach allows us to solve the induction problem with non-uniform and time-dependent conductances, in contrast to previous studies. Because we do not model the propagation of Alfvén waves we assume instantaneous changes along the magnetic fields lines. Therefore, the vertical variations of the fields are not included, but this is not a problem, as we are mainly interested in the electric fields and currents at the ionospheric level.
The structure of this paper is such that in Sect. 2 we develop the theoretical part of our calculation method. In Sect. 3 we test the accuracy of the method in the case of Alfvén wave reflection, using the analytical solution given by Yoshikawa and Itonaga (1996) as reference. The results are summarised and future developments are discussed in Sect. 4.

The new method
We use a Cartesian coordinate system where the ionospheric current sheet is taken to be the xy-plane and the z-axis points vertically downwards. The Earth's magnetic field is assumed to be parallel to the z-axis, which is a reasonable approximation in the northern polar region. We also use the thin-sheet approximation, i.e. we assume that the ionospheric Hall and Pedersen conductivities are where δ(z) is the Dirac delta function. This means that all horizontal currents flow as a thin sheet, at altitude z=0. We assume that the ionospheric Hall and Pedersen conductances H and P are known. In the general case, these conductances may be non-uniform and also time-dependent. We also require that the time series of the potential part of the ionospheric electric field, E pot with (∇×E pot ) z =0, is given as input. It may have been obtained from measurements or may be given as output by some magnetospheric model after mapping the magnetospheric electric field down to the ionosphere along the magnetic field lines. The output of the calculation method is the induced rotational part of the electric field, E rot with ∇·E rot =0. The general outline of the calculation method is the following: -Express the potential (E pot ) and rotational (E rot ) parts of the electric field using the Cartesian Elementary Current Systems (CECS, these are discussed in Sect. 2.1).
-E pot is associated with a current system j 1 and E rot with j 2 . Express these currents with CECS and use Ohm's law to relate the CECS representations of E pot and j 1 , as well as E rot and j 2 .
-Calculate the magnetic field B created by the currents.
-Faraday's law gives an equation that relates the unknown scaling factors of the CECS representation of E rot to the scaling factors of the input field E pot and conductances H , P .
Detailed description of the above steps is given in Sect. 2.2 for a very simple situation and in Sect. 2.3 for the general case.

Cartesian Elementary Current Systems
We represent the ionospheric electric fields and currents by using a special non-local vector basis function, the Cartesian Elementary Current Systems (CECS). CECS were introduced by Amm (1997) and although for historical reasons the name CECS refers to current systems, they can be used to represent any smooth enough (continuously differentiable), 2-dimensional vector field in Cartesian geometry. There are two different kinds of CECS, one type is divergence-free (DF) and the other curl-free (CF). The elementary systems, illustrated in Fig. 1, are defined as Here we use a cylindrical coordinate system that is centered at the pole of the CECS. The scaling factors of the elementary systems are denoted by I cf and I df , while U is the Heaviside unit step function. The above definition is suitable for representing ionospheric current densities, as the vertical currents are restricted to the plane z=0 and also the FAC are included in Eq.
(2). The ionospheric electric field, however, is not limited to one plane. In fact, the horizontal electric field is nearly constant in the z-direction, while the z-component is negligible due to the very high field-aligned conductivity. Consequently, a suitable form of CECS for representing electric fields is In Eqs.
(1-4) we have already used the convention that scaling factors of CECS which are used to represent currents and electric fields, are marked with I and V (and have units of amperes and volts), respectively. The CECS of Eqs.
(1-4) have the property that the curl and divergence of the horizontal part vanish everywhere except at the poles of the elementary systems. At the pole the horizontal part of CF CECS has a δ-function source and DF CECS has a δ-function curl.
An arbitrary ionospheric current system j (or any similar vector field) can be uniquely expressed by placing a sufficient (in principle infinite) number of CF and DF CECS at different positions in the ionospheric plane. In practical calculations both the current vectors and CECS scaling factors are given at some discrete grid points. In that case we can construct a matrix relation where j and I are vectors containing the vector components of the current density and CECS scaling factors at different grid points, respectively. Matrix K gives the relation between the two representations of the current density. K depends only on the geometry of the grids which are used and can be constructed using Eqs.
(1-2). For a given current system Eq. (5) can be inverted to give the CECS representation of the currents. This is a convenient way of decomposing the currents into CF and DF parts. Ionospheric electric fields can be handled analogously.
The magnetic fields associated with the elementary currents of Eqs. (1-2) are needed in the following calculations. A straightforward calculation (Amm and Viljanen, 1999) using the vector potential gives

Simplest case
In order to illustrate the new technique in a case that can be treated analytically, we consider the simplest possible situation where the ionospheric conductances, H , P are uniform and the input potential electric field consists of one CF CECS with time-dependence e iωt . The induced rotational electric field is represented by a number of DF CECS placed at some grid points, as illustrated in Fig. 2. It should be noted that the grid can be quite arbitrary, although regular square grids are used in numerical calculations in Sect. 3. The total electric field is The first term is the input field and the summation over all the grid cells gives the induced field. In the summation ρ l andê φ l are the distance and unit vector in the φ-direction in the coordinate system centered at CECS pole l.
Because the conductances are uniform, the current system associated with the electric field in Eq. (8) is easy to write down. The ionospheric Ohm's law relating the horizontal electric field and horizontal sheet current density J is Each electric field CF (DF) CECS is associated with one CF (DF) CECS of the Pedersen current and one DF (-CF) CECS of the Hall current, Here we have not explicitly written out the FAC. The magnetic field associated with the above current density can be calculated using Eqs. (6-7). The z-component of the magnetic field at z=0 is Using Eqs. (8) and (11) we can integrate the z-component of Faraday's law, over an arbitrary grid cell k. The time-derivate gives just a factor iω while the curl of E gives a δ-function at every position where the poles of DF CECS are located. The integral over grid cell k is just The remaining integrals in the above equation give just geometrical factors. Let's call the first G k and the second H k,l . If we gather the scaling factors V rot k to a vector V rot and similarly µ 0 G k /(4π ) to a vector G and µ 0 H k,l /(4π ) to a matrix H, we can write the above equation in the form The unknown DF CECS scaling factors V rot can be solved as The solution given in Eq. (15) can be used with arbitrary input fields E pot just by expressing the input in terms of CF CECS and using linear superposition. With non-harmonic time-dependence we do not have an algrebraic solution, but Faraday's law, Eq. (12), gives a first-order differential equation in time for the unknown DF CECS scaling factors V rot .
In the next section, we discuss the more complicated case of non-uniform conductances.

General case
The ionospheric Hall and Pedersen conductances can be very non-uniform and also change in time scales of few minutes or seconds, especially during high auroral activity. In this section we discuss the most general situation, where the input potential field E pot and conductances P , H may be arbitrary (yet physically reasonable) functions of position and time.
The input electric field E pot can be expressed in terms of CF CECS as in Eq. (5), Here the matrix N 1 depends only on the geometry of the grids that are used with the vector field E pot and the CF CECS scaling factors V pot , and it can be constructed using Eq. (4).
The current system j 1 associated with the electric field E pot can be calculated using Ohm's law, Eq. (9). Because Ohm's law is a linear relation between the electric field and current, it is possible to define matrix M 1 so that However, j 1 can also be expressed in terms of CECS as in Eq. (5), j 1 =K 1 ·I 1 . Inverting matrix K 1 we have a relation between the CECS representations of E pot and j 1 , A similar relation also holds for the unknown rotational electric field E rot and currents j 2 associated with it, It should be emphasized that matrices K 1 , K 2 (which depend on geometry of the calculation grids) and M 1 , M 2 (which depend on geometry and conductances) can be constructed using Eqs.
(1-4) and (9) without knowing the electric field. The next step is to calculate the z-component of the magnetic field associated with the currents. According to Eqs. (6-7), only the divergence-free part of currents j 1 and j 2 is needed. We can define new matrices L 1 and L 2 so that This can be done by picking the appropiate rows of the matrices inv(K 1 )·M 1 and inv(K 2 )·M 2 in Eqs. (18) and (19). Strictly speaking, the above matrix relations are only valid if the electric fields and currents have no sources and curls outside the grid area. In practise it is enough to choose the grid to be suitably larger than the area of interest. The zcomponent of the magnetic field can be written out as Here the first summation is over the currents associated with the input potential field and the second is over the induced system. As in the previous section, integration of Faraday's law, Eq. (12), over an arbitrary grid cell k gives We again define factors G k,m and H k,l so that they contain the geometrical factors from the remaining integrals, together with the constant µ 0 /(4π ). The matrix form of the above equation is where we have gathered the geometrical factors G k,m and H k,l into matrices G and H, respectively, and made use of Eqs. (20) and (21). The matrices L 1 and L 2 are easy to calculate numerically by following the steps outlined above. If the time series of the input potential field V pot and conductances P , H are known, Eq. (24) can be step-by-step integrated to give the scaling factors of the induced rotational electric field, V rot .

Comparison with previous results
Ionospheric induction effects have previously been studied, e.g. by Yoshikawa and Itonaga (1996). They presented analytical results for Alfvén wave reflection from a uniformly conducting ionosphere. The calculation method presented in Sect. 2 does not describe wave reflection per se, so a direct comparison with results by Yoshikawa and Itonaga (1996) is not possible. However, comparison can be made using the following steps: -Define the potential electric field of the incident wave.
-Calculate the reflected potential and rotational fields using the reflection coefficients given by Yoshikawa and Itonaga (1996).
-Use the total potential field (sum of incident and reflected parts) as input to the calculation method of Sect. 2.
-Compare the output rotational field given by the new calculation method with the reflected rotational field calculated using the reflection coefficients of Yoshikawa and Itonaga (1996).
The purpose of the comparison is to demonstrate that the new method gives the correct results in this special case of the Alfvén wave reflection from a uniformly conducting ionosphere. As mentioned above, the new calculation method is not based on the concept of Alfvén wave reflections, but instead uses the potential part of the total ionospheric electric field as input. Therefore, in this comparison we have to first determine the incident and reflected potential fields using the reflection coefficient given by Yoshikawa and Itonaga (1996). After that is done we may compare the induced rotational fields calculated using the two different methods.
In the test cases we assume the electric field of the incident Alfvén wave to consist of one CF CECS with a timedependence e iωt and an amplitude of 10 4 V. We take the ionospheric conductances to be P =2 S, H =4 S and the Alfvén speed above the ionosphere is V A =500 km/s. The reflected potential and rotational fields are also represented in terms of CF and DF CECS placed at a regular grid. The reflection coefficients given by Yoshikawa and Itonaga (1996) are functions of the horizontal wave number of the incident field, so calculations must be done in Fourier space. Some details of the calculation are given in the Appendix.
The ionospheric potential field that is obtained from the incident and reflected waves is used as input in Eq. (15). In this case the input field consists of several CF CECS, so we solve the problem for each input CF CECS separately and sum the results. In each test case the calculation grid we use in Eq. (15) is twice the size of the grid where the input is given. So, if the input obtained from the Alfvén wave reflection is given in an N×N grid, we use a (2N-1)×(2N-1) grid when solving Eq. (15).
The results of the four test cases are given in Tables 1-4. In each case the single CF CECS of the incident Alfvén wave is located in the middle of the calculation grid. The situation is rotationally symmetric, so in Tables 1-4 we show the resulting DF CECS scaling factors only in the upper left quarter of the grid. In each table four quantities are given: the amplitude |V rot | and phase arg(V rot ) of the CECS, together with the relative error in the amplitude (|V rot |−|V rot 0 |)/|V rot 0 | and phase arg(V rot )−arg(V rot 0 ) with respect to the reference results V rot 0 by Yoshikawa and Itonaga (1996). In test case 1 we use 11×11 grid for the input CF CECS (so in Eq. (15) we use 21×21 grid). The grid spacing is 50 km and the angular frequency is ω=2π/(60 s). Test case 2 is otherwise similar, but with a higher frequency ω=2π/(1 s). In case 3 the grid is still 11×11, but the grid spacing is reduced to 10 km while the frequency is the same as in case 2. Test case 4 is similar to case 3, but the input grid is now 27×27 although only the center part is given in Table 4.
In case 1 the errors in the amplitude and phase are reasonably small near the center of the grid. At the edges the relative errors increase, probably due to some boundary effects, but in absolute terms the errors are about the same size as in the center.
With the higher frequency used in case 2 the errors increase and only the few centermost grid cells are calculated Therefore, we may expect that with high frequences and/or a highly conducting ionosphere we need to use smaller grid spacing, in order to maintain a reasonable accuracy. Table 3 gives the results of test case 3, where the grid spacing is 10 km. The results are clearly better, although errors are still quite large near the boundaries. In order to minimize the boundary effects, the input in test case 4 is given in 27×27 grid (so in Eq. 15 we use a 53×53 grid), although only the same centermost grid cells as in the previous examples are shown in Table 4. Now the errors are reduced to a quite reasonable range throughout the center of the calculation grid. A similar calculation with the lower fre- quency ω=2π/(60 s) (not shown here) gives even smaller errors, less that 4% in amplitude and 3 degrees in phase.

Summary and discussion
We have presented a new calculation method for solving the inductive electric fields in the ionosphere. In contrast to many previous studies the new method does not solve the induction problem in terms of incident and reflected Alfvén or fast magnetosonic waves. Instead, the input quantities in this method are the potential part of the ionospheric electric field, together with the Hall and Pedersen conductances as functions of time. The output quantity is the induced rotational part of the electric field. The calculation method works in the time-domain and can handle non-homogenous and also time-dependent conductances.
The new method makes use of special non-local vector basis functions, CECS, that are used to represent the ionospheric electric fields and currents. Already the basis functions are divided into CF and DF types, which is very convenient in many situations encountered in ionospheric electrodynamics. The entire spatial structure of the vector fields can be "hidden" into the scaling factors of the CECS that are used to represent the fields. For example, the solution of the general problem, Eq. (24), is a differential equation only in time, while all the spatial relations are in the matrices G·L 1 and H·L 2 .
Another interesting feature is that we did not have to specify any explicit boundary conditions when we derived the solution in Eq. (24). However, as mentioned in Sect. 3,  8 −179.7 −171.5 −164.2 −157.6 −152.6 −150.7 −175.8 −166.5 −157.6 −148.5 −140.5 −137.1 −173.1 −163.1 −152.6 −140.5 −126.7 −118.6 −171.6 −161.8 −150.7 −137.1 −118.6 −101.8 Error in amplitude, % when we use the CECS representation, we implicitly assume that the electric fields and currents do not have any curls or sources outside the analysis region. In practise this does not seem to be a problem, at least if the analysis region is chosen to be suitably larger than the area of interest. This is illustrated in test cases 3 and 4 of Sect. 3, where the boundary effects decreased significantly with increased grid area. The test cases of Sect. 3 also showed that the new calculation method is reasonably accurate, if the calculation grid is chosen correctly. Some part of the differences between the two compared methods, especially near the boundaries of the grid, are probably explained by the very different calculation techniques. In our own method the calculation area and the spatial resolution are determined by the grid used. On the other hand, Yoshikawa and Itonaga (1996) used Fourier space, where the calculation area is the whole xy-plane and the spatial resolution is, in principle, infinite. The improvement of accuracy from test case 2 to 3 and 4 indicates that the results of the new method would converge to those of Yoshikawa and Itonaga (1996) in the limit of an arbitrarily large and dense calculation grid. It should also be noted that in absolute terms the errors in test case 4 are already very small, since the largest relative errors occur near the boundaries where the CECS scaling factors are quite small. The test cases were limited to the rather special situation of uniform conductances, but the results prove our new calculation method to be correctly formulated and there seems to be no reason why it should not also perform equally well in the general situation with non-uniform conductances.
In future papers we will apply the presented calculation method to different phenomena that are observed in the ionosphere. These will certainly include the Westward Travelling Surge (WTS) and -band models already studied in a more approximate way by Vanhamäki et al. (2005). These models provide a realistic potential electric field and conductance distributions that are based on measurements. With the presented calculation method we can calculate the associated rotational electric field and the currents driven by it. This gives us a direct method for studying inductive processes in the ionosphere. In principle, it would be possible to perform similar studies using the concept of Alfvén wave reflections, but in practise it would be very difficult for the following reasons: -Previous considerations of the Alfvén wave reflection process, e.g. by Yoshikawa and Itonaga (1996); Buchert (1998), andSciffer et al. (2004), are based on the assumption of uniform ionospheric conductances. In many cases this is not a valid approximation, especially in the auroral regions.
-The total ionospheric electric field is relatively easy to obtain from measurements, but it is very difficult to decompose it into incident and reflected waves. Consequently, there are no models that describe the spatial structure of the incident Alfvén waves that are associated with specific ionospheric phenomena, like the WTS.
Thus, we see that the new method allows us to perform more general studies of ionospheric induction phenomena using realistic electric field and conductance configurations. A major topic of future studies will be the quantitative estimation of the role of inductive phenomena in ionospheric electrodynamics and especially in the ionospheremagnetosphere coupling. We also endeavor to develop a 3-dimensional generalization of the presented calculation method, so that we could make use of altitude-dependent conductances and electric fields. Additionally, one interesting possibility would be to use the presented calculation scheme as an ionospheric solver in a global MHD simulation. Yoshikawa and Itonaga (1996) include a perfect conductor at depth d below the ionospheric plane for modelling the solid Earth. In order to obtain results that can be compared with our own calculation method, we have modified the above formulas so that this perfect conductor is not present. The required modification is to simply change atm coth(k ⊥ d)→ atm in all the formulas given by Yoshikawa and Itonaga (1996). This can be verified by solving the ionospheric boundary value problem as done by Yoshikawa and Itonaga (1996) with the new condition that below the ionosphere there are only downward propagating waves.
We assume that the incident electric field consists of one CF CECS, The divergence of the incident field is a δ-function and its Fourier transform is just a constant, V pot /(2π ). According to Eqs. (A4) and (A6) the curl of the reflected field is where a = µ 0 ω P + 2 H A + P .
Here we have used the fact that the reflection coefficient R A→F is a function of the amplitude of the spatial wave number k ⊥ . The angular integral in the above equation can be calculated using formula 3.387.2 of Gradshteyn and Ryzhik (1965) and J 0 is the zeroth order Bessel function of the first kind. According to Eqs. (A3) and (A5) the divergence of the reflected field, α r , is The δ-function at the position of the incident CF CECS's pole in the above equation comes from the constant term in Eq. (A5). The integral in Eq. (A10) has to be calculated numerically.
The scaling factors of the reflected DF CECS are obtained by integrating β r over the grid cells, This integral can be evaluated numerically using, e.g. Gaussian integration. The only exception is the grid cell k=0 that contains the incident CF CECS, because β r diverges at ρ=0. This can be handled by approximating the cell k=0 by a circle with some radius r, for in this case the area integral can be evaluated analytically using formula 5.52.1 of Gradshteyn and Ryzhik (1965), The scaling factors of the reflected CF CECS are