Numerical Simulation of Current Sheet Formation in a Quasi-Separatrix Layer using Adaptive Mesh Refinement

The spontaneous formation of thin current sheets is believed to play an important role in the dynamics of astrophysical plasmas like the solar corona and, in particular, for the onset of magnetic reconnection. In fact, impulsive events like solar flares and coronal mass ejections are often associated with magnetic reconnection as a mechanism that effectively releases magnetic energy, which in turn calls for highly concentrated electrical currents to explain the breakdown of ideal plasma behavior by means of e.g. micro-instabilities in the non-collisional coronal plasma.

Recently, the study of magnetic field configurations as candidates for reconnection has shifted from separatrix surfaces towards quasi-separatrix layers (QSLs) Démoulin et al. (1996a, b); Milano et al. (1999). In contrast to genuine separatrix configurations, QSL fields do not necessarily contain magnetic null points in 3D, which makes them relevant in many more situations than strict separatrices. QSLs describe a magnetic field mapping between two boundaries that is continuous, but changes rapidly in space. This change has been quantified in terms of a flux tube “squashing factor” $Q$ Titov and Hornig (2002) and more recently generalized to remove boundary projection effects Titov (2007). In accordance with the significance of QSLs for magnetic reconnection, the problem of current sheet formation has been investigated. Here, Titov et al. Titov et al. (2003) and Galsgaard et al. Galsgaard et al. (2003) studied the formation of current sheets in a straight hyperbolic flux tube (HFT) configuration both analytically and numerically, and found exponential growth of the current density as a response to shearing magnetic footpoint motion. One major conclusion of that work was that a shear in the applied boundary motion was important, if not essential, for the formation of thin current sheets, and that the creation of a stagnation point flow in the HFT was the key element in this process. However, this has been later questioned by Aulanier et al. Aulanier et al. (2005): They argued that the initial squashing factor $Q$ in the simulations of Galsgaard et al., probably together with additional symmetry properties, was too small to account for highly concentrated currents, and that the shear boundary motion and stagnation point flow in Galsgaard et al. (2003) served as to dynamically create thin QSLs only during the simulation. This argument was underpinned in Aulanier et al. (2005) through extensive comparison between MHD simulations of less symmetric potential magnetic field configurations which initially result from four magnetic point sources submerged below a photospheric boundary. They contain QSLs with squashing factors of up to $Q\approx 10^{5}$ and were exposed to different boundary driver patterns, namely a shear and a translational motion. Aulanier et al. observed the formation of thin, intense current sheets, almost irrespective of the type of applied boundary driving, which stresses the relevance of the initial field geometry and QSL strength, rather than the boundary motion, for current sheet formation. Later, that work has been extended by finite resistivity to simulate magnetic reconnection at that thin current sheet Aulanier et al. (2006) with its strong temporal change in magnetic connectivity.

A general limitation of these previous numerical studies was the lack of reasonable numerical grid resolution at the sites of current sheet formation Aulanier et al. (2006): During the formation process, the currents get more and more concentrated locally so that the sheet scales quickly reach the numerical grid spacing. This left open the question whether the current concentration would continue to length scales small enough to account for the onset of micro-instabilities, or if this process would get arrested at some larger length scale. In order to address this question, but also get more insight into the local dynamics in the immediate vicinities of the current sheets, we have carried out numerical simulations in the spirit of Aulanier et al. (2005), using one of their initial configurations and boundary drivers. However, we used the technique of adaptive mesh refinement (AMR), implemented in the code racoon Dreher and Grauer (2005), in order to obtain a much higher grid resolution in the current sheet vicinity than was reported earlier. This allowed us to obtain insight into the formation process that reaches beyond the previous results.

In the next section, we will briefly describe the model that we employed in our studies, while results that complement the work of Aulanier et al. (2005)   are given in Sec. III, followed by a summary and conclusion in the last section.

Our simulations are based on a reduced subset of the MHD equations, appropriate for the quasi-static, low-$\beta$ regime, where $\beta$ is the ratio of thermal to magnetic pressure. In normalized quantities, it reads

Pressure terms are omitted from the momentum equation (1a) because of the low-$\beta$ approximation. As we are interested in a quasi-stationary evolution of the system, i.e. the limit of high wave speeds, we replace the continuity equation by a direct prescription of a relaxation mass density, namely $\rho:=B^{2}$. This approach results in an homogeneous Alfvén velocity $c_{A}:=|B|/\sqrt{\rho}=1$ and fast communication of unbalanced forces throughout the system, which shortens relaxation times and has been successfully used in other studies Arnold et al. (2008). Constant kinematic viscosity $\nu=5\cdot 10^{-4}$ and resistivity $\eta=5\cdot 10^{-6}$ are included only to guarantee numerical stability on the grid scale and have little effect on the overall evolution. The artificial scalar function $\Psi$ and its dynamic equation (1d) serve as a convenient means to constrain any finite $\nabla\cdot\bm{B}$, resulting from discretization errors, to negligible values Dedner et al. (2002): Combining $\partial_{t}\nabla\cdot(\ref{eq:induction})$ and $\Delta(\ref{eq:cleaning})$ results in the mixed equation

for $\nabla\cdot\bm{B}$. The crucial term on the right side is the first Laplacian, which describes a hyperbolic transport of $\nabla\cdot\bm{B}$ with velocity $c_{h}$ and leads to radiative distribution of $\nabla\cdot\bm{B}$ throughout the computational domain, while it gets dissipated by phase mixing and diffusion according to the other two terms. Note that there is freedom in a specific choice for the $\Psi$-equation (1d) as it is of the order of the discretization error anyway. For instance, while Dedner et al. Dedner et al. (2002) discuss the term $-\frac{c_{h}^{2}}{c_{p}^{2}}\Psi$ (see their Eq. (19) resp. (24e)) to arrive at a telegraph equation for $\Psi$ and $\nabla\cdot\bm{B}$ in the continuous case, we found this unnecessary, although possible, for our computations. The main reason for this seems to be that in our case, the sources of $\nabla\cdot\bm{B}$-errors are highly localized regions in space, namely the regions of intense currents, so that the hyperbolic transport is the dominating cleaning effect. On the other hand, we added the term $c_{l}\Delta\Psi$ in Eq. (1d), which is not discussed in this specific form in Dedner et al. (2002). Its motivation is, however, not so much a change in the $\nabla\cdot\bm{B}$-cleaning property itself, but the observation that the purely hyperbolic choice, i.e. Eq. (1d) with $c_{l}=0$, tends to introduce odd-even-decoupling on the centered finite difference grid that we used. This decoupling, which was not an issue in Dedner et al. (2002) since they employed finite volume schemes in finite element discretisations, could be healed satisfactorily through the additional Laplacian term which couples $\Psi$-values on neighboring grid points with each other. We finally found overall good $\nabla\cdot\bm{B}$-cleaning properties when choosing the numerical parameter values to $c_{l}=5\cdot 10^{-4}$ and $c_{h}^{2}=5\cdot 10^{-2}$.

The equations are discretized in a domain of $(x,y,z)\in[-0.7,0.7]\times[-0.5,0.5]\times[0,0.5]$, where we identify the plane $z=0$ with the solar photosphere. Integration is performed with a strongly stable third-order Runge-Kutta scheme Shu and Osher (1988), spatial differentiation is realized with standard second-order finite differences. In order to resolve the expected small-scale features, we employed the block-structured adaptive mesh refinement (AMR) framework racoon Dreher and Grauer (2005). Here, we used the norm of the magnetic field gradient, $\nabla\bm{B}$, to derive a local length scale for $\bm{B}$ that serves as an indicator criterion for local mesh refinement. Effective local grid resolutions obtained with this technique were $4096^{3}$ in the present work.

Initial conditions are adapted from Aulanier et al. (2005), where we address the magnetic field configuration resulting from two pairs of opposite polarity photospheric flux patches whose respective connecting axes intersect at an angle of 150 degrees. Specifically, the initial magnetic field stems from 4 virtual magnetic point sources, indexed by $i$, below the photosphere

with respective source strengths $F_{1}=-F_{2}=1$ and $F_{3}=-F_{4}=0.4$ and locations $\bm{r}_{1,2}=(\pm\frac{1}{2},0,-\frac{1}{5})$ and $\bm{r}_{3,4}=(\mp\frac{\sqrt{3}}{20},\pm\frac{1}{20},-\frac{1}{10})$, respectively. This field geometry is known to contain two symmetric dome-shaped QSLs with squashing factors of $Q\approx 10^{5}$, intersecting in a hyperbolic flux tube (see Aulanier et al. (2005) for details).

In the course of the simulation, the field is exposed to a horizontal photospheric vortex flow around the magnetic source $i=3$ in Eq. (2), realized by prescribing the boundary condition for $\bm{v}$ at $z=0$. The maximum flow velocity is $\max(|\bm{v}_{Bnd.}|)\approx 2\cdot 10^{-2}$ and it gets ramped up in time according to $\propto\frac{1}{2}\left[\tanh(\frac{5}{2}(t-1))+1\right]$.

The detailed treatment of the lower boundary is as follows: As all quantities are discretized as cell-centered variables, boundary conditions have to provide values for $\bm{v}$, $\bm{B}$ and $\Psi$ at $z_{1/2}:=\Delta_{z}/2$ in a way that is consistent with the evolution equations and the overall second order accuracy in the grid spacing $\Delta_{z}$. Denoting the boundary values of individual quantities at $z=0$ with a hat, the velocity $\hat{\bm{v}}=\bm{v}_{Bnd.}$ is explicitly given from the prescribed horizontal vortex flow, and in particular, $\hat{v}_{z}=0$. All components of $\bm{v}$ at $z_{1/2}$ can now simply be interpolated in $z$-direction with second order accuracy between $z=0$ and the interior grid points above $z_{1/2}$. Using this direct forcing, there is no need to compute $\bm{j}$ or the Lorentz force at $z_{1/2}$ at all, hence we do not need the horizontal components of $\hat{\bm{B}}$ at this stage. Now, to integrate $\bm{B}$, we note that the convective part of the $B_{z}$-equation becomes autonomous in the boundary plane, involving no other undetermined quantities nor any derivatives normal to the boundary: Writing the convective electric field as $\bm{\mathcal{E}}:=-\bm{v}\times\bm{B}$ gives $\hat{\mathcal{E}}_{x}=-\hat{v}_{y}\hat{B}_{z}$ and $\hat{\mathcal{E}}_{y}=\hat{v}_{x}\hat{B}_{z}$ so that $\partial_{t}\hat{B}_{z}=-\nabla_{h}\cdot(\bm{\hat{v}}_{h}\hat{B}_{z})$ where the index $h$ indicates horizontal vector projections, e.g. the prescribed $\bm{\hat{v}}_{h}=(\hat{v}_{x},\hat{v}_{y})$, and $\nabla_{h}:=(\partial_{x},\partial_{y})$. In other words, when ignoring the numerically motivated resistive diffusion and $\nabla\Psi$-terms, we are able to integrate the “proper” $\hat{B}_{z}$ completely for its own, which in turn allows us again to interpolate $B_{z}$ at $z_{1/2}$ up to $O(\Delta_{z}^{2})$, similar to $\bm{v}$ above. Now, this “proper” $\hat{B}_{z}$-equation also determines $\bm{\hat{\mathcal{E}}}_{h}$, which allows to update $\bm{B}_{h}$ at $z_{1/2}$. At this point, the magnetic diffusion and $\nabla\cdot\bm{B}$-cleaning terms are taken into account as usual, where the former requires an extrapolation of $\bm{B}_{h}$ across $z=0$. Finally, the right side of Eq. (1d) is easily evaluated at $z_{1/2}$, using $\hat{B}_{z}$ from above and the Dirichlet condition $\hat{\Psi}=0$.

While the upper and lateral boundaries can in principle be handled in the same fashion, using the no-slip and no-penetration condition $\bm{v}=0$, we used a simpler approach there: Setting $\bm{v}$ to zero ahead of the boundaries, keeping the tangential components of $\bm{B}$ fixed in ghost cells and applying solenoidal and homogeneous Dirichlet conditions to the normal components of $\bm{B}$ and $\Psi$, respectively, proved to be sufficient for those passive boundaries. Note, in particular, that no artificial Lorentz forces act there either.

The boundary treatment described above presumes a homogeneous numerical grid. It has been applied in the same spirit to the AMR simulations that we present here, where additional complications occur at the interfaces between neighboring grid blocks of different mesh resolution that abut the physical boundaries. Without going into details, we only mention that additional coarse-fine and fine-coarse interpolations are needed for the computation of $\hat{B}_{z}$ and $\hat{\bm{\mathcal{E}}}$ at those junctions, but that they do not pose any fundamental new challenge apart from the programming complexity.

After applying the photospheric boundary driving in $\bm{v}$, dynamic shear modes travel into the domain and mix there. On the scale of a few Alfvén transit times, the magnitude of the current density grows significantly and a quasi-stationary current system as shown in Fig. 1 builds up.

It consists partly of relatively weak currents which are distributed on a large scale in a dome-like structure that is pre-determined by magnetic field lines connecting the driver region with the opposite polarity regions of the photosphere. On top of this, a highly localized thin current sheet can be identified in the vicinity of $(x,z)\approx(0,0.18)$ in Fig. 1. This thin current sheet actually lies inside the pre-existing QSL of the initial magnetic field. We interpret the striation patterns at $x\approx-0.4$ in Fig. 1 as signatures of MHD waves on field lines which connect the strong photospheric field region with the current sheet during the evolution.

These results are basically in good agreement with those published earlier by Aulanier et al. Aulanier et al. (2005). However, the sheet thickness at the stage shown in Fig. 1 is already on the numerical grid scale of Aulanier et al. (2005), so that their studies were unable to investigate into the further evolution or features on smaller scales.

The temporal evolution of the sheet is displayed in Fig. 2 which shows the vertical profiles of $|\bm{j}|$ at $x=y=0$ for four different times. It is evident that the sheet thickness decreases with time while the value of maximum current density increases accordingly to reach values of $\approx 3\cdot 10^{3}$ in the latest stage. At the same time, the sheet moves upwards, as indicated by the inset graphs. In fact, we find that the entire magnetic structure expands gradually as a response to the boundary driving, and the current sheet as a substructure is embedded in this motion. An other detail that emerges from Fig. 2 is the fact that up to $t\approx 4.4$, the most intense currents are not yet found in the thin current sheet itself, but in the large-scale system at $z\approx 0.02$, i.e. close to the photospheric driver (see also Fig. 1). It is only at later times, that the thin sheet dominates in the current intensity.

This phenomenon is also visible in the temporal evolution of $\max(|\bm{j}|)$, which is plotted logarithmically against time in Fig. 3: The early phase, with rather fast growth of $\max(|\bm{j}|)$, corresponds to the ramp-up of the boundary driver, which essentially reaches its maximum magnitude around $t\approx 1.7$. This is followed by a slower growth rate of the current maximum up to $t\approx 4.7$. During this stage, the maximum values stem from the extended currents close to the photosphere (compare with Fig. 2), which eventually get overtaken by the faster growth of the thin embedded current sheet. Further intensification continues, with amplification of $\max(|\bm{j}|)$ by roughly a factor of 5, until the growth slows down at $t\approx 5.5$. At this time, the sheet thickness is only a few times the numerical grid scale and thereby poorly resolved with artificial diffusion effects becoming competitive.

Fig. 4 shows details of the sheet in the cut plane $y=0$ for three different times. Again, the overall upward motion and the thinning and intensification of the sheet are well visible. In addition, we have plotted the plasma velocity as arrows, projected into the displayed plane, to give an impression of the flow in the vicinity of the current sheet. While the upper row shows the velocity relative to the fixed simulation frame, the plots in the lower row show the flow transformed to a frame which is co-moving with the expansion velocity of the structure. To this end, we identified the locations of maximum current density in the plane $y=0$ from plots of successive data output sets, and computed a pattern velocity from their difference. This velocity was then subtracted from the plasma flow velocity in the lower plots of Fig. 4. There have been controversial discussions as to whether the current sheet formation at the hyperbolic flux tube embedded in the QSL is related to a stagnation-point flow Galsgaard et al. (2003); Aulanier et al. (2005). In particular, Aulanier et al. claim that no stagnation point exists at the intense current sheet. This is certainly confirmed by our simulation, however we remark that the focus on a strict definition of stagnation point, i.e. $\bm{v}=0$, maybe somewhat misleading because i) the velocity is sensitive to the chosen frame of reference, and ii) the component along the current sheet should be discarded from these considerations anyway, because it largely decouples from the mechanism of magnetic shearing discussed in Titov et al. (2003) and Galsgaard et al. (2003).

The flow pattern projected into the $y$-$z$-plane is shown in Fig. 5, again transformed into a frame that moves upward with the current sheet and, in addition, results in $v_{y}=0$ in the current sheet center. This figure also demonstrates that the current sheet is indeed elongated in the $y$-direction. Hence, at least in the latest stage displayed in Figs. 4 and 5, it can be treated as a quasi-one-dimensional sheet.

Finally, we remark that the assumption of a quasi-stationary evolution loses its validity in the late stage of the sheet evolution: Obviously, the collapse becomes a localized, dynamic process associated with significant magnetic forces. This can be seen from the field line plot shown in Fig. 6, where magnetic field lines connecting the thin current sheet with the photosphere have been colored with the quantity $\alpha:=|\bm{j}|/|\bm{B}|$. For a force-free field, $\bm{j}\times\bm{B}=0$, the value of $\alpha$ is constant along magnetic field lines. This condition is obviously not met in the QSL, which means that the currents close locally across field lines.

We carried out numerical simulations of current sheet formation in a quasi-separatrix layer using a simplified MHD model appropriate for the quasi-static evolution of a low-$\beta$ plasma. The setting under consideration has been investigated before Aulanier et al. (2005) and our results agree well with that work as long as the current sheet structure is well resolved in both studies. With the use of local adaptive mesh refinement (AMR), we were able to follow the thinning of the current sheet further down to a scale which is about one order of magnitude smaller than previously investigated. In particular, our simulations reached a stage in which the maxima of $|\bm{j}|$ in the upper part of the QSL grow significantly beyond the values close to the photospheric boundaries, which gives clear evidence that the most intense current densities actually can be expected in the upper part of the QSL. This late stage is characterized by a relatively fast collapse of the locally almost one-dimensional sheet with an approximately exponential increase of $\text{max}(|\bm{j}|)$ in time, and the evolution is no more quasi-static at this point.

When magnetic forces become significant in this late stage of the current sheet formation, full nonlinear MHD dynamics will take place. Previous studies have addressed details of the local dynamics of such current sheets using appropriate initial conditions and periodic systems (e.g. Grauer and Marliani (2000) and references therein). There, one particular question has been whether the current density might form a singularity in finite time, or whether its growth is limited to merely e.g. exponential behavior. On theoretical grounds, it could be shown that a dynamical alignment between the velocity and the magnetic field would bound $|\bm{j}|$ to exponential growth. Analyzing our data further, we actually found indications of such an alignment (not shown here), so that we expect to see a collapse of the sheet with exponential growth, i.e. a continuation of the phase observed between $t\approx 4.7$ and $\approx 5.5$ in Fig. 3, given that it could be resolved numerically. At present however, even our AMR simulations are limited by the lack of further resolution and by numerical side-effects like artificial stabilizing diffusion.

The plasma flow pattern has been analyzed in a cut plane that is approximately perpendicular to the local direction of the current density in the sheet: If transformed into a frame that moves with the expanding structure, it exhibits a shear pattern which arises from a vortex below and a large-scale flow above. The vortex maps to the vortical boundary driver, while the large-scale flow is related to the slow overall expansion of the magnetic structure.

In this paper, we have only addressed the case of one specific boundary perturbation, namely a twisting motion at the footpoints. We have also conducted a number of simulations with a translational motion analogous to that used in Aulanier et al. (2005), and observed comparable behavior in these cases as well. In particular, the achieved current densities continued to rise at similar rates until they were restricted by finite grid resolution.