Initiation of ballooning instability in the near-Earth plasma sheet prior to the 23 March 2007 THEMIS substorm expansion onset

In this work, we analyze the ballooning stability of the near-Earth plasma sheet prior to the initial expansion onset of the 23 March 2007 THEMIS substorm event. Using solar wind data from WIND satellite observation for the substorm event as an input at dayside, we reconstructed a sequence of global magnetospheric configurations around the expansion onset by means of OpenGGCM simulation. These simulations have reproduced most of the salient features, including the onset timing, observed in the THEMIS substorm event (Raeder et al., 2008). The ballooning instability criterion is evaluated using eigenvalue analyses for the near-Earth plasma sheet region when the configuration attains a quasistatic equilibrium condition prior to the onset. Our analysis of the evolution of the near-Earth magnetotail region during the substorm growth phase (10:10 UT to 10:43 UT) reveals a correlation between the breaching of the ballooning stability condition and the initial expansion onset in the temporal domain. The analysis indicates that the ballooning instability is at first initiated by the stretching of closed field lines near the region with a local nonzero minimum in normal component of magnetic field ( Bz>0) associated with the formation of a bipolar-flow pattern in the near-Earth plasma sheet around 10RE . The subsequent initiation of ballooning instability with enhanced growth rate occurs (about 28 min later) on those most stretched closed field lines Earthward of the neutral line following the formation and the tailward recession of a plasmoid.


Introduction
Conventional scenarios of the substorm onset trigger mechanism invoke either current sheet instabilities in the near-Earth magnetotail region or magnetic reconnection processes in the middle magnetotail region of the plasma sheet.Both scenarios have drawn upon suggestive but sometimes ambiguous observational and theoretical evidence (see, for example, Lui, 2003).This ambiguity is in part due to the lack of sufficient resolution in the observational data and numerical simulations.On the other hand, it appears that the trigger mechanism may actually involve multiple entangled processes (Voronkov et al., 2004).In this work, we examine one such possibility of a trigger mechanism that involves multistage processes.A ballooning stability analysis of the near-Earth plasma sheet prior to the initial expansion onset of the 23 March 2007 THEMIS substorm event is performed.
Most previous investigations of ballooning instabilities in the near-Earth magnetotail have been carried out under the assumption of a magnetostatic magnetosphere that is absent of either flow or reconnection (Miura et al., 1989;Hameiri et al., 1991;Lee and Wolf, 1992;Pu et al., 1992;Ohtani and Tamao, 1993;Hurricane et al., 1995Hurricane et al., , 1996Hurricane et al., , 1997;;Bhattacharjee et al., 1998;Cheng and Lui, 1998;Lee, 1999;Wong et al., 2001;Crabtree et al., 2003;Zhu et al., 2003;Cheng and Zaharia, 2004;Schindler and Birn, 2004;Zhu et al., 2004Zhu et al., , 2007)).In reality, the magnetosphere exhibits persistent convection and reconnection activities in the magnetotail region, which is often intermittent, as evidenced by the presence of bursty bulk flows (BBF) in both observations and simulations.The major processes in the near-Earth plasma sheet, ballooning, reconnection, and plasma bulk flow, each may have contributed indispensably to the onset of a substorm.The roles of each, and the interactions among them, may hold a key to understanding the trigger mechanism of the substorm onset.
P. Zhu et al.: Ballooning instability prior to THEMIS substorm onset Recent observations from the multi-satellite mission THEMIS of the 23 March 2007 substorm event (hereafter referred to as the THEMIS substorm event) provide a unique opportunity for us to examine more closely the roles of ballooning instability and magnetic reconnection in the triggering of a substorm onset in a dynamically evolving magnetospheric configuration (Angelopoulos et al., 2008;Keiling et al., 2008).Using WIND solar wind data for the THEMIS substorm event, a sequence of global magnetospheric configurations have been reconstructed around the substorm onset by means of OpenGGCM simulations (Raeder, 2003).These simulations have reproduced many salient features of the substorm growth phase, such as the field and plasma evolutions detected by THEMIS satellites in the magnetotail (Raeder et al., 2008).The configurations of the near-Earth plasma sheet region obtained from these simulations are found to be in quasi-static equilibrium prior to the initial expansion onset, which allows us to evaluate the conventional ideal MHD criterion for the ballooning instability of the configurations in an eigenvalue analysis.Our ballooning analysis of the near-Earth magnetotail configurations reconstructed from OpenGGCM simulations has demonstrated a close correlation between the breaching of the ballooning stability criterion and the timing of the initial expansion onset.The simulations and analysis indicate that the magnetic field line stretching associated with the formation of a bipolarflow pattern, and the subsequent formation of a neutral line region Earthward of a tailward receding plasmoid in the near-Earth plasma sheet, are able to set up and enhance the favorable condition for the initiation of a ballooning instability in that region.The breaching of ballooning instability boundary might directly contribute to the initial disruption of the current sheet.
In rest of this paper, we first describe the simulation setup (Sect.2) and the evolution of the near-Earth plasma sheet configurations prior to the initial expansion onset as obtained from the simulation (Sect.3).We then report the methods and results of our ballooning analysis of the pre-onset configurations (Sect.4).Finally, the major findings and implications are summarized and discussed (Sect.5).

Simulation setup
The OpenGGCM simulation used in this study covers a spatial domain of x : [20, −500], y : [−36, 36], z : [−36, 36] where the GSM coordinates are used and normalized by the Earth radius (denoted as "R E ").The size of the nonuniform Cartesian grid is 630×200×300.In the near-Earth plasma sheet region (x : [−5, −40], y : [−5, 5], z : [−5, 5], unit: Earth radius), which is of primary interest for this paper, the typical grid spacings are in the ranges of 0.145<( x)<0.147,y∼0.25, and z∼0.16.The grid spacings minimize around y=0 and z=0 planes, and increase monotonically in the tailward direction.The major simula-tion results remain qualitatively the same at higher spatial resolutions.To reconstruct the THEMIS 23 March 2007 substorm event, the WIND solar wind data at the upstream of bow shock is imposed as the dayside boundary conditions.The ionosphere boundary conditions are determined by the Coupled Thermosphere Ionosphere Model (CTIM) (Fuller-Rowell et al., 1996).

Simulation results
The maximum electron precipitation energy flux for discrete aurora from simulation indicates that the start of the first aurora intensification occurs shortly before 10:45 UT, consistent with the time series of the Knight potential and ground magnetic perturbation (Fig. 1).Whereas the timing of the first aurora intensification in simulation (which is around 10:40 UT∼10:45 UT) is somewhat earlier than the minor activation at 10:54 UT and the major activation at 11:19 UT as recorded by THEMIS ground and in-situ observations (Angelopoulos et al., 2008;Keiling et al., 2008), comparisons of many other features of the event, including the field and plasma evolutions during the substorm growth phase, have shown considerable agreement (Raeder et al., 2008).In this paper, we refer the timing of the first aurora intensification in simulation (10:40 UT∼10:45 UT) as the timing of the initial substorm expansion onset.As shown in Fig. 1, the solar wind condition upstream of the bow shock suggests the northward turning of the magnetic field at the day side magnetopause occurs later at about 11:10 UT (Raeder et al., 2008).
Four representative snapshots of the near-Earth magnetotail configurations from 10:10 UT to 10:43 UT prior to the first aurora intensification are described in the following.In particular we focus on the evolving pressure (p) and tailward flow (−u x ) patterns from their 2-D contours within the domain x : [−5, −40], z : [−5, 5] in the meridional plane (Fig. 2) and 1-D profiles along the y=z=0 line (Fig. 3).
Around 10:10 UT, the ionosphere and ground magnetic activity levels are low, and the near-Earth region (−x∼10−20 R E ) of the plasma sheet is moderately stretched, as seen from the pressure pattern.The plasma sheet stretching and thinning are further enhanced at about 10:15 UT, with the thinnest point located around 12R E .Soon afterwards, the pressure contours at 10:35 UT and 10:43 UT indicate the formation and recession of a plasmoid-like structure in the tailward direction.Before the formation of the plasmoid (10:10 UT and 10:15 UT), the pressure profile in the near-Earth region is flat along the x direction, but the peak value of the pressure starts to increase from 700 pP to 900 pP.After the plasmoid formation (10:35 UT and 10:43 UT), there is a dramatic steepening of the pressure profile in the region −x 10 R E , with the peak pressure value reaching 2300 pP.During the time period from 10:35 UT to 10:43 UT, however, the magnetic configuration in that region (−x 10 R E ) has become less stretched and more dipole-like (Fig. 2).
In the meanwhile, a weak bipolar flow, in both Earthward and tailward directions, starts to appear at 10:10 UT and becomes prominent at 10:15 UT, with the stagnation site located around x=−10 R E in the equatorial plane, as seen in the tailward flow contours (Fig. 2).The initial formation of the bipolar flow coincides with the thinning of the plasma sheet, and the initial stagnation site locates at the thinnest section of the plasma sheet near x=−10 R E .The stagnation site, along with the flow pattern, steadily recedes tailward, reaching x=−17 R E at 10:43 UT.The development of the bipolar flow pattern and the formation of plasmoid-like structure suggest the occurrence of a reconnection event in the near-Earth plasma sheet before the first substorm intensification.
The profile of the normal component of magnetic field (B z ) along x-axis (y=z=0), together with that of the tailward flow (−u x ), may be used to determine the timing and location of the reconnection event.As shown from the B z profiles in Fig. 4, a neutral line (as defined by B z =0) first appears around 10:35 UT and 16 R E , marking the start of a magnetic reconnection event in the near-Earth tail region.The neutral line and the stagnation sites do not share the same location, due to the asymmetry in outflows from the reconnection site as a consequence of the asymmetric distributions of magnetic field and pressure in x-direction.The emerging neutral line at 10:35 UT is accompanied by the formation of a plasmoid (Fig. 2).As the plasmoid recedes tailward, the stagnation site traces the tailward moving neutral line.However, the growing formation of the bipolar flow pattern around 10 R E started as early as 10:10 UT, well before the appearance of the neutral line and the plasmoid at 10:35 UT.As the bipolar flow grows from 10:10 UT to 10:15 UT, a dip forms in the B z profile, with the minimum B z located near the stagnation site.
The continuing decrease of the B z minimum eventually leads to the formation of neutral line and plasmoid at 10:35 UT.
One direct consequence of the formation of the bipolar flow pattern in the near-Earth region (−x 10 R E ) is a dramatic and subsequent pressure build-up in that same region from 10:35 UT to 10:43 UT, as induced by the Earthward flow from that bipolar flow pattern.The causal relation may be more evident from the time delay between the development of the two processes as seen in Fig. 3.At 10:15 UT, the Earthward flow is strong but the pressure build-up is weak.Later at 10:35 UT and 10:43 UT, the Earthward flow is weakened whereas the pressure build-up becomes stronger as a consequence of the continuous action from previous Earthward flow.The Earthward flow here is somewhat different from the steady magnetospheric convection described in Erickson (1992) and Hau and Voigt (1992).As a part of the bipolar flow that has spontaneously emerged in the near-Earth region, this Earthward flow seems to be a more localized and dynamic process, which occurs in addition to the adiabatic magnetospheric convection.
However, the pressure build-up in region −x 10 R E is also accompanied by a decrease in magnetic curvature due to the dipolarization induced by the Earthward flow.The more dipole-like configuration becomes more susceptible to a nearly complete cancellation between the interchange drive and the compressional stabilization.As a result, the pressure build-up does not directly contribute to the initiation of ballooning instability.As also shown in our analysis in the next section, it is the magnetic field line stretching associated with the formation of the bipolar flow pattern that gives rise to interchange drive and its dominance over the compressional stabilization, which results in the initiation of ballooning instability.
As the initial bipolar flow pattern moves tailward, it eventually develops into the outflows of a near-Earth reconnection

Ballooning analysis
We begin our analysis with an ideal MHD model of ballooning mode near marginal stability.In a 2-D magnetotail configuration described by B=ŷ×∇ (x, z), the ballooning mode along a closed field line can be described by the eigenmode equation (Hameiri et al., 1991;Lee and Wolf, 1992;Zhu et al., 2003) 2 where the line-tied boundary condition is imposed.Here ŷ is the unit vector in y-direction, (x, z) the magnetic flux function defining the configuration, l the coordinate along a field line, κ =κ•∇ /B denotes the normal component of magnetic curvature κ, B the magnitude of magnetic field, γ the adiabatic index, β=µ 0 p/B 2 , u A =B/ √ µ 0 ρ the local Alfvén speed, ρ the mass density, ξ the plasma displacement across magnetic surfaces, and the eigenvalue which is a measure of instability and approximates the linear growth rate of the ballooning mode.A =( dlA/B)/( dl/B) is the average along the field line l, with dl≡ N S dl, where "S" and "N" denote the locations of the field line foot points near the ionosphere.On the right hand side of Eq. ( 1), the first term represents the line bending strength, the second term measures the pressure gradient drive due to unfavorable magnetic curvature (which we refer to as "interchange drive" in the rest of this paper), and the last term is the contribution from plasma compression due to the line-tied boundary condition.Alternatively, Eq. ( 1) can be written in the following form: where Here 2 int , 2 ben , and 2 com represent the energy contributions to the ballooning growth 2 from interchange drive, linebending force, and compression, respectively.
The ballooning stability criterion in Eq. ( 1) or ( 2) is valid for near-Earth plasma sheet configurations that are quasistatic.An appropriate form of the quasi-static condition can be given by with the left hand side of the above inequality measuring (the square of) the ratio of the instability growth time −1 to the configuration evolution time τ eq .Here J is the plasma current, L eq ∼R E is the configuration scale length, and τ −2 eq ∼|J×B−∇p|/(ρL eq ).For ideal MHD instabilities, the growth time scale is estimated using −1 ∼τ A ∼L eq /u A .In the substorm growth phase of the near-tail region as reconstructed from the OpenGGCM simulation, the MHD ballooning instability time scale as estimated from τ A , as well as shown later in the ballooning analysis based on Eq. ( 1), is an order of magnitude faster than the configuration evolution time scale τ eq : −2 ∼τ 2 A ∼ 10 2 τ 2 eq ∼10 3 −10 4 .The quasi-static condition is well satisfied.
We numerically solve the eigenvalue Eq. ( 1) for the semiclosed flux tubes within the meridian plane (y=0) crossing the region 7 −x 20 in the equatorial plane (z=0) at each of the four time slices prior to the first aurora intensification, subject to the line-tied boundary condition at the ends of the flux tubes.In the eigenmode calculation, the end location of each field line is chosen to be around −x=2, beyond which (10 Fig. 5. Left column: Magnetic field lines within the meridian plane (y=0) crossing the region 7 −x 20 in the equatorial plane (z=0) at the four selected time slices prior to the first aurora intensification; Right column: Square of growth rate 2 (red sold line with the symbol "×"), the line-bending contribution 2 ben (blue short dashed line), the interchange contribution 2 int (magenta dotted line), and the compression contribution 2 com (black dot dashed line), for the ballooning modes along the magnetic field lines shown in the left column.Each field line is identified by the x-coordinate of its crossing location in the equatorial plane (z=0).In both columns, the green dashed lines represent the zero horizontal axes.
(−x<2) the eigenvalue does not vary much.Also for the convenience of the calculation, the mode is assumed to be symmetric about the peak point where B x =0 and the magnetic curvature κ is the maximum, since the south-north asymmetry in the field line about the peak point is expected to be rather small.A Newton iteration method is used in a shooting code to search for the most unstable eigenmodes.In Fig. 5, we plot those magnetic field lines within the meridian plane, along with the growth rates of the ballooning modes along the field lines, as a function of the x-coordinate of field line crossing location in the equatorial plane (z=0).Applying the ballooning mode eigenfunction to Eqs. (3-5), we also evaluate the energy contributions to ballooning instability from 2 int , 2 ben , and 2 com for each flux tube at each time slice.At all times, the near-Earth magnetotail is highly stable to ballooning modes in the region that is Earthward of the pressure profile peak due to the tailward-increasing pres-sure profile that is stable to interchange drive ( 2 int <0) and a strong stabilization from compression ( 2 com >0).The tail-side boundary of that region, which is the location of the maximum pressure profile, shifts from −x∼8.75 R E at 10:10 UT to −x∼7.5 R E by 10:43 UT.In this slowly and slightly shrinking region, the configuration supports stable ballooning mode oscillations in the ULF range of ∼10−100 s throughout the time slices of interest.The initiation of ballooning instability and the induced current disruption can only occur outside this region.
At 10:10 UT, the outside region 8.75 R E −x 20 R E is essentially marginally stable to ballooning mode as the result of a relatively flat pressure profile and an almost equally compensating stabilization from compression ( 2 com ∼ 2 int ).The mode profile broadly spreads along the field line, and the line-bending force is rather weak ( 2 ben ∼0).
P. Zhu et al.: Ballooning instability prior to THEMIS substorm onset At 10:15 UT, the region between 8.5 R E and 10 R E becomes weakly ballooning unstable ( ∼0.04 s −1 ).The increase in pressure gradient between 8.5 R E and 10 R E does not result in a proportional increase in 2 int of that region, because the magnetic field lines in that region are more dipolarized (with reduced magnetic curvature) as they move closer to Earth along with the Earthward flow.In contrast, the magnetic field lines in the region 10 R E −x 13 R E are more stretched, due to the continuing development of the bipolar flow pattern, consistent with the observation from the pressure contour (Fig. 2).As a consequence there is a dramatic increase in 2 int in that region.However, because 2 com ∝κ p, the change in 2 com profile due to the field line stretching and the bipolar flow is similar to that of 2 int .The two remain balanced in most of the region.Only in a narrow region between 10 R E and 11.5 R E , the interchange drive 2 int dominates the compression stabilization 2 com , and the flux tubes becomes considerably more unstable, with the peak value of ballooning growth rate ∼0.1 s −1 .The increase in magnetic curvature κ and hence the interchange drive 2 int due to the field line stretching is mostly responsible for the enhanced ballooning instability of that narrow region.
As the field line stretching continues following the tailward recession of the stagnation site of the bipolar flow, the peak of 2 int profile also moves tailward, reaching the location x=−14 at 10:35 UT.In the region −x<13 R E , interchange and compression energy nearly cancel each other resulting in a marginally ballooning stable state.In the region −x>17 R E , the ballooning mode is weakly unstable due to a weak interchange drive and nearly zero contributions from line-bending and compression.In the region between 13 R E and 17 R E , the interchange drive remains dominant with a peaked profile.Unlike the situation at 10:15 UT, the compression contribution 2 com becomes much less than the interchange drive 2 int due to an enhanced β level near the neutral line and more localized mode structure; the linebending contribution 2 ben becomes much stronger in this region, which compensates a large portion of 2 int yielding a moderate, peaked ballooning growth rate comparable to that at 10:15 UT.At 10:43 UT, near the initial onset, the ballooning instability is significantly enhanced on flux tubes in a very narrow region around 15 R E .For the flux tubes whose equatorial crossing locations are in the region 8 −x 14, the magnetic configuration has become more dipole-like, the interchange drive and the compression remain nearly cancelling each other, and the line-bending contribution is negligibly small.As a result this region is marginally (un)stable to ballooning mode.However, there is a sharp increase of the interchange drive for flux tubes in the narrow region 14 −x 15.2, easily surpassing the line-bending and compression contributions.This is the consequence of the highly stretched magnetic configuration in the narrow region Earthward of tailward moving plasmoid.The ballooning mode profiles on these flux tubes are highly localized around the middle point of the field line where the magnetic curvature κ is maximum, with a dominant interchange drive and a weak compression contribution, characteristic of the incompressible ballooning mode.They are in contrast to the ballooning mode structures on those dipole-like flux tubes where the mode profile broadly spreads along the field line; in those latter cases the interchange contribution is of the same magnitude as the compression, characteristic of the compressible ballooning mode.For all those closed flux tubes considered in the analysis at this time (10:43 UT), the maximum ballooning growth occurs on the flux tube crossing 15 R E , with ∼0.24 s −1 .Such a growth is also the maximum value obtained in the analysis during 10:10 UT to 10:43 UT.The e-folding time of the ballooning growth (∼4 s) is the same order as the system Alfvén time (τ A ∼10 s), which is consistent with the time scale of a typical substorm expansion onset.

Conclusions
In summary, an eigenvalue analysis of ballooning instabilities has been performed for the near-Earth plasma sheet configurations prior to the 23 March 2007 THEMIS substorm expansion onset.The configurations are obtained from an OpenGGCM simulation of the event.The development of the pre-onset plasma sheet configuration and their associated ballooning properties demonstrate a close correlation between the breaching of ballooning instability criterion and the initial onset of the substorm event in timing.The initial robust ballooning growth starts around 11 R E in a narrow near-Earth region.That ballooning unstable narrow region moves tailward following the recessions of a bipolar flow pattern and a subsequently formed plasmoid.The narrow region reaches the most ballooning unstable state near 15 R E when the configuration evolution approaches the initial expansion onset.The maximum ballooning growth rate is estimated to be Alfvénic, with a characteristic scale of about 0.2 s −1 .
The simulation and analysis indicate that the ballooning instability is at first initiated by the closed field line stretching near a minimum B z location that is associated with the formation of a bipolar flow along x around 10 to 11 R E in the near-Earth plasma sheet.The most unstable region coincides with the equatorial crossing locations of the most stretched flux tubes.The ballooning growth of the most unstable region is further enhanced on those most stretched closed field lines Earthward of the neutral line following the formation and the tailward recession of a plasmoid.
The development of ballooning instability on a thin current sheet was previously noted by Bhattacharjee et al. (1998) in a similar eigenvalue ballooning analysis of a series of near-Earth plasma sheet configurations obtained from a 2-D Hall MHD simulations under model solar wind conditions.There the current sheet thinning in the near-Earth plasma sheet was preceded by a neutral line in the equatorial plane approaching Earthward from the middle magnetotail region of plasma sheet.In contrast, in this work, the initiation of ballooning instability first occurs in a current sheet thinning region in near-Earth plasma sheet without the presence of an Earthward approaching neutral line from middle tail, even though the stretched thin current sheet is associated with a bipolar flow and a minimum B z that has originated directly from the near-Earth region (−x∼10 R E ).In addition, this work reports a second stage initiation of ballooning instability on highly stretched closed field lines following the formation of a neutral line and plasmoid that starts in the near-Earth region and moving tailward.Both major differences from the previous work in Bhattacharjee et al. (1998) may be due to differences in the model, dimensionality, domain, and boundary conditions that define the simulations.
The significance of the formation of a B z minimum in the near-Earth plasma sheet was previously noted by Erickson (1992).The physics mechanism underlying the initial formations of the bipolar flow and the B z minimum associated with the closed field line stretching in the near-Earth plasma sheet (−x∼10 R E ) remains unclear.The bipolar flow pattern was interpreted as a consequence of an intrinsic resistive tearing instability in response to the enhancement in local anomalous resistivity due to field line stretching by Zhu et al. (2008a,b).Such an interpretation is not conclusive, since it is not determined whether the minimum B z is still strong enough to prohibit the tearing instability.Siscoe et al. (2008) also observed the formation of the bipolar flow in the near-Earth plasma sheet prior to a substorm onset in their global MHD simulations, and they interpret the process as a certain macroscopic instability or a loss of equilibrium in analogy to the model of corona mass ejection, similar to the scenario previously proposed by Hau and Voigt (1992).Further study of this process is planned for future work.

PFig. 2 .
Fig. 2. Pressure (left column) and tailward flow (right column) contours in the meridional plane at four time slices prior to the first aurora intensification.

Fig. 3 .
Fig. 3. Tailward flow (left column) and pressure (right column) profiles along line y=z=0 at four time slices prior to the first aurora intensification.In the left column, the green dashed lines represent the zero horizontal axes.

Fig. 4 .
Fig. 4. Tailward flow (left column) and normal (z-) component of the magnetic field (right column) profiles along line y=z=0 at four time slices prior to the first aurora intensification.In both columns, the green dashed lines represent the zero horizontal axes.