Three-dimensional oscillatory magnetic reconnection

Here we detail the dynamic evolution of localised reconnection regions about three-dimensional (3D) magnetic null points by using numerical simulation. We demonstrate for the first time that reconnection triggered by the localised collapse of a 3D null point due to an external MHD wave involves a self-generated oscillation, whereby the current sheet and outflow jets undergo a reconnection reversal process during which back-pressure formation at the jet heads acts to prise open the collapsed field before overshooting the equilibrium into an opposite-polarity configuration. The discovery that reconnection at fully 3D nulls can proceed naturally in a time-dependent and periodic fashion is suggestive that oscillatory reconnection mechanisms may play a role in explaining periodicity in astrophysical phenomena associated with magnetic reconnection, such as the observed quasi-periodicity of solar and stellar flare emission. Furthermore, we find a consequence of oscillatory reconnection is the generation of a plethora of freely-propagating MHD waves which escape the vicinity of the reconnection region


Introduction
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, to 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 categorised as; (i) collisionless/kinetic effects [Yamada et al., 2010], (ii) the extension to three dimensions (3D) of a long history of two-dimensional (2D) theory (where profound differences make this extension highly non-trivial [Priest et al., 2003, Pontin, 2012), (iii) timedependent 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 continuouslydriven (and hence at least quasi-steady) systems, whilst 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 outside of 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 due to 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. What has so far not been studied is the subsequent reconnection dynamics triggered by the waves in a configuration where the reconnection region is connected naturally to an external, nonreconnecting region. That 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 behaviour, termed Oscillatory Reconnection [Craig & McClymont, 1991, Ofman et al., 1993, McLaughlin et al., 2009, Murray et al., 2009, Pucci et al., 2014. However, the extension to 3D is not trivial; the process by which waves trigger reconnection (null collapse, see e.g. Chapter 7 of Priest & Forbes 2000) is little-studied in three-dimensions, 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 also demonstrate that 3D reconnection can itself produce MHD waves. The results show that reconnection is inherently periodic and highlight the deeply-interconnected nature of MHD waves and reconnection.

Simulation Setup
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 nondimensionalised. We consider a background magnetic field of the form: which is known as the linear, proper potential null point [Parnell et al., 1996]. The magnetic null point itself (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 fieldline, running along the z-axis towards the null point, and a fan plane z = 0, consisting of fieldlines pointing radially outwards. Other fieldlines, separated by the fan, have a hyperbolic structure ( Figure 1a). Upon this magnetic field we impose a finite-amplitude perturbation B 1 = ∇ × A 1 where Figure 1: (a) The equilibrium magnetic field B 0 . Black fieldlines illustrate the behaviour of the spine and fan fieldlines, where the fan fieldlines point radial from the null in the z = 0 plane and the spine line which points towards the null along the z-axis. Red fieldlines are traced from above the fan, and blue from below. (b) The perturbing flux ring B 1 , superimposed upon the background field. Coloured isosurfaces profile the increasing perturbation field strength from zero (transparent) through weak (blue) to strong (red). The circulation of flux ring fieldlines in planes of fixed y is illustrated by the yellow line.
which is a ring of magnetic flux centred about the null point, circulating in planes of fixed y (Figure 1b). 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 initialisation site and propagates outwards. Special care is taken so that it and any subsequently generated waves do not reflect at the boundaries and return inwards to influence our solution at the null (see the appendix). The background plasma is taken to be initially stationary (v = 0) and of uniform density (ρ = 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 B 0 is scale-free the relative sizes of these three parameters control the dynamics of the initial implosion. Here we present results for ψ = 0.05, σ = 0.21, η = 10 −3 , and p = 0.005. This choice of parameters is such that the implosion forms a localised current sheet which is isolated from the boundaries and free to evolve in a self-consistent manner.

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 towards the null point as a fast magnetoacoustic mode. In linearised MHD it would evolve in a self-similar manner, focusing its energy towards the null, increasing in amplitude. Due to the linearly-decreasing Alfvén speed profile (c A ∝ |B 0 |) length scales decrease exponentially, and so current density at the null point increases exponentially in time in the linear regime whilst 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, Hassam, 1992, Craig & McClymont, 1993, Craig & Watson, 1992, McClymont & Craig, 1996. However, given that our simulation considers a full nonlinear solution, there exists the opportunity for increasing 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 towards 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. Past 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 towards 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 ycomponent of 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 current is that which is flowing 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 which is localised 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 however fully localised in 3D, as can be seen in Figures 2a and b. As we will see 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  The reconnection of two thin flux tubes which are initially symmetric about either side of the initial current sheet is shown. The top row spans the first reconnection event (positive current at the null), where fieldlines anchored above and below the null (i.e. at larger z) flip around the spine, and fieldlines anchored closer to the fan plane (i.e at lower z) flip up the spine, thus being transferred through the fan plane. This connectivity change is spine-fan reconnection. By the completion of the first reversal the selected flux tubes have bifurcated into four separate tubes. The bottom row follows them through the second reversal. The exterior flux tubes do not participate in the second reconnection event and are so recoloured as a transparent grey for clarity. The interior flux ropes are again reconnected in the spine-fan mode, but in the opposite direction to the initial reconnection due to the current reversal. See also the corresponding supplementary animation. than driven manner -once the current sheet is formed, no further energy is injected to sustain the sheet or otherwise alter it.

Reconnection reversals
The longer term evolution of the current flowing through the plane of collapse (j y in the y = 0 plane) is shown in Figure 3. In the first few frames, we see the formation of a current sheet by t ≈ 0.6 due to the nonlinear implosion collapsing the null point's spine and fan towards each other (note the spine has locally 'collapsed' onto the fan in Figure 2b). Shortly after the initial implosion halts, strong opposite (negative) polarity currents form at the ends of the current sheet and begin to propagate inwards towards the null point. These current concentrations -called deflection currents due to their close association with the formation of a shock in the reconnection outflow (see Forbes 1986)-act to prise open the collapsed spine and fan fieldlines and so relieve the force imbalance about the null point which supports the pre-existing (positive) current sheet. The pre-existing current sheet thus shortens along what was its length-wise axis and broadens along its width-wise axis, with a decrease in the current density at the null point itself. Eventually, this opposite polarity current completely reverses the spine-fan collapse, and overshoots, causing a collapse of the spine and fan magnetic fieldline structures in the opposite direction. The deflection current concentrations thus eventually coalesce to form a new current sheet by around t ≈ 2. The process then repeats, restoring the sheet to its original polarity before the simulation end time of t = 6. The dynamics of these reconnection reversals are considered later.
We confirm magnetic reconnection occurs across these oscillating current sheets and investigate its nature by tracking the changes to fieldline connectivity within selected magnetic flux tubes during the simulation (Figure 4). Initially, we see that the two selected fieldline bundles are carried in towards the null by the wave-generated inflow. As they enter the vicinity of the current sheet (diffusion region), the pairwise magnetic connectivity between the fluid elements at ends of the ropes is lost as fieldlines are continuously reconnected across both the spine line and fan plane. In this fashion, the reconnection resulting from 3D null collapse is distinct from the 2D case. Unlike in 2D reconnection, where connectivity changes occur always in a pair-wise sense at the null point only, here connectivity is changed instant-to-instant throughout the whole of the non-ideal region [Priest et al., 2003]. The observed connectivity change is indicative of the spine-fan mode of 3D reconnection , Priest & Pontin, 2009, where magnetic flux is redistributed across both the spine and fan plane. By the point at which the first reversal has completed (t ≈ 2, cf. Figures 3 and 4) the selected flux ropes have been fully redistributed into opposite lobes about the spine and fan. As all magnetic connectivity is lost between the four groups of seed particles at the ends of the selected flux ropes during the first reconnection event (see Appendix E), the selected flux ropes are bifurcated into four distinct ropes of flux by the completion of the first reversal. As time further progresses during the period of negative current at the null point, the inner-most of the selected flux ropes are reconnected through the diffusion region again, but in the opposite direction due to the reversal of current. The figure therefore demonstrates that the transfer of flux about the spine and fan due to spine-fan reconnection is oscillatory and time-dependent.

Reconnection jets and reversal dynamics
The spine-fan reconnection not only redistributes 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 'disc-shaped' outflow, relieving the compression of the implosion symmetrically, but rather plasma is directed predominantly along the current sheet axis as this is the dominant direction of the Lorentz force associated with sharply curved reconnected fieldlines. 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 which separate the fieldlines in the inflow regions and the outflow regions and are analogous to those found in Petschek's 2D model [Petschek, 1964].
As the outflow evolves, we observe the pooling of hot, over-dense 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 leads in turn to the formation of fast MHD shocks (termination shocks). These refract magnetic field towards the shock's 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 which 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 which 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 which drive the next reversal via the same processes.

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. Figure 6a,b shows time-distance diagrams which trace disturbances seen in |v| and ρ propagating along the positive z-axis. A number of different wavefronts are visible. The features labelled I and O respectively correspond to the incoming and outgoing waves generated by our initial perturbation, i.e., they are not generated by reconnection.
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 characteristic of Alfvén waves, incompressible and are as such not visible in the figure for ρ (Figure 6b). These waves are also manifest in the generation of oppositely signed vorticity tubes in the reconnection region which 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 which 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 neighbouring fieldlines 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.
The features marked S are identified as large-amplitude, anharmonic slow-mode wavesdue 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., both the slow waves and termination shock have the same source). Owing to their strong nonlinearity these anharmonic waves not only propagate plasma compression but actually 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 ≈ 0) and fan (z ≈ 0) (which indicates these waves propagate in a field-aligned sense, as expected of slow waves -see also the corresponding movie), and the under-dense region that remains in the vicinity of the null point as a consequence of the material transport.

Discussion
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 localised spine-fan reconnection, (iii) associated reconnection reversals due to back-pressuredriven 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 both produces periodic reconnection reversals and periodically excited propagating MHD waves due to a self-generated oscillation. Both the inherent periodicity of 3D reconnection events and the generation of waves may 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 magneto-seismology, 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 process such as magneticdripping [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 involve truly 3D null point geometries as opposed to quasi-2D X-lines (such as seen in e.g. Zuccarello et al. 2009, Masson 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 which is sufficiently energetic (relative to the plasma pressure and resistivity) that it implodes in a nonlinear fashion as per McClymont & Craig [1996] and forms a current sheet which then goes on 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 modelling 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., 2013Leake et al., , 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 such models which place observational restrictions on permissible null geometries, driving energies, and plasma pressures, whilst 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; however the nature of their genesis is an outstanding question. The generation of masstransporting slow modes 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 fieldlines 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 non-thermal 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, localised 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 large current densities yielding favourable 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 ∼ | ln η| 2 in the linear regime and becoming virtually independent of η in the nonlinear case (so-called super-fast reconnection as per 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 leading to secondary, faster phases of reconnection, as hypothesised by McClymont & Craig [1996] and discussed in Priest & Forbes [2000]. There also exists the possibility 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. 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 localised current sheet which precipitates Oscillatory Reconnection, which itself is our main focus. We do however note that we have found that three-dimensional null collapse is at least qualitatively similar to the known 2D behaviour. A detailed, quantitative study of the scaling of the 3D null collapse mechanism, and it's relation to the 2D case, will be important to undertake in future.

B Initial Conditions
A summary of the (nondimensional) initial conditions on the primitive variables, discussed in the main document, is as follows: with resistivity taken uniformly as η = 10 −3 throughout. The ratio of specific heats is taken as γ = 5/3.

C Boundary conditions and reflectivity testing
The boundary conditions are taken as no-slip (v = 0) with zero-gradient conditions on 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 inwards 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. As such the only robust method of avoiding reflections influencing the solution is to simply place the boundaries sufficiently far away from the null point that the shortest possible signal travel time from the initialisation 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 whilst 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 ≈ ±0.25 out to the boundary is t out ≈ 0.5 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, ≈ 0.18) after a further t in ≈ 0.5 ln(40/0.18) = 2.70, yielding the shortest possible theoretical reflection time of t reflect ≈ 5.23. This was also practically tested with a 'boundary convergence test' where simulations were run with boundaries which 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 per Figure 6a, which further confirmed these theoretical times. The maximum boundary length used (±40) was a compromise between maximising 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 ln(2). Thus, by simple wave theory, the fastest possible theoretical reflection time is t reflect ≈ 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. Further, we note in the later stages of the second reconnection reversal (t > t reflect ) up to the end time of t = 6 that no incoming fronts of significant intensity are detected in Figure 6a, which is in fact taken along this shortest possible signal travel path.

D Grid geometry and convergence testing
To accommodate the vast distances to the boundary whilst 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: where ξ i is a uniformly distributed computational coordinate ξ ∈ [0, 1] subdivided amongst 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 λ = 15.5, and use 512 3 computational cells (512 in each direction). This yields maximal resolution of ∆x ≈ 0.0011 close to the null point itself (i.e., this simulation, with the far boundaries, would require approximately (8 × 10 4 ) 3 cells on a uniform grid). The sufficiency of the resolution in the regions surrounding the null which 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 the current is not concentrated at the grid scale, meaning that unphysical numerical diffusion is negligible. Our second approach was to perform lower resolution simulations using the grid stretching (128 3 and 256 3 cells), which produced similar results. We note that below 256 3 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 which originally required a vast uniform grid, on a small stretched grid.
The final simulations required approximately 5 days of run time on 16 nodes of hyperthreaded 28 core 2.4GHz Intel Xeon processors (896 threads), with high-speed omnipath interconnect between the nodes and 64Gb RAM per node. The use of 512 3 element arrays for a simulation of 8 primitive variables (plus arrays for ancillary calculations) does not present a significant RAM overhead for modern nodes, rather the large computation time is due to the requirement of a large number of simulation cycles due to the physical parameters and grid dictating very small hydrodynamic and resistive timesteps.
E Tracking connectivity change with tracer particles 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 timestep and tri-linearly interpolated velocity at the particle's previous position (x t+1 = x t + ∆t v). The tracked particle positions are used then as instantaneous seed points for fieldline 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 resultant four flux ropes which are corresponding to fluid elements which are initially connected at t = 2.5. These seeds are used for tracing in the same fashion for t > 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.