A&A 444, 961-976 (2005)
DOI: 10.1051/0004-6361:20053600

Current sheet formation in quasi-separatrix layers and hyperbolic flux tubes

G. Aulanier - E. Pariat - P. Démoulin

Observatoire de Paris, LESIA, 92195 Meudon Cedex, France

Received 9 June 2005 / Accepted 26 July 2005

Abstract
In 3D magnetic field configurations, quasi-separatrix layers (QSLs) are defined as volumes in which field lines locally display strong gradients of connectivity. Considering QSLs both as the preferential locations for current sheet development and magnetic reconnection, in general, and as a natural model for solar flares and coronal heating, in particular, has been strongly debated issues over the past decade. In this paper, we perform zero-$\beta$ resistive MHD simulations of the development of electric currents in smooth magnetic configurations which are, strictly speaking, bipolar though they are formed by four flux concentrations, and whose potential fields contain QSLs. The configurations are driven by smooth and large-scale sub-Alfvénic footpoint motions. Extended electric currents form naturally in the configurations, which evolve through a sequence of quasi non-linear force-free equilibria. Narrow current layers also develop. They spontaneously form at small scales all around the QSLs, whatever the footpoint motions are. For long enough motions, the strongest currents develop where the QSLs are the thinnest, namely at the Hyperbolic Flux Tube (HFT), which generalizes the concept of separator. These currents progressively take the shape of an elongated sheet, whose formation is associated with a gradual steepening of the magnetic field gradients over tens of Alfvén times, due to the different motions applied to the field lines which pass on each side of the HFT. Our model then self-consistently accounts for the long-duration energy storage prior to a flare, followed by a switch-on of reconnection when the currents reach the dissipative scale at the HFT. In configurations whose potential fields contain broader QSLs, when the magnetic field gradients reach the dissipative scale, the currents at the HFT reach higher magnitudes. This implies that major solar flares which are not related to an early large-scale ideal instability, must occur in regions whose corresponding potential fields have broader QSLs. Our results lead us to conjecture that physically, current layers must always form on the scale of the QSLs. This implies that electric currents around QSLs may be gradually amplified in time only if the QSLs are broader than the dissipative length-scale. We also discuss the potential role of QSLs in coronal heating in bipolar configurations made of a continuous distribution of flux concentrations.

Key words: magnetohydrodynamics (MHD) - methods: numerical - Sun: magnetic fields - Sun: flares

   
1 Introduction

The energy needed to power solar flares and to sustain coronal heating is thought to come from the coronal magnetic field, since its energy dominates over all other forms of stored energy. However, the coronal plasma, like most natural magnetized plasmas, has typical Lundquist numbers far larger that unity. So the resistive term in the induction equation can become large enough only if small-scale magnetic field gradients (i.e. narrow current layers) are created. Regions in which either the magnetic field, or the velocity field, or the Alfvén speed initially have small scale gradients can naturally result in such current layers. However, those are not necessarily typical of the solar corona, and other situations can also exist. Magnetic configurations with a complex topology, i.e. with separatrices, are the most obvious configurations where current sheets can form, when no steep gradient is initially present in the system. Separatrices are magnetic surfaces where the magnetic field line linkage is discontinuous. A particularly important location for current-sheet formation, then for reconnection in a classical view, is the intersection of two separatrices, which is a null point (a point where the magnetic field vanishes) or a separator. In most cases, a separator is a singular field line joining two null points. More generally, current sheets are thought to form along the separatrices when arbitrary footpoint motions are imposed at a line-tied boundary at the separatrices (e.g. Lau 1993; Low & Wolfson 1988; Aly 1990).

The initial studies of the topology in 3D magnetic configurations have been realized by defining a magnetic field created by discrete sub-photospheric sources (Baum & Bratenahl 1980). Hénoux & Somov (1987) proposed that reconnection along the separator can interrupt currents flowing there, thus permitting to release the energy stored in these currents. Gorbachev & Somov (1988) further applied the theory to an observed solar flare and showed that field lines passing close to the separator connect to observed chromospheric bright ribbons. Furthermore, numerous analyses of flares have shown that H$\alpha$ and UV flare brightenings are typically located along the intersection of separatrices with the chromosphere: they are connected by field lines which are expected to have formed through magnetic reconnection in the given configuration (e.g. Mandrini et al. 1991; Démoulin et al. 1994b; van Driel-Gesztelyi et al. 1994; Mandrini et al. 1995).

The description of the magnetic field with sub-photospheric sources, as well as its related topological analysis irrespective of the location of the line-tied photosphere, is only an approximation to describe the organization of the magnetic field in flux tubes: it implicitly assumes that the origin of a flare is rooted below the line-tied boundary, where the magnetic configuration has no reason to be that of what is prescribed by the sources. Assuming that all the sources are located in the line-tied plane permits this difficulty to be bypassed. However it also leads to some undesirable effects which may not be relevant of the solar photosphere, especially within active regions: wide areas in this plane have purely tangential magnetic fields, except in the vicinity of the sources and of null points located at the boundary. Also, these boundary nulls fully constrain the topology above, which may be considered to be at least restrictive, if not artifical. Finally, the very definition of point charges prevents modeling twisting motions. This approach is nevertheless very interesting, since it allows to use powerful mathematical (analytical) tools, which permit many aspects of the complex problem to be explored without the need of heavy numerical simulations (e.g. Priest et al. 2005; Longcope & Klapper 2002).

The above results have demonstrated that the location of energy release in solar flares is defined by the magnetic topology and that the physical mechanism is magnetic reconnection. However, it has been shown that the energy release did not involve all of the separatrix; e.g. H$\alpha$ flare brightenings were always present only on a restricted part of the chromospheric footprint of the computed separatrices. Moreover, for some observed events, a coronal magnetic null related to the flare was not always present in the configurations associated with the observed photospheric magnetic field (Démoulin et al. 1994a). Another well-known possibility getting separatrices is when field lines are tangent to the photospheric boundary (called "bald patches'', Titov et al. 1993). But just as with coronal nulls, bald patches have been found only in a small fraction of observed events (e.g. Aulanier et al. 1998). In fact, in many flaring configurations, the computed separatrices were only associated with the magnetic nulls being also located below the photosphere (located between the assumed sub-photospheric sources). These studies teach us that coronal magnetic reconnection must occur in a broader variety of magnetic configurations than traditionally thought, as derived from studies of 2D configurations. It is also worth noticing that the separatrices of a magnetic configurations invariant by translation along one direction disappear in most cases when the configuration is fully extended to 3D (Schindler et al. 1988). This structural instability of separatrices point also to the need of a broader concept.

In order to address these difficulties, Démoulin et al. (1996b,a) proposed that "quasi-separatrix layers'' (QSLs) generalize the definition of separatrices to cases where there is no coronal magnetic null. QSLs are regions where there is a drastic change in field line linkage, while the linkage is truly discontinuous at separatrices. For each observed flare studied with this approach, the brightenings were always found along, or just nearby, the intersection of QSLs with the chromosphere (Mandrini et al. 1997; Démoulin et al. 1997; Bagalá et al. 2000, and references therein). These results demonstrate that flares are coronal events, where the release of free magnetic energy is due to the presence of regions where the magnetic field line linkage changes drastically, and not necessarily discontinuously.

Physically, the magnetic energy available for flaring must be associated with non-potential magnetic fields, so in the presence of extended and/or narrow electric current distributions. Indeed, when photospheric vector magnetic field measurements were available in the studies quoted above, two photospheric current concentrations of opposite sign were always found in the close vicinity of the computed QSLs, both linked by modeled coronal field lines (Démoulin et al. 1997, and references therein). Theoretically, the formation of a strong current layer in any QSL is expected with almost any kind of boundary motion which crosses a QSL, as conjectured analytically by Démoulin et al. (1996a). The main reason is that the magnetic stress of very distant regions can be brought close to one another, typically over the QSL thickness. Unfortunately, analytical arguments for current layer formation in QSLs cannot go too far in configurations without symmetry: the derivation of the currents is both a non-linear and a non-local problem requiring in particular integration over field lines. So Titov et al. (2003) considered a straightened magnetic configuration between two plates, with a Hyperbolic Flux tube (or HFT) at the center of QSLs. They calculated analytically that currents form and increase exponentially with time in a HFT, only when the boundary shearing motions create a stagnation point in the middle of the configuration.

MHD numerical simulations are required to analyze the evolution of general magnetic configurations having QSLs. A numerical difficulty is that the currents are expected to form on the scale of the QSLs, a scale that can be many orders of magnitude lower than the scale of the whole studied magnetic configuration. The simulation of Milano et al. (1999) was the first to show that currents did form along QSLs. But the latter were not present in the initial uniform field configuration. They were dynamically formed by the prescribed boundary motions, which consisted of two vortices with a stagnation point in between. The numerical simulation of Galsgaard et al. (2003) was aimed at addressing the problem of pre-existing QSLs. Unfortunately, it considered very broad initial QSLs, with a thickness of about one tenth of the numerical domain, so that only weak currents formed there. The selected boundary motions produced a stagnation point inside the domain, which resulted in the formation of a strong current sheet, as predicted analytically by Titov et al. (2003). However, Démoulin (2005) extensively explained that these strong currents were not associated with the initial broad QSLs, but rather with a new set of much thinner QSLs that dynamically formed thanks both to the stagnation point and to the large-scale boundary displacements. It is then unclear how the results of Milano et al. (1999) and of Galsgaard et al. (2003) can be generalized to magnetic configurations which initially possess narrow QSLs.

In order to really address this issue, the one conjectured in previous studies of observed solar flares, in this paper we instead perform MHD simulations of the slow current build-up associated with pre-existing narrow QSLs. In Sect. 2, we consider and analyze the properties of two potential magnetic configurations that have large differences in QSL thickness. In Sect. 3, we describe the numerical method which we used to evolve both configurations with two different forms of line-tied motions. In Sect. 4, we discuss our results in terms of the formation conditions of narrow current layers within QSLs, in general, and at the HFT, in particular, as a function of the boundary flow. The results are summarized and discussed in the frame of solar flare and coronal heating modeling in Sect. 5. Numerical and physical issues on the width of current sheets are discussed in Appendix A.

2 Definition and topology of the magnetic configurations 

2.1 Initial magnetic configurations

 

We considered a 3D Cartesian domain $x \in [-0.65 ,~ 0.65]$, $y \in [-0.41 ,~ 0.41]$, $z \in [0 ,~ 0.65]$, where z is altitude. z=0 is considered as the photospheric plane, which is treated as a line-tied boundary in which kinematic motions were prescribed in the MHD simulations. In this domain, we calculated two potential magnetic field configurations ( $\triangle \vec{b} =0$) made up of four polarities: P1 (resp. N1) is the positive (resp. negative) polarity of an outer bipole, and P2 (resp. N2) is the positive (resp. negative) polarity of an inner bipole which contains less magnetic flux than the outer one, but which has stronger magnetic field concentrations.

Each of the four polarities results from point-sources located at various depths beneath the photosphere. Throughout this paper, both configurations are labeled by $\Phi =120^\circ $ and $150^\circ $, which are the angle between the axes of both bipoles when the field is potential at t=0. The magnetic field $\vec{b}$of these configurations is given by:

                               bx(x,y,z) = $\displaystyle \Sigma_{i=1}^{4}~ F_{i} ~ (x-x_{i})~ r_{i}^{-3} ,$  
by(x,y,z) = $\displaystyle \Sigma_{i=1}^{4}~ F_{i} ~ (y-y_{i})~ r_{i}^{-3} ,$  
bz(x,y,z) = $\displaystyle \Sigma_{i=1}^{4}~ F_{i} ~ (z-z_{i})~ r_{i}^{-3} ,$  
ri = $\displaystyle \sqrt{ (x-x_{i})^2 + (y-y_{i})^2 + (z-z_{i})^2 } .$ (1)

The values of the free parameters for each polarity ( xi,yi,zi,Fi) are given in Table 1. Typical field lines for both configurations are shown in the upper row of Figs. 1 and 2. Note that the point sources are used only to define the initial magnetic field. Later, only the magnetic field for $z\geq 0$ is considered (Sect. 3).

Table 1: Parameters of the magnetic configurations.

The chosen settings imply that there is neither a magnetic null point in $ z\ge 0$ nor field lines tangential to the photospheric boundary at z=0, so there are no separatrices in the domain. Topologically speaking, the configurations are bipolar, equivalent to an arcade, even though they display four contrasted magnetic field concentrations. The following numbers permit to estimate the degree of contrast at z=0 for $\Phi =120^\circ $ and $150^\circ $: $b_z^{ {{\rm {max}}}}({{\rm in P2}}) \simeq 35$, $b_z^{ {{\rm {max}}}}({{\rm in P1}}) \simeq 25$, and $b^{ {{\rm {min}}}} \simeq 3$.


   \begin{figure} \par\mbox{\psfig{figure=3600f1a.ps,width=5.75cm,angle=-90} \psfig... ...,angle=-90} \psfig{figure=3600f1f.ps,width=5.75cm,angle=-90} } \par \end{figure} Figure 1: Top views of the magnetic configurations for two orientations of the central bipole with respect to the large one ( $\Phi =120^\circ $ and $150^\circ $): potential field ( top row) and configurations resulting from the translational ( middle row) and twisting ( bottom row) motions applied at z=0 (see Fig. 5). The horizontal (resp. vertical) axis corresponds to the x (resp. y) coordinate. Pink (resp. blue) contours show positive (resp. negative) values of $b_z(z=0)=\pm 5,10,15,20,25,30$. The inversion line bz(z=0)=0 is plotted in yellow. The other lines are magnetic field lines. In each panel, they are plotted starting from the same positions at z=0 in the negative polarities, so these lines are comparable from one panel to another since their fixed footpoint is not displaced by the flows (see Fig. 5). The time unit corresponds to the transit time of an Alfvén wave on a distance of 0.2 at the initial homogeneous Alfvén speed.
Open with DEXTER


  \begin{figure} \par\mbox{\psfig{figure=3600f2a.ps,width=5.75cm,angle=-90} \psfig... ...75cm,angle=-90} \psfig{figure=3600f2f.ps,width=5.75cm,angle=-90} }\end{figure} Figure 2: Same as Fig. 1, but projection views in the full numerical domain are shown. The viewing angle is rotated approximatively $180^\circ $ around the z axis relative to Fig. 1.
Open with DEXTER

2.2 Definition of quasi-separatrix layers 

QSLs are defined as regions where there is a drastic change in field-line linkage (Démoulin et al. 1996b,a). More precisely, let us consider the mapping from one photospheric polarity to the opposite one: $\vec{r}_{+}(x_{+},y_{+}) \mapsto \vec{r}_{-}(x_{-},y_{-})$ and the reversed one $\vec{r}_{-}(x_{-},y_{-}) \mapsto \vec{r}_{+}(x_{+},y_{+})$. These mappings can be represented by some vector functions [X-(x+,y+), Y-(x+,y+)] and  [X+(x-,y-), Y+(x-,y-)], respectively. The norms  $N(\vec{r}_{+})$ and  $N(\vec{r}_{-})$ of the respective Jacobian matrices in Cartesian coordinates are:

 
                               $\displaystyle N_{\pm}$ $\textstyle \equiv$ $\displaystyle N(x_{\pm},y_{\pm})$ (2)
  = $\displaystyle \sqrt{ % \left( {\partial X_{\mp} \over \partial x_{\pm}} \right)... ... \right)^2 + \left( {\partial Y_{\mp} \over \partial y_{\pm}} \right)^2 } \cdot$  

A QSL was first defined by the condition $N(x_{\pm},y_{\pm}) \gg 1 $in both photospheric polarities (Démoulin et al. 1996b).

Let us now consider a field line linking photospheric locations  (x+,y+) and  (x-,y-), which both have different normal field components bz+ and bz-. In this case, a difficulty with the definition of QSLs by Eq. (2) is that  $N(x_{+},y_{+}) \neq N(x_{-},y_{-})$ if  $b_{z+} \neq b_{z-}$, so QSLs do not fulfill a unique condition, in general, when defined by Eq. (2). Recently, Titov et al. (2002) defined another characteristic function for QSLs which is independent of the mapping direction: the squashing degree Q. It is calculated as follows:

 \begin{displaymath}Q_{+} = \frac{N_{+}^2}{\vert b_{z+} /b_{z-}^{*}\vert} ~ ~ \... ...ac{N_{-}^{*2}}{\vert b_{z-}^{*} /b_{z+}\vert} ~ ~ \equiv Q~, \end{displaymath} (3)

where asterisking the functions indicates that their arguments x-and y- are substituted in  X-(x+,y+) and  Y-(x+,y+), respectively. With this new prescription, a QSL is defined by $Q\gg2$, the value Q=2 being the lowest value possible (Titov et al. 2002). By definition, Q is uniquely defined along a field line by:

 \begin{displaymath} (\vec{b} \cdot \vec{\nabla}) Q = 0 . \end{displaymath} (4)

The physical meaning of this apparently complex definition can, in fact, be simply explained as follows. Let us consider an elementary flux tube rooted in an infinitesimal circular region in one polarity. Q simply measures the aspect ratio of the distorted ellipse defined by the footpoint mapping of this flux tube in the other polarity. In other words, Q measures how much the initial elementary region is squashed by the field-line mapping.

Other quantities which define the mapping properties of QSLs do exist, one of them being the ratio |bz+ /bz-*|. Their full description and meaning are analyzed in Titov et al. (2002). In this paper, we only use Q, since it is sufficient to localize the QSLs and to define their thickness. A global view of QSL properties and their application to coronal physics is presented in Démoulin (2005) and Titov (2005).

   
2.3 Numerical calculation of quasi-separatrix layers

3D magnetic configurations, where the maximal value of Q is large, are challenging for numerical computation of the associated QSLs, because their related widths are often orders of magnitudes below the spatial resolution of any numerical mesh. For local quantities such as the magnetic field, computing gradients below the mesh size has no meaning. But this is not true for N and Q, because their values are dominantly determined by the large-scale properties of the magnetic configuration. This was thoroughly explained in Démoulin et al. (1997,1996a), where the effect of spatial discretization defining the analytical magnetic fields was tested. There it was demonstrated that calculating QSLs below the scale of the discretization is relevant and reliable as long as the large scale-lengths of the magnetic field are well resolved.

In order to accurately compute N (or Q) in a plane (e.g. z=0), one needs to determine a 2D grid (e.g. in x,y) which is locally adapted to the QSL width. The latter must then be different than the grid used later in the MHD simulations. One must also be able to compute the connectivities with high precision. In the case of a uniform grid, with a mesh interval $\delta$ small enough to resolve all the connectivity gradients, the computation of N is only limited both by the numerical precision of the field line integration and by the numerical derivatives in Eq. (2). With the best integration algorithms, the position  $\vec{r}_{i,j}$ of the second footpoint of the field line passing by the grid point of index (i,j) can be known with relative precision of 10-11. This permits us to calculate the field line connectivities very accurately. Since local small-scales of the magnetic field have a low effect on the QSLs, the magnetic field in between the mesh points is computed simply by a linear interpolation of the nearest points values. Then from Eq. (2), the computed value Ni,j of N at the grid point (i,j) is given by the knowledge of the field line connectivity of the four surrounding points on the grid:

 \begin{displaymath} N_{i,j} = \frac{\sqrt{(\vec{r}_{i+1,j}-\vec{r}_{i-1,j})^2 + (\vec{r}_{i,j+1}-\vec{r}_{i,j-1})^2}}{2 \delta} \cdot \end{displaymath} (5)

Qi,j is then derived from Ni,j using Eq. (3).


  \begin{figure} \par\mbox{\psfig{figure=3600f3a.eps,width=6.0cm,angle=-90} \psfig... ...space*{0.13cm} \psfig{figure=3600f3d.eps,width=5.78cm,angle=-90} }\end{figure} Figure 3: Top views of the quasi-separatrix layers (QSLs) for the potential field configurations for  $\Phi =120^\circ $ and $150^\circ $. The greyscale images show the distribution of the squashing degree Q (Eq. (3)) at z=0 in logarithmic scale. For $\Phi =120^\circ $ (resp.  $150^\circ $), white corresponds to Q=5 (resp. 6) and black to $Q=1.3 \times 10^{5}$ (resp. 1011). Typical field lines that are rooted in the vicinity of the QSLs, and the same bz(z=0) contours as in Fig. 1, are overplotted in the bottom row.
Open with DEXTER

In practice and in order to save computational time, we never use a uniform grid. Q is first computed on a coarse grid, whose resolution progressively improves adaptively in a multi-step procedure. In a first step, Q is calculated everywhere, but only those regions where Q is the highest are kept. In these regions, the spatial resolution is doubled. In a second step, Q is re-computed in the previously selected points, as well as in the new four neighboring point of the improved grid. This procedure is repeated until the number of points where Q is computed reaches some previously fixed value. In the present paper, this number of point was 3500, which corresponds to a spatial resolution of $\sim$ $3 \times 10^{-3}$, twice larger than the smallest cell which is considered in the MHD simulations (see Sect. 3.1).

With such grid resolution, however, Q is still approximatively computed where the gradients of connectivity are strong, i.e. where the QSL are the thinnest. For the refined calculation of Q in 2D, we consider a local square around each saved position of the latter grid. In a second multi-step procedure, Q is then successively re-computed in the central point of each square, as the resolution progressively increases by a factor two at each step. The recurrence is stopped when the values of Q computed at two consecutive steps converge, i.e. when their ratio exceeds a fixed value (which we choose to be equal to 0.9 in this paper). At that stage, Q is calculated well at the 3500 selected points. Note than in separatrices, this second iteration never converges since Q tends to infinity.

The results of these two iterative procedures were used to generate all Q maps at z=0 shown in this paper and to calculate the corresponding $Q^{ {{\rm {max}}}}$. Figure 3 shows such maps for the potential magnetic field configurations defined in Sect. 2.1. The shape of the QSLs and their significance are discussed below in Sect. 2.4.

We define the QSL width as the full width at half maximum of the Q profile. In order to calculate the latter, we recompute Q along several 1D segments that cross the QSL at various angles, in the plane z=0. The 2D Q maps are used to choose the position of these segments. Q can there be computed using various spatial resolutions  $\delta ^\prime $. The true QSL width along a given direction is reached when $\delta ^\prime $ is small enough to ensure that any further refinement does not change the Q profile. The minimum full width at half maximum of the Q profiles along each of the segments finally results in the QSL width.

Table 2: Maximum amplitude of the footpoint motions, of the squashing degree and minimum width of its profile across the QSLs, in each configuration studied. The amplitude of the translation motions, $2\; \delta y$, is measured by the ratio of the maximum displacement ($\delta y$), divided by characteristic size along y of the system (0.5). The amplitude of the twisting motions is measured by the maximum numbers of turns $\Delta N$ within P2. The importance of the mapping distortion is given by the maximum value of squashing degree (Eq. (3)), $Q^{ {{\rm {max}}}}$. Finally the thicknesses of the QSLs are given by the ratio of the full width at half maximum of Q, $\delta _{Q}$ with the smallest cell size, d (Sect. 3.1). Both $Q^{ {{\rm {max}}}}$ and $\delta_{Q}~ /d$ are computed with a low spatial resolution of $d = 1.5\times 10^{-3}$ (resp. at a much higher spatial resolution, see Sect. 2.3) and the values are noted with a subscript "L'' (resp. "H'').

Figure 4 shows QSL profiles for both configurations defined in Sect. 2.1, using two different segment lengths, thus with two different spatial resolutions  $\delta ^\prime $. With the lower resolution which is of the order of the MHD mesh resolution (see Sect. 3.1), one can neither obtain the correct value of $Q^{ {{\rm {max}}}}$, nor the real width of the central peak of the QSL profile. The broad wings of the QSL are, however, well visible. The properties of the QSL, as calculated with both resolutions, are given in Table 2. It shows that the issue of resolution is the most sensitive for $\Phi =150^\circ $. The QSL for $\Phi =120^\circ $ at t=0, however, is almost resolved by the numerical mesh used in the MHD simulations: "almost'', because the mesh is non-uniform and d is only the smallest grid size (see Sect. 3.1).

   
2.4 Topology of the potential fields


  \begin{figure} \mbox{\psfig{figure=3600f4a.ps,width=4.2cm,angle=-90}\hspace*{0.1... ...\hspace*{0.10cm} \psfig{figure=3600f4d.ps,width=4.2cm,angle=-90} }\end{figure} Figure 4: Profiles of $\log Q(z=0)$ perpendicular to the QSL for the initial configurations $\Phi =120^\circ $ ( left column) and $\Phi =150^\circ $ ( right column). The cuts are centered at (x,y)=(xQ,yQ), with (xQ,yQ)=(-0.155,-0.145) for $\Phi =120^\circ $ and (xQ,yQ)=(-0.121,-0.203) for $\Phi =150^\circ $. The profiles are calculated using a spatial resolution corresponding to the smallest cell in the numerical mesh $\delta ^\prime = d = 1.5\times 10^{-3}$ ( upper row) and using an optimized higher spatial resolution $\delta ^\prime $ ( lower row), as explained in Sect. 2.3.
Open with DEXTER

In the bipolar potential configurations defined in Sect. 2.1, the magnetic field line linkage has four basic sets of magnetic connectivities (see Figs. 1 and 2), just as in 2D quadrupolar magnetic configurations, but without separatrix between them. For both configurations, the intersections of the QSLs with the z=0 boundary have only two extended thin strips, one over each magnetic polarity (see Fig. 3, top row). These potential configurations are thus very similar to the one analyzed by Démoulin et al. (1996a) and Titov et al. (2002).

Two close field lines rooted at z=0 on both sides of one strip rapidly diverge in the volume to connect, on the other strip, regions which are very far from each other, as shown in Fig. 3, bottom row.

The thin volume, where Q has the highest values, is of particular interest: the way field lines diverge there suggests to call the magnetic structure of QSLs a HFT Titov et al. (2002). The 3D shape of this HFT is better understood as one follows its 2D cross-section from one polarity to the other one on the boundary. Let us define the edge of the QSL by the value $Q_{\rm e}$, which is a fraction of the maximal value of Q. The HFT starts as an elongated strip over one polarity, then it is transformed progressively in a cross shape in the volume, and it ends in the form of another elongated strip on the other polarity. Each strip at z=0 involve one branch of the cross at z>0. A cartoon of the cross-section from one polarity to the other is then:

 \begin{displaymath}\setlength{\unitlength}{0.75mm} \begin{picture} (110,8)(0,0... ...} % \multiput(100,0)(-0.3,0){6}{\line(-1,1){10}} \end{picture}\end{displaymath} (6)

This shape is similar for any values $Q_{\rm e}\gg 2$. The volume defined by $Q^\prime_{\rm e}$ fully encloses the volume defined by $Q_{\rm e}$ if  $Q^\prime_{\rm e}>Q_{\rm e}$, so that defining an increasing series of $Q_{\rm e}$ values defines a series of volumes that are somehow organized like Russian dolls.

The QSL shape is robust to the transformation of the magnetic configuration, as the locations of the highest values of Q for the two configurations of Fig. 3 have similarly curved shapes. The slight modifications of these shapes mostly follow the displacement of the polarities.

The maximum value of Q however, is extremely sensitive to modifications of the magnetic configurations: when Q is calculated at a spatial resolution much higher than the numerical discretization chosen for the MHD simulation so as to obtain its true profile (see Sect. 2.3), $Q^{ { {\rm {max}}}} \sim 4 \times 10^3$ for $\Phi =120^\circ $, while $Q^{ {{\rm {max}}}} \sim 6 \times 10^8$ for  $\Phi =150^\circ $, five orders of magnitude higher. Asymptotically, $Q^{ {{\rm {max}}}}$ tends to infinity as $\Phi$ tends to $180^\circ $.

   
3 Method for MHD evolution

   
3.1 Equations and mesh

We use a simplified version of our zero-$\beta$ (pressureless) time-dependent 3D MHD code, which is extensively described in Aulanier et al. (2005). The present version solves the following equations:

  
                          $\displaystyle \rho~ \frac{\partial \vec{u}}{\partial t}$ = $\displaystyle - \rho~ (\vec{u} \cdot \vec{\nabla}) \vec{u} + \vec{\jmath} \times \vec{b} + \rho~ \mathfrak{D} \vec{u}$ (7)
$\displaystyle \frac{\partial \vec{b}}{\partial t}$ = $\displaystyle \vec{\nabla} \times (\vec{u} \times \vec{b}) + \eta \triangle \vec{b}$ (8)
$\displaystyle \vec{\nabla} \times \vec{b}$ = $\displaystyle \vec{\jmath}$ (9)
$\displaystyle \vec{\nabla} \cdot \vec{b}$ = 0 , (10)

where $\rho$ is the mass density, $\vec{u}$ the plasma velocity, $\vec{b}$ the magnetic field, $\vec{\jmath}$ the electric current density and $\eta$ the magnetic resistivity. The calculations are achieved in a dimensionless form, so that the magnetic permeability is set to 1. $\mathfrak{D}$ is a diffusion operator for the velocity (see Sect. 3.3).

Since we are only interested in quasi-static evolutions and to save computer time, we fix $\rho$ in time to its initial value given below:

 \begin{displaymath}\rho(x,y,z) = c_\circ^{-1} b^2(x,y,z,t=0) , \end{displaymath} (11)

so that $c_A(x,y,z,t=0)=c_\circ=0.2$. This leads to define the time unit as the transit time of Alfvén waves over a distance of 0.2, which corresponds to the physical spacing between both central polarities P2 and N2. This setting does not lead to any singularity, since the studied configurations contain no magnetic null point in the domain.

The boundary conditions at z=0 are line-tied, and those of the five other faces are open. Their numerical implementation with the use of ghost cells is described in details in Aulanier et al. (2005).

The simulations are done in the domain defined in Sect. 2.1, using a non-uniform mesh $n_x \times n_y \times n_z = 191 \times 161 \times 170$ points. The mesh intervals vary in the range d $x \in [1.5\times 10^{-3}~ ,~ 1.8\times 10^{-2}]$, d $y \in [1.5\times 10^{-3}~ ,~ 1.2\times 10^{-2}]$, d $z \in [1.5\times 10^{-3}~ ,~ 0.8\times 10^{-2}]$, expanding from x=y=z=0 following dxi+1/dxi=dyj+1/dyj=1.027 and dzk+1/dzk=1.01.

   
3.2 Boundary motions


  \begin{figure} \par\mbox{\psfig{figure=3600f5a.eps,width=5.58cm,height=4.05cm,an... ...fig{figure=3600f5d.eps,width=5.58cm,height=4.05cm,angle=-90} } \par \end{figure} Figure 5: Boundary flow patterns applied for $\Phi =120^\circ $ and $150^\circ $, shown by dark arrows superposed to the same bz(z=0) isocontours as shown in the top row of Fig. 1. The flows are fixed in time for t>13 (see Eq. (15)).
Open with DEXTER

The magnetic configurations evolve in response to large-scale kinematic motions  $\vec{u}_{x,y}$ which we prescribe in the line-tied plane at z=0. Since we want to study the dependence of the generation of electric currents at QSLs with respect to the precise nature of the footpoint motions, we consider various types of motions which move only a part of the QSLs. Firstly, we only apply motions within the positive central polarity P2. Secondly, we choose two types of motions which contain neither X-type stagnation point, nor small scales. The first type of motion is a nearly solid translation of P2 along y. The second type of motion is twisting of the strongest fields in P2, which has a nearly solid rotation over more than half of the vortex radius. Both types of boundary motions are shown in Fig. 5, superposed on contours of bz(z=0). They both do not directly affect the field lines which have a footpoint in P1.

The translation motion is defined by:

 
                         ux(z=0) = 0 ,  
uy(z=0) = $\displaystyle \frac{u^\circ}{4} \bigg{[}\tanh\bigg{(} \frac{y-y^\circ(x)}{\delta y^\circ}\bigg{)}+1\bigg{]}$  
    $\displaystyle \bigg{[}\tanh\bigg{(}\frac{x-x_1^\circ}{\delta x^\circ}\bigg{)} -\tanh\bigg{(}\frac{x-x_2^\circ}{\delta x^\circ}\bigg{)} \bigg{]} ,$  
$\displaystyle y^\circ(x)$ = $\displaystyle 4 (x-x_3^\circ)^2+y_1^\circ ,$ (12)

where $u^\circ$ is the maximum velocity at the boundary. In all our simulations, we set $u^\circ = 1.5\times 10^{-3} = 0.75 \%$ of $c_\circ$ (defined with Eq. (11)) so as to ensure a relatively slow driving of the system.

The twisting motion in P2 is defined by:

  
    $\displaystyle {u}_{x}(z=0) = \frac{\partial \psi}{\partial_y}~ ~ ,~ ~ {u}_{y}(z=0) = -\frac{\partial \psi}{\partial_x} ~ ,$ (13)
    $\displaystyle \psi = \psi^\circ \tanh\big{[}\alpha_1^\circ~ b_z^2(z=0)\big{]} \tanh^4\big{[}\alpha_2^\circ~ b_z^2(z=0)\big{]},$ (14)

where $\psi^\circ$ is a parameter which is adjusted so as to prescribe a maximum twisting velocity of $u^\circ$ (same value as above). With the velocity written as in Eq. (13), bz(z=0) is only advected with time (without modifying its Lagrangian value), and this is re-enforced numerically at each numerical iteration. The choice of such complex $\psi$ functions (Eq. (14)) was motivated to prescribe a nearly uniform twisting motion in the strong field regions, surrounded by a region of fast velocity decrease and a last outer region in which both the velocities and their horizontal derivatives tend to zero close to the inversion line around P2. In this way, the central part of P2 does not incorporate small scales, and its outer regions do not lead to numerical instabilities, since no small unresolved field line is advected and since the velocity is numerically derivable everywhere.

The values for the remaining free parameters in Eqs. (12) and (14) are given in Table 3. In the simulations, all the velocities given above are multiplied by $\gamma(t)$:

 \begin{displaymath} \gamma(t) = \frac{1}{2}~ \tanh\bigg{[}\frac{2(t-10)}{3}\bigg{]} + \frac{1}{2} , \end{displaymath} (15)

which allows the system to first relax to a numerical equilibrium for 0<t<7, followed by an early acceleration phase for 7<t<13, towards a constant boundary driving for t>13.

Table 3: Parameters for line-tied boundary motions.


  \begin{figure} \par\mbox{\psfig{figure=3600f6a.ps,width=8.5cm,angle=-90}\hspace*... ...}\hspace*{0.0cm} \psfig{figure=3600f6d.ps,width=8.5cm,angle=-90} }\end{figure} Figure 6: Greyscale images of the coronal currents j(x,z) at y=0.07 in linear scale. In all panels, dark grey corresponds to j(x,z)=0. White corresponds to j(x,z)=100,300,100,150 respectively, for the ( upper-left), ( lower-left), ( upper-right), and ( lower-right) panels. Each image shows the co-existence of "extended'' currents which result from the line-tied footpoint motions, of "narrow'' currents layers within the whole QSLs, and of an intense current layer at a Hyperbolic Flux Tube ("HFT'') located where the narrow current layers intersect. The plots are drawn a few Alfvén times before the magnetic field gradients reach the scale of the mesh in the HFT.
Open with DEXTER

   
3.3 Diffusion operators

Some strong Lorentz forces develop during the calculations on small-scales. They lead to strong vorticity layers on the scale of a few cells, which are typically located around the QSLs at z>0. Their proper numerical treatment leads us to use the following diffusion operator for velocity:

 \begin{displaymath} \mathfrak{D}~ u_i = \frac{u_{\nu^\star}}{d} \left( \delta_x^2 u_i + \delta_y^2 u_i + \delta_z^2 u_i \right) \\ \end{displaymath} (16)

where ui is the velocity component along either axis (x,y,z) and d is the smallest cell size in the domain. $u_{\nu^\star}$ is the characteristic diffusion speed. We set to $u_{\nu^\star} = 0.03 = 15\%~c_\circ$ in all our simulations. This leads to strong viscous effects, which are unfortunately required: setting $u_{\nu^\star}=u^\circ$, which is the standard value in turbulent simulations, leads to numerical instabilities in the HFT typically at $t \sim 20$, which do not permit us to follow the development of strong electric currents over long time-scales. $\delta_x^2$ is a second-derivative operator with respect to the mesh rather than to spatial units. For any quantity f, this operator is equal to:

\begin{displaymath}\delta_x^2 f = f(x^{i+1},y^{j},z^{k}) -2 f(x^{i},y^{j},z^{k}) + f(x^{i-1},y^{j},z^{k}). \end{displaymath} (17)

So as to reach a compromise which ensures that the effects of resistivity are small, but enough to ensure numerical stability for a long time, we set $\eta=1.5\times 10^{-6}$ in all our simulations. This leads to a low characteristic resistive speed of $\eta/d = 10^{-3} = 0.5\%~c_\circ$.

Considering $u^\circ$ as the characteristic velocity of the system implies the following magnetic Reynolds and Lundquist numbers: $R_{\rm m}=1.5$ and Lu=200 at the scale of the smallest cell, $R_{\rm m}=200$and $Lu=2.6\times 10^{4}$ at the scale of the central bipole P2N2, $R_{\rm m}=10^{3}$ and $Lu=1.3\times 10^{5}$ at the scale of the full magnetic configuration.

   
4 Development of electric currents

In spite of the very small scales in the QSLs, which are intrinsic to the studied configurations, our MHD simulations do not result in numerical instabilities for several tens of Alfvén times. During this long time interval, electric currents develop in various regions (Sect. 4.1). In particular narrow current layers develop all along the QSLs (Sect. 4.2), while the strongest current layer is formed at the HFT where the QSLs are the narrowest (Sect. 4.3).

   
4.1 Extended and narrow current layers

In all our simulations, the footpoint motions naturally lead to the formation of nearly-field-aligned currents distributed over wide volumes which are defined by the envelope of the field lines which are transported. We call them "extended current layers''. These currents are stronger for the twisting motions than for the translation motions (see Fig. 6). Even though these currents are relatively strong, they do not tend to dissipate easily, since they result from large-scale magnetic field gradients in the domain induced by the line-tied motions. Indeed, the resistive dissipation term is $\propto$ $\triangle \vec{b} \sim \jmath / L$, which clearly shows that for equal electric current densities, the narrowest current layers will dissipate more quickly. Since the boundary motions are less than 1% of the Alfvén speed, the electric current remains nearly aligned with the magnetic field in these extended regions, so the whole configuration is always very close to a force-free state. This behavior is typical of every MHD simulation with slow line-tied boundary motions (e.g. Török & Kliem 2003; DeVore & Antiochos 2000; Aulanier et al. 2005).

It must be noted that in the time interval during which our motions are prescribed, the footpoint displacements remain relatively small, at most of the order of $\sim$1/4th (resp. 1/7th) of the characteristic size of the system for $\Phi =120^\circ $ (resp. $150^\circ $). This is shown in Figs. 1 and 2. Also, as explained in Sect. 3.2, none of the prescribed line-tine motions possess very small scales. In spite of all this, the footpoint displacements in our simulations lead to the development of "narrow current layers'' at $ z\ge 0$. They begin to form on small scales as soon as the motions start, so they mostly do not come from some time-varying steepening effect. Another property is that, for a given magnetic field configuration, these narrow currents layers form in the same specific locations (see Fig. 6), whatever the prescribed motions, translation or twisting.

All the above properties lead to the conclusion that these narrow current layers are not a direct consequence of the prescribed velocity gradients at z=0, as is usually the case in line-tied MHD simulations in which large-scale and long-duration braiding or twisting or shearing motions are applied (see e.g. DeVore & Antiochos 2000; Galsgaard & Nordlund 1996; Aulanier et al. 2005; Mikic et al. 1989; van Ballegooijen 1986; Galsgaard et al. 2003). 2D slabs (in x,z) of the 3D currents layers are shown in Fig. 6, a few Alfvén times before magnetic field gradients reach the scale of the mesh and halt the simulations. These currents display a shape which is reminiscent of separatrices with a null point or with a separator, though none of the latter exist in the 3D magnetic configurations that are analyzed.

In the translation cases, the electric current densities in the narrow layers are almost everywhere larger than the extended currents. They are also associated with Lorentz forces, so that they are not fully force-free. In the twisting cases, the extended currents are the highest at low altitude above the polarity P2 for $z \le 0.05$, but they have similar magnitudes than those in the narrow layers almost everywhere else. For both types of motions, the magnitude of the currents in the narrow layers and in the extended regions increase in time at similar rates, except in the region where two narrow layers intersect (Fig. 6). In all runs, the smallest-scale currents eventually form in this latter region. Their time-evolution is described in Sect. 4.3.

   
4.2 Current layers at QSLs


  \begin{figure} \par\mbox{\psfig{figure=3600f7a.ps,width=6.00cm,angle=-90} \psfig... ...} \psfig{figure=3600f7h.ps,width=6.00cm,height=4.29cm,angle=-90} }\end{figure} Figure 7: ( Left column): greyscale images of the electric currents at the lower boundary j(z=0) in logarithmic scale. In all panels, white (resp. dark grey) corresponds to j(z=0)=850 (resp. 10-3). ( Right column): greyscale images of the squashing degree at the lower boundary Q(z=0) in logarithmic scale. As in Fig. 3, for $\Phi =120^\circ $ (resp. $150^\circ $), white corresponds to Q=5 (resp. 6) and black to $Q=1.3 \times 10^{5}$ (resp. 1011).
Open with DEXTER


  \begin{figure} \par\mbox{\psfig{figure=3600f8a.eps,width=5.2cm}\hspace*{0.2cm} \... ...5.2cm}\hspace*{0.2cm} \psfig{figure=3600f8c.eps,width=5.2cm} }\par\end{figure} Figure 8: Color images of 2D slabs (in x,z) at y=0.07 of the HFT for the configuration $\Phi =150^\circ $ at t=38 evolved with twisting motions (see Fig. 6, lower-right panel). ( Left panel): logarithm of the squashing degree $=\log Q$. ( Middle panel): inverse of the scale-length of the magnetic field gradients in current layers $=\jmath /b$. ( Right panel): magnitude of the electric currents $=\jmath $.
Open with DEXTER

In order to investigate the relation between the current layers and the QSLs, we calculate the distribution of the squashing degree Q(z=0) with exactly the same procedure as described in Sect. 2.3 for the potential fields. For all configurations, the maximum values of the squashing degree  $Q^{ {{\rm {max}}}}$ and the associated widths  $\delta _{Q}$ of the QSLs are given in Table 2, as calculated with spatial resolutions that are typical of the numerical mesh and with resolutions that are much finer than the mesh. With both resolutions, we note that $Q^{ {{\rm {max}}}}$ is larger than 2 by several orders of magnitudes. The precise characteristics of the QSLs, however, strongly depend on the resolution at which they are computed. For a resolution equal to the smallest grid size d of the MHD simulations, $\delta_{Q} \simeq d-3d$ for both configurations. But for the higher resolution, the minimum value of  $\delta _{Q}$ is always much smaller than d (except for the potential field of the configuration $\Phi =120^\circ $ where $\delta _{Q} \sim d$). Thus the numerical grid of the MHD simulations does not resolve most of the QSLs which result from the footpoint motions.

The 2D maps of Q and j at z=0 are drawn in Fig. 7 at the same times as in Fig. 6. Apart from the regions that were directly affected by the boundary motions (i.e. within and around P2), Q and j show a striking resemblance in all four cases. The similarity is most evident for the translation motions, since the latter produce less extended field-aligned currents in the displaced flux tubes. But the similarity is also very visible when twisting motions are applied. Also, apart from the region covered by the flows, the QSLs are weakly deformed by the line-tied motions (compare Figs. 3 and 7). The currents then spontaneously form where the QSLs are located in the potential fields. Since the selected motions have no relationship with the QSLs, we then reach the interesting conclusion that any boundary line-tied motion invariably results in the formation of current layers all along narrow QSLs. In our simulations, most of the spatial locations of these current layers are defined by the intrinsic properties of the magnetic configurations that already exist for the corresponding potential field. They are not defined by the topological properties of the boundary motions. Then these current layers are formed just like current sheets in configurations which have separatrices. We then reach an opposite conclusion from Titov et al. (2003) and Galsgaard et al. (2003), who pretend that the nature of the boundary motions is a determining factor in the formation of current layers in QSLs.

Figure 7, especially for $\Phi =120^\circ $ for which larger twists could be applied, clearly shows how the rotational motions deform the QSL in the middle of P2, and how they tend to develop new wide and well resolved QSLs (with weaker Q) around the envelope of the twisted area. These new QSLs result directly from the twisting profile, which rapidly decreases to zero away from the center of P2. It is worth noticing that these new QSLs are also matched by electric currents, but they are neither as intense nor as narrow as the current layers which form in the main QSLs (see Fig. 6).

It is finally interesting to note that in our four simulations, the widths of the narrow current layers that form around QSLs tend to be larger for initially broader QSLs, as seen in Figs. 6 and 7. Also, the width of the current layers is well resolved, of the order of  $\delta _{Q}$ as calculated with the resolution of the numerical simulations. This issue and its consequences are discussed further both in the context of QSLs and separatrices in Appendix A.

   
4.3 The hyperbolic flux tube

Apart from the regions right above sheared/twisted polarities, the strongest electric currents ,which eventually form in all our simulations, are always located at high z, even though the magnetic fields are the strongest at low z (see Fig. 6). The location of these currents corresponds to the region in the QSL that has the highest squashing degree Q. It is the core of the QSL. In the limit of infinitely thin QSLs, this region corresponds to the intersection of two separatrices, which is called a separator. Contrary to a separator, which is a singular line, the HFT is a complex layer-like volume that takes the very elongated shape of the QSLs at the boundary, as shown in Eq. (6)).

For the specific configuration $\Phi =150^\circ $ at t=38 evolved with twisting motions, Fig. 8 shows the comparison of the squashing number in the vicinity of the HFT, calculated from global field line integrations with both the magnitude and the width of the current layers calculated from the local magnetic field derivatives. The Q map was calculated as explained in Sect. 2.3, except that the grid was defined on the plane y=0.07 (instead of z=0) and that the field lines were integrated in both directions from this plane. A similar behavior was found in all our runs. It is obvious that even though the strongest currents are the distributed ones at low altitude, the current layers are narrower within QSLs. These narrow currents reach their minimum thickness within the core of the QSLs, i.e. in the HFT. As mentioned in Sect. 4.2, it is also clear from Fig. 8 that the currents which form in the MHD  simulations are broader than the unresolved central peaks of the QSLs.

At early times, the current layer, which develops in the vicinity of the HFT at high z, first has a nearly circular shape in the (x,z) plane around y=0, with four extensions along the QSLs. Its diameter is $\sim$ $3 \times 10^{-2} \sim 20 d$. In the case of twisting motions, it is a combination of the outer parts of the extended currents and of the currents which form right in the middle of the HFT. This combination explains the spatial shift between the center of the current sheet and of the HFT visible in Fig. 8. In the case of translation motions, however, the maximum currents are almost co-spatial with the center of the HFT. Then in all our runs, as time progresses, this current layer flattens vertically along z and slowly expands horizontally, mostly along x. We thus find that, whatever the precise form of the boundary motions, HFTs are preferential places for the formation of an intense current layer.

We checked that no stagnation point for the velocity ever forms in the vicinity of the HFT in any of our simulations. This is again contradictory to the restricted conditions that Galsgaard et al. (2003) found for current sheet formation in HFTs. This quantitative difference is probably due to our much thinner pre-existing QSLs combined with the absence of special symmetry properties in our models. All our magnetic field lines are rooted in one single line-tied plane (whereas Galsgaard et al. 2003, considered a straightened configuration between two opposite line-tied plates), and only one of our four polarities is located in the boundary flow region (whereas all four polarities are displaced in Galsgaard et al. 2003). The precise dynamics and geometry of the HFT current sheets in our simulations are still controlled by the form of the line-tied motions, as shown in Fig. 6. The steepening of the current layer is mostly due to local compressive shearing motions, which result from a combination of the different vertical expansions and of the horizontal rotations of the field lines across the HFT, since they have different sizes and are not rooted in the same regions at z=0. It is then natural that Galsgaard et al. (2003) could not obtain this behavior and thus needed to create a stagnation point so as to create a current sheet at the HFT, considering the absence of both short and long field lines in their straightened magnetic field configuration.


  \begin{figure} \par\mbox{\psfig{figure=3600f9a.ps,width=4.1cm,angle=-90,bbllx=73... ...le=-90,bbllx=73pt,bblly=187pt,bburx=523pt,bbury=695pt,clip=} }\par\end{figure} Figure 9: Plots of the three components of $\vec{b}$ and $\vec{j}$ along z at (x, y)=(-0.02, 0.07). The whole domain along z is shown. In the first and second rows, (continuous-red, dashed-blue, long-dashed-green) lines respectively correspond to (bx, by, bz). In the third row, (continuous-red, dashed-blue, long-dashed-green, thick-continuous-black) lines respectively correspond to ($\jmath _x$, $\jmath _y$, $\jmath _z$, $\jmath $). In both cases, the location of the smallest scale for $\vec{b}$ associated with the narrowest peak of $\jmath $ corresponds to the HFT. Note that in this region, when the same small scale is achieved, the jump in $\vec{b}$, so the current magnitude, is larger for initially wider QSLs (i.e. for $\Phi =120^\circ $ shown in the left column).
Open with DEXTER

Electric current and magnetic field profiles along z for fixed (x,y) positions passing through the middle of the HFT are shown in Fig. 9; from these plots, one can estimate more quantitatively what the greyscale levels correspond to in Figs. 6 and 7. The potential field profiles are also drawn for comparison. These plots are comparable in the sense that they correspond to the formation of similar small scales in the HFT. These plots suggest that for a given magnetic configuration, the broader the QSLs are for its potential field, the longer the twist can be applied on the boundary and the higher the electric currents can be generated in the HFT, before the latter reach the scale of the mesh, i.e. the dissipative scale.

In all our simulations, provided that the viscous term was well adapted, the steep magnetic field gradients which progressively form in the HFT invariably caused numerical instabilities after several tens of Alfvén times, which eventually halted the simulations. We verified that increasing $\eta$ permits to further evolve the systems for longer times. But we did not pursue in this direction, since the aim of this paper was to study the formation of current layers and their possible collapse at the scale of the mesh, with reduced diffusive effects.

The Lorentz forces are the strongest at the HFT. First analyses show that once the scale-lengths are small enough, the Lorentz forces lead to an undriven collapse of the current layer, and they accelerate the plasma at its outer edges for non-zero resistivities. This results in a magnetic reconnection-like process, but the related change of connectivity during the diffusion is not discontinuous. Instead the field lines tend to slip along each other on both sides of the HFT, while their footpoints at z=0 quickly shift along the QSLs over long distances. Field line slippage was in fact first envisioned in the general context of magnetic reconnection with no null point by Priest & Démoulin (1995). Theoretical arguments for it were developed by Priest et al. (2003). It was only recently identified in non-zero $\beta$ MHD simulations of reconnection within a thick HFT (Pontin et al. 2005) and between sheared arcades in the frame of prominence modeling (Aulanier et al. in preparation; DeVore et al. 2005). We thus believe that this specific behavior is the generalization magnetic reconnection from 2D to 3D, when neither null points nor separators are present in the system. Detailed analysis of this process in zero-$\beta$ for our modeled configurations will be the object of a forthcoming paper following the present one.

   
5 Discussion

5.1 Summary 

We considered two quadrupolar configurations (Figs. 1 and 2). They only differed by the respective angle made between the axes of the large outer and the small inner bipoles within one configuration. In spite of their quadrupolar nature, these configurations were, strictly speaking, bipolar. They did not possess separatrices. They still had strong gradients of field line connectivity in regions called QSLs (Figs. 3 and 4).

We considered two types of line-tied boundary motions, with zero-$\beta$ resistive MHD simulations, using a $191 \times 161 \times 170$ non-uniform mesh. These motions were prescribed so as to displace only the field line footpoints within one of the polarities of the inner bipole, either by translation, or by twisting motions (see Fig. 5). Their maximum velocity was very sub-Alfvénic, allowing tens of wave reflections from one footpoint to another. They led to the advection of field lines over distances that were small compared to the characteristic scale-lengths of the configurations (see Fig. 1). Their gradients had typical scale-lengths which were between that of one single polarity and that of the whole quadrupolar configuration. These flows neither possessed any stagnation point at the line-tied boundary nor did they result in the formation of stagnation points in the domain as a result of the MHD evolutions.

The prescribed motions firstly resulted in the development of extended quasi force-free currents. The location and amplitude of these currents were directly related to the form of the motions, as is the case in all line-tied magnetic field simulations.

The key result is that these motions also invariably resulted in the formation of very narrow current layers, even though no true magnetic separatrices were present in the systems (see Fig. 6). These narrow current layers were always cospatial with the QSLs for various footpoint motions (see Fig. 7). Most QSLs already existed in the potential fields, and the evolution of their shapes mainly resulted from the advection by the boundary motions. Some secondary QSLs also formed where the boundary motions had the steepest shear gradients. Current layers naturally developed in these QSLs as well.

The thin volume corresponding to the highest squashing degree Q of the QSLs had a specific shape which led us to call it a HFT. For long enough motions, the strongest and narrowest current layer developed around the HFT (see Figs. 6 and 8), even though no stagnation point ever formed in this region. In typical magnetic configurations that possess separatrices, a current sheet is known to form with almost all kinds of boundary evolution at the separator, or at the null point if the 3D separatrix is only made of a fan surface and a singular spine field line. The current layer forming at the HFT is a generalization of the latter for configurations without separatrices, but with QSLs instead. When the magnetic field gradients reached the scale of the mesh, numerical instabilities developed as a natural result of the formation of unresolved gradients in this region. This instability could only be prevented by increasing the resistivity. Comparisons of several configurations have shown that the wider the QSLs were in the potential field, the stronger the currents became in the HFT before they reached the dissipative scale (see Fig. 9).

Since we varied both magnetic field configurations and footpoint motions, we started exploring the parameter space. The generic characteristics of our results on the development of electric currents suggests that they must also be valid in any magnetic configurations that have thin QSLs.

   
5.2 A model for solar flares: topology, energy build-up and switch-on of reconnection

Our results have strong implications for the physics of solar flares in general. Flare models that are based on magnetic reconnection in narrow current layers can be divided into two main classes. The first class involves a large-scale MHD instability (e.g. Amari & Luciani 1999) or a global non-equilibrium (e.g. Forbes 2000), which results in a fast flux tube deformation on time-scales that can be Alfvénic. The latter leads to strong vortical and/or compressive motions, which naturally results in the dynamic formation of narrow current layers and in the triggering of reconnection, with or without complex topologies in the pre-flare configuration. Since we have considered only modest footpoint motions, our simulations are not directly relevant to these models. The second class of flare models considers the slow buildup of large current sheets in separatrices (Somov 1992, and references therein). They predict that the width of the current sheets that spontaneously form in separatrices tends to zero in the limit of infinite Lundquist number (e.g. Lau 1993; Aly 1990). So they require that the current sheets do not diffuse for a long time (which is problematic, as discussed by Low & Wolfson 1988), before reconnection is switched on due to the triggering of plasma (or MHD) instabilities within the current sheets when some threshold is reached. Our simulations extend the latter models, and provide natural solutions to their difficulties.

In magnetic configurations which initially contained broader QSLs, the electric currents in the HFT increased to higher magnitudes when the magnetic field gradients reached the dissipative scale (Sect. 4.3). Then if the fast energy release is not the result of a global instability, as in our simulations, the narrower the initial QSLs are, the shorter the time it takes to reach the dissipative scale and the less energy is likely to be accumulated before (and released during) a flare-like event.

We then argue that the most energetic solar flares that are not triggered by an early large-scale ideal instability must occur in magnetic configurations whose corresponding potential field have broad QSLs. This is rather counter-intuitive if one considers the long history of the separatrix-related flare models mentioned above, which involve the formation of long current sheets, which are spontaneously infinitely thin, during the energy build-up phase.

Let us now rescale our models to solar units for an active region. Typically, distances between P2 and N2 should be $\sim$20 Mm, photospheric velocities should be $\sim$0.1 km s-1, and Alfvén speeds should be $\sim$103 km s-1. In this context, the Alfvén time is $\sim$20 s and the photospheric velocity is $\sim$10-4 of the Alfvén speed. In our simulations for $\Phi =120^\circ $, the currents in the HFT reached the scale of the mesh in $\sim$102 Alfvén times, and we used a line-tied velocity of $\sim$10-2 of the Alfvén speed. The energy build-up phase converted to solar units should then be of the order of $\sim$104 Alfven times, i.e. $\sim$2.3 days. This is of the order of the observed time-scales. We then propose that the above estimations, combined with the slow driven gradual steepening of the magnetic field gradients in the vicinity of a HFT, until they reach small dissipative scales, permits to solve the long standing paradigm for both the long-duration energy storage before a flare takes place, and for the switch-on of magnetic reconnection during the impulsive phase of the flare. This has been, in fact, one of the main problems in line-tied separatrix-related models, as discussed by Low & Wolfson (1988).

Our results then support and extend past works that associate temporal and spatial properties of observed solar flares with QSLs computed from magnetic field extrapolations. When the magnetic configuration has a low free magnetic energy stored or when the configuration is strongly quadrupolar, the potential field extrapolations of observed photospheric magnetograms and calculation of the resulting QSLs are sufficient to predict where a flare can potentially take place (Démoulin et al. 1997; Gaizauskas et al. 1998). When the distributed currents are important (so the free energy is high) and the configuration is more bipolar, force-free field extrapolations are needed to determine the location of the QSLs with more accuracy, so the flare location (Bagalá et al. 2000; Mandrini et al. 1996; Schmieder et al. 1997). Extrapolations and QSLs should then also be useful for predictingthe acceleration sites and trajectories of solar energetic particles in flares.

   
5.3 The role of QSLs in coronal heating

Many models exist for heating the corona by the dissipation of thin current layers, as extensively reviewed in Mandrini et al. (2000). The related currents can either be of the AC (alternative current) or of the DC (direct current) type. Only the latter are related to low-frequency perturbations such as sub-Alfvénic line-tied footpoint motions. Presently, observational constraints permit to select the most relevant models (Démoulin et al. 2003; Schrijver & Title 2005). DC type models are among the ones which fulfill observational requirements.

Recently, Gudiksen & Nordlund (2005) have performed MHD simulations of turbulent flux braiding in a potential field that was extrapolated from an observed magnetogram. The development of narrow current layers in their simulations was due to the line-tied driving. Since their boundary flows follow a turbulent power-spectrum, it should naturally create small scales and stagnation points in the velocity profiles, as directly prescribed in the simulations of van Ballegooijen (1986) and Galsgaard et al. (2003). Thus, one may argue that the topology of the flow was directly at the origin of the current layers which develop within the coronal volume.

Here we propose that another effect might play a non-negligible role in this particular simulation, and on the Sun more generally, based on the idea that well developed active regions are typically composed of numerous flux concentrations. Even with strictly bipolar and potential configurations, Démoulin & Priest (1997) found very thin QSLs, when several flux concentrations were embedded in a (non-zero) weaker vertical field background. They show that the QSL thickness strongly depends on the intensity of this background. In the present paper, we have shown numerically that narrow current layers spontaneously develop in such QSLs, even though we considered much simpler configurations with only two flux concentrations on each side of the inversion line. So we believe that at least some of the current layers in Gudiksen & Nordlund (2005) must be associated with QSLs defined by the magnetic flux concentrations at the boundary, rather than with the topology of the boundary flows, which must anyway create others QSLs, even if the magnetic field is initially homogeneous.

This conjecture and its associated time-scales should be tested in the future. If QSLs associated to the flux concentrations dominate in general, their calculation in potential (or force-free) field extrapolations of any high resolution magnetogram could be a good proxy not only for the occurrence of solar flares, but also for the locations where coronal heating results in the illumination of discrete loops in EUV images of the corona. Wang et al. (2000) and Fletcher et al. (2001) provided first observational evidences of this.

   
Appendix A: Numerical and physical issues regarding resistive and viscous effects in separatrices and QSLs

For the bipolar magnetic configurations that we studied, especially for $\Phi =150^\circ $, but also for more general configurations, the profile of the squashing degree Q can be very strongly peaked (Titov et al. 2002; Démoulin et al. 1996a,b). This peak is, in fact, infinitely high and thin in the case of separatrices. The related full width at half maximum of the Q profiles can then often be orders of magnitudes below the grid resolutions presently achievable in MHD simulations. We argue that in an ideal plasma, the continuous characteristics of the ideal MHD equations (with neither resistive nor viscous effect) should physically advect the field lines in the volume according to the line-tied boundary motions. We believe that this physical advection should then typically form thin current layers, with a thickness comparable to the QSL thickness (zero in the case of a separatrix, see e.g. Lau 1993; Low & Wolfson 1988; Aly 1990). The width of the current sheets that form is of major importance for flare and reconnection modeling, since it clearly determines the diffusive time, thus the duration of energy storage (see Sect. 5.2) and the reconnection rate (see e.g. Petsheck 1964; Sweet 1958).

If the physical behaviour described above is also true when discretized equations are used in a numerical simulation, the line-tied motions should try to form sharp current layers at the QSLs on a time-scale of the order of the travel time of Alfvén waves along a QSL field line, and at a scale below the grid resolution, thus further leading to a quick numerical instability over a few time steps. This fast instability clearly does not happen in our simulations. It does happen in the HFT, but only after tens of Alfvén time units. Outside the HFT, but still within the QSLs, all our simulations result in current layers that are indeed narrow, but that are still resolved and far wider than the QSL itself (see Fig. 8 and Table 2). The same behavior is found in many separatrix line-tied 2.5D MHD simulations (e.g. Ma et al. 1995; Longcope & Magara 2004), which should physically result in spontaneous zero-width current sheet formation (as discussed by Low & Wolfson 1988). In the case of MHD numerical simulations, three questions therefore arise:

(1)
How is the discretized system of MHD equations sensitive to the QSL small-scales (or separatrix zero-scale) that already exist below the mesh resolution at t=0?
(2)
How do the resistive and viscous terms control the width of the current layers that form around separatrices and QSLs?
(3)
How reliable are line-tied simulations of unresolved QSLs and separatrices in the calculation of time-scales associated with reconnection (duration of the energy storage and reconnection rate)?
In order to discuss these major issues, we performed two other simulations (hereafter labeled "S1'' and "S2'') of our sharpest case: $\Phi =150^\circ $ with twisting motions. Both simulations were done with exactly the same settings as the original (let us call it "S0''), except that $\eta=0$ in S1 and $\eta=0$ and  $u_{\nu^\star}(t>13)=0$ in S2. S1 surprisingly gave very similar results to S0. In particular, the narrow current layers that formed in the QSL, although they were slightly more intense, spontaneously developed on the same (resolved) scale, for tens of time units. The main difference is that S1 was halted earlier than S0, at about $t \sim 31$ instead of 40, because of numerical instabilities developing in the HFT, not elsewhere in the QSL. This earlier instability was expected off course, since the current build up at the HFT was not weakened by diffusion. However the long duration of the ideal simulation was not expected according to the idea that magnetic diffusion was not there to ensure that the current layers within the QSLs remained resolved. S2 remained stable for about one Aflvén time after the viscous term was set to zero at t=13, and fast growing numerical instabilities stopped the run by $t \sim 15$. Even though common sense tells us that setting all diffusive terms in an MHD simulations cannot result in anything else than a numerical instability, we still performed this run to see where and how the instability first developed. We found that it did not occur in the HFT, but on both sides of the magnetic configuration, around $(x,z)=(\pm 0.25,0.18)$. These regions did not correspond to the strongest squashing degree Q, but rather to the regions where the ratio between the (non-uniform) grid size to the width of the unresolved QSL peak was the highest.

Comparing the three simulations permits us to give some answers to the questions addressed above.

(1)
Since with no diffusive term, numerical instabilities first and very quickly develop where the QSLs are the least resolved, we argue that our MHD code is somehow aware of the field line connectivities, even if the latter have strong unresolved gradients. Since our numerical scheme (described in detail in Aulanier et al. 2005) has no explicit feature to ensure this, we believe that it is a general property of all explicit codes that do not incorporate any important scheme-based diffusive effects.
(2)
Since a simulation with only viscosity is barely distinguishable from another one with both viscous and resistive effects, it appears that not only resistivity, but also viscosity, is a determining factor in the broadening of current sheets. The effect of the latter is in fact to smoothe the velocity profiles that result from the Aflvén waves that transport the information from one field line footpoint to another. So the shearing profile of the magnetic field on both sides of the QSL at z>0 is artificially broadened, regardless of any resistive effect. This issue is not mentioned often in numerical simulations of separatrices and null points, which tend to focus their parameter study on the value of the Lundquist number.
(3)
Since viscosity and resistivity both limit the width of the current layers (to achieve long runs), it seems that the slow amplification of electric currents in magnetic field gradients that are on the scale of the mesh (i.e. the pre-flare energy build-up phase) is fully conditioned by the diffusive terms, which are always much higher than in a real plasma. Hence, the current sheet gradual stretching that is calculated in unresolved HFTs (see Fig. 8, middle panel) or at the intersection of separatrices (e.g. Ma et al. 1995; Longcope & Magara 2004; Biskamp 1986) may then be unphysical. This could have strong consequences in the evaluation of reconnection rates.
In a real physical system such as the solar corona, however, the continuous MHD equations are applicable to a much broader range of spatial scales than in present numerical simulations. We conjecture that the true width of the current layers within QSLs should then be defined by the true profile of the squashing degree Q (as in Fig. 8, left panel), while the boundary motions should define the precise current magnitude within these thin layers. The central peak of Q can have very high values, tending to infinity for a true separatrix, and the peak width can be extremely narrow, tending to zero for a true separatrix). Physical limitations to achieve these small scales should then be micro-physics (diffusive) effects and finite time evolutions, if the characteristic wave propagation speed is locally slow, e.g. in the vicinity of a magnetic null point. More generally, this implies that the gradual amplification of electric currents in QSLs can only be possible as long as the later are broader than the physical dissipative scale, whatever its precise nature (resistive or collisionless). This suggests that slow line-tying motions around separatrices cannot account for anything more than steady coronal heating, since there, the currents readily form at the dissipative scale.

All these conjectures are difficult to test quantitatively, since resolving thin QSLs numerically is difficult. For very thin QSLs, MHD simulations should be performed with very small grid scales (see the values of  $\delta _{Q, {\rm {H}}}$ given in Table 2), which are hardly reachable with present computers facilities. Adaptive mesh refinement techniques are plausible ways to study a much larger range of scales. For wider QSLs, the amplitude of footpoint motions required to generate enough currents in the initial QSLs may be so large that new QSLs may form earlier on smaller scales, due to the prescribed velocity-field topology at the boundary. The newly formed QSLs can become much thinner than the initial QSLs, thereby dominate the current buildup. This is what happened in the simulations of Galsgaard et al. (2003), as discussed in Démoulin (2005). Since our configuration for $\Phi =120^\circ $ has QSLs that remain resolved during about one fourth of the duration of our simulations, a compromise between both possibilities should be possible to investigate in a close future.

Acknowledgements
The authors thank C. R. DeVore and R. Grappin for fruitful discussions. The calculations in this paper were done on the Digital-UNIX ES40 and on the Compaq-HP Quadri-Opteron computers of the Service Informatique de l'Observatoire de Paris (SIO).

References

 
Copyright ESO 2005