Effect of In-Plane Shear Flow on the Magnetic Island Coalescence Instability

Magnetic reconnection (MR) is a fundamental phenomenon in dynamically evolving plasma systems that is responsible for relaxing the magnetic field topology by converting the magnetic energy into kinetic energy. This phenomenon heats up the plasma Biskamp_2000_MRbook ; priest_forbes_2000_MRbook and is frequently observed in both laboratory and space plasma. Examples include saw-tooth crash driven by tearing mode instability (TMI) in tokamaks Wesson1987 , solar flares, coronal mass ejection, magnetospheric substorms, etc. (see Zweibel2009 and references therein). The most favorable location for two-dimensional MR to occur is the region of thin current sheet that develops at the $X$-type magnetic null points ($X$-type null lines in 3D are topologically unstable) generated through MHD instabilities, e.g. the TMI Lapenta2008 , plasmoid instability Hsseinpour2018 and island coalescence instability Stainner2013 ; Stainner2017 ; rmhd1 ; rmhd2 ; rmhd3 ; Pritchett-Wo-1979 ; Pritchett1992 ; Karimabadi2011 ; finn_Kaw1977 ; Knoll_chacon2006 ; Bard2018 ; Makwana2018 , to name a few. Any change in plasma parameters around these local $X$-points is known to affect the entire process of MR. A large number of numerical simulations have been conducted to understand the role of various parameters on MR, such as the effect of non-uniform resistivity Baty2007 , presence of strong guide field Stainner2017 , presence of asymmetric magnetic field and density on the upstream side of reconnection current sheet Doss2015 ; Eastwood2013 and more. Additionally, as the plasma bulk flow is common in fusion plasma experiments (e.g., generated indirectly by neutral beam injection or wave heating/current driving mechanisms) as well as in space plasma environment (e.g., caused due to plasma jets from astronomical bodies, stellar wind, etc.), it can also affect MR in many important ways by altering the upstream and downstream flow pattern.
As is well known, the magnetic field can also affect the flow dynamics. For example, the presence of magnetic field parallel to shear flow with velocity discontinuity generates Kelvin-Helmholtz instability (KHI) if the difference in velocity across the discontinuity layer is greater than twice the Alfvénic velocity ($v_{A}$) Keppens1999 ; SChandrasekhar_KHI_book . In the past, most of the research work on MR in the presence of shear flows have used two kinds of MHD equilibrium: the Harris current sheet equilibrium Harris1962 and the Fadeev or magnetic current island equilibrium Fadeev1965 (used interchangeably throughout this work). The Harris type equilibrium has a long thin current sheet that separates two regions of anti-parallel magnetic field. In the absence of shear flow, when the current sheet becomes unstable to TMI, it breaks down into large magnetic islands through MR. However, in the presence of in-plane parallel Chen-Morrison1990 ; Ofman-Morrison1993 ; Cassak2011 ; Hsseinpour2018 or anti-parallel Ofman-Morrison1993 ; Chen1997 shear flows, the rate of island formation, equivalently the reconnection rate, reduces monotonically with an increase of shear flow amplitude up to a critical value, after which the TM mode transits to the KH mode. This critical shear flow amplitude for a resistive MHD model is $v_{A}$. The inclusion of Hall physics reduces the critical flow amplitude to sub-Alfvénic values Chacon2003 enabling the KHI and TMI to couple. However, as stronger shear flows suppress the growth and size of islands, it becomes difficult to study the TM generated magnetic islands in the presence of Alfvénic and super-Alfvénic flows. One of the drawbacks of using the Harris type equilibrium is that the initiation of TMI driven MR requires either an external driving force (e.g. Newton’s challenge problem Birn_Hesse2007 ) or the current sheet aspect ratio needs to be very large (thin and long current sheet). Moreover, in natural reconnecting systems, thin current sheets develop dynamically at the $X$-points. However, the pre-formed current sheet structure of Harris equilibrium is less suitable to address reconnection physics (see Ref. Stainner2017 for further details). Interestingly, Fadeev equilibrium Fadeev1965 has such features inherently built into it. As Fadeev equilibrium contains a 1D chain of current filaments separated by $X$-type magnetic null points, in the presence of finite dissipation, the force of attraction between the parallel current strands brings them closer leading to the formation of a thin reconnection current sheet at the $X$-point and drives the MR. This self-driven mechanism for MR is also one of the advantages of this equilibrium Stainner2017 . Mutually attracting magnetic islands finally coalesce to form a bigger island.
In the past, several attempts have been made to understand the island coalescence problem. For example, stability analysis of a 2D island configuration has been studied analytically by Finn and Kaw finn_Kaw1977 . To understand the observed fast MR time scale, 2D island coalescence problem has been studied numerically using resistive MHD rmhd1 ; rmhd2 ; rmhd3 ; Pritchett-Wo-1979 , Hall-MHD Knoll_chacon2006 ; Bard2018 , kinetic Pritchett1992 ; Karimabadi2011 , ten moment two-fluid Stainner2017 and hybrid Makwana2018 simulation models. In the above-mentioned body of work, the role of several external parameters such as system size, guide field strength on the properties of MR has been addressed. However, the effect of shear flow on island coalescence has not been studied so far. Furthermore, Fadeev equilibrium has well developed current islands to couple with both sub-Alfvénic and super-Alfvénic shear flows. Previously, several analytical studies have reported the shear flow effects on the Kelvin-Stuart finn_Kaw1977 island configuration (same as Fadeev equilibrium). Throumoulopoulos et. al. Throumoulopoulos2009 have constructed a class of magnetic island equilibrium in presence of shear flows by solving the Grad-Shafranov equation. Considering only sub-Alfvénic shear flows that are relevant to fusion experiments, they have shown stabilization of island equilibrium and formation of pressure islands. Analytical work by Waelbroeck et. al. Waelbroeck2007 , using a resistive MHD model, verifies modest influence on the stability of magnetic islands when subjected to a sub-Alfvénic shear flow. The aim of our present work is to understand the effect of shear flow on the island coalescence problem for a wide range of values of flow amplitude ($v_{0}$) from sub-Alfvenic to super-Alfvenic and for a wide range of velocity shear scale length ($a_{v}$). Using the VR-RMHD model and very high resolution computer simulations, we have investigated Fadeev equilibrium in the presence of a tan-hyperbolic flow profile (shear flow is symmetric and anti-parallel on both sides of magnetic islands). We find that in line with previous studies, the sub-Alfvénic flows have a negligible effect on coalescence instability as the time required to attain the peak value as well the magnitude of reconnection rate is close to that in the absence of shear flow. However, for super-Alfvénic shear flows and irrespective of shear length scales (compared to magnetic island width), we show the existence of coalescence instability and suppression of the magnetohydrodynamic Kelvin-Helmholtz instability (MHD-KHI). In the absence of in-plane shear flow, previously reported numerical work Birn_Hesse2007 using fully compressible resistive MHD model on island coalescence problem has shown that the effect of compressibility has no role on the reconnection rate and island dynamics. Moreover, in the case of MHD-Kelvin-Helmholtz instability (MHD-KHI), the condition for the fastest growing mode remains unchanged for both incompressible and compressible models Miura1982 . Therefore to start with the simplest, nontrivial model, here we have considered an incompressible resistive MHD model (VR-RMHD model). We also present extensive results on its quasi-linear and non-linear fate for the entire range of parameters addressed here.

The rest of our paper is organized as follows: the VR-RMHD model and the BOUT++ numerical framework used for its study are discussed in Section II. Then, a series of benchmark results for the resistive island coalescence instability in the absence of in-plane flow are discussed in Section III. In Section IV, the dynamics of magnetic island - flow shear system for different $v_{0}$ and $a_{v}$ values are discussed. Finally, our conclusions and potential future work have been discussed in Section V.

We solve the VR-RMHD model in a 2D Cartesian geometry (see Ref. rmhd1 and references there in) using the BOUT++ framework bout1 ; bout2 . As mentioned in the preceding section, this model assumes plasma as an incompressible magnetized fluid i.e. $\nabla\cdot\mathbf{u}=0$, using which the full compressible MHD equations are simplified to the vorticity ($\omega$)-vector potential ($\Psi$) formalism. The governing equations of the $\omega$-$\Psi$ formalism read as rmhd1

In the above equations, the out-of-plane or $y-$component of vorticity $\omega=\hat{y}\cdot\left(\vec{\nabla}\times\mathbf{u}\right)=-\nabla_{\perp}^{2}\varphi$; $u_{x}=-\partial_{z}\varphi$, $u_{z}=\partial_{x}\varphi$, where $\mathbf{u}~{}(=u_{x}\hat{x}+u_{z}\hat{z})$ is the in-plane (“poloidal”) velocity of the plasma and $\varphi$ is the corresponding stream function. Similarly, the out-of-plane current density $J_{y}=\hat{y}\cdot\left(\vec{\nabla}\times\mathbf{B}\right)=-\nabla_{\perp}^{2}\Psi$; $B_{x}=-\partial_{z}\Psi$, $B_{z}=\partial_{x}\Psi$, where $\mathbf{B}~{}(=B_{x}\hat{x}+B_{z}\hat{z})$ is the in-plane magnetic field and $\Psi$ is the corresponding magnetic vector potential. The above equations are normalized as follows: length $L$ to system length $L_{x}$, velocity to Alfvén velocity $v_{A}=B/\sqrt{\mu_{0}\rho}$ and time to the Alfvénic time $t_{A}=L/v_{A}$. Here, normalized density of plasma $\rho=1$ and magnetic permeability $\mu_{0}=1$. The normalized viscosity $\hat{\nu}$ and resistivity $\hat{\eta}$ are defined as $\hat{\nu}=\nu/(Lv_{A})$ and $\hat{\eta}=\eta/(\mu_{0}Lv_{A})$, where $\nu$ is the kinematic viscosity and $\eta$ is a constant resistivity. The main variables $\omega_{y}$ and $\Psi$ are normalized as $\omega_{y}L/v_{A}$ and $\Psi/(BL)$. The Poisson bracket is defined as $\left[f,g\right]_{z,x}=(\partial_{z}f)(\partial_{x}g)-(\partial_{x}f)(\partial_{z}g)$.
Equations (2) and (3) are solved in a 2D uniform Cartesian grid ($0\leq x\leq L_{x}$ and $0\leq z\leq L_{z}$) where $L_{x}=L_{z}$. For the rest of the paper, we use dimensionless, normalized quantities unless otherwise specified. Initial equilibrium current density $J_{y0}$, vorticity $\omega_{y0}$ and perturbed vector potential $\Psi_{1}$ profiles arermhd1 ,

Here, the parameter $\epsilon=0.2$ determines the island width, $a_{B}=1/2\pi$ determines the simulation domain size as $L_{x}=L_{z}=4\pi a_{B}=2$, $\Psi_{1}$ is the perturbed vector potential with amplitude $A=0.01$. For all our simulations, we have set $\hat{\eta}=\hat{\nu}=10^{-4}$. In the limit of $\epsilon=0$, Fadeev equilibrium reduces to the Harris current sheet of shear width $a_{B}$. Fig. 6(a) shows spatial profiles of $J_{y0}$ (left panel colormap) and $\Psi_{0}$ (left panel contour) and $\omega_{y0}$ (right panel colormap, for $a_{v}=2a_{B}$ case). Periodic and Dirichlet boundary conditions are respectively implemented along the $z$ and $x$ directions. In BOUT++, we have used FFT based calculations along $z$ direction and finite difference based calculations along $x$ direction. Value of $J_{y},~{}\Psi,~{}\omega_{y}$ and $\varphi$ at $x=$ 0 and 2 is zero (Dirichlet boundary). All the runs are carried out for two different grid sizes $\mathrm{d}x=4\times 10^{-3}$, $\mathrm{d}z=2\times 10^{-3}$ $(N_{x}=512,N_{z}=1024)$ and $\mathrm{d}x=2\times 10^{-3}$, $\mathrm{d}z=10^{-3}$ (corresponding $N_{x}=1024,N_{z}=2048)$ As given in Refs. rmhd2 ; rmhd3 , the width of the reconnection current sheet (say $\delta$) generated during island coalescence $\sim$ $\hat{\eta}^{2/3}$ when $\hat{\eta}\geq 10^{-4}$. In our case, $\hat{\eta}=1\times 10^{-4}$ (corresponding $\delta\simeq 0.00215)$, hence our $z$-directional grid size is sufficient enough to resolve the reconnection current sheet.
Equations (2)-(3) are solved using the BOUT++ framework which was originally developed for tokamak edge turbulence studies. The BOUT++ framework uses a finite difference technique based 3D nonlinear solver. Recent developments and the use of “Object Oriented Programming (OOP)” concepts of C++ language have enabled the users to solve an arbitrary number of coupled, general partial differential equations bout1 ; bout2 . In solving our set of VR-RMHD equations, we have used an energy conserving Arakawa bracket method for calculating the nonlinear terms and the implicit CVODE time solver from the SUNDIALS CVODE package. For all the runs presented here, the time solver takes a time step of 0.05$t_{A}$.
To check the correctness of the newly implemented model, we first test for energy conservation. In ideal-RMHD, total energy ($E_{ideal}=1/2\int\mathrm{d}\tau(|\nabla_{\perp}\varphi|^{2}+|\nabla_{\perp}\Psi|^{2}$)) of the system remains conserved. Taking $J_{y0}$ (see Eq. 4) as the initial profile and $\hat{\eta}=\hat{\nu}=0$, $E_{ideal}(t)$ as a function of time is plotted in Fig. 1. Constant $E_{ideal}(t)$ values at our presently working grid sizes (1024$\times$1024 and 2048$\times$2048) shows energy conservation by our solver in the ideal limit. Also, this verifies negligible numerical dissipation even at large grid sizes. As can be expected, the solver numerically sustains Fadeev equilibrium. Here, for our 2D model, the integration is over $\mathrm{d}\tau=\mathrm{d}x\mathrm{d}z$ and the operator $\nabla_{\perp}=\hat{x}\partial_{x}+\hat{z}\partial_{z}$.

One of the important parameters to investigate is the reconnection rate. In literature rmhd1 , this reconnection rate ($M$) is calculated from the reconnection flux $\psi_{r}$ as $M=\partial_{t}\psi_{r}$ ($\psi_{r}=\Psi_{X}-\Psi_{O}$, with $\Psi_{X}$ and $\Psi_{O}$ the magnetic flux at $X$-point and O-point, respectively). In the absence of in-plane shear flows, the line joining the O-points always aligns itself with the line $x=1$ (see for instance Fig. 2), hence it becomes easy to calculate $M$ from the time derivative of $\psi_{r}$. However, in the presence of shear flows, the line joining the O-points gets twisted because of flow dynamics and it becomes numerically expensive to compute the reconnection rate from reconnection flux. Reconnection rate can also be calculated by measuring the reconnection electric field $E_{y}$ at the $X$-point using Eq. 3, as $E_{y}=-\eta J_{y}|_{\scriptscriptstyle X}$ rmhd1 ; Stainner2013 . In the absence of shear flows, the nonlinear term (first term in right side of Eq. 3) has negligible contribution for the calculation of $E_{y}$ inside the reconnecting current sheet rmhd1 . In the presence of shear flows, this is also found to be valid. Hence, in this work, the reconnection rate is calculated as $M=E_{y}=-\eta J_{y}|_{\scriptscriptstyle X}$ for our benchmark studies in the absence of shear flows as well as studies with shear flows.

For the benchmark purpose, we have considered the resistive island coalescence problem rmhd1 ; rmhd2 ; rmhd3 without any shear flow. Benchmark study uses the initial profile of equilibrium $J_{y0}$ and perturbed $\Psi_{1}$ (as given by Eq. 4) with $v_{0}=0$. We have used the grid size $\mathrm{d}z=5\times 10^{-4}$, $\mathrm{d}x=1\times 10^{-3}$. Six different resistivity values have been considered between $5\times 10^{-6}$ and $2\times 10^{-4}$ (see Fig. 3(b)). Here our length scale $L=2$, whereas it is $L=1$ in Ref. rmhd1 because of their quarter domain simulation; our time scale is therefore twice of that in Ref. rmhd1 .

In the absence of in-plane shear flow, the time evolution of a Fadeev equilibrium is shown in Fig. 2, similar to Fig. 2 of Ref. rmhd1 . As discussed earlier, Fadeev equilibrium is an ideal MHD equilibrium, hence when perturbed with finite resistivity, the plasma is able to break the frozen-in condition and cross the $X$-point finn_Kaw1977 . This allows the attraction force between the current filaments to become dominant. The typical perturbation profile $\Psi_{1}$, having a maximum of magnetic flux at the $X$-point, further accelerates the instability. This movement of islands towards the $X$-point is shown in Fig. 2(a). Fig. 2(b) shows that the coalescence process has started and the reconnection sheet has formed at the $X$-point. In Fig 2(c), the current sheet is fully developed and the reconnection rate is at its maximum. After this, the magnetic flux piles up on both sides of the current sheet causing a slow down of the coalescence process, and hence the reconnection rate. The reconnection rate at these time frames are marked in Fig. 3(a). In Fig. 1 of Ref rmhd1 , the reconnection rate for $\hat{\eta}=2\times 10^{-5}$ attains a peak value at $t\simeq 7t_{A}$. Likewise in our case, as seen in Fig. 3(a), $E_{y}$ peaks around $t\simeq 3.3t_{A}$ (recall that our $t_{A}$ is twice that defined in Ref. rmhd1 ).

We now turn on the shear flow through an initial vorticity profile as given in Eq. 4. The shear flow profile and simulation parameters are chosen such that in the absence of current filaments ($\epsilon\rightarrow 0$), the KHI gets destabilized. For our presently chosen domain size ($L_{z}$) and parameters ($v_{0}=1.4v_{A}$ and different values of $a_{v}$), we use the analytical formula given by Miura Miura1982 to calculate the fastest growing MHD-KHI mode (of mode number $m$) given as, $m=L_{z}/2\pi a_{v}$. In Table 1, we have listed the values of $m$ for two different $L_{z}$ value.

To verify the stability of these modes indicated, we perturb the velocity shear profile with different mode numbers. Left panel of Fig. 4 shows the variation of growth rate for different modes for various $a_{v}$ values and domain sizes. As expected, we do not get an unstable mode for $a_{v}=2a_{B}$ for domain size $L_{z}=2$. To get an unstable mode for this velocity shear width, we double both $L_{x}$ and $L_{z}$ i.e. from $L_{z}=L_{x}=2$ to $L_{z}=L_{x}=4$ as well as grid numbers $N_{z}$ and $N_{x}$ in order to keep the grid resolution same. As shown in the right panel of Fig. 4, doubling the length $L_{x}$ (for $L_{z}=4$ case) takes the x-boundary farther from the vorticity sheet, resulting a marginal change in growth rate of the MHD-KHI modes. The fastest growing mode for $a_{v}=2a_{B}$ and $L_{z}=4$ is found to be $m=1$.
Effects of in-plane shear flow in the island coalescence problem are studied by changing three important parameters in the vorticity (or velocity) profile: (1) flow shear strength $v_{0}$ (2) shear width $a_{v}$ and (3) the direction of shear flow parallel or anti-parallel to the magnetic field. Parameters used for these simulation are given in Section II.

The simulations are performed on the Antya cluster at the Institute for Plasma Research (IPR). We would like to thank Dr N. Bisai, IPR for his valuable inputs.