Three-dimensional Oscillatory Magnetic Reconnection

Magnetic fields play a key role in determining the dynamics of plasmas at all scales: from fusion experiments and laboratory plasmas to planetary magnetospheres, the Sun and stars, and galaxies and accretion disks. Magnetic reconnection is a fundamental plasma process associated with dynamic energy release in these systems, and is believed to explain a broad range of phenomena including solar and stellar flares, coronal mass ejections, astrophysical jets, and planetary aurorae. Current theoretical frontiers in fundamental reconnection research can be broadly categorized as (i) collisionless/kinetic effects (Yamada et al. 2010), (ii) the extension to 3D of a long history of 2D theory (where profound differences make this extension highly non-trivial; Priest et al. 2003; Pontin 2012), (iii) time-dependent and transient effects, and (iv) the interaction of the local reconnection dynamics with those of the global system. Here, we present significant results concerning the latter three within the framework of nonlinear magnetohydrodynamics (MHD).

Much of our present understanding of reconnection has focused primarily on continuously driven (and hence at least quasi-steady) systems, while transient effects that are crucial in many plasma environments are typically neglected (we note an important exception in the case of current sheet instabilities, which are beyond the scope of this paper, where studies are necessarily time-dependent, e.g., Loureiro et al. 2007, Comisso et al. 2016). Combined with questions regarding the influence of chosen boundary conditions, the applicability of any given reconnection model to a particular situation is often unclear. Here we study instead the self-consistent formation and evolution of a reconnecting current sheet in response to a finite external perturbation. The reconnection is triggered in the vicinity of a null point, a natural weakness in the magnetic field, by an external driver in the form of an MHD wave. It is well established that such null points trap these waves through refraction (McLaughlin et al. 2011), leading to the formation of electric current layers around the null in which reconnection may take place (Pontin & Craig 2005; Pontin et al. 2007; McLaughlin et al. 2009). The subsequent reconnection dynamics triggered by the waves in a configuration where the reconnection region is connected naturally to an external non-reconnecting region has so far not been studied. This is the aim of this paper.

Within the reduced framework of 2D MHD, where the assumption of an invariant direction simplifies the system greatly, it has recently been demonstrated that in such a scenario reconnection does not proceed in a quasi-steady manner but rather exhibits a time-dependent behavior, termed oscillatory reconnection (Craig & McClymont 1991; Ofman et al. 1993; McLaughlin et al. 2009; Murray et al. 2009; McLaughlin et al. 2012; Pucci et al. 2014). However, the extension to 3D is not trivial; the process by which waves trigger reconnection (null collapse, see e.g., Priest & Forbes 2000) is little studied in 3D, and also reconnection itself can be fundamentally changed (Priest et al. 2003). Additionally, mode conversion between magnetoacoustic and Alfvén waves becomes permissible.

Here we demonstrate for the first time that impulsively initiated reconnection at realistic 3D magnetic null points naturally proceeds in a time-dependent and oscillatory fashion, and we also demonstrate that 3D reconnection can itself produce MHD waves. The results show that reconnection is inherently periodic, and they highlight the deeply interconnected nature of MHD waves and reconnection.

The simulation involves the numerical solution of the MHD equations using the Lare3d code (Arber et al. 2001). Here we outline the initial conditions, with full technical details deferred to the appendix. All variables in this paper are nondimensionalized. We consider a background magnetic field of the form

$$\begin{eqnarray}&&{{\boldsymbol{B}}}_{0}=[x,y,-2z],\end{eqnarray} \tag{ 1 }$$

which is known as the linear proper potential null point (Parnell et al. 1996). The magnetic null point itself (${\boldsymbol{B}}=0$) is located at the origin. The field is free from electrical currents and so constitutes a minimum magnetic energy state. Topologically, it consists of a spine field line, running along the z-axis toward the null point, and a fan plane z = 0, consisting of field lines pointing radially outward. Other field lines, separated by the fan, have a hyperbolic structure (Figure 1(a)).

Zoom In Zoom Out Reset image size

Upon this magnetic field we impose a finite-amplitude perturbation ${{\boldsymbol{B}}}_{1}={\boldsymbol{\nabla }}\times {{\boldsymbol{A}}}_{1}$, where

$$\begin{eqnarray}&&{{\boldsymbol{A}}}_{1}=\psi \,\exp \left[\displaystyle \frac{-({x}^{2}+{y}^{2}+{z}^{2})}{2{\sigma }^{2}}\right]\hat{{\boldsymbol{y}}},\end{eqnarray} \tag{ 2 }$$

which is a ring of magnetic flux centered about the null point, circulating in planes of fixed y (Figure 1(b)). When the simulation begins, the perturbation immediately splits into inward and outward propagating magnetoacoustic waves. The incoming part is our intended driver for initiating magnetic reconnection and acts implosively to collapse the magnetic field into the first reconnection configuration. The outgoing part of the perturbation leaves the initialization site and propagates outward. Special care is taken so that it and any subsequently generated waves do not reflect at the boundaries and return inward to influence our solution at the null (see the appendix).

The background plasma is taken to be initially stationary (${\boldsymbol{v}}=0$) and of uniform density ($\rho =1$). We are now left with a choice of a value to assign to the perturbation energy (set by ψ and σ), and the uniform pressure p and resistivity η. As ${{\boldsymbol{B}}}_{0}$ is scale free, the relative sizes of these three parameters control the dynamics of the initial implosion. Here we present results for $\psi =0.05$, $\sigma =0.21$, $\eta ={10}^{-3}$, and p = 0.005. This choice of parameters is such that the implosion forms a localized current sheet that is isolated from the boundaries and free to evolve in a self-consistent manner.

3.1. Initial Implosion

The incoming wave is initially symmetric as per the geometry of the flux ring and generates a cylindrical ring current flowing through the null. It propagates toward the null point as a fast magnetoacoustic mode. In linearized MHD it would evolve in a self-similar manner, focusing its energy toward the null, increasing in amplitude. As a result of the linearly decreasing Alfvén speed profile (${c}_{A}\propto | {B}_{0}|$), the length scales decrease exponentially, and so the current density at the null point increases exponentially in time in the linear regime while retaining its cylindrical symmetry (McLaughlin & Hood 2004; McLaughlin et al. 2008). This focusing would continue until either resistive diffusion or plasma back-pressure act to halt further collapse, the scalings for which are well understood only in 2D (Craig & McClymont 1991, 1993; Hassam 1992; Craig & Watson 1992; McClymont & Craig 1996). However, given that our simulation considers a full nonlinear solution, the opportunity to increase the amplitudes to allow the excess flux carried by the wave to overwhelm the background field and so form shocks.

The approaching pulse imbalances the distribution of magnetic flux about the null point, enhancing and diminishing the overall field strength in different lobes about the null point's spine and fan. The anisotropic nature of the Lorentz force leads to a net force that is directed alternately toward and away from the null point, in the regions of magnetic enhancement and rarefaction, respectively. It follows that the implosion carries fronts of magnetic and plasma compression and rarefaction about the spine and fan. For a sufficiently energetic perturbation (such as the one considered here), nonlinear steepening occurs in these compressive fronts, leading to a breaking of the cylindrical symmetry. Previous studies of 2D null collapse demonstrated that the implosion becomes quasi-1D as its fronts accelerate and stall in regions of compression and rarefaction, leading to collapse of the separatricies toward one another (Craig & Watson 1992; McClymont & Craig 1996; Gruszecki et al. 2011).

An analogous breaking of the symmetry is observed in our simulations, with the disturbance flattening into a planar structure. Figure 2 shows selected transparent isosurfaces of the y-component of the current at an early time (t = 0.1) and at the approximate time at which the collapse is halted (t = 0.6). The y-component of the current is component that flows through the null, generating the parallel electric field required for 3D reconnection. We see that at t = 0.1, the current profile has already started to depart from the initial cylindrical symmetry in planes of fixed y, and by t = 0.6 a quasi-planar current sheet has been formed by the aforementioned processes of nonlinear steepening in the compressive fronts. The increase in amplitude is due to the wave focusing. Thus, the perturbation creates a current sheet that is localized in the vicinity of the null point, and so we find that 3D null collapse, for our choice of parameters, is qualitatively similar to 2D collapse in planes of fixed y. It differs in that the current distribution is fully localized in 3D, however, as can be seen in Figures 2(a) and b. As we show in the following section, there is a crucial difference in that the 3D collapse triggers a 3D reconnection mode, which is significantly different to 2D. We highlight that this sheet is detached from the influence of the boundaries and has been generated in an impulsive rather than a driven manner—once the current sheet is formed, no further energy is injected to sustain the sheet or otherwise alter it.

Zoom In Zoom Out Reset image size

3.2. Reconnection Reversals

3.3. Reconnection Jets and Reversal Dynamics

The spine-fan reconnection redistributes not only magnetic flux about the null point, but also plasma, with the establishment of a system of flow through the current sheet (Figure 5). It is not a simple "disk-shaped" outflow, relieving the compression of the implosion symmetrically, but plasma is instead directed predominantly along the current sheet axis as this is the dominant direction of the Lorentz force associated with sharply curved reconnected field lines. Along the outflow, plasma is accelerated to super-magnetosonic speeds in two reconnection jets. The jets are bounded by standing slow-mode MHD shock fronts that separate the field lines into the inflow and outflow regions and are analogous to those found in Petschek's 2D model (Petschek 1964).

Zoom In Zoom Out Reset image size

As the outflow evolves, we observe the pooling of hot overdense plasma at the heads of the reconnection jets, leading to back-pressures that choke off the outflow. This is accompanied by the relief of (plasma) compression in the current sheet itself via the expulsion of the excess plasma swept in by the initial implosion. A plasma rarefaction forms in the inflow lobes, with plasma accelerated through the current sheet unreplenished in the absence of continual driving flows. The consequence of this plasma and flux redistribution is that the balance of forces acting across the spine and fan planes (the collapse of which is necessary to sustain the reconnection and current sheet) must change as reconnection proceeds.

The competition of these effects is the source of dynamism in oscillatory reconnection. Shortly after the initial collapse halts, strong plasma pressure gradients build up at the jet heads, exerting a force against the outflow. In conjunction with deleterious effects in the inflow region (such as the diffusion/redistribution of the driving magnetic field and current sheet decompression), this leads to a deceleration of the jets, which in turn leads to the formation of fast MHD shocks (termination shocks). These refract magnetic field toward the shock tangent with increases in overall plasma pressure and field strength. This bending of the field away from the normal gives rise to the opposite-polarity deflection currents that were earlier identified in Figure 3. These fast nonlinear waves formed in the jet head travel in toward the null, and act implosively with similar physics to the initial collapse. Thus, a nonlinear fast wave collapses in a quasi-1D fashion, altering the magnetic field in the vicinity of the null, first undoing the spine-fan collapse of the initial implosion, then overshooting into a secondary collapse that pushes the spine and fan together in the opposite sense. This reverses the sign of current at the null, and thus reverses the sense of the reconnection. With the completion of the reversal and the secondary implosion, we see the establishment of new (weaker) reconnection flows that drive the next reversal via the same processes.

3.4. Wave Generation

In these simulations the reconnection only occurs in a small region near the null point. However, the global effects of the magnetic field restructuring are not artificially confined to this area by computational boundaries, but rather may escape the vicinity of the current layer as freely propagating waves. We find that many escaping MHD waves are generated and highlight a selection here.

Figures 6(a), (b) shows time-distance diagrams that trace disturbances seen in $| {\boldsymbol{v}}|$ and ρ propagating along the positive z-axis. A number of different wavefronts are visible. The features labeled I and O correspond to the incoming and outgoing waves, respectively, that are generated by our initial perturbation, i.e., they are not generated by reconnection.

Zoom In Zoom Out Reset image size

We identify the features marked A as Alfvén waves. They propagate at the linearly increasing background Alfvén speed, which leads to the exponential profile along its fronts (by the solution of ${dz}/{dt}=| {B}_{z}| =2z$). They are, as is characteristic of Alfvén waves, incompressible, and as such are not visible in the figure for ρ (Figure 6(b)). These waves are also manifest in the generation of oppositely signed vorticity tubes in the reconnection region that propagate up along either side of the spine line (Figure 7). Vorticity propagation is a marker of Alfvén waves in 3D MHD, as it corresponds to a rotational disturbances of the magnetic field that allows magnetic tension to act as the predominant restoring force giving rise to the key properties of the wave (Goossens et al. 2011; Thurgood & McLaughlin 2012; Shelyag et al. 2013). These waves are generated by the swaying of the spine line and neighboring field lines from side to side along the x-axis during each reversal event, which produces two counter-rotational vortices in the flow pattern surrounding the spine (somewhat reminiscent of the "m = 1" kink wave in coronal loops). In this way, the vorticity produced by the magnetoacoustic modes driving the swaying of the spine to the Alfvén mode in a process of mode conversion.

Download figure:

Video Standard image High-resolution image

The features marked S are identified as large-amplitude, anharmonic slow-mode waves—due to their magnetoacoustic propagation speed (always lower than the Alfvén speed)—and have signal present in both the velocity and density diagrams. They are excited in the outflow regions of the reconnection jets (i.e., the slow waves and termination shock both have the same source). Owing to their strong nonlinearity, these anharmonic waves not only propagate plasma compression, but transport material from the reconnection outflow away from the null and out along the spine axis (and fan plane). The profile of material carried away from the null in y = 0 at t = 4 by such waves is also shown in Figure 8; note the multiple large-amplitude fronts along the spine ($x\approx 0$) and fan ($z\approx 0$) (which indicates that these waves propagate in a field-aligned sense, as expected of slow waves—see also the corresponding movie), and the underdense region that remains in the vicinity of the null point as a consequence of the material transport.

Download figure:

Video Standard image High-resolution image

We have demonstrated unambiguously that magnetic reconnection at a 3D null point, driven by a perturbation of finite duration, naturally proceeds in a time-dependent oscillatory manner. We find that the process consists of four key elements, namely (i) the creation of current sheets by MHD waves in a process of 3D null collapse, (ii) flux and plasma redistribution via localized spine-fan reconnection, (iii) associated reconnection reversals due to back-pressure-driven overshoots, and (iv) energy escape processes in the form of wave generation.

The identification of 3D oscillatory reconnection is a milestone in our understanding of energy release in high magnetic Reynolds number plasmas, demonstrating that reconnection triggered in an aperiodic manner produces both periodic reconnection reversals and periodically excited propagating MHD waves through a self-generated oscillation. The inherent periodicity of 3D reconnection events and the generation of waves may both play a role in explaining the observed quasi-periodicity in solar and stellar flare emission (Nakariakov & Melnikov 2009; Nakariakov et al. 2016; Van Doorsselaere et al. 2016). Understanding the physical mechanism underpinning this quasi-periodicity is crucial for developing the next generation of flare models and for exploiting these observations as a diagnostic tool (via magnetoseismology, see De Moortel & Nakariakov 2012). At present, there are a number of competing physical mechanisms that may explain quasi-periodic pulsations (QPPs), including oscillatory reconnection, the magnetic tuning fork (Takasao & Shibata 2016), self-oscillatory processes such as magnetic-dripping (Nakariakov et al. 2010), and dispersive wave trains (Nakariakov et al. 2004) (see also the aforementioned review papers). Our work finds that oscillatory reconnection may proceed at fully 3D nulls and therefore is one potential mechanism that may provide an explanation for QPPs, since a subset of flaring events involves truly 3D null-point geometries as opposed to quasi-2D X-lines (such as seen in, e.g., Masson et al. 2009; Zuccarello et al. 2009; Zhang et al. 2012).

The initial conditions presented in this paper were chosen by experimentation to correspond physically to an MHD pulse impinging on the null that is sufficiently energetic (relative to the plasma pressure and resistivity) that it implodes in a nonlinear fashion according to McClymont & Craig (1996) and forms a current sheet that then proceeds to participate in oscillatory reconnection. Pulses that are not sufficiently energetic, e.g., lower amplitudes (for the same pressure and resistivity), remain linear and do not perturb the spine and fan asymmetrically. They therefore do not form current sheets but rather ring currents, which evolve in a qualitatively different fashion to the system presented here. As our experiment is free of global scales, its applicability to a given scenario, such as determining its period for a typical solar event, must be addressed by combining situation-specific global modeling with analysis of observational and experimental data. We expect that in the solar atmosphere, given the abundance of energetic waves (e.g., Long et al. 2017), flux emergence events (e.g., Leake et al. 2013, 2014), low pressures, and low resistivity, the nonlinear regime of null-point collapse is expected to be accessible and so oscillatory reconnection should proceed in a qualitative sense as presented here. Determining the specific periods accessible by the mechanism in a solar scenario will require models that place observational restrictions on permissible null geometries, driving energies, and plasma pressures, while also accounting for gravitational stratification.

Furthermore, MHD waves are observed to be ubiquitous in the solar atmosphere (Tomczyk et al. 2007; McIntosh et al. 2011; Morton et al. 2015). Understanding these waves is essential since they are believed to play a central role in the atmospheric mass and energy balance; the nature of their genesis is an outstanding question, however. The generation of mass-transporting slow modes that are due to reconnection may explain the origin of propagating intensity disturbances observed in the open-field corona (DeForest & Gurman 1998; Banerjee et al. 2011). Additionally, models predict that the dissipation of torsional Alfvén waves could provide significant chromospheric and coronal heating (Antolin & Shibata 2010). Thus, the natural excitation of torsional Alfvén waves by oscillatory reconnection is especially interesting given that twist and torsional motions have recently been shown to be highly prevalent in the low solar atmosphere (De Pontieu et al. 2014).

There are also significant implications for the restructuring of the global field in response to the reconnection process. For example, reconnection at a 3D null point in the Sun's atmosphere can facilitate the exchange of plasma between magnetic field lines that are closed (connect at both ends to the solar surface) and those that are open to interplanetary space (Pontin et al. 2013), driving material outflows (Bradshaw et al. 2011). Our results thus provide a potential mechanism for (quasi-)periodic ejections observed in solar jets (Morton et al. 2012; Chandrashekhar et al. 2014; Zhang & Ji 2014; Zhang et al. 2016).

Additionally, to our knowledge, Figure 5 is the first detailed consideration of Petschek-like reconnection jets in a spine-fan system, replete with standing slow shocks and fast termination shocks, with particular implications for nonthermal particle acceleration. We note that the steady-state Petschek reconnection in uniform-resistivity MHD is widely considered unattainable, with the requirement of an anomalous resistivity, localized at the null, to act as an obstacle in the flow (e.g., Biskamp & Schwarz 2001). Here we stress that the finite-duration "Petschek-like" system obtained in this work, with uniform resistivity, is intimately connected to the time-dependent nature of our problem.

Finally, we note that previous studies involving (2D) null collapse were chiefly motivated in terms of examining the possibility of the collapse to small scales and high current densities, yielding favorable scaling laws of the reconnection rate with resistivity. Early results in the zero-β limit suggested a "fast" scaling, with decay rates of the perturbations scaling as $\sim | \mathrm{ln}\eta {| }^{2}$ in the linear regime and becoming virtually independent of η in the nonlinear case (so-called super-fast reconnection according to Forbes 1982; McClymont & Craig 1996). However, studies at finite-β suggest that in the corona back-pressure in the current layer may halt the collapse before it enters a regime of fast reconnection. There are a number of possible resolutions—in MHD this includes the possibility of strong nonlinear effects associated with highly energetic perturbations mitigating the effects of pressure, which leads to secondary faster phases of reconnection, as hypothesized by McClymont & Craig (1996) and discussed in Priest & Forbes (2000). The possibility also exists that even if the collapse is pressure limited, provided the sheet width can reach kinetic length scales, then the resulting microphysical reconnection rates will be sufficiently fast. In the context of 2D null-point collapse, the fast rates of collisionless reconnection have been demonstrated by particle-in-cell simulations in a series of papers by Tsiklauri and collaborators (Tsiklauri & Haruki 2007, 2008; Tsiklauri 2008; Graf von der Pahlen & Tsiklauri 2014, 2016). As such, an important outstanding issue is the efficiency of 3D null collapse as a reconnection and dissipation mechanism. In this paper, 3D null-collapse features as the means of establishing a localized current sheet that precipitates oscillatory reconnection, which itself is our main focus. We do note, however, that we have found that 3D null collapse is at least qualitatively similar to the known 2D behavior. It will be important to conduct a detailed quantitative study of the scaling of the 3D null-collapse mechanism and its relation to the 2D case.

The authors acknowledge generous support from the Leverhulme Trust, and this work was funded by a Leverhulme Trust Research Project Grant: RPG-2015-075. The authors acknowledge IDL support provided by STFC. The computational work for this paper was carried out on HPC facilities provided by the Faculty of Engineering and Environment, Northumbria University, UK.

The simulation is the numerical solution of the nondimensional resistive MHD equations:

$$\begin{eqnarray}&&\displaystyle \frac{D\rho }{{Dt}}=-\rho {\rm{\nabla }}\cdot {\boldsymbol{v}}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{D{\boldsymbol{v}}}{{Dt}}=\displaystyle \frac{1}{\rho }({\rm{\nabla }}\times {\boldsymbol{B}})\times {\boldsymbol{B}}-\displaystyle \frac{1}{\rho }{\rm{\nabla }}p+{{\boldsymbol{F}}}_{\mathrm{shock}}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\displaystyle \frac{D{\boldsymbol{B}}}{{Dt}}=({\boldsymbol{B}}\cdot {\rm{\nabla }}){\boldsymbol{v}}-{\boldsymbol{B}}({\rm{\nabla }}\cdot {\boldsymbol{v}})-{\rm{\nabla }}\times (\eta {\rm{\nabla }}\times {\boldsymbol{B}})\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\displaystyle \frac{D\epsilon }{{Dt}}=-\displaystyle \frac{p}{\rho }{\rm{\nabla }}\cdot {\boldsymbol{v}}+\displaystyle \frac{\eta }{\rho }{j}^{2}+\displaystyle \frac{{{\boldsymbol{H}}}_{\mathrm{visc}}}{\rho }\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\boldsymbol{j}}={\boldsymbol{\nabla }}\times {\boldsymbol{B}}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\boldsymbol{E}}=-{\boldsymbol{v}}\times {\boldsymbol{B}}+\eta \,{\boldsymbol{j}}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&p=\epsilon \rho (\gamma -1),\end{eqnarray} \tag{ 9 }$$

which are solved using the Lare3d code. All results presented are in nondimensional units. Algorithmically, the code solves the ideal MHD equations explicitly using the Lagrangian remap approach and includes the resistive terms using explicit subcycling (Arber et al. 2001, 2016). The solution is fully nonlinear and captures shocks via an edge-centered artificial viscosity approach (Caramana et al. 1998), where shock viscosity is applied to the momentum equation through ${{\boldsymbol{F}}}_{\mathrm{shock}}$ and heats the system through ${{\boldsymbol{H}}}_{\mathrm{visc}}$. Extended MHD options available within the code, such as the inclusion of Hall terms, were not used in these experiments. Full details of the code can be found in the original paper (Arber et al. 2001) and the users manual. Requests for access to the code should be addressed to its developers.

A summary of the (nondimensional) initial conditions on the primitive variables, discussed in the main document, is as follows:

$$\begin{eqnarray}&&p=0.005\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\rho =1\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{\boldsymbol{v}}=0\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{\boldsymbol{B}}={{\boldsymbol{B}}}_{0}+{\boldsymbol{\nabla }}\times {{\boldsymbol{A}}}_{1}\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{{\boldsymbol{B}}}_{0}=[x,y,-2z]\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{{\boldsymbol{A}}}_{1}=\psi \,\exp \left[\displaystyle \frac{-({x}^{2}+{y}^{2}+{z}^{2})}{2{\sigma }^{2}}\right]\hat{{\boldsymbol{y}}}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\psi =0.05\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\sigma =0.21,\end{eqnarray} \tag{ 17 }$$

with resistivity taken uniformly as $\eta ={10}^{-3}$ throughout. The ratio of specific heats is taken as $\gamma =5/3$.

The boundary conditions are taken as no-slip (${\boldsymbol{v}}=0$) with zero-gradient conditions on ${\boldsymbol{B}}$, ρ and p. A simulation was run with no perturbation to check that the boundaries do not launch spurious waves into the domain, and to check the overall stability of the setup.

Special care must be taken so that outgoing waves, whether from the initial perturbation or those subsequently generated, do not reflect off the domain boundaries and return inward to influence our solution at the null. During our initial testing, we found that open boundaries and artificial damping zones are imperfect strategies for removing outgoing waves in linear null-point geometries. The only robust method of avoiding reflections that influence the solution is therefore to simply place the boundaries sufficiently far away from the null point that the shortest possible signal travel time from the initialization site to the boundary and back is in excess of the period we wish to study. Computationally, this is not a trivial task. The shortest signal travel time is governed by the Alfvén speed, which increases linearly with distance from the null point (at the center of the domain). Thus, the overall travel time to the boundary increases only logarithmically with increasing distances to the boundaries. Boundaries must be placed at vast distances while maintaining sufficient resolution in the dynamic and diffusive regions close to the null point.

For the results shown here, we placed the x, y, and z boundaries at ±40. The shortest possible signal travel time permitted is that of a wave propagating at the Alfvén speed directly along the spine (where B_z increases as $2z$), and so by separation of variables the travel time from the outer edges of the initial perturbation $z\approx \pm 0.25$ out to the boundary is ${t}_{{out}}\approx 0.5\mathrm{ln}(40/0.25)=2.53$. The fastest possible reflection will return to the diffusion-dominated region about the null (estimated at the length of the initial current sheet, $\approx 0.18$) after a further ${t}_{{in}}\approx 0.5\mathrm{ln}(40/0.18)=2.70$, yielding the shortest possible theoretical reflection time of ${t}_{\mathrm{reflect}}\approx 5.23$. This was also practically tested with a "boundary convergence test", where simulations were run with boundaries that were closer to the null (hence having shorter reflection times). We found that results only differed at times in excess of the theoretical reflection time calculated as above for each case. In these test cases, it was often possible to track reflections along paths using time-distance diagrams as in Figure 6(a), which further confirmed these theoretical times. The maximum boundary length used (±40) was a compromise between maximizing the reflection time and maintaining the resolution available close to the null—note that doubling the distance to the boundary will only increase the time by $\mathrm{ln}(2)$. By simple wave theory, the fastest possible theoretical reflection time is therefore ${t}_{\mathrm{reflect}}\approx 5.23$, and so the non-interference of reflections during the initial implosion, the first reversal, and the initial stages of the second reversal is guaranteed. Therefore, the self-generated oscillations of the null point are indeed due to reconnection itself. Furthermore, we note in the later stages of the second reconnection reversal ($t\gt {t}_{\mathrm{reflect}}$) up to the end time of t = 6 that no incoming fronts of significant intensity are detected in Figure 6(a), which is in fact taken along this shortest possible signal travel path.

To accommodate the vast distances to the boundary while sufficiently resolving the diffusion at the null point and the wave dynamics in the surrounding region, a grid-stretching scheme is employed.

The grid is stretched according to a variation on the scheme of Roberts (1971), namely the cell boundary positions x_b along the x-direction are distributed according the transformation:

$$\begin{eqnarray}&&{x}_{b}=40\left\{1+\displaystyle \frac{\sinh [\lambda ({\xi }_{i}-{\rm{\Gamma }})]}{\sinh [\lambda {\rm{\Gamma }}]}\right\}-40\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{\rm{\Gamma }}=\displaystyle \frac{1}{2\lambda }\mathrm{ln}\left[\displaystyle \frac{1+({e}^{\lambda }-1)0.5}{1-(1-{e}^{-\lambda })0.5}\right],\end{eqnarray} \tag{ 19 }$$

where ${\xi }_{i}$ is a uniformly distributed computational coordinate $\xi \in [0,1]$ subdivided among the number of cells used in the x-direction. The degree of grid clustering at the origin is controlled by the stretching parameter λ. Likewise, the same form and parameters are used for the distribution of cells in y and z. This choice thus clusters cells near the origin (null point).

In our final simulations we take $\lambda =15.5$, and use 512³ computational cells (512 in each direction). This yields a maximum resolution of ${\rm{\Delta }}x\approx 0.0011$ close to the null point itself (i.e., this simulation, with the far boundaries, would require approximately ${(8\times {10}^{4})}^{3}$ cells on a uniform grid). The sufficiency of the resolution in the regions surrounding the null that contribute to the solution was tested in two ways. First, we ran quantitative convergence tests on the initial implosive phase (up to approximately t = 0.6) on a version of the simulation using uniform grids of successively finer resolutions (facilitated with much closer boundaries, as the reflection time is not prohibitive for such purposes). Using this approach, we were able to ensure that the collapse is adequately resolved and the solution converges. These simulations were also compared to the same stages of the final simulation using the stretched grid (showing good agreement). The resolution close to the null gives approximately 50 grid points across the width of the current sheet (diffusion region) itself, ensuring that the current is not concentrated at the grid scale, meaning that unphysical numerical diffusion is negligible. Our second approach was to perform simulations with lower resolution using the grid stretching (128³ and 256³ cells), which produced similar results. We note that below 256³ cells the solution in the wave-dominated region surrounding the null begins to deteriorate (e.g., the structure across the jets). Additionally, the ability of the stretched geometry to maintain the self-similar nature of implosions and explosions was tested by comparison of the stretched grid to uniform grid problems in runs with Sedov and MHD blasts, and with the successful reconstruction of the implosion/explosion of cylindrical MHD fast waves problem of McLaughlin et al. (2009), a 2D problem that originally required a vast uniform grid on a small stretched grid.

The final simulations required approximately five days of run time on 16 nodes of hyperthreaded 28 core 2.4 GHz Intel Xeon processors (896 threads), with high-speed omnipath interconnect between the nodes and 64 Gb RAM per node. The use of 512³ element arrays for a simulation of eight primitive variables (plus arrays for ancillary calculations) does not present a significant RAM overhead for modern nodes; the long computation time is instead due to the requirement of a large number of simulation cycles because of the physical parameters and grid, dictating very small hydrodynamic and resistive time steps.

Figure 4 was created as follows: (initially) magnetically connected fluid elements at either side of the flux tubes selected were tracked by writing a custom tracer module for Lare3d (i.e., fluid elements are followed in a Lagrangian sense as passive particles in the flow). This module is called after each hydrodynamic step of the main code and simply updates each particle position with an Eulurian push based on the hydrodynamic time-step and tri-linearly interpolated velocity at the particle's previous position (${{\boldsymbol{x}}}_{t+1}={{\boldsymbol{x}}}_{t}+{\rm{\Delta }}t\,{\boldsymbol{v}}$). The tracked particle positions are then used as instantaneous seed points for field line tracing using the inbuilt features of VisIt (Childs et al. 2012) at each time frame shown in Figure 4 and its supplementary movie. After the bifurcation, new seeds are traced from either end of the resulting four flux ropes, which are corresponding to fluid elements that are initially connected at t = 2.5. These seeds are used for tracing in the same fashion for $t\gt 2.5$. The particle module was tested exactly in linear flow fields. The tracers were also found to behave in the expected fashion for a number of test problems, including MHD wave propagation, the Kelvin–Helmholtz instability, and blast waves.