Geophysicae Simulation of the interchange instability in a magnetospheric substorm site

We perform modeling of the interchange instability driven by longitudinal pressure asymmetry in the region of the pressure buildup that forms in the inner magnetosphere at the substorm growth phase. The simulation refers to the dawnward side of the Harang discontinuity and times after Bz IMF turning northward. The solution for the equilibrium state indicates tailward flows associated with vortices, which is in agreement with a previous finding of Ashour-Abdalla et al. (1999, 2002). We show that in the regions of equilibrium field-aligned currents (FACs), small initial perturbations in pV γ (p is the isotropic plasma pressure, V is the unit magnetic flux tube volume, γ=5/3 the adiabatic exponent), set up as ripples inclined to azimuth, grow in time. For the background FAC of∼10−6 A/m2, the linear growth rate of the instability is ∼ 6 min. Starting from the 12th min of evolution, the perturbations exhibit nonlinear deformations, develop undulations and front steepening. An interesting peculiarity in the distribution of the associated small-scale FACs is that they become asymmetric with time. Specifically, the downward currents are more localised, reaching densities up to 15×10−6 A/m2 at the nonlinear stage. The upward FACs are more dispersed. When large enough, these currents are likely to produce the aurora. We also run our simulation for the initial perturbations of large transverse scales in order to demonstrate that the interchange instability can be responsible for pressure and cross-tail current spatial variations of great extent.


Introduction
Efforts of numerous studies are directed toward finding an internal magnetospheric instability that could lead to onset arc formation followed by auroral breakup (Lui, 1991;Ohtani Correspondence to: I. Golovchanskaya (golovchanskaya@pgi.kolasc.net.ru) et al., 1999;Lyons et al., 2002;Voronkov et al., 2002).It is generally accepted that enhanced plasma transport to the Earth during the substorm growth phase results in a significant pressure buildup in the near tail at premidnight-to-postmidnight MLTs (e.g.Wang et al., 2004).Therefore, it is not surprising that the possibility of near-Earth tail destabilisation by the ballooning/interchange motions, which tend to destroy the large pressure gradients that form in this region, has been repeatedly addressed.Different modifications of the ballooning and interchange instabilities are considered (Swift, 1967;Roux et al., 1991;Ivanov et al., 1992;Ohtani and Tamao, 1993;Liu, 1997;Golovchanskaya et al., 2004).The stability of the radial pressure profile near geosynchronous orbit and farther down the tail was tested statistically by Ohtani and Tamao (1993).They showed that at the substorm growth phase this profile does not develop earthward pressure gradients large enough to satisfy the k p −γ (k b +k c ) >0 criterion for the ballooning instability (k p and k b are, respectively, the reverse scales of pressure and the magnetic field e-fold increase toward the Earth, k c is the magnetic curvature) and refuted the basic ballooning as a possible substorm trigger.The radial pressure profile also appears to be stable to interchange perturbations, which is signified by pV γ increasing tailward for both quiet and stretched magnetospheric configurations (Lee and Wolf, 1992).This implies that the manner of pressure buildup disruption shown in Fig. 1 (left) is hardly the case.
Later, Liu (1997) revised the ballooning mechanism with regard to substorm triggering.He performed a more rigorous analysis of the MHD equations for hot inhomogeneous plasma in a curvilinear magnetic field.One more ballooningtype mode was revealed, which, unlike the basic ballooning, is essentially related to the field-aligned velocity perturbations and thus has more in common with a slow magnetosonic wave than with a shear Alfvén wave.This mode has a looser destabilisation criterion and suggests rapid magnetic reconfigurations (on a time scale of ≤1 min) that may be relevant to the explosive growth phase of the substorm (Liu, 1997).
Published by Copernicus GmbH on behalf of the European Geosciences Union.Alternatively, in Ivanov et al. (1992), Golovchanskaya and Maltsev (2003) an interchange modification driven by longitudinal asymmetry in plasma pressure distribution (Fig. 1, right) has been proposed (see also Liu, 1996).This mode, whose evolution in the substorm site is a subject of the present study, turns unstable whenever the magnetic tension force, acting on a hot plasma in the curvilinear magnetic field, and plasma pressure gradient are not strictly antiparallel.The measure of the directional displacement between these two vectors is a background field-aligned current (FAC) generated by the Vasyliunas (1970) mechanism.For the substorm growth phase the presence of the large-scale FACs for night MLTs is well established (Iijima et al., 1993;Watanabe and Iijima, 1993) (Fig. 2).In the near tail these currents have a sense of the Region 2 FACs, with the average intensities in the ionosphere increased up to 10 −6 A/m 2 .
The plasma configuration with displaced ∇p and ∇V vectors is always unstable to the perturbations in pV γ , set up as ripples lying between the contours of V and pV γ (Golovchanskaya and Maltsev, 2003).Growing in amplitude interchange features extend at a certain (typically not very large) angle to azimuth.The accompanied small-scale FACs can reach the densities sufficient to produce the field-aligned electron acceleration into the ionosphere and arc formation.
The interchange instability driven by a displacement of ∇p and ∇V vectors has a close similarity to the atmospheric baroclinic instability.The latter develops whenever the atmospheric pressure gradient becomes inclined to the gravity direction, owing, for example, to horizontally nonuniform heating.
Further in Sect. 2 we describe the equilibrium configuration, which then will be subject to interchange destabilizing.The temporal evolution of ripple-type perturbations, including the nonlinear stage, is considered in Sect.3. The results of the simulation with regard to substorm physics are discussed in Sect. 4.

Finding equilibrium state
It is well known that in a general case there is no stationary solution that would describe self-consistently the magnetospheric convection, pressure distribution, large-scale FACs and magnetic field (Erickson and Wolf, 1980;Wolf et al., 1991).At the same time, setting the non-stationary or nonanalytical equilibrium is an obvious inconvenience for instability simulation.Therefore, we restrict our treatment by a specific equilibrium state, still relevant to substorm onset, for which it is possible, at least locally, to bring in a reasonable consistence the pressure gradients, magnetic volume, FACs, and plasma flow in the near tail.Namely, taking into account that substorm onsets are typically preceded by the B z IMF turning northward (Lyons, 1996), we address the times when the Region 1 Birkeland FACs, being directly driven, as well as the associated earthward convection, are essentially reduced (or ceased) in response to the B z IMF turning northward.On the other hand, the Region 2-type FACs, which arise in the nighttime near-Earth tail during the growth phase (Fig. 2), are related to the pressure gradients developing here and thus have much a larger relaxation time of ∼1 h (Rezhenov et al., 1980;De Zeeuw et al., 2004).During this period they solely control the convection in the inner magnetosphere.We suggest that these currents are generated by the Vasyliunas (1970) mechanism and further in this section try to adjust the distributions of other equilibrium quantities to the observed FAC densities and directions.
Having chosen a simulation domain as a 4 R E ×4 R E square site in the magnetospheric equatorial plane that extends approximately from geosynchronous orbit to 11 R E down the tail, we direct the x-axis opposite to the ∇V vector, the z-axis vertically upward and the y-axis so that it completes the right-hand coordinate system, that is, the simulation domain is = {−11 < x < −7 , −4 < y < 0 , z = 0 } , where x, y, z are in R E .We set 1-D model distributions of the equatorial magnetic field B z (x) and magnetic volume V (x) (Fig. 3a), typical for the growth phase of a moderate substorm.Our simulation is non-self-consistent with respect to the magnetic field, i.e.B z and V remain unchanged in the course of instability evolution (formally this corresponds to the low plasma parameter β=2µ 0 p B 2 , where µ 0 is the magnetic permeability of a vacuum).A possible explanation for such a limitation is that, in our logics, the final cause of pressure/tail current structuring, leading to onset arc formation, is the instability, whereas magnetic variations is a consequence which can influence but cannot govern instability evolution.
With B z (x) and V (x) given, it is easy to express all other equilibrium quantities, i.e. the electric potential ϕ (x , y) (whose contours are coincident with those of pV γ ), plasma pressure p (x , y), field-aligned current j || (x) and convection velocity v (x , y) in terms of V (x), B z (x) and eight input constants: p , B 0 , V 0 , J 0 , α , α B , α p , θ 0 .Here p is the height-integrated ionospheric Pedersen conductance, B 0 , V 0 , J 0 are the scales of the equilibrium magnetic field, unit magnetic flux tube volume and FAC, respectively.
The equilibrium distributions have been found as a selfconsistent analytical solution of the following set of equations where Eq. ( 1) is the Vasyliunas (1970) expression for the FAC, Eq. ( 2) is the current continuity requirement, Eq. ( 3) describes the adiabatic transport in a stationary equilibrium.The potential electric field E is determined by Eq. ( 4), with the convection velocity given by Eq. ( 5).In Eq. ( 5) e z is a unit vector along the z-axis.We note that the gradients of p, V , ϕ, and pV γ , entering the governing Eqs. ( 1)-( 5), are presumed to be across the magnetic field, since none of these quantities change along the magnetic lines.Then, the problem is reduced to a 2-D one and can be treated in any plane perpendicular to the magnetic field (further we choose the magnetospheric equatorial plane), accounting for the convergence of the magnetic field tubes.We refer to the MLT sector dawnward from the Harang discontinuity and thus choose the upward (away from the ionosphere) direction for j || (x) (see Fig. 2).
The model distributions for B z (x), B z in nT, and V (x), V in R E nT , R E being a numerical value equal to 6400 km, are taken as  Figure 3a shows V (x) normalised over V 0 for α = 6.The solution for the large-scale FAC (Fig. 3b) has the form The FAC is negative, i.e. flowing into the equatorial plane (away from the ionosphere).Its density is maximum at the near-Earth boundary and subsides with distance down the tail.
The equilibrium electric potential is described as As mentioned above, this is the potential of the remnant polarization electric field, or the "overshielding" field, as termed by De Zeeuw et al. (2004), persisting in the inner magnetosphere after a northward swing of B z IMF and a reduction in the earthward convection.The contours of ϕ(x, y) shown in Fig. 3d outline a part of the convection vortex that forms dawnward of the Harang discontinuity during this time period.It is seen that in the considered region the plasma flow is tailward and toward the dawn flank.
One should expect that duskward of the Harang discontinuity the flow is tailward and toward the dusk flank, since there the background FAC is of the opposite sign (Fig. 2).The tailward direction of the flow near local midnight is consistent with a previous finding of Ashour-Abdalla et al. (1999,2002), who reported on the reversed flow in the inner magnetosphere at substorm onset and related this convection state to the high near-Earth plasma pressure.A mirror symmetry of the convection with respect to the Harang discontinuity is well documented for ionospheric auroral features during breakup (Akasofu, 1964).
The solution for equilibrium pressure distribution (Fig. 3c) is given by From comparing Figs.3a and c, it is seen that because of the pressure increase in the y direction, the gradients of the magnetic volume and plasma pressure are not strictly antiparallel.In accordance with Eq. ( 1), this is consistent with the presence of the field-aligned current.
The distribution of pV γ (Fig. 3e) indicates that the equilibrium is stable against the basic interchange, as pV γ increases in the direction of ∇V .Therefore, only the interchange modification, which is shown in Fig. 1 (right), can lead to destabilization.
To estimate the characteristic quantities, we adopt the following values of the constants for a moderate disturbance: B 0 =40 nT, V 0 =0.18R E /nT, J 0 =1.08 nA/m 2 (referred to the magnetospheric equatorial plane), p =6 S, α=6, α B =3, α p =10, θ 0 =π /3.Then Eqs. ( 11) and ( 12) yield 0 =66 kV, P 0 =2.3 nPa, A 0 ≈0.5, and the plasma parameter appears to be β≈3.6.The obtained characteristic pressure P 0 is by a factor of 3.3 larger than the pressure at the same distances under quiet conditions, the latter being ≈0.7 nPa (Lui and Hamilton, 1992).We note that for strong substorms, the βparameter in the near tail can be greater than that adopted in our simulation by a factor of 1.5-2.The characteristic convection velocity v in the substorm site turns out to be about 130 km/s (or ∼ 4 km/s in the ionosphere) in the negative x direction, and ∼ 70 km/s in the negative y direction.These velocities can be significantly reduced by particle precipitation into the ionosphere, which leads to an enhancement in the ionospheric conductance which we ignore in the present treatment.
The equilibrium that we have found is an idealization in many respects.One point is that we neglect the influence of the FAC of opposite polarity, flowing in the dusk MLT sector, on the convection pattern in the simulation site.However, the lack of an analytical solution for the equilibrium in case this effect is included prevents us from setting a more complicated pressure and FAC distributions for the background.Generally, finding a more accurate local equilibrium for the substorm site is an independent problem, which has to be addressed in the future.
We complete this section with the reference to the study of De Zeeuw et al. (2004), in which the physics of duskto-dawn electric field formation in the inner magnetosphere in response to B z IMF turning northward (a so-called overshielding effect) is considered in detail.

Temporal evolution of interchange perturbations
For modeling the temporal evolution of baroclinic-type interchange perturbations (Fig. 1, right), we use Eqs.( 1), ( 2), ( 4), ( 5) and a nonstationary adiabatic transport equation It should be noted that the interchange motions, by definition, keep the form of the magnetic tube unchanged, that is, allowing for the metric factor, the displacement of the ionospheric end of the tube should be the same as that of any other point of the tube.Consequently, the field-aligned derivative of the potential ϕ should be zero everywhere along the tube, including its ionospheric end.This takes the role of the ionospheric boundary condition for our problem.
Representing the electric potential ϕ(x, y, t), plasma pressure p(x, y, t), field-aligned current j || (x, y, t) and convection velocity v(x, y, t) as a sum of the stationary equilibrium solution (hereafter designated by subscript 0) and a perturbed quantity (designated by subscript 1), we substitute these to Eqs. ( 1), ( 2), ( 4), ( 5), ( 13) and obtain the following equations for the perturbed quantities: in which the perturbed potential ϕ 1 satisfies the boundary condition In the equilibrium found in Sect. 2 a small single depletion in pV γ is set up which has a 2-D Gaussian form (with the amplitude of 5% at maximum) and goes between the flow contours and the volume contours (Fig. 4, left column, t=0).Such an orientation is known to be in favor for the growth of a baroclinic-type interchange (Golovchanskaya and Maltsev, 2003).We reproduce the temporal evolution of the perturbation by solving numerically the set of 2-D equations, Eqs. ( 14)-( 16), with boundary condition Eq. ( 17).The mesh width is h=25 km and the difference approximation of the fourth order accuracy is used.
The main difficulty that we encountered in the simulation is that the nonuniform, large-scale convection tends to remove the perturbation from the simulation site, whose size cannot be made much larger, considering the high resolution required and limited computational abilities.It can be easily estimated that the perturbation remains in the simulation region shown in Fig. 3 for only about 6 minutes.Therefore, we roughly suggested that the perturbation evolves in the background, which characteristics are spatially averaged characteristics of the equilibrium found above, and centered it in the middle of the considered region (formally v 0 in Eq. ( 15) was taken to be zero).
Figure 4 shows the temporal evolution of an elongated bubble (left column) with the associated perturbed electric potential ϕ 1 in kV (middle column) and FACs in nA/m 2 (right column, the light (dark) grey corresponding to the FACs into (away) from the ionosphere).The transverse scale of the per- turbation is rather small, roughly corresponding to that of a discrete auroral arc when reduced to the ionosphere.One can see that by the 18th min of evolution the perturbation p 1 V γ has developed nonlinear deformations, exhibiting undulations and front steepening.The small-scale FACs rapidly increase in intensity, becoming essentially asymmetric and displaying distortions at the nonlinear stage.The downward FAC is more localized and has larger local densities compared to the upward FAC.
Figure 5 illustrates the temporal evolution of a bipolar perturbation in pV γ , having a much larger transverse scale.By this run we intend to demonstrate that the interchange instability is capable of producing pressure (and cross-tail current) inhomogeneities which are large in extent.Finally, we compare the time scales of perturbation evolution inferred from the simulation with the linear growth times of the instability under study.The linear analysis of Eqs. ( 1), ( 2), ( 4), ( 5), (13) yields the following expression for the growth rate (Golovchanskaya and Maltsev, 2003 and references therein) From Eq. ( 18) it is seen that increases with equilibrium FAC increasing and drops in case the ionospheric conductance is enhanced.For the characteristic values of j || , P and B z adopted in Sect.2, the growth rate can be estimated as ∼2.6•10 −3 s −1 , with the corresponding growth time being ∼6 min.In our simulation, the time scale of perturbation growth at the beginning of the evolution is roughly consistent with this estimate.

Discussion
Developing a pre-substorm equilibrium plasma state that has the proper properties in the region of impending onset is of primary concern in modeling instabilities relevant to the substorm.In the present study we focus on the overshielding electric fields and associated convection pattern in the substorm site after the northward turning of B z IMF.
It is known that in the transition from weak to enhanced earthward convection during the growth phase, the inner plasma sheet develops pressure gradients, whose polarization by the gradient-curvature current produces field-aligned currents and a shielding electric field.This field prevents the plasma from penetrating into the regions of strong magnetic field.The formation of the appropriate pressure gradients takes ∼30-60 min, with the shielding being ineffective during this period (Wang et al., 2004).Then, it is very reasonable to suggest that the relaxation of these pressure gradients after the northward turning of B z IMF and the reduction in the sunward convection, occurs on the same time scale, inferring the persistence of the shielding (dusk-to-dawn) electric field for about an hour.The validity of this suggestion has been proved in the computer simulation of De Zeeuw et al. (2004).Using the RCM MHD code, the authors showed that about 10 min after the northward turning of the IMF reaches the magnetopause the inner magnetosphere exhibits an overshielding electric field which lasts just over an hour.
The equilibrium solution that we found in Sect. 2 refers to this time interval.The observed distribution of the fieldaligned currents in the near tail at the growth phase (Fig. 2) presupposes the formation of a large-scale vortical nightside convection pattern after a reduction in the dawn-to-dusk electric field.This feature is reproduced in our solution.The equilibrium streamlines shown in Fig. 3d indicate the dawnward diversion of the tailward flow on the dawnside.The duskward diversion of the flow is expected on the duskside due to the opposite polarity of the FAC there.Since the overshielding electric fields are static, they transfer to the ionosphere, manifesting in the dynamics of the auroral substorm.The convection pattern in the magnetosphere can be much more complicated due to the induced electric fields, which are ignored in our treatment.However, previously, Ashour-Abdalla et al. (1999, 2002), using a combination of spacecraft observations and global MHD simulations, noted the presence of tailward flows associated with vortices in the inner magnetosphere at substorm onset.The authors related this finding to the high near-Earth plasma pressure which limits or partially blocks the sunward convection.
The major importance of the interchange modification considered in Sect. 3 for substorm physics is that it is capable of redistributing the cross-tail current along the tail (see also Liu, 1996).On small transverse scales this leads to the formation of multiple arc-like structures.On large scales the instability can coherently produce giant domains of pressure enhancements alternating with pressure depletions.Multiple current circuits that set up due to the instability may be relevant to multiple substorm onsets (e.g.Wiens and Rostoker, 1975;Ohtani et al., 2004).
Our theory predicts electric field variations of large amplitude in a substorm site.This prediction is consistent with the empirical scenario of substorm onset developed by Maynard et al. (1996) from the CRRES observations of magnetic and electric fields in the inner magnetosphere near times of substorm onsets.According to those authors, the process grows from ripples at the inner edge of the plasma sheet, associated with dusk-dawn excursions of the electric field, prior to the beginning of dipolarization.
The background FAC is nonuniform along the tail (Fig. 3b).The onset arc, while having exactly the same genesis as any other arc, i.e. being a result of an intense largescale FAC disruption into smaller-scale FAC sheets, is most likely coincident with this FAC maximum and thereby first intensifies, brightens and deforms.
Our modelling gives evidence of a nearly identical interchange evolution on different transverse scales.This can explain the close similarity in dynamics of auroral forms with transverse sizes ranging from 1 km to 10 3 km (Oguti, 1975).
To reproduce the auroral electrojet in further treatments, we should include precipitation and strongly enhanced ionospheric conductance associated with the upward FAC after it reaches significant densities (Fig. 5, right column).Electrojet filamentary structure, which is well known from observations (e.g.Wiens and Rostoker, 1975;Elphinstone, 1996), can be achieved in the simulation as a direct consequence of interchange destabilizing on different scales.Another nonadiabatic process that should be included in the simulation is magnetic tube heating in the regions of strong upward FACs.All these improvements are the subject of further work.
Throughout the paper it has been implied that it is the interchange instability that causes plasma sheet and auroral structuring.Other mechanisms should also be considered.In particular, the stability of flow velocity shears that arise in the inner magnetosphere in response to B z IMF northward turning should be carefully examined.

Fig. 1 .
Fig. 1.Disruption of the near-Earth pressure buildup at the late growth phase by the basic ballooning/interchange instability (left) and baroclinic-type interchange instability (right).Hot and dense plasma filaments are shown in grey.Colder and more tenuous plasma features are blank.The arrows indicate directions of velocity perturbations inside the interchange structures.The solid lines indicate pressure contours in the magnetospheric equatorial plane.

Fig. 2 .
Fig. 2. Large-scale FACs enhanced during the substorm growth phase in the pre-midnight-to-post-midnight MLT sector as observed by Iijima et al. (1993).

Fig. 4 .
Fig. 4. Time evolution of a single depletion in pV γ with a small transverse scale (left), associated perturbed electric potential ϕ 1 (middle) and FACs (right).
Mingalev et al.:Simulation of the interchange instability expect alternating regions of pressure depletion and enhancement, with the formation of multiple local current circuits.