The magnetotail reconnection region in a global MHD simulation

. This work investigates the nature and the role of magnetic reconnection in a global magnetohydrodynamic simulation of the magnetosphere. We use the Gumics-4 simulation to study reconnection that occurs in the near-Earth region of the current sheet in the magnetotail. We locate the current sheet surface and the magnetic x-line that appears when reconnection starts. We illustrate the difference between quiet and active states of the reconnection region: variations in such quantities as the current sheet thickness, plasma ﬂow velocities, and Poynting vector divergence are strong. A characteristic feature is strong asymmetry caused by non-perpendicular inﬂows. We determine the reconnection efﬁciency by the net rate of Poynting ﬂux into the reconnection region. The reconnection efﬁciency in the simulation is directly proportional to the energy ﬂux into the magnetosphere through the magnetopause: about half of all energy ﬂowing through the magnetosphere is converted from an electromagnetic into a mechanical form in the reconnection region. Thus, the tail reconnection that is central to the magnetospheric circulation is directly driven; the tail does not exhibit a cycle of storage and rapid release of magnetic energy. We ﬁnd similar behaviour of the tail in both synthetic and real event runs.


Introduction
Satellite observations during the second half of the past century provided cornerstone information for constructing a general picture on the basic structure and dynamics of the magnetosphere. While satellites remain indispensable as the only source of in-situ measurements, new tools are needed to complement the sparse measurements of the few Correspondence to: T. V. Laitinen (tiera.laitinen@helsinki.fi) spacecraft in the vast magnetosphere. An important class of these new tools are global simulations (e.g. Janhunen, 1996;Fedder and Lyon, 1995;Raeder, 1999). They can model the whole solar wind-magnetosphere-ionosphere system self-consistently, and they allow data acquisition everywhere simultaneously and continuously.
All global simulations of the magnetosphere are presently based on magnetohydrodynamics (MHD), since a kinetic description of the whole plasma system is still beyond the capacity of computers. MHD provides a good description of the plasmas of the solar wind and outer magnetosphere, for example, the location of the magnetopause; on the other hand, it is less succesful in the near-Earth region where many features depend on non-MHD physics, for example, multiple temperatures . Despite the shortcomings, the dynamical behaviour produced by global MHD simulations in many respects reproduces the observed one (e.g. Pulkkinen et al., 1998;. One of the less-explored questions in the correspondence of global MHD simulations with nature is the role played by magnetic reconnection. Reconnection itself, as a theoretical physical process, has been widely studied (e.g. Vasyliunas, 1975;Birn et al., 2001), and many applications have been proposed; an extensive introduction is given by Priest and Forbes (2000). As they define, "reconnection is essentially a topological restructuring of a magnetic field caused by a change in the connectivity of its field lines. This change allows the release of stored magnetic energy, which in many situations is the dominant source of free energy in a plasma".
The basic magnetospheric circulation model proposed by Dungey (1961) contains two neutral points (or neutral lines in three dimensions), one at the magnetopause and one in the tail. Reconnection is also an important ingredient in many modern models for magnetospheric dynamics, especially so in the near-Earth neutral line model for substorms, as reviewed by Baker et al. (1996). There is ample observational evidence for reconnection in the magnetosphere, both at the magnetopause (e.g. Chisham et al., 2004, and references therein), in the near-Earth and mid-tail (e.g. Nagai et al., 1998;Øieroset et al., 2001), and in the distant tail (e.g. Ho et al., 1994). The most directly observed signatures include simultaneous plasma velocity and magnetic field direction changes coincident with a density depletion, and Hall magnetic fields in the surroundings. Perhaps the most common signatures are fast plasma jets.
A global simulation based on ideal MHD is certainly not the right tool for investigating details of reconnection physics; apart from resolution issues, Birn et al. (2001) have shown that conventional MHD produces much slower reconnection than simulations based on the Hall MHD, hybrid or full particle approach. Furthermore, observations point to the existence of notable Hall effects (e.g. Øieroset et al., 2001), which are not included in ideal MHD. However, as MHD simulations are widely used in magnetospheric research, it is essential to understand how they handle such an important process as reconnection, and how that affects their overall performance. This will be our aim in this article.
We investigate near-Earth tail reconnection in the Gumics-4 global MHD simulation (Janhunen, 1996). First, we develop methods for locating and illustrating reconnection. We show that Gumics-4 does produce a phenomenon that we recognize as reconnection. We then define a quantitative measure for reconnection efficiency and compare its time development to solar wind parameters and energy flux through the magnetopause. Applying these methods to several simulation runs, we study the role of reconnection in the dynamics of the simulated magnetosphere. In addition to finding out what reconnection in a global MHD simulation looks like and how it behaves, we discuss how it affects the accuracy and capabilities of global MHD in modelling the behaviour of the magnetosphere.

Model description
In this investigation we use the Gumics-4 simulation code developed by Janhunen (1996). It is a three-dimensional global simulation of the coupled magnetosphere-ionosphere system, based on ideal magnetohydrodynamics in the magnetosphere and electrostatic equations in the ionosphere. Figure 1 shows an image of the magnetosphere as rendered by  Gumics-4 uses the geocentric ecliptic (GSE) coordinate system, which will be our basic coordinate system in this article. The simulation box is rectangular, extending from +32 to −224 R E in the x direction, and to ±64 R E in the y and z directions. As boundary conditions at the sunward edge of the simulation box the solar wind density, temperature and velocity, as well as the interplanetary magnetic field (IMF), are given as functions of time. On the other walls of the simulation box, outflow conditions are applied. The inner boundary the MHD simulation is a spherical shell with a radius of 3.7 R E centred at the Earth. Here the boundary conditions are determined by the coupling with the ionosphere, and by the Earth's magnetic dipole field.  1. The magnetosphere as rendered by Gumics-4. Colour coding on the yz-plane depicts plasma density; yellow lines are magnetic field lines and blue ones plasma flow lines. Behaviour of the grid can also be seen on the coloured plane.
The grid in Gumics-4 is an adaptive Cartesian octogrid with semiautomatic refinement. This means that the base grid is composed of cubic cells 8 R E across. Where better resolution is needed, the grid spacing is halved by splitting a cell into eight smaller cubes, and the splitting can be repeated recursively. While running, the simulation code regulates the degree of grid refinement automatically cell-by-cell, using steepness of gradients as the spatial resolution refinement criterion. However, there is a built-in preference to refine the grid more readily near the Earth. In all simulation runs of this study, the size of the smallest grid cells is 0.25 R E . The time-accurate code also uses temporal subcycling, which means that the time step is not the same in all cells, but is smaller near the Earth and in small cells, in order to save computation time while ensuring the fullfilment of the Courant-Friedrichs-Lewy condition (Courant et al., 1928), which states that the time step must be smaller than the information travel time across the simulation grid cell; in practice, this is the travel time of an Alfvén wave.
Gumics-4 solves the ideal MHD equations in their conservative form. No explicit resistivity term is added, and thus all diffusion in the simulation is of numerical origin. The code applies the finite volume method where quantities are treated as cell averages, and uses the Roe solver to solve the Riemann problem of calculating fluxes of conserved quantities through the cell interfaces (LeVeque, 1992). In the rare cases when the Roe method produces intermediate states with negative pressure, the Harten-Lax-van Leer (HLL) solver is used instead (Janhunen, 2000). Elliptic cleaning (Brackbill and Barnes, 1980) is used to maintain ∇·B=0, since this is not otherwise guaranteed by the finite volume method.
The ionosphere in Gumics-4 is modelled as a spherical shell at an altitude of 110 km. The space between the ionosphere and the inner boundary of the MHD domain at 3.7 R E is treated as a passive medium which only transmits the electric potential and field-aligned currents and where no currents flow perpendicular to the magnetic field. The field-aligned currents and electron precipitation originating from the magnetosphere are mapped to the ionosphere along dipole field lines. The precipitation power is used together with an empirical formula for solar UV radiation contribution (Moen and Brekke, 1993), to calculate height-integrated ionospheric Pedersen and Hall conductivities in a triangular finite element grid that covers the whole sphere and has a spacing of about 100 km in the auroral oval regions. The horizontal current distribution is then solved using the magnetosphereoriginated field-aligned currents as a source term. This procedure also gives the ionospheric potential, which is mapped back to the 3.7 R E shell and differentiated there to calculate the boundary condition and the transverse velocity through E×B, for the MHD magnetosphere.
In this investigation we use mainly synthetic events, i.e. simulation runs with artificial solar wind data. They allow us to discern the effects of different solar wind variables separately. We also include one real event, a moderate substorm on 15 August 2001, to establish the degree of applicability of our results to realistic situations.
Our four synthetic events form a set that we call the IMF rotation runs. In each of them, the IMF is initially directed northward; it then rotates through a full 360 • clock angle so that its direction changes by 10 • every 10 min. Hence, the duration of the runs after initialisation is 6 hours, plus a short period of northward IMF, to allow time for the last effects to propagate. All other solar wind parameters other than the IMF clock angle are constant throughout the runs. The four runs differ from each other by the values of the IMF magnitude and solar wind dynamic pressure. The dynamic pressure is either 2 nPa or 8 nPa, with the larger value being a compound effect of both larger velocity and larger density. (2 nP=m p ·7.47 cm −3 ·(400 km/s) 2 and 8 nP=m p ·13.3 cm −3 ·(600 km/s) 2 , where m p is the proton mass.) The IMF magnitude is either 5 nT or 10 nT. Each of the four possible combinations of these values is used in one of the IMF rotation runs. Figure 2 illustrates the classical two-dimensional x-point reconnection geometry, adapted to the near-Earth magnetospheric tail. (The figure is drawn from one of the IMF rotation runs.) The yellow colouring depicts the current sheet that separates the two tail lobes. The black lines are magnetic field lines and the blue dashed lines are plasma flow lines. The distortion in the magnetic field on the right is due to the terrestrial dipole. However, the most important difference in comparison with textbook pictures of reconnection is the non-perpendicularity of the inflows. The lobe plasma flows towards the reconnection region at quite a small angle relative to the current sheet. Therefore, the process lacks the symmetry in the x direction that is often assumed in theoretical reconnection studies. We will return to the consequences of this asymmetry later; first, we discuss methods for locating The familiar x-point geometry of classical reconnection theories can be recognized, but note the non-perpendicular inflows and the effect of the terrestrial dipole field on the right. the reconnection from simulation results, and for analyzing and quantifying the process.

The current sheet
As tail reconnection takes place on the current sheet between the tail lobes, our first step in finding the reconnection region is to locate the current sheet centre. For practical reasons we define it as a two-dimensional surface. Two alternative and physically relevant criteria exist: the surface of maximal current density or the surface of magnetic field direction change. We have applied both of them for a consistency check and conclude that they agree extremely well: in the thin current sheet region that we are interested in, the separation of the two surfaces produced by these two different criteria is generally less than the simulation grid cell size (0.25 R E in this region), i.e. below the limit of meaningful resolution. We have chosen to use the magnetic criterion in this analysis, because we will use the surface found as an intermediate step in locating the reconnection x-line, which is a feature of the magnetic configuration only. Moreover, the magnetic criterion turns out to produce a numerically smoother result. Mathematically, the current sheet surface is described as z(x, y). Its first approximation is located by finding the sign change of B x (z). The surface is then refined using components of B adapted to the surface, but the effect of this refinement turns out to be practically negligible. The procedure is described in more detail in Laitinen et al. (2004). From now on, we will refer to the surface thus found as the "the current sheet".
Because the current sheet is a physically relevant reference surface, it is useful to decompose vectors into one normal and two tangential components. For this purpose we introduce the adapted (x , y , n) coordinate base: the x axis is tangent to the current sheet and perpendicular to the y GSE axis, the y axis tangent to the current sheet and perpendicular to the x GSE axis, and n is normal to the current sheet. Note that the x and y axes are not necessarily perpendicular to each other, although in our application they are approximately so. In the uppermost panels of Figs. 4, 5 and 9 we show the current sheet as a projection onto the xy plane, but the arrows always depict the tangential part of the vector quantity, i.e. the x and y components. The middle and lowermost panels in the figures depict a plane perpendicular to the current sheet, thus providing a conventional side view of the reconnection region. The location and relative orientation of these surfaces are illustrated in Figure 3.

The neutral line
The focal point of classical reconnection models is the neutral line, where the inflowing plasma is divided into two oppositely directed outflows along the current sheet and the normal component of the magnetic field changes its sign. This characterization of the neutral line actually includes two definitions, that of the flow reversal line and that of the magnetic x-line. We locate these lines on the current sheet using the criterion B n =0 for the x-line and v x =0 for the flow reversal line. Both neutral line candidates are drawn in Figs. 4a and 5a. These figures are from the synthetic IMF rotation run with small pressure and large IMF magnitude, since this combination produces the clearest distinction between quiet and reconnecting times.
In Fig. 4a, which depicts the current sheet at a quiet time when the IMF is oriented northward, there is only the flow reversal line present and no magnetic x-line at all. Furthermore, the flow velocities in the region are so small, a few tens of kilometres per second, that even the flow reversal line is not very relevant. Figure 5a shows the same region when the IMF is oriented southward. There is now a smooth x-line passing through the region of highest current density, i.e., the thinnest part of the current sheet. The peak current density is  Fig. 4. The reconnection region in quiet time, i.e. before the reconnection starts. All panels depict the same instant, 2:00 in the synthetic IMF rotation run with small solar wind dynamic pressure and large IMF magnitude. The colour and vector scales are the same as in Figs. 5 and 9. (a) The current sheet seen from above, projected onto the xy plane in GSE coordinates. Coordinate units are Earth radiae. Colouring depicts current density, and the arrows are plasma flow velocity vectors. The blue dashed line is the flow reversal line. (b) The perpendicular plane (see Fig. 3), providing a side view of the same region. Colouring is again current density and the arrows are now magnetic field vectors on the plane. The magnetic field earthwards of x= − 10 is not drawn, since the arrows would be too long and clutter the image. (c) The perpendicular plane again. Now the colouring is the divergence of the Poyntig vector, ∇·S, and the arrows are plasma flow velocity vectors (nearly invisible due to the smallness of velocities at this time).
roughly tenfold compared to that of the quiet time, indicating that the change in the current sheet has been significant. The flow reversal line has conformed to the shape of the xline, and in the region of enhanced current density, there is a   constant separation of about 1.5 R E between them. Furthermore, the plasma outflow velocities have increased dramatically. Together these changes strongly suggest that we have found the reconnection site.
There is one detail, though, that does not fit the common expectations about reconnection: the separation between the magnetic x-line and the flow reversal line. In classical reconnection theory these lines always coincide. While it could be argued that the coincidence of the neutral lines in most reconnection studies is dictated just by the complete symmetry of the configuration, which is normally assumed in those studies but absent in ours; this feature does require closer inspection.
We interpret that the cause of neutral line separation in our simulations is the unsuitability of a discrete finite volume grid for representing accurately the flow pattern occurring here. As can be seen in Figs. 2 and 5c, the inflow to the reconnection region is fast and strongly inclined, having a large component in the −x direction. The outflow towards the Earth is relatively slow. The flow pattern is thus composed of a thin wedge of plasma moving slowly earthward, surrounded by plasma that is moving in the opposite direction at a very high speed. However, a simulation running in a cubic finite volume grid cannot produce a wedge thinner than the cell size. If we first created an ideal flow pattern like the one described, and then tried to save the velocity field values in a grid by calculating cell averages, the tip of the wedge would move radically due to the averaging. This is illustrated in Fig. 6. Of course, there is no underlying ideal flow pattern in a simulation, only the cell averages; but this example illustrates the limitations associated with a discrete grid. When the simulation cannot produce the kind of configuration that the boundary conditions of the region would ideally lead to, it seems to produce an approximation where one feature, the tip of the wedge, becomes misplaced.
One might now suppose that better resolution would bring the neutral lines closer to each other. Nevertheless, we have not deemed it worthwhile to run events with increased resolution. There are three important reasons for this. First, our grid cell size in the reconnection region is now 0.25 R E , and a 3 keV proton in a 5 nT field has a Larmor radius of the same size. These are not atypical numbers in the central plasma sheet. We are therefore already at the validity limit of the MHD approximation, where better resolution would no longer be physically meaningful. Second, at best, the neutral line separation would decrease proportional to the cell size; the discrepancy cannot be removed by scaling it down. Third, for our purposes it is enough that we can confidently recognize reconnection and determine where and when it is taking place, realizing that the accuracy in location along the current sheet is limited by the numerics.
Unlike the velocity field, the magnetic field around the xline varies in a simple, regular manner with no wedge-like gradient patterns. There should thus be no discretisationoriginating problems associated with the determination of the magnetic x-line location. We therefore conclude that the reconnection neutral line is best approximated by the magnetic x-line when it appears in a region of enhanced current density with a flow reversal line in proximity.
The behaviour of the magnetic neutral line reported here resembles very much that observed by Birn and Hesse (1991) in their resistive MHD simulation of the magnetotail. In both cases, the x-line forms at quite a large distance tailward from the flow reversal and then moves close to it when reconnection intensifies. The location of the neutral line (after the initial movement) is about x≈−15 R E in both simulations. Walker et al. (1993) also report tail reconnection at approximately the same location in their simulation. It is interesting to contrast this with observations, which generally place the near-Earth neutral line at −20 R E ≤x≤−30 R E (e.g. Nagai et al., 1998), but sometimes even at x≈−16 R E (Runov et al., 2003).

The noon-midnight meridional views
To illustrate the reconnection in a way that is as familiar as possible, we use a "side view" on a perpendicular plane (Fig. 3). This plane is defined so that it goes through that point on the magnetic x-line where y GSE =0. Furthermore, the plane is perpendicular to the current sheet at that point and parallel to the x GSE axis. One may imagine it as an xz plane tilted in the same way as the current sheet is tilted with respect to the xy plane.
The magnetic field and current density are shown on the perpendicular plane during a quiet time in Fig. 4b and during strong reconnection in Fig. 5b. It is evident from the figures that thinning of the current has taken place between the times depicted, and a notable change in the magnetic field geometry is also seen.
One of the most important properties of reconnection is its ability to release magnetic energy in the form of thermal and kinetic energies. Figures 4c and 5c illustrate this property by showing the divergence of the Poynting vector ∇·S colourcoded, together with the plasma velocity field depicted as arrows. Here the change from a quiet to a reconnecting state is profound. While the quiet state shows practically no energy conversion in the tail, during reconnection the whole thin current sheet region converts electromagnetic energy at a rate of roughly 10 −11 W/m 3 , or 10 −4 W/m 2 when integrated in the z direction. Plasma acceleration in the reconnection region is evidenced by a fast tailward flow reaching ∼3000 km/s. Maps of plasma temperature and pressure (not shown) reveal that considerable heating is also taking place in the region.
The earthward outflow is much weaker than the tailward one, and has a maximum velocity of about 700 km/s. A single value for inflow velocity cannot be given. Typical plasma inflow velocities at a perpendicular distance of 2 R E from the centre of the current sheet are of the order of 500 km/s, at an angle of about 66 degrees to the normal of the current sheet. The perpendicular component is thus about 200 km/s. However, to the earthward end of the reconnection region comes a high-speed stream exceeding 1000 km/s. The Alfvén speed in the region varies very strongly with the location, from a few hundred to several thousand km/s. However, all inflows are sub-Alfvén and both outflows super-Alfvén. The largest Alfvén Mach numbers are 8 for the tailward outflow and 1.3 for the earthward one.
Plasma density at the x-line is very low, about 0.02 protons/cm 3 . The lobe magnetic field below and above the reconnection site is about 20 nT. Finally, we emphasize that the numbers given above are to be taken as crude estimates only, as the values of all mentioned quantities are strongly varying functions of location and the choice of the most representative location is ambigious.
A curious feature in Figs. 5b, c is the bifurcation of the current sheet behind the reconnection region, beyond x=−20 R E . It might be tempting to associate this structure with the shocks of, for example, Petschek reconnection. But inspection of plasma densities reveals that these fronts are decompressive, i.e. plasma density decreases when passing through, while MHD shocks are compressive fronts. Moreover, this structure does not appear together with reconnection, but only after reconnection has attained considerable strength. If reconnection remains relatively weak, the tail does not bifurcate; this is seen in Fig. 9. This indicates that the structure in question is not an integral part of reconnection, but rather a consequence of it.
The bifurcated tail structure resembles strikingly a supersonic wake. Its three-dimensional structure, however, is more complicated than the pictures selected for this article can reveal. Only when solar wind B y ≈0 does the current sheet split symmetrically into two layers; otherwise, a more accurate description would be that the current sheet becomes folded such that its cross section in the yz plane resembles the letter S. The direction of the folding (S or Z) is determined by the sign of solar wind B y . The magnetic field thus plays a decisive role in this phenomenon, wherefore the analogy of a hydrodynamic wake is incomplete.
Current sheet bifurcation has been reported earlier both in observations (Asano et al., 2005;Runov et al., 2003) and in simulations (Birn and Hesse, 1996). The bifurcation in our simulation has a much larger size and is probably not the same phenomenon. Thorough treatment of this issue is, however, outside the scope of this article.  reault and Akasofu, 1978), and the most intense loading is expected to occurr when the cos θ curve reaches its minimum.
The middle panel in Fig. 7 shows the distances of the xline and the flow reversal line from the Earth, measured at y=0. These curves confirm what Figs. 4a and 5a suggest: that the flow reversal can always be found at roughly the same location, whereas the x-line only appears as a signature of reconnection. The flow pattern in the near-Earth tail thus stays qualitatively unchanged, and the location of the flow reversal is only slightly modified while other quantities undergo strong variations.
Unlike the flow pattern, the magnetic field structure undergoes a qualitative change when the x-line appears at the 02:40 simulation time and then disappears at 06:00. In addition, the x-line first moves towards the Earth as reconnection intensifies and then tailward when reconnection weakens. Similar behaviour is observed in all the simulation runs that we have looked at. However, the x-line distance from the Earth varies along the neutral line, and is thus a function of the y coordinate. In some runs the x-line is considerably more undulating on the scale of ∼1 R E than in the figures shown here. Furthermore, it may first appear in segments that then combine to form one cross-tail x-line. It also tends to disappear segmentwise. Although the behaviour of the xline is interesting in itself, it does not provide a sound basis for any numerical quantification of reconnection.
The uppermost panel in Fig. 7 shows several different quantities that might characterize the strength of reconnection. The quantities are scaled such that their time development can be compared. One prerequisite for reconnection is the formation of a thin current sheet, which is measured by the blue curves: the dark blue curve is the maximal current density in the region, and the light blue one is the area of that part of the current sheet where the current density exceeds 4 nA/m 2 . The total cross-tail current is determined by the lobe magnetic fields via Ampère's law, and does not vary greatly. The ten-folding of the current density is therefore mostly due to thinning of the current sheet. This thinning is restricted by the grid: the half-thickness of the current sheet rapidly shrinks to 1-2 grid cell sizes, i.e. the current sheet becomes as thin as it can. The value of maximum j thus "saturates", which manifests itself in the top of the dark blue curve being wide and flat. The size of the thin current sheet, measured by the light blue curve, is not as sensitive to grid properties but would also show saturation in other runs with stronger reconnection. It also involves the arbitrary choice of a treshold value, here 4 nA/m 2 .
An important consequence of reconnection is the accelerated plasma outflows. In the top panel of Fig. 7, the black curve shows the largest tailward flow speed on the current sheet, and the green one is the earthward flow speed at y=0, 2 R E , earthward of the flow reversal line. Both show at least a ten-fold increase, witnessing acceleration by reconnection. However, both are single-point measurements, and as proxies they also suffer from grid-dependent uncertainty, somewhat similar to that described in Section 3.2. As a numerical measure of the reconnection rate, we wish to find a number that represents globally the whole reconnection process, is not grid-dependent, and involves as little arbitrariness in the definition as possible.

The reconnection power: definition
Besides the actual reconnection of magnetic field lines, the most intrinsic and important property of reconnection is the conversion of magnetic energy into mechanical energy of the plasma. Locally, this is measured by ∇·S. Figure 5c reveals the spatial distribution of this energy conversion, and comparison with Fig. 4c shows the drastic difference between non-reconnective and reconnective states of the tail. The region of most intense energy transformation is quite sharply defined in that figure as a dark blue slab, which coincides well with the thin current sheet region.
A volume integral of the local energy conversion rate over the reconnection region gives the total energy conversion rate in the process. This is the quantity that we shall use as the numerical measure of reconnection, and we will call it "the reconnection power". It is the red curve in Fig. 7.
The only arbitrariness involved here lies in the choice of the integration volume. The strongest energy conversion is concentrated compactly in a bounded region, and the inclusion of surrounding regions with a relatively small conversion rate has no significant effect on the outcome, so we choose a simple rectangular box as the integration volume. The dimensions of the integration box are: −30 R E <x<−10 R E , −10 R E <y<10 R E and −6 R E <z<6 R E . It thus encompasses the whole thin current sheet region while staying clear of the magnetopause. The earthward, dawnward and duskward faces of the box are set to coincide with the respective edges of the thin current sheet. The tailward face is so far away that the box also includes a considerable portion of the wake-like structure behind the thin current sheet; as seen in Fig. 5c, this bifurcated part also exhibits energy conversion. The z range of the box is chosen so wide that the aforementioned structures stay within the box even when the tail or current sheet are tilted.
We have carried out the integration with several different boxes of different sizes. The results differ only by a constant factor; scaled appropriately, the resulting time development curves are virtually identical. Thus, except for the overall level, the reconnection power is not dependent on the integration volume chosen (as long as the volume includes the central part of the thin current sheet). By choosing a large box encompassing the whole energy conversion region, we obtain not only an indicator of the time development but also an estimate of the total amount of energy processed by reconnection.
By Gauß's theorem, the reconnection power may equally well be calculated as a surface integral. This method is numerically more robust, so we define Here ∂V is the boundary of the box described above, and a minus sign is included to obtain a positive number as a result.  The main result in Fig. 8 is that the reconnection power is directly proportional to the energy flux through the magnetopause. The constant of proportionality is about 0.5 in all the runs, although small differences of the order of 10% can be seen between the runs. The exact value of the constant of proportionality is irrelevant, as it depends on the size of the integration box in the definition of the reconnection power. But its order of magnitude is interesting, because it shows that reconnection has a central role in the energy budget: about half of the energy flowing through the magnetopause is converted from a magnetic into a mechanical form in the reconnection region. Even more of the energy passes through that region, as the reconnection power only measures energy conversion, i.e. it does not include the energy that goes through the reconnection region unconverted.
Both the tail reconnection power and the magnetopause energy flux are governed by the IMF direction within each run, but comparison of the runs shows that the IMF magnitude and the solar wind dynamic pressure also play an important role. The runs with a four-fold dynamic pressure (green and black curves) exhibit about a three-fold magnetopause energy flux and tail reconnection power, as compared to the other two runs. A larger IMF magnitude raises the overall level of these quantities only slightly, but it enhances their variation considerably: compare the black and blue curves (|B|=10 nT) to the green and red ones (|B|=5 nT), respectively.

Simulation of a substorm event
So far we have only looked at synthetic input data. Now we will apply the same methods to a simulation of a real event, a moderate substorm on 15 August 2001. The solar wind measurements were recorded by the Geotail spacecraft. The solar wind B z was around zero at the beginning of the simulated time period and turned weakly southward at about 03:39 UT. The solar wind dynamic pressure was small, below 1 nPa, during the event. CANOPUS magnetograms recorded the onset of a modest substorm (∼−500 nT at 04:27 UT), and at about the same time the parameter reached a peak value of 1·10 11 W. In the simulation, the onset occurred half an hour later, but the simulated substorm was faster, so that it recovered at about the same time as the observed one. A more detailed description of the event, including several data plots, can be found in Palmroth et al. (2004). Figure 9 shows the reconnection region at 05:20 UT when the reconnection power peaks. At the current sheet (panel a), the x-line has the same basic shape as in Fig. 5a, although it is more undulating on the small scale. This undulation is typical for the x-line in events with real solar wind data, and it can be considerably stronger than seen here. While our numerical search method necessarily produces a continuous line, in some cases a better physical interpretation would be to say that there are several separate stretches of the x-line that can be grouped together as one approximate x-line.
Due to a z component in the solar wind velocity, the tail as well as the current sheet are tilted with respect to the xy plane. The north-south asymmetry also manifests itself in different inflow velocities from the lobes to the reconnection region, as can be seen in Fig. 9c. However, this asymmetry does not affect the qualitative structure of the reconnection region. Figure 10 shows the reconnection power during the event, together with the magnetopause energy flux and the total energy consumption in the ionosphere. All three quantities behave very similarly. The magnetopause energy flux leads the other two by about 5 min, which only confirms that it is the driver of the system. The time resolution of 5 min does not allow one to discern a systematic temporal difference between reconnection and the ionosphere, so the question of their causal interdependence is left open. However, the maximum of the ionospheric energy consumption occurs as much as 15 min after that of the magnetopause and reconnection.
The magnetopause energy input exhibits sharp, short-time scale fluctuations. They may be caused by the motion of the magnetopause surface, as its location is not as stable in the real event runs as in the synthetic ones. On the other hand, these fluctuations are also mostly discernible also in tail reconnection and the ionosphere, but there they are strongly smoothed. This may be a sign of a buffer effect in the system, and since it affects the tail reconnection, as well as the ionosphere, it is a global rather than just as innermostmagnetospheric phenomenon. The energy buffer is not very large, as the big, longer-time scale fluctuation -the substorm -has the same relative intensity in all three places. Furthermore, no sign of loading-unloading behaviour of the tail can be seen.

Summary and discussion
We set out to study the role of magnetotail reconnection in the global MHD simulation Gumics-4. We began by locating a thin current sheet which is traversed by a magnetic x-line with flow reversal in proximity. On a plane perpendicular to  the current sheet, the magnetic field geometry is that of the x-line reconnection. The thin current sheet region is a sink of the Poynting vector and a site of acceleration and heating of plasma. Thus, this region not only looks like a site of reconnection with the right magnetic field geometry and plasma flow pattern, it also acts like one, converting magnetic energy into thermal and kinetic forms. All these facts taken together evidence conclusively that although Gumics-4 deals with ideal MHD where classical reconnection by definition cannot occur, numerical diffusion allows it to produce a phenomenon in the magnetotail that may be characterized as magnetic reconnection.
The role of numerical diffusion in MHD simulation has been the subject of a lively debate in recent years (see, e.g. Gombosi et al., 2000). Discretization of the ideal MHD equations always leads to numerical errors proportional to second derivatives of the magnetic field. Sharp boundary structures, such as current sheets, thus tend to become smeared out, i.e. the effect resembles that of physical resistivity, at least qualitatively. Numerical diffusivity is hard to quantify because it depends on several details of the numerical solution, for example, on the speed at which the structure propagates with respect to the Eulerian grid. A thorough treatment of this important question is unfortunately outside the scope of the present study. However, in different test runs with a smallest grid spacing of either 0.5 or 0.25 R E , the thickness of the current sheet is proportional to the cell size. Better resolution thus allows sharper gradients, which counteracts the smaller numerical diffusivity associated with smaller cell size. This suggests that while some of the local characteristics, such as maximal current density in the reconnection region, depend on the grid resolution and the amount of numerical diffusion, the global picture is not drastically affected. The reconnection site in Gumics-4 appears as a simple diffusion region with unsurprising x-point magnetic field geometry and a strong consumption of Poynting flux. There is no evidence of shocks, such as those in the Petschek reconnection model or those observed by Ho et al. (1994) in the magnetotail. Discrete grid produces some small-scale effects, such as the separation between the x-line and flow reversal, or a narrow minimum in |∇·S| at the centre of the diffusion region. This is not surprising, as the reconnection site is a place of steep gradients in several quantities; the code uses a discrete grid with limited resolution, and ideal MHD does not include all physics of real reconnection in space plasma.
Conventional reconnection theory often assumes perpendicular inflows and also otherwise complete symmetry in the x direction (see, e.g. Birn et al., 2001, and companion articles). In our example of reconnection this symmetry is absent due to the fast and strongly inclined inflows from the lobes. When the solar wind v z = 0, these inflows also have different speeds in the north and south lobes. The nearby Earth dipole field adds to the asymmetry. We suggest that in particular the non-perpendicularity of the inflows merits a closer study by other methods, for example, theoretical considerations or devoted reconnection simulations. Could it cause qualitative asymmetry in the configuration, such as shocks, on only one side of the diffusion region?
As a quantitative measure of reconnection efficiency, we chose the reconnection power, i.e. the energy conversion rate, which we calculate as the net Poynting flux into the reconnection region. A common measure of the reconnection rate is the electric field along the neutral line (Vasyliunas, 1975). It is not applicable in an ideal MHD simulation, because the electric field is a secondary quantity calculated as E=−v×B. The neutral line is defined as the magnetic x-line, and thus the magnetic field either vanishes or is parallel to it. The electric field then either vanishes or is perpendicular to the neutral line. This reflects the fact that the diffusion allowing reconnection in the simulation is numerical, not resistive. The electric field outside the diffusion region could be used to measure the rate of flux transfer towards reconnection, but the choice of the location where to determine the value would be ambiguous. Reconnection power, as defined in Section 4.2, provides a global measure of the process that is spatially quite extensive.
We then compared the reconnection power to the net energy flow into the magnetosphere through the magnetopause. Figure 11 summarizes our main result: the reconnection power is proportional to the energy input through the magnetopause. Half of the energy coming to the magnetosphere is converted from the electromagnetic into a mechanical form by tail reconnection. Thus tail reconnection has a central role in the magnetospheric energy circulation in Gumics-4. However, the tail reconnection in our simulation is a completely passive phenomenon, in the sense that it is directly driven by what comes through the magnetopause. It does not exhibit any spontaneous behaviour, nor is there any sign of a loading-unloading cycle. The tail seems rather to damp short-time scale fluctuations in the energy flow from the magnetopause, while long-time scale fluctuations are simply transmitted through.
Storage and explosive release of energy in the magnetotail is an essential and observationally supported feature of substorm models (e.g. Baker et al., 1996). On the other hand, the results of Gumics-4 are fairly consistent with observations in many other respects , despite the lack of loading and unloading. This indicates that the magnetospheric dynamics also includes notable directly driven aspects, which is consistent with the analysis of Sun et al. (1998).
We note that tail reconnection in Gumics-4 responds to driving in the same way in synthetic and real event runs. This testifies to the consistency of the simulation and increases confidence in other results obtained from synthetic runs.
To summarize, we have shown that in the global MHD simulation Gumics-4, an increase in energy flow into the magnetosphere leads to the formation of a diffusion region between the tail lobes, with a plasma flow pattern and a magnetic field configuration that are characteristic of magnetic reconnection. In the diffusion region Poynting flux is converted into mechanical energy of the plasma at a rate which is half of the total energy flux entering through the magnetopause.