© ESO 2019

1. Introduction

Solar flares release up to ∼10²⁵ J of magnetic energy over the course of minutes to hours (Benz 2017; Fletcher et al. 2011). Of this energy, up to 50% can be transported away by ∼10³⁹ energetic charged particles (Aschwanden et al. 2016; Vilmer 2012) and the mechanism responsible is still an area of contention. The primary means by which magnetic energy stored in the corona is converted into kinetic energy is widely accepted to be magnetic reconnection, due to the broad range of observed phenomena that it accounts for (Janvier et al. 2015; Shibata & Magara 2011; Priest 2014).

The reconnection of magnetic field lines comes in many flavours, each of which is relevant in some particular coronal magnetic field topology. In the case of solar flares, a significant factor in the acceleration of particles is the direct electric field at the site of the reconnection (Zharkova et al. 2011) – in the standard flare model this is described as being part of a “monolithic” current sheet above the flaring loop, which is borne out by observations of a hard X-ray source at the loop-top (Su et al. 2013; Vilmer 2012).

There are numerous outstanding problems with the standard model interpretation, such as the minimal likelihood of a monolithic current sheet forming in the coronal environment and then persisting for long enough to accelerate large quantities of particles (Aschwanden 2006). The corona is a dynamic place: any equilibria which arise can be perturbed by any number of phenomena, photospheric footpoint motion (Parker 1972) or magnetohydrodynamic (MHD) waves for example. It is therefore expected that current sheets will form and will be perturbed (Longcope & Cowley 1996) – hence the relevance of forced magnetic reconnection which consideres reconnection as an externally driven and transient process.

Bobrova & Syrovatskii (1979) found that by deforming the boundary of some region of plasma in a sheared force-free field, a current sheet may be formed at the a resonant surface (where k ⋅ B = 0). Due to the non-zero resistivity of hot coronal plasma, the “frozen-in” conditition of ideal MHD may be violated inside current sheets; the topology of the magnetic field is no longer necessarily conserved, and the field may relax, releasing stored magnetic energy. Hahm & Kulsrud (1985) considered Taylor’s problem: a Harris neutral sheet equilibrium, stable to the tearing mode, to which a brief boundary perturbation of small amplitude was applied.

During forced magnetic reconnection, the equilibrium current sheet tears into an infinite chain of magnetic islands through tearing mode-like evolution, accompanied by the rapid release of stored magnetic energy. The analysis of Hahm & Kulsrud (1985) described forced reconnection as a process involving a linear resistive phase followed by a transition into a non-linear Rutherford regime (Rutherford 1973). This transition was subsequently revealed to be more complex depending on the conditions – either passing into a non-linear Sweet–Parker regime (Waelbroeck 1989) or bypassing the Rutherford phase entirely (Wang & Bhattacharjee 1992). Numerical results of the early 2000s closely followed those of the analytical work, with the visco-resistive MHD simulations of Fitzpatrick (2003) showing that the findings of Hahm & Kulsrud (1985) were accurate and that at large perturbation amplitudes the non-linear behaviour described by Wang & Bhattacharjee (1992) could be reproduced.

Forced magnetic reconnection was later studied in force-free fields, which are more relevant to the corona, by Vekstein & Jain (1998), finding that the released magnetic energy can be much greater than that which was supplied by the boundary perturbation – a stable, coronal plasma equilibrium could be induced to relax to a lower energy state by applying small amplitude, spatially periodic boundary perturbations. Thus, the driving perturbation is a trigger for energy release, but not the main energy source. The non-linear extension of this linear, analytical theory was explored using numerical simulations by Browning et al. (2001).

There is some commonality between the plasmoid instability and the chain of islands produced by forced reconnection. In the latter, however, the size of the islands is determined by the external perturbation. Thin current sheets are unstable to the plasmoid instability in sufficiently large systems when the Lunquist number S_L is greater than some critical value. The current sheet spontaneously tears into chains of islands across multiple length scales in a process of cascading reconnection (Bárta et al. 2011). These systems achieve a reconnection rate of around one order of magnitude greater than that of Sweet–Parker reconnection (Bhattacharjee et al. 2009). By virtue of providing a mechanism for fast magnetic reconnection in the corona where the Lundquist number is large (S_L ∼ 10¹²), the plasmoid instability has attracted a great deal of attention and indeed has been observed many times in solar flares and magnetotail reconnection (Loureiro et al. 2007).

Recent studies have shown that the evolution of systems undergoing forced reconnection may take a different path due to the formation of a plasmoid unstable non-linear current sheet. The first numerical evidence of this was obtained in the work of Comisso et al. (2015), where it was shown that plasmoids may form in the current sheet thereby triggering a fast-reconnection regime. Subsequent analysis by Vekstein & Kusano (2015) supports this and frames these results as an interruption of the previously accepted non-linear evolution prior to the Rutherford phase due to the characteristic timescale of the plasmoid instability being much shorter.

Whatever the formation process, each island in the chain formed by the initial reconnection event is attracted to its neighbour due to the parallel currents at their O-points, which can result in islands coalescing in a further reconnection event (Schumacher & Kliem 1997). This, coupled with the self-similar structure of the thin current sheet connecting coalescing islands, is the mechanism by which the plasmoid instability accelerates reconnection to a fast regime (Bhattacharjee et al. 2009; Huang & Bhattacharjee 2010) and could contribute significantly to the energy release in the forced reconnection scenario.

Most work on the plasmoid instability has focussed on 2D simulations but the validity of some of these results has recently been called into question, since there is evidence that plasmoid-unstable configurations in 2D fail to produce plasmoids in 3D simulations (Kulpa-Dybeł et al. 2010; Kowal et al. 2017; Takamoto 2018). In these cases it seems that the absence of plasmoids is caused by turbulence, which has different properties in two or three dimensions since in 3D the extra degree of freedom allows field lines to slip past each other. This can result in a broader reconnection region, with a reconnection rate independent of any small scale diffusion process such as resistance or viscosity i.e. fast reconnection, consistent with the model proposed by Lazarian & Vishniac (1999). This process has been shown to occur in the absence of any driving, due to turbulent reconnection outflows seeding the formation of turbulent reconnection regions (Oishi et al. 2015; Beresnyak 2016; Kowal et al. 2017). In light of this it has been proposed that failure to reproduce the plasmoid instability in 3D MHD simulations could be caused by the onset of the instability driving a turbulent mode which becomes dominant. These works are in tension with the findings of Huang & Bhattacharjee (2016), where 3D MHD simulations of merging flux ropes revealed a plasmoid-mediated turbulent reconnection with a rate one order-of-magnitude lower than that reported by Kowal et al. (2009) in their study of externally-driven turbulent reconnection. Furthermore, recent analytical work by Comisso et al. (2018) has described a theory of MHD turbulence in the presence of the plasmoid instability, showing that the magnetic Reynolds number above which plasmoid formation is sufficiently statistically-significant to alter the turbulent cascade is quite small: implying that plasmoids can be expected to modify the turbulent behaviour of many astrophysical systems.

Regardless of how this tension between the turbulent and plasmoid models of reconnection is resolved, forced reconnection is a model of a fundamentally different physical situation. Rather than having a large, unstable current sheet or applying small-scale driving, the forced reconnection scenario considers a tearing-stable current sheet subjected to large-scale, coherent, transient driving. The size and number of plasmoids produced in the initial reconnection depends only on the nature of the applied perturbation rather than the Lundquist number, as in the case of tearing instabilities. We build here on the previous body of work on forced reconnection which considers only 2D geometries, and the focus of our study here is to extend this work to consider more realistic perturbations and coalescence of plasmoids. Therefore, in order to be able to compare with previous work and determine the significance of these novel effects, we restrict to 2.5D simulations – having field components in all 3 coordinate directions but with an invariant direction (here, z). There is no doubt that the extra physics introduced with the addition of a 3rd dimension will introduce novel and interesting phenomena, and may significantly alter the dynamics, but we leave these questions for future work. For example, Ripperda et al. (2017a) investigated reconnection and particle acceleration in a pair of twisted flux ropes with anti-parallel currents using a 2.5 model, and subsequently considered the effects of fully 3D geometry Ripperda et al. (2017b). Much of the behaviour was common between the two models, but instabilities such as kink and tilt play a role in the 3D case, and, with respect to test particles, they showed that the 2.5D model unrealistically allows large accelerations due to the effectively infinite acceleration length in the invariant direction (although this effect can be corrected for).

The MHD simulations of forced magnetic reconnection developed by Browning et al. (2001), Jain et al. (2005) and Gordovskyy et al. (2010a) ended at the initial reconnected equilibrium – a chain of magnetic islands. This was due to the limitations of a small domain – as only one island in an infinite chain was simulated, there was no possibility of the island merging with another in the chain. Zhang et al. (2014) employed a similar method of enquiry when investigating island coalescence in a perturbed Harris-sheet configuration (see GEM Challenge, Birn et al. 2001). By using a larger simulation domain, the formation of plasmoids and their merging was observed. Resistivity in this study was uniform as opposed to the anomalous resistivity method employed by Gordovskyy et al. (2010a), which will have caused an overestimation in the resistive diffusion observed.

We investigate the dynamics of forced magnetic reconnection with respect to the form of the perturbation used, with the intention to better understand the phenomenon of secondary reconnection in the form of island coalesence. A number of questions arise which we consider in this paper. Our main focus is on how the process of forced reconnection depends on the spatial distribution of the driving perturbation, and particularly how this affects the release of magnetic energy. The analytical theories of forced reconnection – and most simulations based on this theory – assume a spatially-sinusoidal perturbation, which extends over an infinite length of the boundary. This is clearly not realistic, and driving in the solar atmosphere (or any other context) will clearly have a range of spatial forms, and will be to some extent localised in space. Are the main principles of forced reconnection valid for more general forms of spatial driving, included localised perturbations? And how do the reconnection dynamics and the energy release depend on the form of the driving disturbance? We further consider how allowing for island coalescence affects forced reconnection, and in particular, how the reconnection during coalescence differs from the initial reconnection during island formation. A particular focus, relevant to flares and coronal heating, is to consider the energy release through the different stages of reconnection. This has important implications for the acceleration of non-thermal particles, which will be considered in a subsequent paper: we will use the computed MHD fields as input to a test particle code with a view to investigating particle acceleration in each stage of reconnection. Thus, it is important to determine the structures of parallel electric fields.

To this end we extend the model of forced reconnection developed by Gordovskyy et al. (2010a,b). In the next section, we describe the MHD model and parameters. Then follows results in Sect. 3, with a discussion of the different types of perturbation used and their effects on the evolution of the system. After the summary is an appendix with an abridged linear analysis of forced reconnection for the Harris current sheet used in the numerical work.

2. MHD simulation set-up

In our model, the MHD code Lare2d (Arber et al. 2001) is used to investigate the 2.5D spatial and temporal evolution of a transient (forced) magnetic reconnection event, triggered by perturbing a stable current layer in a force-free field. In previous works on forced reconnection, the perturbation applied is a sinusoidal flow pattern at the upper and lower boundaries. Here, we consider instead multiple sinusoidal perturbations, varying the amplitude and wavelength of both sinusoids, as well as localised (narrow Gaussian) perturbations. In this sense, we build up from the analytical work, which requires a simple sinusoid across the whole boundary, to more localised perturbations.

To simulate forced reconnection events, a 2D equilibrium which is initially force-free is perturbed by inflows and outflows along its upper and lower edges (Fig. 1) and its evolution is calculated by solving the resistive MHD equations (Priest 2014).

Fig. 1. Left panel: diagram representing the initial conditions. The domain is periodic in the x-direction, with rigid boundaries in the y-direction after the perturbation is applied. The guide field Bz (out of the plane) is shown in blue, overlaid with field lines showing the x-component of the magnetic field Bx. The red dotted line shows the form of the perturbation applied at the boundary: a sinusoidal inflow and outflow pattern with red arrows indicating direction, equal and opposite at the upper and lower y boundaries for a short period of time (see Eq. (9)). Right panel: the initial force free magnetic field components.

The initial equilibrium field and perturbation are based on the forced reconnection model by Gordovskyy et al. (2010a,b), which in turn is based on the analytical force-free forced reconnection model by Vekstein & Jain (1998). The equilibrium is a force-free Harris layer; comprised of a field reversal in the x-direction about y = 0 with an out of plane guide field (Fig. 1):

$$\begin{array}{c}\hfill {\mathbf{B}}_{\mathrm{ini}}=B_0[\text{tan}h\frac{y}{y_0};0;\text{sec}h\frac{y}{y_0}],\end{array}$$(1)

where y₀ is the initial characteristic width of the current layer centred on y = 0. This field is stable to the tearing mode for sufficiently large values of y₀.

The linear forced reconnection solution for this particular field, valid for small amplitude driving perturbations, is presented in Appendix A.

2.1. Lare2d

LareXd is a Lagrangian remap code used to solve the MHD equations in 2 or 3D (Arber et al. 2001) (the 2D version is used here). The timestep is split into two substeps: the Lagrangian step (where the grid can be distorted) and the remap step where the solution is conservatively remapped back onto the original grid. The original grid is staggered: in 2D, the magnetic field is calculated on cell edges and is updated such that the zero-divergence condition ∇ ⋅ B = 0 is satisfied to machine precision. Velocity is defined on cell vertices and updates such that mass is conserved.

We solve the Lagrangian form of the MHD equations given by

$$\begin{array}{c}\hfill \frac{D\rho }{\mathit{Dt}}=-\rho \mathbf{\nabla }·\mathbf{v}\end{array}$$(2)

$$\begin{array}{c}\hfill \frac{D\mathbf{v}}{\mathit{Dt}}=\frac{1}{\rho }[(\mathbf{\nabla }\times \mathbf{B})\times \mathbf{B}-\mathbf{\nabla }p]\end{array}$$(3)

$$\begin{array}{c}\hfill \frac{D\mathbf{B}}{\mathit{Dt}}=(\mathbf{B}·\mathbf{\nabla })\mathbf{v}-\mathbf{B}(\mathbf{\nabla }·\mathbf{v})-\mathbf{\nabla }\times (\eta \mathbf{\nabla }\times \mathbf{B})\end{array}$$(4)

$$\begin{array}{c}\hfill \frac{Dϵ}{\mathit{Dt}}=\frac{1}{\rho }(p\mathbf{\nabla }·\mathbf{v}+\eta {\mathbf{j}}^2),\end{array}$$(5)

$$\begin{array}{c}\hfill \mathrm{\nabla }\times \mathbf{B}={\mu }_0\mathbf{j}\end{array}$$(6)

where B and j are the magnetic field and current density, respectively. ρ, v, p and η are the plasma mass density, flow velocity, pressure and resistivity. Note that D/Dt is the advective derivative of a quantity in a velocity field. The specific internal energy density is defined as

$$\begin{array}{c}\hfill ϵ=\frac{P}{\rho (\gamma -1)}=\frac{k_{\text{B}}T}{{\mu }_{\mathrm{m}}(\gamma -1)},\end{array}$$(7)

where μ_m is the reduced mass, calculated by averaging the mass of all particles in the plasma.

Anomalous resistivity was employed, whereby the resistivity was enhanced in regions where the current was sufficiently strong:

$$\begin{array}{c}\hfill \eta =\{\begin{array}{cc}{\eta }_{\mathrm{bg}}\hfill & \text{if}j<j_{\text{crit}}\hfill \\ {\eta }_0\hfill & \text{if}j\ge j_{\text{crit}},\hfill \end{array}\end{array}$$(8)

where η_bg and η₀ are the background and anomalous resistivities, respectively, and j_crit is the critical current.

While anomalous resistivity is a commonly used computational convenience, there are some physical arguments for its inclusion which are well documented in the work of Wilkins (1980). In short, it can be treated as a numerical method to approximate particle effects in an MHD code – that is, there are current-dependent micro-instabilities that operate on unresolved kinetic scales which have the effect of enhancing resistive diffusion in regions of high current (Büchner & Elkina 2006; Wu et al. 2010; Bárta et al. 2010; Bareford & Hood 2015). Since the scale on which these effects operate is not captured in fluid models, we cannot know the “correct” values of the parameters. However, since resistivity sets the timescale on which the reconnection unfolds, we can fit the parameters to the timescales observed in the corona (∼10 s–100 s) (Che 2017). As a practical convenience, it sets a lower limit on the thickness of any current sheets which form, hence the threshold current j_crit is chosen to be higher than that which would occur in a grid-scale current sheet. This allows dissipation to occur mainly through reconnection, rather than global resistive diffusion.

2.2. Perturbations

To instigate forced reconnection, the y = ±L₀ boundary must be driven. This is achieved by applying plasma flows across the boundary, which drag the field with them via the frozen-in condition, for a short period of time (O(10τ_A)).

In classic forced magnetic reconnection, this would be sinusoidal in x. Here we investigate different forms of perturbation with a view to understanding how the analytical theory extends to more realistic scenarios. Note that non-linear effects can be important even for small perturbation amplitudes (Δ), since the island width scales are ${\mathrm{\Delta }}^{\frac{1}{2}}$ (Browning et al. 2001). Since the magnetic field evolution is non-linear, it is non-trivial to determine how the effects of different perturbations may combine. We begin with sinusoial perturbations, then pairs of sinusoids of different wavelengths and amplitudes, then localised perturbations which have a Gaussian profile. Equations (9)–(11) show the y-component of the velocity used in each case. In all cases, υ_x = 0, υ_z = 0 and υ_y is set to zero when t >  τ_p.

The sinusoidal case:

$$\begin{array}{c}\hfill {\upsilon }_y=\pm \frac{\mathrm{\Delta }}{{\tau }_{\mathrm{p}}}[1-\text{cos}\left(\frac{2\pi t}{{\tau }_{\mathrm{p}}}\right)]\text{cos}\left(\frac{2\pi x}{\lambda }\right),\end{array}$$(9)

where Δ is the perturbation amplitude and λ is the wavelength.

For the multiple sinusoids:

$$\begin{array}{c}\hfill {\upsilon }_y=\pm \frac{1}{{\tau }_{\mathrm{p}}}[1-\text{cos}\left(\frac{2\pi t}{{\tau }_{\mathrm{p}}}\right)][{\mathrm{\Delta }}_1\text{cos}\left(\frac{2\pi x}{{\lambda }_1}\right)+{\mathrm{\Delta }}_2\text{cos}\left(\frac{2\pi x}{{\lambda }_2}\right)]·\end{array}$$(10)

For the localised perturbation:

$$\begin{array}{c}\hfill {\upsilon }_y=\frac{\mathrm{\Delta }}{{\tau }_{\mathrm{p}}}[1-\text{cos}\left(\frac{2\pi t}{{\tau }_{\mathrm{p}}}\right)]\text{exp}(-\frac{x^2}{l_x^2}),\end{array}$$(11)

where l_x is the characteristic width of the Gaussian.

2.3. Simulations

Table 1 contains the parameters used for each simulation to be presented. For each multiple sinusoidal perturbation (MSP), IDs beginning with 2, each component sinusoid is also simulated as a lone perturbation. In this way, we build from single sinusoidal perturbations to multiple interacting perturbations to a general perturbations (which can be expressed as an infinite Fourier summation). The summation of the effects is highly non-linear, with strong dependence on the wavelength and amplitude of each perturbation.

Parameters common to all simulations are the (normalised) current sheet thickness y₀ = 0.25 and the anomalous resistivity terms defined in Eq. (8): threshold current j_crit = 6.0, background resistivity η_bg = 10⁻⁷ and the enhanced resistivity η₀ = 10⁻³. Changing y₀ would affect the island size, for a diven driving disturbance, but this variation has been investigated previously (Browning et al. 2001).

The code parameters are non-dimensionalised with respect to the peak magnetic field strength (B₀), box half-width (L₀) and initial (uniform) density (ρ₀). Thus velocities are scaled with respect to the Alfven speed ${\upsilon }_{\mathrm{A}}=B_0/\sqrt{(}m{u}_0\rho )$ and times with respect to the Alfven time τ_A = L₀/υ_A.

3. Results and discussion

First, “white noise” perturbations were investigated, in which the amplitude of the velocity at each location along the boundaries parallel to the x-axis were randomly assigned, but resulted in no reconnection, confirming tearing mode stability; and the importance of the form of the driving perturbation as a trigger.

3.1. Sinusoidal perturbations

Case 1a is a reproduction of the work of Gordovskyy et al. (2010b) (Fig. 2), while case 1b (Fig. 4) is the same model with a simulation domain that is 4 times wider in the x-direction. Note that in case 1a, and in the early stages of case 1b, a chain of magnetic islands is formed, as expected in the forced reconnection scenario (Hahm & Kulsrud 1985; Vekstein & Jain 1998). The periodicity matches the driving perturbation. Note that, because we have a guide field (B_z), these islands are actually 3D twisted magnetic flux ropes (of infinite length). Since the domain is periodic in the x-direction, the two cases might be expected to be identical in physical terms. However, island coalescence requires symmetry to be broken in a way that is not possible if only one island is captured by the domain. Figures 3 and 4 show the energy profiles over time and the MHD snapshots, respectively, for case 1b. The primary difference in the outcome is that rather than completing the initial reconnection and remaining in a new saturated equilibrium, the magnetic islands attract each other and coalesce. Strong electric fields develop in current sheets oriented in the y-direction where the islands are merging. At t = 350τ_A, the system had been in steady state for 100τ_A with no further merging taking place.

Fig. 2. B and Bz for case 1a (single wavelength λ = 4.0) at final snapshot t = 100τA; a replication of the final state of the model by Gordovskyy et al. (2010b).

Fig. 3. Energetics for case 1b (Table 1), with corresponding MHD fields shows in Fig. 4. At around 100tA, the first stage of reconnection is slowing to a halt as islands reach saturation. At around 150tA the large drop in magnetic energy (panel a) and spike in kinetic energy (panel b) correspond to the island coalesence event seen in Fig. 4 (panels e–f). The total energy is shown as an indication of numerical accuracy over the course of the simulation: the small initial rise is due to the driving at the boundary, remaining constant thereafter (observed in all simulations, witihin 5%).

Fig. 4. MHD snapshots for case 1b (Table 1) – a single sinusoidal perturbation in a wide box. Left to right panel: magnetic field lines B with Jz contours, magnetic field lines with Bz contours and the plasma flow v. Arrows on v plots indicate the in-plane velocity direction. Snapshots (panel a) through (panel f) correspond to times t = 0, 8, 64, 112, 144 and 300τA. Note that the aspect ratio is not to scale – the width of the domain 16 in x and 2 in y. The colour scheme is adaptive to enhance contrast. Note that colour representing zero is consistent across similar plots throughout.

Figure 3 shows that this second reconnection event results in a rapid decrease in magnetic energy – where the primary reconnection takes around 100τ_A, the coalesence is an order of magnitude faster. It is also accompanied by a much greater increase in kinetic energy. The reconnection associated with coalescence is thus more violent than the initial island-forming event. This has implications for particle acceleration which will be explored in a Paper II of this series – a multi-stage acceleration process can be expected, with further acceleration of particles taking place during magnetic island merger events.

The first, island-forming phase of reconnection releases about 8% of the initial magnetic energy while the second, coalescence phase releases a further 4%. Although there is significant kinetic energy during the reconnection periods, the final state is close to a new magnetostatic equilibrium and almost all the dissipated magnetic energy is converted into internal energy.

The flow patterns can also be seen in Fig. 4: in b, inflows in the y-direction are indicating where the x-points will form and outflows indicate where the o-points will form. By snapshot c, the vortices responsible for growing and shaping the island can be seen. By snapshot f, we see a strong increase in the y-component of velocity as the islands attract each other and begin reconnecting.

3.2. Multiple wavelength perturbations

For cases 1c, 1d and 2a, combinations of perturbations of wavelength λ₁ = 2 and λ₂ = 16(=L_X) were investigated. As expected from linear theory, the longer wavelength perturbation (case 1d) releases a larger fraction of the initial magnetic energy (compared with cases 1a and 1b): thus, about 16% of the total magneic is eventually dissipated. Figure 5 shows the final state of the λ₁ = 2 lone perturbation simulation (case 1c) – little reconnection occurs, with small islands forming but failing to coalesce. The magnetic field loses energy monotonically, which drives some Ohmic heating. As predicted by linear analytical theory (see Vekstein & Jain 1998 and Appendix A), the island width (for a given driving perturbation) is smaller for shorter wavelength perturbations.

Fig. 5. B and Jz for case 1c (single wavelength λ = 2.0) at final snapshot t = 350τA.

Figure 6 shows key moments in the magnetic field evolution as well as the magnetic and kinetic energy profiles for the λ = 16 lone perturbation (case 1d). Since the perturbation is of wavelength equal to the width of the box, initially only one island begins to form (Fig. 6a) – apparently similar to the shorter wavelength case shown in Fig. 2. As the island grows, a thin current sheet connecting the halves of the island begins to form and we observe the onset of the plasmoid instability as it tears into a shorter-length magnetic island. This is associated with quite prolonged and strong increases in kinetic energy. The formation of this large island in the centre – despite the fact that the boundary disturbance is creating inflows at this location – may be interpreted in terms of the linear theory (see Appendix A), since the related equilibrium indeed has an island in the box centre (at the location of strongest boundary inflow). As the smaller island grows it has the effect of increasing the rate at which magnetic energy is depleted (Fig. 6). This leads to energy profiles that follow a similar pattern to the reconnection into coalesence seen in Fig. 3. Eventually, the smaller island is absorbed into the larger island and the kinetic energy exhibits a similar spike to that of case 1b.

Fig. 6. Panels a–d: B and Jz for case 1d (single wavelength λ = 16.0) at time t = 40, 92, 196 and 212τA. Panels e–f: corresponding energetics: magnetic energy (left panel) and kinetic energy (right panel). The snapshot times shown correspond to the different stages of reconnection: initial perturbation of equilibrium, first saturated state, island attraction and island coalescence.

The results of applying both perturbations combined differ from both previous cases, as shown in Fig. 7. At t = 24τ_A, the field topology is similar to that found at the same time in the case of the lone pertubation of λ = 16 (case 1d, Fig. 6), except with short-wavelength curvature patterns superimposed by the λ − 2.0 perturbation – this is around the end of the linear stage of the growth. By time t = 40τ_A (snapshot b) however, the difference between cases 1d and 2a become more pronounced in that no central island has formed. Instead the periodic inflow from the small wavelength perturbation supresses the growth of any secondary islands by pushing them out towards the large island. This results in a series of smaller plasmoid absorption events. These plasmoids coalesce with the large island at a much earlier time, so the simulation reaches a similar end-state (to the lone long-wavelength perturbation case) around 100τ_A sooner. The difference in the distribution of islands and reconnection events also means a different distribution in the electric field, which will have implications for particle acceleration.

Fig. 7. Panels a–d: B and Jz for case 2a (multiple wavelengths λ1 = 2.0, λ2 = 16.0) at time t = 24, 40, 92 and 350τA. Panels e–f: corresponding energetics: magnetic energy (left panel) and kinetic energy (right panel). Case 2a uses the perturbations from cases 1c and 1d in tandem, resulting in different reconnection dynamics and energetics, such as the lack of formation of the central island as seen in Fig. 6.

It should also be noted that the evolution of the magnetic energy is rather different from the cases previously discussed. The eventual loss of magnetic energy is almost the same as in case 1d (the long wavelength lone perturbation) – but as well as reaching the final equilibrium much quicker, Fig. 7e shows that the island-forming and island-merging stages of reconnection are not now clearly distinguishable. Similarly, there is a bursty series of releases of kinetic energy, rather than the two peaks seen previously. The reconnection here is even more dynamic, with levels of kinetic energy reaching as high as 0.3 around t = 40 (in dimensionless units), accounting at that stage in the evolution for the bulk of the released magnetic energy. Such a large kinetic energy results from flows close to the Alfven speed over quite large fractions of the domain.

Cases 1c to 2a demonstrate that longer wavelength perturbations dominate the evolution of the system in terms of location and size of final islands, and while shorter wavelength perturbations do not greatly effect the system alone, they can significantly modify the behaviour of systems subjected to both large and small wavelength perturbations. Crucially, the result is not merely a superposition of the two resulting sets of islands, but rather a significant departure from the dynamics of either individual case.

Case 2b (Fig. 8) is the same combination of perturbing wavelengths as case 2a but the amplitude of the small wavelength component has been increased so that its contribution to the overall dynamics would be comparable to that of the long wavelength component. The reconnection onset occurs sooner and proceeds more rapidly than in case 2a: for 2a, the final reconnected state is reached at around 150τ_A; for 2b, this happens around 30τ_A. This is a result of the larger overall perturbation, which makes the two cases difficult to compare in these terms. However, it is clear from Fig. 8 that in the early stages of reconnection the small wavelength perturbation dominates: plots a and b show the formation and saturation of a chain of islands of width matching the wavelength of 2. The stronger the amplitude, the stronger the effect. The long wavelength flows cause early coalescence, leading to wider islands, fewer in number.

Fig. 8. Panels a–d: B and Bz for case 2b (multiple wavelengths λ1 = 2.0, λ2 = 16.0, with increased amplitude on the smaller wavelength component such that δ1 = 0.35) at time t = 8, 20, 28 and 350τA. The simulation reaches its final state much more rapidly than in case 2a due to the increased total flux. Bz was chosen to be displayed rather than Jz since, in this case, the island structure was more visible.

3.3. Localised perturbations

An infinitely-long periodic disturbance clearly cannot occur in reality, and we now condifer the effects of applying a spatially localised driving disturbance.

In order to investigate localised perturbations, a narrow Gaussian disturbance was applied. Figure 9 shows an example of a perturbation applied from the y = +L₀ boundary at position x = −2. Interestingly, similar magnetic field evolution is exhibited to that which was found in the large wavelength lone perturbation (Fig. 6, case 1d), which is also borne out in the energy profiles. Similar features of the magnetic field evolution include the formation of a current sheet between the islands in the infinite chain which tears into a smaller island. This smaller island is absorbed, and small plasmoids are seen forming in the regions between islands in both cases 4a and 2a. The magnetic and kinetic energy profiles are remarkably similar to those in Fig. 7, characterised by decaying periodic spikes in kinetic energy accompanied by a rapid decrease in stored magnetic energy. This is due to the multiple island growth and coalesence events common to both this localised case the multiple wavelength case.

Fig. 9. B and Jz for case 4a (localised perturbation) at time t = 24, 40, 60, 84, 112 and 152τa.

Since, as discussed above, larger wavelengths contribute more to the evolution of the system, and a Gaussian curve can be expressed as an infinite series of weighted sinusoids, it is likely that the longer wavelength components of the Gaussian are dominating the shorter wavelengths.

4. Summary and conclusions

The solar corona is a dynamic array of magnetic structures, interacting in a myriad of ways and storing free magnetic energy which may be released as a large-scale solar flare or a smaller nanoflare associated with coronal heating. The forced reconnection scenario shows us that even stable magnetic fields may be perturbed by a transient external disturbance such that a current sheet forms and undergoes magnetic reconnection, releasing magnetic energy stored in the field, which was initially in a state of force-free equilibrium. It should be emphasised that the driving disturbance is a trigger for releasing stored energy, and the energy input/dissipation associated directly with the disturbance itself are usually quite insignificant. Whilst the traditional models of forced reconnection assume an idealised periodic sinusoidal perturbation, we consider the effects of more complex and realistic driving disturbances. An important outcome of our study is to demonstrate that forced reconnection still “works” for such realistic drivers, and the main features of forced reconnection are not just artefacts of the idealised model. Nevertheless, the relationship between the dynamics and energetics of the reconnection and the form of the driving disturbance are quite subtle, and also depend critically on the effects of the coalescence instability. By considering longer numerical boxes than in previous studies, we found that the chains of magnetic islands formed in forced reconnection are inevitably subject to coalescence instability, with the merger of islands leading to further release of stored magnetic energy. However, the distribution of island sizes, both from the primary driving and from the secondary plasmoid instability, depends strongly on the spatial form of the driving disturbance.

Island coalescence events were found to be more violent than their primary reconnection counterparts (which occurred in the initial current sheet creating the chain of magnetic islands). The reconnection is much faster, and strong bursts of kinetic energy are produced, associated with both the rapid motion of islands towards each other, and the subsequent fast reconnection outflows. Once the thin current sheet forms between two islands that are merging, they rapidly merge in powerful reconnection events. The main reason that these island mergers are much faster than the initial reconnection is that coalescence is an ideal instability, driven by increasing imbalances of the Lorentz force as the islands move together (“like currents attract”); thus, the inflow speed scales with the Alfven velocity, whereas in the island formation phase, the scaling of the reconnection rate depends on the resistivity as well. One interesting consequence of this sudden increase in reconnection rate could be found in the work of Bárta et al. (2011). Through the use of adaptive mesh refinement in MHD simulations, it was revealed that inside reconnection current sheets there exist increasingly small, self-similar, fragmenting-coalescing plasmoid structures. Lare2d has no adaptive grid refinement, but future work could examine the fields presented here at much smaller scales.

For perturbations with a dominant low-wavenumber mode, islands of larger width formed as predicted in the linear analysis of Vekstein & Jain (1998). With these larger islands come longer secondary current sheets connecting them – rather than the simple x-points which might have been expected (and which are observed for simulations in shorter numerical boxes, or for shorter wavelength disturbances). These long current sheets are observed to exhibit the plasmoid instability which increases the reconnection rate and allows access to a state of lower magnetic energy, thus allowing more effective conversion of magnetic energy into (eventually) thermal energy. Interestingly, the large plasmoid or magnetic island seen with long wavelength disturbances tends to form at the location of maximum inflow from the boundary, where an x-point rather than a o-point would intuitively be expected; this may be partly explained in terms of the linear theory of forced reconnection, which predicts that the lower-energy reconnected equilibrium has such a structure, if the initial current layer is sufficiently narrow. The fomation of this large secondary island is associated with a significant release of magnetic energy, with further energy release occurring as this island merges with the original primary island. In both these phases, a large fraction (about 10%) of the released magnetic energy is initially converted into kinetic energy, but this is subsequently dissipated, giving increases in internal energy.

For a driving perturbation with multiple wavelength components, the effects combine in a highly non-linear way, even for quite weak driving disturbances. This is partly because, according to linear theory, the island width scales as Δ^−1/2 (where Δ) is the amplitude of the boundary displacement), and hence island sizes are typically comparable or larger than the width of the initial current layer. When long-wavelength and short-wavelength disturbances of comparable amplitude are combined, the resultant magnetic field evolution and energy release are dominated by the long-wavelength element. This is as expected from linear theory, due to the increasing island width with wavelength. Thus, the initial magnetic islands form predominantly on the scale of the longest wavelength in the driver. However, the short-wavelength component has significant effects on the reconnection dynamics. It can entirely suppress the formation of large secondary plasmoids, and instead, a series of much smaller plasmoids are formed which subsequently coalesce into the primary island. The final relaxed state is reached much more quickly, and there is no clear separation of the different phases of evolution in the energy release. This has implications for particle acceleration, chiefly that in cases of faster reconnection the distribution of current sheets (and thereby of parallel electric field) is less concentrated. Many, smaller, more widely distributed acceleration sites may be expected to produce different particle energy spectra. A strongly localised perturbation behaves in remarkable similar way to a long-wavelength sinusoid, with the primary island formation dominated by the length-scale of the system.

One major limitation of this study is the use of a 2.5D scheme, since there are many important physical phenomena that cannot be captured, and complex solar magnetic fields are not well-represented by any model with simple symmetry. A fully 3D investigation of forced reconnection could yield significantly different results, in the same manner that 3D studies of the plasmoid instability have demonstrated (Kulpa-Dybeł et al. 2010; Kowal et al. 2012). Future work will consider 3D simulations, allowing consideration of important effects such as 3D turbulence. The major new factors that would be introduced in a fully 3D model are firstly, that the twisted magnetic flux ropes could be of finite-length (rather than effectively infinite), and secondly, that 3D instabilities such as the kink mode of the twisted magnetic flux ropes (plasmoids) may appear (Browning et al. 2008; Hood et al. 2009; Ripperda et al. 2017b). Restrictng the length of the flux ropes will introduce new effects due to boundary conditions, such as line-tying at the ends of the ropes, and also have significant consequences for test-particle acceleration by restricting the distance they can travel. An interesting study will be to consider the effects of applying boundary perturbations with 2D structure (rather than 1D as here), thereby producing multiple chains of islands at different planes which may interact. This may introduce additional physical effects such as field line stochasticity. Furthermore, the models might be extended in future to fully 3D configurations which more realitically represent solar magnetic fields, including loop curvature. . Nevertheless, the present work is a significant step in the direction of more physically realistic models of forced reconnection, with respect to previous models.

This is a the first paper in a series which is concerned with both the MHD aspects of forced reconnection, and its effects on acceleration of charged particles, with application to explaining the origin of non-thermal energetic particles in solar flares. The second companion paper will use a test particle approach to explore electron and ion trajectories and the evolution of particle populations in forced reconnection, with different driving disturbances. In particular, this will consider the effects of the merger of magnetic islands (twisted flux ropes in 3D) through coalescence instability on the particle energy spectra.

Acknowledgments

M. Potter acknowledges support from a STFC studentship. P. K. Browning and M. Gordovskyy acknowledge support from STFC grant ST/P000428/1. This work used the DIRAC 1, UKMHD Consortium machine at the University of St Andrews and the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by a ST/K00042X/1, STFC capital grant ST/K00087X/1, DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National EInfrastructure.

References

Appendix A: Forced reconnection in a Harris current sheet: approximate analytical solution for the perturbed flux

Here we follow the approach of Vekstein & Jain (1998; hereafter refered to as VJ98) to generate an approximate analytical solution to the perturbed flux function for a system undergoing forced magnetic reconnection. This approach assumes linearity, which is valid for small driving perturbations. In VJ98 the initial magnetic field configuration considered was B_i = {B_x(y);0;B_z(y)}: a force-free field satisfying ∇ × B = α(r)B, where there parameter α (representing shear) was constant. To describe the initial fields presented in this work (Eq. (1)) in these terms, α must be set to 1/cosh(y/y₀). The plasma β is very small so density and pressure should be unimportant for the dynamics of the process.

The model is 2D, periodic in the x-direction (x=[−x₀;x₀]) and finite in the y-direction (x=[−y₀;y₀]). All three components of the magnetic field and plasma velocity field are non-zero. The initially stable system is perturbed by plasma flow through the boundaries.

Let us assume that the plasma flow through the boundary causes a deformation of the magnetic field (“frozen-in” condition) such that the magnetic field line along the boundary becomes

$$\begin{array}{c}\hfill y_{\text{pert}}(-y_0)=-y_0+{\mathrm{\Delta }}_0\text{cos}\left(kx\right),\end{array}$$

where k = π/x₀. It is assumed that Δ/y₀ ≪ 1. Since the deformed magnetic field remains invariant in the z-direction, and defining the poloidal magnetic flux function ψ(x,y) as

$$\begin{array}{cc}\hfill B_x& =\frac{\partial \psi }{\partial y},\hfill \\ \hfill B_y& =-\frac{\partial \psi }{\partial x},\hfill \end{array}$$

then the z-component of the magnetic field is constant along contours of ψ i.e. B_z(x,y) = ℱ(ψ) (Priest 2014). As a result, ψ is a solution to the Grad–Shafranov equation (see VJ98):

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^2\psi +\mathcal{F}\frac{d\mathcal{F}}{d\psi }=0.\end{array}$$(A.1)

Let us assume that the perturbation at the boundary yields the flux function

$$\begin{array}{c}\hfill \psi (x,y)={\psi }_0\left(y\right)+\delta \psi (x,y),\end{array}$$(A.2)

where ψ₀(x) is the initial unperturbed flux function ψ = ∫B_xi(y) dx. This gives a form for the perturbation to the flux function:

$$\begin{array}{c}\hfill \delta \psi (x,y)={\psi }_1\left(y\right)\text{cos}\left(kx\right).\end{array}$$(A.3)

The initial field can be represented as

$$\begin{array}{c}\hfill \mathbf{B}=B_0\{\text{sin}\left(\theta \left(y\right)\right);0;\text{cos}\left(\theta \left(y\right)\right)\},\end{array}$$(A.4)

with the function θ(y) determined as α = dθ/dy. Since ℱ as a function of ψ remains unchanged (see VJ98), then Eq. (A.1) can be rewritten for perturbed flux ψ₁ as

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2{\psi }_1\left(y\right)}{\mathrm{d}{y}^2}+[{\left(\frac{\mathrm{d}\theta \left(y\right)}{\mathrm{d}y}\right)}^2-k^2-\frac{1}{\text{tan}\left(\theta \left(y\right)\right)}\frac{{\mathrm{d}}^2\theta \left(y\right)}{\mathrm{d}{y}^2}]{\psi }_1\left(y\right)=0.\end{array}$$

It is at this point we depart from VJ98: the case they considered had constant α, here we look at the case where α = 1/cosh(y/y₀), which, in conjunction with Eq. (A.4), allows the fields described in Eq. (1) to be recovered.

This gives us θ(y) = asin(tanh(y/y₀)), which reduces our Grad–Shafranov to

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2{\psi }_1\left(y\right)}{\mathrm{d}{y}^2}+(\frac{2}{y_0^2}{\frac{1}{\text{cos}h}}^2(y/y_0)-k^2){\psi }_1=0.\end{array}$$(A.5)

This second-order differential equation is similar to the Poschl–Teller type equation, the analytical solution to which was considered by Barut et al. (1987). We employ a similar procedure to solve Eq. (A.5).

First, let us introduce the new variable z = −sinh²(y/y₀). Then, it follows that

$$\begin{array}{cc}& \frac{\mathrm{d}{\psi }_1}{\mathrm{d}y}=-\frac{2}{y_0}\frac{\mathrm{d}{\psi }_1}{\mathrm{d}z}\text{sin}h(y/y_0)\text{cos}h(y/y_0),\hfill \\ \hfill & \frac{{\mathrm{d}}^2{\psi }_1}{\mathrm{d}{y}^2}=\frac{2}{y_0^2}(2z-1)\frac{\mathrm{d}{\psi }_1}{\mathrm{d}z}+\frac{4}{y_0^2}z(1-z)\frac{{\mathrm{d}}^2{\psi }_1}{\mathrm{d}{z}^2}·\hfill \end{array}$$

By substitution, Eq. (A.5) becomes

$$\begin{array}{c}\hfill 4z(1-z)\frac{{\mathrm{d}}^2{\psi }_1}{\mathrm{d}{z}^2}+2(1-2z)\frac{\mathrm{d}{\psi }_1}{\mathrm{d}z}-(\frac{2}{1-z}-k^2{y}_0^2){\psi }_1=0.\end{array}$$(A.6)

Further, we introduce function F(z):

$$\begin{array}{cc}& {\psi }_1\left(z\right)={(1-z)}^{\frac{\lambda +1}{2}}F\left(z\right)\hfill \\ \hfill & \frac{\mathrm{d}{\psi }_1}{\mathrm{d}z}=-\frac{\lambda +1}{2}{(1-z)}^{\frac{\lambda -1}{2}}F\left(z\right)+{(1-z)}^{\frac{\lambda +1}{2}}F\prime \left(z\right)\hfill \\ \hfill & \frac{{\mathrm{d}}^2{\psi }_1}{\mathrm{d}{z}^2}=\frac{{\lambda }^2-1}{4}{(1-z)}^{\frac{\lambda -3}{2}}F\left(z\right)-\frac{\lambda +1}{2}{(1-z)}^{\frac{\lambda -1}{2}}F\prime \left(z\right)\hfill \\ \hfill & +{(1-z)}^{\frac{\lambda +1}{2}}F″\left(z\right).\hfill \end{array}$$

Substituting the above terms into Eq. (A.6) gives

$$\begin{array}{cc}\hfill z(1-z)F″+& (\frac{1}{2}-(\lambda +2)z)F\prime \hfill \\ \hfill & -\frac{1}{4}(\frac{{(\lambda +1)}^2(1-z)-{\lambda }^2-\lambda +2}{(1-z)}-k^2{y}_0^2)F=0.\hfill \end{array}$$(A.7)

By choosing λ² + λ − 2 = 0, the above equation can be reduced to the form of a hypergeometric equation:

$$\begin{array}{c}\hfill z(1-z)F″+[c-(a+b+1)z]F\prime -abF=0.\end{array}$$(A.8)

The solution to Eq. (A.8) is known in the following form (see Milne-Thomson et al. (1972)):

$$\begin{array}{cc}& F=C_{12}{F}_1(a,b,c;z)\hfill \\ \hfill & +C_2{z}^{1-c}{}_2{F}_1(a+1-c,b+1-c,2-c;z),\hfill \end{array}$$

where C₁ and C₂ are arbitrary constants and ₂F₁ is the Gauss hypergeometric function.

The parameters a, b and c can be determined by comparing Eqs. (A.7) and (A.8), provided λ is one of two roots (λ₁ = 1 and λ₂ = −2). Hence, c = 1/2 and a, b = (1/2)((λ+1)±ky₀).

Therefore, the solution to Eq. (A.5) can be written as

$$\begin{array}{cc}& \psi \left(y\right)=C_1{\text{cos}h}^{\lambda +1}\left(\frac{y}{y_0}\right)\hfill \\ \hfill & {}_2{F}_1[\frac{(\lambda +1+k{y}_0)}{2},\frac{(\lambda +1-k{y}_0)}{2};\frac{1}{2};-{\text{sin}h}^2\left(\frac{y}{y_0}\right)]\hfill \\ \hfill & +C_2\text{sin}h\left(\frac{y}{y_0}\right){\text{cos}h}^{\lambda +1}\left(\frac{y}{y_0}\right)\hfill \\ \hfill & {}_2{F}_1[\frac{(\lambda +2+k{y}_0)}{2},\frac{(\lambda +2-k{y}_0)}{2};\frac{3}{2};-{\text{sin}h}^2\left(\frac{y}{y_0}\right)]·\hfill \end{array}$$

By imposing boundary conditions of ψ(±a) = δ₀B₀, where a is the width of the domain in y, it is trivial to determine the constant C₂ = 0 and

$$\begin{array}{c}\hfill C_1={\delta }_0{B}_0{\left[{\text{cos}h}^2{\left(\frac{a}{y_0}\right)}_2{F}_1(\frac{2+k{y}_0}{2},\frac{2-k{y}_0}{2};\frac{1}{2};-{\text{sin}h}^2\left(\frac{a}{y_0}\right))\right]}^{-1},\end{array}$$(A.9)

which gives the full solution for the perturbed flux as

$$\begin{array}{c}\hfill {\psi }_1=C_1{\text{cos}h}^2{\left(\frac{y}{y_0}\right)}_2{F}_1[\frac{2+k{y}_0}{2},\frac{2-k{y}_0}{2};\frac{1}{2};-{\text{sin}h}^2\left(\frac{y}{y_0}\right)]·\end{array}$$(A.10)

Figures A.2 and A.3 show contours of the total flux function for different values of current layer thickness y₀, while Fig. A.1 shows the corresponding perturbed flux ψ₁ as a function of y. For large values of y₀, the perturbed flux approaches that of the VJ98 solution and that is reflected in the islands forming in at the edges of Fig. A.2 as expected. However, for thin current sheets the perturbed flux becomes negative about the y = 0 line. This results in contours of the total flux revealing islands which form in the centre of the plot: in the inflow region. This is a novel and surprising result, but is partly supported by the numerical work presented in this paper where, in cases where current sheets were sufficiently thin, islands were observed to form in the inflow region.

Fig. A.1. ψ1 against y for large and small values of current layer thickness y0.

Fig. A.2. Contours of the total flux function (Eqs. (A.2) and (A.3)) for y0 = 16.0 i.e. much greater than 1.0, approaching the VJ98 solution (y0 = ∞).

Fig. A.3. Contours of the total flux function for y0 = 1/3. For all thin current layers (y0 ≪ 1) the islands form in the inflow region.

All Tables

All Figures

Fig. 1. Left panel: diagram representing the initial conditions. The domain is periodic in the x-direction, with rigid boundaries in the y-direction after the perturbation is applied. The guide field Bz (out of the plane) is shown in blue, overlaid with field lines showing the x-component of the magnetic field Bx. The red dotted line shows the form of the perturbation applied at the boundary: a sinusoidal inflow and outflow pattern with red arrows indicating direction, equal and opposite at the upper and lower y boundaries for a short period of time (see Eq. (9)). Right panel: the initial force free magnetic field components.

In the text

	Fig. 2. B and Bz for case 1a (single wavelength λ = 4.0) at final snapshot t = 100τA; a replication of the final state of the model by Gordovskyy et al. (2010b).
In the text

Fig. 3. Energetics for case 1b (Table 1), with corresponding MHD fields shows in Fig. 4. At around 100tA, the first stage of reconnection is slowing to a halt as islands reach saturation. At around 150tA the large drop in magnetic energy (panel a) and spike in kinetic energy (panel b) correspond to the island coalesence event seen in Fig. 4 (panels e–f). The total energy is shown as an indication of numerical accuracy over the course of the simulation: the small initial rise is due to the driving at the boundary, remaining constant thereafter (observed in all simulations, witihin 5%).

In the text

Fig. 4. MHD snapshots for case 1b (Table 1) – a single sinusoidal perturbation in a wide box. Left to right panel: magnetic field lines B with Jz contours, magnetic field lines with Bz contours and the plasma flow v. Arrows on v plots indicate the in-plane velocity direction. Snapshots (panel a) through (panel f) correspond to times t = 0, 8, 64, 112, 144 and 300τA. Note that the aspect ratio is not to scale – the width of the domain 16 in x and 2 in y. The colour scheme is adaptive to enhance contrast. Note that colour representing zero is consistent across similar plots throughout.

In the text

	Fig. 5. B and Jz for case 1c (single wavelength λ = 2.0) at final snapshot t = 350τA.
In the text

	Fig. 6. Panels a–d: B and Jz for case 1d (single wavelength λ = 16.0) at time t = 40, 92, 196 and 212τA. Panels e–f: corresponding energetics: magnetic energy (left panel) and kinetic energy (right panel). The snapshot times shown correspond to the different stages of reconnection: initial perturbation of equilibrium, first saturated state, island attraction and island coalescence.
In the text

Fig. 7. Panels a–d: B and Jz for case 2a (multiple wavelengths λ1 = 2.0, λ2 = 16.0) at time t = 24, 40, 92 and 350τA. Panels e–f: corresponding energetics: magnetic energy (left panel) and kinetic energy (right panel). Case 2a uses the perturbations from cases 1c and 1d in tandem, resulting in different reconnection dynamics and energetics, such as the lack of formation of the central island as seen in Fig. 6.

In the text

Fig. 8. Panels a–d: B and Bz for case 2b (multiple wavelengths λ1 = 2.0, λ2 = 16.0, with increased amplitude on the smaller wavelength component such that δ1 = 0.35) at time t = 8, 20, 28 and 350τA. The simulation reaches its final state much more rapidly than in case 2a due to the increased total flux. Bz was chosen to be displayed rather than Jz since, in this case, the island structure was more visible.

In the text

	Fig. 9. B and Jz for case 4a (localised perturbation) at time t = 24, 40, 60, 84, 112 and 152τa.
In the text

	Fig. A.1. ψ1 against y for large and small values of current layer thickness y0.
In the text

	Fig. A.2. Contours of the total flux function (Eqs. (A.2) and (A.3)) for y0 = 16.0 i.e. much greater than 1.0, approaching the VJ98 solution (y0 = ∞).
In the text

	Fig. A.3. Contours of the total flux function for y0 = 1/3. For all thin current layers (y0 ≪ 1) the islands form in the inflow region.
In the text