Three‐Dimensional Closure of Field‐Aligned Currents in the Polar Ionosphere

Using a simplified three‐dimensional Hall‐magnetohydrodynamics simulation, we investigated the current closure of field‐aligned currents (FACs) in the polar ionosphere. Ion‐neutral collision was taken into consideration. To excite a pair of FACs, an electric field perturbation is applied to the upper boundary of the simulation box. The flow shear propagated downward accompanied with the FACs. When the electron density was initially uniform, most of the FACs are connected with the Pedersen current due to electrostatic processes. Some of them are connected with the Hall current due to inductive processes. When the density was initially enhanced in a longitudinally elongated region (high‐density band), overflow of the Hall current takes place near the edge of the high‐density band. Additionally, localized FACs appear to bridge between the Pedersen current layer (high altitude) and the Hall current layer (low altitude). The formation of the additional FACs is closely associated with the field‐aligned gradient of ∇∙E, where E is the electric field. We compared the current lines with those evaluated by the traditional thin‐layer assumption. The current closure obtained by the thin‐layer assumption is fully different from that obtained by the three‐dimensional model. In the full three‐dimensional simulation, a current line flowing in the Hall layer can pass underneath the current flowing in the Pedersen layer. Such “intersection” is not allowed in the thin‐layer assumption. We believe that the three‐dimensional model offers advantages for fully understanding the closure of the FACs together with the current closure on the magnetospheric side.

At high latitudes, the inclination of the geomagnetic field is almost vertical, so that the horizontal electric field may be independent of altitude. In the thin-layer approximation, one may introduce height-integrated Ohm's Law as , Equation 4 implies that the FACs can be connected with the Hall current when the gradient of the Hall conductivity is present. Such connection has been thought to occur near the bright aurora where the large gradient of the Hall conductivity is present (Akasofu et al., 1981;Ebihara & Tanaka, 2015a;Fujii et al., 2011;Kan et al., 1984). Three questions arise. The first one is how the FACs are connected with the ionospheric current three-dimensionally. Since the peak altitude of the Pedersen conductivity is higher than that of the Hall current conductivity (Baumjohann & Treumann, 1996), the Pedersen current flows at higher altitudes than the Hall current does. The three-dimensional current closure is problematic. The second question is regarding energy. The Hall current is orthogonal to the electric field, so that it does not contribute to dissipation. When the gradient of the Hall conductivity is present, space charge will be accumulated. The excess charge (resulting from the discontinuity of the Hall current) generates the secondary electric field that can result in the Joule dissipation due to the Pedersen conductivity. Fujii et al. (2011) raised a question where the energy setting up the additional electric field comes from. They introduced the Pedersen current layer and the Hall current layer, and pointed out that the third term on the right-hand side of Equation 4 can be closed by two routes. One is a connection with FACs to the magnetosphere, and the other one is a connection with the Pedersen current at different altitudes. They proposed three-dimensional current system in the presence of the conductivity gradient ( Figure 5 in Fujii et al., 2011). The third question is the importance of the last term on the right-hand side of Equation 3. Yoshikawa and Itonaga (2000) and Yoshikawa (2002) formulated the equations describing the electromagnetic (inductive) processes, and suggested the importance of the term. Apart from theoretical considerations, self-consistent numerical simulations have been performed to understand the magnetosphere-ionosphere coupling without relying on the thin-layer assumption. Neukirch et al. (1995) developed a stationary two-dimensional fluid model, and showed the current closure. Dreher (1997) performed a two-dimensional Hall magnetohydrodynamics (MHD) simulation that takes into consideration the ion-neutral collision. They showed the current closure and the density perturbations. Lysak (1997) developed the two-dimensional simulation model which used the height-resolved conductivities. This model indicates the coupling of the shear Alfvén wave to the compressional mode wave by the Hall conductivity. Zhu et al. (2001) developed a two-dimensional magnetosphere-ionosphere coupling simulation, incorporating altitude-dependent ionization rates and cooling rates of electrons and ions. The ionospheric plasma parameters were modified, which were fed back into the transportation of plasma. Ion heating primarily resulted from plasma-neutral friction, whereas electron heating occurs due to the precipitation of energetic electrons as well as Joule dissipation. A global two-dimensional multifluid-collisional Hall MHD simulation was developed by Tu and Song (2016). They demonstrated propagation of the Alfvén wave and the coupling of the FACs and the horizontal current. When the wavefront reached the bottom of the ionosphere, the Pedersen current appears so as to connect the FACs. The Hall currents were much weaker than the Pedersen current. Lysak (2004) developed a three-dimensional magnetosphere-ionosphere coupling simulation. The ionosphere was treated as a thin layer. He showed the propagation of ultralow-frequency waves in the dipole magnetic field. Lysak et al. (2013) introduced a new model including a height-resolved ionospheric conductivity. They suggested that the inclusion of the height-resolved ionosphere results in the enhancement of the damping of ionospheric Alfvén resonator (IAR) mode waves.
From the theoretical considerations, three-dimensional current closure has been suggested to present in the polar ionosphere, but the detail is not well known. As far as we know, few simulation studies have focused on the three-dimensional connection between the FACs and the ionospheric current. This study is aimed at clarifying the three-dimensional current closure by using self-consistent three-dimensional simulations. Although our simulation is somewhat simplified, we believe that this is a key step in understanding the whole current closure between the magnetosphere and the ionosphere.

Simulation Model
We applied the governing equations used by Dreher (1997) for the three-dimensional problem. The governing equations were derived from the equation of continuity of plasma, the law of conservation of momentum of ions and electrons, Ohm's Law, Ampère's Law, and Faraday's Law as summarized below. Assuming electrical quasi-neutrality, we can describe the equation of continuity as where n is the number density of electrons and v i is the ion velocity. Following Dreher (1997), we introduced S to take into account the exponential relaxation of the electron density to the initial one n 0 as where r is a constant, which is assumed to be 1 s −1 for all the altitudes. To make the magnetosphere-ionosphere coupling processes clear, ionization due to external sources, including solar extreme ultraviolet radiation and precipitation of energetic particles, was not considered. We described the conservation of momentum of ions as where m i , j, B, and ν in are the mean ion mass as m i = 25 m p = 4.25 × 10 −26 kg, the current density, the magnetic field, and the ion-neutral collision frequency, respectively. m p stands for the mass of a proton. The ion-neutral collision frequency ν in is dependent on the altitude z, and is expressed as where in E  = 2 × 10 3 s −1 , L ν = 40 km, and z min = 70 km altitude. Figure 1 shows ν in as a function of altitude.
Equations 5 and 7 describe the advection of the density and the momentum, respectively. The advection equations were solved by using the implicit scheme and the Lax-Wendroff scheme with the Superbee limiter function (Roe, 1986). To derive the electric field E, Ohm's Law is solved as where e is the elementary charge. The term on the right-hand side of Equation 9 is called a Hall term, which gives rise to the Hall effect. The current density j is described by Ampère's Law, where μ 0 denotes the magnetic constant. To solve the time evolution of the magnetic field B, we used Faraday's Law as We solved the governing equations numerically in the cartesian coordinate system. The x-, y-, and z-directions look eastward, northward, and upward, respectively. The simulation box is bounded from x min = −250 km to x max = 250 km, from y min = −250 km to y max = 250 km, and from z min = 70 km to z max = 1,200 km. To generate a pair of FACs, an electric field perturbation E p was applied to the upper boundary surface at z max as where ˆx E e is a unit vector of the x-direction, and f is frequency = 10 −2 Hz. The perturbed electric field gives rise to flow shear. The scale lengths of L x and L y are set to be 100 km. An example of E p is presented in Figure 2. The compressional wave is not expected to propagate efficiently since the   x [km] wave frequency is below the cutoff frequency (Lee & Lysak, 1990). Only the shear Alfvén waves are expected to propagate selectively along field lines.
Initially, the magnetic field is set to be downward and homogeneous, that is, where 0 E B  is set to be 5 × 10 −5 T. B is fixed to be B 0 at the bottom boundary (at z = z min ), whereas an open boundary condition is applied to the top boundary (at z = z max ). The open boundary condition represents the condition that the spatial derivative is set to be zero. The ion velocity v i is initially set to be 0. We employed three different initial conditions for the electron number density. In Case 1, the number density is homogeneous. In Case 2, the number density is artificially enhanced at |y| ≤ 50 km (narrow high-density band). In Case 3, the number density is enhanced at |y| ≤ 200 km (wide high-density band). Namely, n 01 , n 02 , L 2, and L 3 are set to be 3 × 10 10 m −3 , 6 × 10 10 m −3 , 50, and 200 km, respectively. The open boundary conditions are applied for the ion velocity and the number density.

Results
Due to the perturbation of the electric field imposed at the top boundary, the shear of plasma flow evolves downward with an Alfvén speed (1,250 km s −1 for n 01 ) together with FACs. Figure 3 shows the time evolution of j z . Hereinafter, j z is regarded as j || (FAC, positive upward) because |B z | >> |B x | and |B z | >> |B y |. A pair of the FACs, downward one in the positive x region, and upward one in the negative x region, appear to develop as time proceeds. At t = 4.0 s, the leading edge of the Alfvén wave has completed two round trips between the upper boundary and the lower region. In Case 1 (where the number density is uniform), the FACs are symmetrical with respect to the x-axis. In Case 2 (where the number density is inhomogeneous), the FACs are asymmetric, in particular, at 150 km altitude and lower. One can find the region where j || is locally developed at low altitudes, like the toe of a boot as shown in Figure 3f. The altitude-dependent distribution of j || is more clearly shown in Figure 4. We will explain the formation and closure of the localized FACs later.
At high altitudes where the ion-neutral collision frequency is low, polarization currents are responsible to ensure the current continuity. In contrast, at low altitudes where the collision frequency is high, the conductive currents, the Hall and the Pedersen currents, flow perpendicularly to the magnetic field. In this simulation, the Pedersen current flows dominantly around 140 km (which is called a Pedersen layer), and the Hall current flows dominantly around 100 km (which is called a Hall layer). The generation of the Hall current is associated with the term on the right-hand side of Equation 9. Figure 5 shows the perpendicular electric field and the perpendicular current at 140 and 100 km altitudes. At 140 km altitude, the current (blue arrows) is almost parallel to the electric field (red arrows). That is, the Pedersen current is dominant. At 100 km altitude, the current is almost orthogonal to the electric field. That is, the Hall current is dominant. The current tends to flow clockwise (counterclockwise) in the downward (upward) FAC region, which is associated with the divergent (convergent) electric field. The direction of the electric field is generally outward (inward) in the x-y plane from the center of the downward (upward) FACs. In addition to the radial component of the electric field, the counterclockwise (clockwise) component of the electric field can be seen in the downward (upward) FAC region. The azimuthal component of the electric field is closely associated with the inductive processes (Yoshikawa & Itonaga, 2000). Namely, the radial component of the electric field causes the rotational Hall current, which induces the magnetic disturbances in the z-direction. The change in the magnetic disturbances induces the rotational electric field in the x-y plane. This implies that the FACs can be connected with the Hall current in the uniform conductivity due to the inductive process, which is mentioned in more detail below. A noticeable difference between Case 1 and Case 2 is found in the electric field near the origin (x ∼ 0 and y ∼ 0). The electric field directs almost westward (negative y-direction) in Case 1 (Figure 5b), whereas it looks at the southwest direction (negative x-direction and negative y-direction) in Case 2 ( Figure 5d). The deflection of the electric field is primarily caused by space charge (polarization) associated with the gradient of the electron density (conductivity). At x ∼ 0, positive (negative) space charge is accumulated near the northern (southern) edge of the high-density band. Thus, the resultant electric field related to the polarization directs southward (negative y-direction). We think that the boundary conditions do not significantly affect the electric field shown here because of the following two reasons. First, we did not solve the evolution of the electric field directly. Second, the electric field propagated along magnetic field lines.
We present Figure 6 to show ∇•E at 100 km altitude. In Case 1, ∇•E is symmetrically distributed with respect to the x-axis. Near the x-axis (|y| < 50 km), ∇•E is positive at x > 0 whereas ∇•E is negative at x < 0. The positive value is associated with the positive charge deposited by the downward FAC, while the negative one  is associated with the negative charge deposited by the upward FAC. In Case 2, the distribution of ∇•E is no longer symmetric. At x ∼ 0, positive (negative) space charge is deposited at y ∼ 50 km (∼−50 km). When the current (which is dominated by the Hall current at 100 km altitude) intersects the high-density band, space charge is deposited due to the discontinuity of the Hall current. In the high-density band, the positive and negative charges give rise to the southward electric field (negative y-direction), resulting in the deflection of the total electric field. Figure 7 shows the current lines obtained by line-integrating the instantaneous current density vector j.
Roughly speaking, the current is almost field-aligned above 150 km altitude, whereas it is almost field-perpendicular below 150 km altitude. At a glance, the three-dimensional current lines depend on the original location of the FACs. The downward FACs labeled by 1a and 2a are connected with the horizontal current directing westward (negative x-direction) at relatively high altitude, followed by upward FACs. The downward FACs labeled by 1b and 2b are connected with the horizontal current directing north-west direction at relatively low altitude, followed by upward FACs. Note that the current lines 1b and 2b pass under the current lines 1a and 2a, respectively, near the origin (x ∼ 0 and y ∼ 0).
We compared the three-dimensional current lines with those calculated under the traditional thin-layer approximation. We chose the lower boundary of the region where the current is field-aligned, which is located at 203 km altitude. From the requirement of the current continuity, we solved the following elliptic partial where Σ, j || and L are the height-integrated conductivity tensor, the FAC (positive upward), and the thickness of the thin layer, which is set to 10 km, respectively. The line integral of j tla is drawn by the red line in Figure 8. The current lines under the thin-layer approximation are largely different from those solved by the three-dimensional simulation. The underpass cannot be seen in the thin-layer approximation.
This simulation automatically satisfies the condition that ∇•j = 0 because we neglect the displacement current. Thus, we can decompose the current j into the FACs j z (= j zˆz E e ), the Hall current j H and the Pedersen current j P . Considering the condition that ∇•(j z + j H + j P ) = 0, we can describe how the FACs are connected with the Pedersen and Hall currents. We focus only on the low altitudes where the polarization current is negligibly small in comparison with the conduction current. Figure 9 shows the isosurfaces of ∇•j z , ∇•j H , and ∇•j P . In Case 1, in most region, ∇•j z + ∇•j P ≅ 0, which implies that the FACs are mostly connected with the Pedersen current. This is what is expected from the uniform ionospheric conductivity with the electrostatic processes. Interestingly, a small part of the FACs is connected with the Hall current, which is caused by the inductive processes (Yoshikawa & Itonaga, 2000). In Case 2, FACs are connected with the Pedersen current as well as the Hall current. The connection between the FACs and the Hall current is significant near x = 0 and near the edge of the high-density band (y ∼ ±50 km). This is consistent with the expectation under the nonuniform ionospheric conductivity described by Equation 4. Figure 10 shows the z-component of the current density j z . The blue surface represents the downward current (negative j z ), and the red one represents the upward current (positive j z ). The most significant difference between Case 1 and Case 2 is found at low altitudes (below 150 km altitude), where additional FACs are formed as shown in Figure 3f. The structure of the FAC distribution is like the toe of a boot. To understand the formation of these localized FACs at low altitudes, we introduced the following equation that describes the change in the FACs (Song & Lysak, 2001).  (Note that the unit vector of the parallel direction || is originally given by B/B (Song & Lysak, 2001). In this study, we define the positive value of j || to be the upward FAC, and the unit vector of the parallel direction to be given by −B/B. Equation 19 is valid regardless of the polarity of j || .) Equation 19 shows that FACs change when the gradient of ∇•E has a gradient in the field-aligned direction. The yellow and green surfaces in Figure 10 represent the regions where the value of −∇ z (∇•E) is negative and positive, respectively. (We assumed that ∇ z ≅ ∇ || .) The yellow surface (where the value −∇ z (∇•E) is negative) is almost perfectly overlapped with the blue surface (where the value j z is negative). Figure 11 shows j z and −∇ z (∇•E)/μ 0 at x = 0 and y = −40 km as a function of z. FAC is almost absent in Case 1 (blue line), whereas FAC is present in Case 2 (red line). In Case 2, the FAC has a peak around 120 km altitude, which corresponds to the localized FAC as was viewed as the toe of a boot in Figures 3f and 4d. The height profile of j z is similar to that of −∇ || (∇•E), in particular, below 200 km altitude. This implies that the localized FACs peaking around 120 km altitude are locally formed. The field-aligned gradient of ∇•E efficiently arises from space charge deposited by the discontinuity of the Hall current at low altitudes. Figure 6 clearly demonstrates that the negative space charge is deposited around x ∼ 0 and y ∼ −50 km in Case 2. The deposition of the negative space charge gives rise to the suppression of ∇•E, resulting in the field-aligned gradient of ∇•E, and the formation of the localized FACs. Figure 12 shows isosurfaces of j z for Cases 1, 2, and 3. Case 3 is the same as Case 2 except that the width of the high-density band is doubled as expressed in Equation 15. The scale length of the flow shear in the y-direction L y is 100 km for Case 2 and Case 3. The localized FACs that are identified in Case 2 (as shown in Fig-Figure 10.  ure 12b) are absent in Case 3 (as shown in Figure 12c). Meanwhile, the isosurfaces of j z in Case 3 resemble those in Case 1. This result implies that the localized FACs appear when the horizontal scale length of the density gradient is smaller than that of the FACs.
j•E must be negative to generate FACs. With Equation 9, j•E yields as Equation 20 Figure 13 shows the distributions ] z in the horizontal plane at 100 km altitude. In Figure 13a, the regions where the polarity of rotational plasma velocity is opposite to that of the rotational Lorentz force are clearly overlapped with the localized FACs. By comparing Figure 13b with Figure 13c, we confirmed that the Lorentz force associated with the Hall current acts as the generation of the magnetic energy like a dynamo.

Discussion
Based on the two typical current lines shown in Figure 14a, we can depict the three-dimensional current vectors including the localized FAC as shown in Figure 14b. The blue vectors are composed of a current closure of the FACs by way of the Pedersen layer. The red vectors are composed of a closure by way of the Pedersen and the Hall layers. The current is the sum of the FAC and the horizontal current, so that the actual current lines are slant and smooth as shown in Figure 14a. It is interesting to note that the reddish line passes underneath the bluish line. This cannot be seen in the traditional thin-layer approximation.
In Figure 3f, we showed the additional, localized FACs, in particular, below 150 km altitude. The localized FACs are closely associated with the field-aligned gradient of ∇•E, which results from space charge deposited by the discontinuity of the Hall current. The generation of the localized FACs is a natural consequence of the discontinuity of the Hall current in a limited altitude range, and is found to bridge the Hall and the Pedersen layers. To generate FACs, the condition that j•E < 0 is basically necessary. At the edge of the high-density region, the convergent (divergent) Hall current creates the divergent (convergent) polarization field. The rotational plasma flow associated with the polarization field performs negative work against the Lorentz force associated with the Hall current. This process can be regarded as a dynamo. Since the localized FAC is almost confined in a limited altitude range, most of the magnetic energy generated in the Hall layer may be carried via the localized FACs, and may be dissipated in the Pedersen layer.
Our simulation results are basically consistent with the conventional Cowling channel model. In the Cowling channel, the Hall current associated with the primary electric field can deposit space charge near the edge of the channel, resulting in the polarization electric field. The Hall current associated with the polarization field intensifies the Pedersen current associated with the primary electric field. We believe that the three-dimensional model offers advantages for understanding the closure of the FACs. In the traditional thin-layer approximation, the FACs are, in general, connected with the Pedersen and the Hall currents flowing in the same plane. By simply summing up the Pedersen and the Hall currents flowing in the same plane, one may obtain the net horizontal current, and use it to trace the closure of the FACs. However, this approach is misleading. As shown in Figure 9, the three-dimensional closure of the current is fully different from the two-dimensional closure. To establish the three-dimensional closure, the additional FACs appear to bridge between the Pedersen and Hall layers. The presence of this type of the localized FACs was previously suggested by Fujii et al. (2011). We will investigate the role of the additional FACs on the overall development of the Cowling channel by introducing more realistic configuration in the future.
In Figure 9, the FACs are shown to be connected with the Hall current, even though the ionospheric density is uniform. The inductive processes should be involved to establish the connection. The inductive processes  are known to give rise to the coupling between shear Alfvén waves and compressional waves (Lysak, 1999;Lysak & Yoshikawa, 2006;Lysak et al., 2013). We confirmed the change in the magnetic energy near the Hall layer. However, the compressional waves are not expected to propagate efficiently since the wave frequency is below the cutoff frequency (Lee & Lysak, 1990). The contribution from the compressional wave to the current closure should be one of the important future works.
The additional FACs are expected to occur near a westward traveling surge (WTS) during the substorm expansion phase Inhester et al., 1981;Opgenoorth et al., 1983). Ebihara and Tanaka (2015a, 2015b, 2018 used a global MHD simulation coupled with the thin-layered ionosphere, and succeeded to reproduce the WTS in the global simulation. The additional upward FACs that manifest the WTS are associated with negative ∇•j H and positive ∇•j P as shown in Figure 13 of Ebihara and Tanaka (2015a). Namely, space charge deposited by the Hall current generates the additional FACs, which manifest the WTS. The result shown in Figure 13 is consistent with that obtained by the thin-layer model. With the three-dimensional simulation, the formation of the additional FACs becomes clearer than before. Figure 12 shows that the additional FACs extend to higher altitudes, but their fate is unknown because of the limited simulation box. Further studies are needed to understand the complete closure of the current. Tanaka (2015a, 2018) reproduced the WTS in the global MHD simulation by introducing the ionospheric conductivity depending on the FACs. We assume in this study that the rate of the production of plasma is independent of the FACs. This assumption makes the FACs distribution fairly stable. It is most likely that a surge-like structure of the upward FACs will be present when we improve the source term of the electron density S to be a function of the upward FACs.
It is evident from Figure 14 that the localized FACs are formed when the scale length of the flow shear L y is sufficiently longer than that of the high-dense region. That is, when the steep conductivity gradient is embedded in the FACs, such localized FACs are generated. Simultaneous observations of precipitating electrons and FACs show that some downwardly accelerated electron regions (inverted-V structures) are embedded in an upward FAC region (Fujii et al., 1994). According to simultaneous observations of precipitating electrons, FACs, and optical auroras, auroral structures are also embedded in the inverted-V structure (Fukuda et al., 2014). These observations suggest that the scale length of the ionospheric conductivity can be smaller than that of the FAC. Under such condition, the localized FACs are expected to appear frequently. The direct observations of the localized FAC may be difficult because significant air drag force impacting on a satellite makes it difficult to keep the satellite altitude below 150 km altitude.
An auroral emission sometimes shows double peaks in a vertical luminosity distribution of auroral arcs. This anomaly is called a two-tiered auroral band (Oguti, 1975) and enhanced aurora (Hallinan et al., 1985). Hallinan et al. (1997) proposed that the electrons are locally energized through wave-particle interactions or DC electric fields. The enhanced aurora may be related to the localized FACs seen in the conductivity gradient if the DC electric field is efficiently formed in the localized FACs. The simulation taking into account kinetics of the electrons, like Vlasov, or particle-in-cell simulation, is needed to investigate the formation of the parallel electric field and the consequent acceleration of the electrons in the localized FACs.

Conclusions
By performing a simplified three-dimensional Hall-MHD simulation, we investigated three-dimensional current closure between FACs and the ionospheric current. We obtained the following conclusions.
1. When the ionospheric density (conductivity) is uniform in the horizontal plane, most of the FACs are connected with the Pedersen current primarily due to the electrostatic process. A few FACs are connected with the Hall current due to inductive processes. 2. When the ionospheric density (conductivity) is nonuniform in the horizontal plane, most of the FACs are connected with the Pedersen and Hall currents. Additionally, localized FACs appear near the edge of the high-density region at low altitudes. The localized FACs are closely associated with the field-aligned gradient of ∇•E, which is caused by excess charge deposited by the Hall current. The localized FACs play a crucial role in connecting the Pedersen layer and the Hall layer. The polarity of the rotational component of the flow velocity is opposite to that of the Lorentz force associated with the Hall current in the localized FACs. This implies that the plasma can perform negative work against the Lorentz force (namely, acting as a dynamo), so as to generate the localized FACs. 3. The three-dimensional current lines are fully different from those described by the traditional thin-layer approximation. A current line can pass underneath another one, which cannot be described by the traditional thin-layer approximation.