Three-dimensional oscillatory magnetic reconnection
TThree-dimensional oscillatory magnetic reconnection
J. O. Thurgood *1,2 , D. I. Pontin , and J. A. McLaughlin Department of Mathematics, Physics and Electrical Engineering, Northumbria University, UK. Division of Mathematics, University of Dundee, UK. [email protected] Here we detail the dynamic evolution of localised reconnection regions aboutthree-dimensional (3D) magnetic null points by using numerical simulation. Wedemonstrate for the first time that reconnection triggered by the localised col-lapse of a 3D null point due to an external MHD wave involves a self-generatedoscillation, whereby the current sheet and outflow jets undergo a reconnection re-versal process during which back-pressure formation at the jet heads acts to priseopen the collapsed field before overshooting the equilibrium into an opposite-polarity configuration. The discovery that reconnection at fully 3D nulls canproceed naturally in a time-dependent and periodic fashion is suggestive that os-cillatory reconnection mechanisms may play a role in explaining periodicity inastrophysical phenomena associated with magnetic reconnection, such as the ob-served quasi-periodicity of solar and stellar flare emission. Furthermore, we finda consequence of oscillatory reconnection is the generation of a plethora of freely-propagating MHD waves which escape the vicinity of the reconnection region
Manuscript in press, accepted for publication by ApJ in June 2017. The final published ver-sion will be available with ’gold’ open access, see the main journal for access to supplementaryanimations.
Magnetic fields play a key role in determining the dynamics of plasmas at all scales: fromfusion experiments and laboratory plasmas, to planetary magnetospheres, the Sun and stars, togalaxies and accretion disks. Magnetic reconnection is a fundamental plasma process associatedwith dynamic energy release in these systems, and is believed to explain a broad range ofphenomena including solar and stellar flares, coronal mass ejections, astrophysical jets andplanetary aurorae. Current theoretical frontiers in fundamental reconnection research can bebroadly categorised as; (i) collisionless/kinetic effects [Yamada et al., 2010], (ii) the extensionto three dimensions (3D) of a long history of two-dimensional (2D) theory (where profounddifferences 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 dynamicswith those of the global system. Here, we present significant results concerning the latter threewithin the framework of nonlinear magnetohydrodynamics (MHD). a r X i v : . [ a s t r o - ph . S R ] J un uch of our present understanding of reconnection has focused primarily on continuously-driven (and hence at least quasi-steady) systems, whilst transient effects that are crucial inmany plasma environments are typically neglected (we note an important exception in thecase of current sheet instabilities, which are outside of the scope of this paper, where studiesare necessarily time-dependent, e.g. Loureiro et al. [2007], Comisso et al. [2016]). Combinedwith questions regarding the influence of chosen boundary conditions, the applicability of anygiven reconnection model to a particular situation is often unclear. Here we study instead theself-consistent formation and evolution of a reconnecting current sheet in response to a finiteexternal 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 iswell 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 reconnectionmay take place [Pontin & Craig, 2005, Pontin et al., 2007, McLaughlin et al., 2009]. Whathas so far not been studied is the subsequent reconnection dynamics triggered by the wavesin a configuration where the reconnection region is connected naturally to an external, non-reconnecting region. That is the aim of this paper.Within the reduced framework of 2D MHD, where the assumption of an invariant direc-tion simplifies the system greatly, it has recently been demonstrated that in such a scenarioreconnection does not proceed in a quasi-steady manner but rather exhibits a time-dependentbehaviour, 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. 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´en waves becomes permissible.Here we demonstrate for the first time that impulsively initiated reconnection at realistic3D magnetic null points naturally proceeds in a time-dependent and oscillatory fashion, andalso demonstrate that 3D reconnection can itself produce MHD waves. The results show thatreconnection is inherently periodic and highlight the deeply-interconnected nature of MHDwaves 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 deferredto the appendix. All variables in this paper are nondimensionalised. We consider a backgroundmagnetic field of the form: B = [ x, y, − z ] , (1)which is known as the linear, proper potential null point [Parnell et al., 1996]. The magneticnull point itself ( B = ) is located at the origin. The field is free from electrical currents andso 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 fieldlinespointing 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 = ∇ × A where A = ψ exp (cid:34) − (cid:0) x + y + z (cid:1) σ (cid:35) ˆ y , (2) igure 1: (a) The equilibrium magnetic field B . Black fieldlines illustrate the behaviour of the spine and fanfieldlines, where the fan fieldlines point radial from the null in the z = 0 plane and the spine line which points towardsthe null along the z -axis. Red fieldlines are traced from above the fan, and blue from below. (b) The perturbingflux ring B , superimposed upon the background field. Coloured isosurfaces profile the increasing perturbation fieldstrength from zero (transparent) through weak (blue) to strong (red). The circulation of flux ring fieldlines in planesof 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 andoutward propagating magnetoacoustic waves. The incoming part is our intended driver forinitiating magnetic reconnection and acts implosively to collapse the magnetic field into thefirst reconnection configuration. The outgoing part of the perturbation leaves the initialisationsite and propagates outwards. Special care is taken so that it and any subsequently generatedwaves 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 = ) 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 is scale-free the relative sizes ofthese three parameters control the dynamics of the initial implosion. Here we present resultsfor ψ = 0 . σ = 0 . η = 10 − , and p = 0 . Initial implosion
The incoming wave is initially symmetric as per the geometry of the flux ring and generates acylindrical ring current flowing through the null. It propagates towards the null point as a fastmagnetoacoustic mode. In linearised MHD it would evolve in a self-similar manner, focusingits energy towards the null, increasing in amplitude. Due to the linearly-decreasing Alfv´enspeed profile ( c A ∝ | B | ) length scales decrease exponentially, and so current density at thenull point increases exponentially in time in the linear regime whilst retaining its cylindricalsymmetry [McLaughlin & Hood, 2004, McLaughlin et al., 2008]. This focusing would continue igure 2: Transparent isosurfaces illustrating 3D distribution of j y about the collapsed field close to the null at t = 0 . (a) and (b) t = 0 .
5. Note the changing colour scales. until either resistive diffusion or plasma back-pressure act to halt further collapse, the scalingsfor which are well-understood only in 2D [Craig & McClymont, 1991, Hassam, 1992, Craig &McClymont, 1993, Craig & Watson, 1992, McClymont & Craig, 1996]. However, given thatour simulation considers a full nonlinear solution, there exists the opportunity for increasingamplitudes to allow the excess flux carried by the wave to overwhelm the background field andso 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’sspine and fan. The anisotropic nature of the Lorentz force leads to a net force that is directedalternately towards and away from the null point, in the regions of magnetic enhancement andrarefaction respectively. It follows that the implosion carries fronts of magnetic and plasmacompression 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 demonstratedthat the implosion becomes quasi-1D as its fronts accelerate and stall in regions of compressionand 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 disturbanceflattening into a planar structure. Figure 2 shows selected transparent isosurfaces of the y -component of current at an early time ( t = 0 .
1) and at the approximate time at which thecollapse is halted ( t = 0 . y -component of current is that which is flowing throughthe 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 symmetryin planes of fixed y , and by t = 0 . y . It differs in that thecurrent distribution is however fully localised in 3D, as can be seen in Figures 2a and b. Aswe will see in the following section, there is a crucial difference in that the 3D collapse triggersa 3D reconnection mode, which is significantly different to 2D. We highlight that this sheet isdetached from the influence of the boundaries and has been generated in an impulsive rather igure 3: Longer term evolution of j y in the y = 0 plane, showing the oscillation of the sheet. Note the changingcolour scales. See also the supplementary animation. igure 4: The reconnection of two thin flux tubes which are initially symmetric about either side of the initialcurrent sheet is shown. The top row spans the first reconnection event (positive current at the null), where fieldlinesanchored above and below the null (i.e. at larger z ) flip around the spine, and fieldlines anchored closer to thefan plane (i.e at lower z ) flip up the spine, thus being transferred through the fan plane. This connectivity changeis spine-fan reconnection . By the completion of the first reversal the selected flux tubes have bifurcated into fourseparate tubes. The bottom row follows them through the second reversal. The exterior flux tubes do not participatein the second reconnection event and are so recoloured as a transparent grey for clarity. The interior flux ropes areagain reconnected in the spine-fan mode, but in the opposite direction to the initial reconnection due to the currentreversal. See also the corresponding supplementary animation. than driven manner – once the current sheet is formed, no further energy is injected to sustainthe 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 = 0plane) is shown in Figure 3. In the first few frames, we see the formation of a current sheetby t ≈ . deflection currents due to their close association with the formation ofa shock in the reconnection outflow (see Forbes 1986)– act to prise open the collapsed spineand fan fieldlines and so relieve the force imbalance about the null point which supports thepre-existing (positive) current sheet. The pre-existing current sheet thus shortens along whatwas its length-wise axis and broadens along its width-wise axis, with a decrease in the currentdensity at the null point itself. Eventually, this opposite polarity current completely reversesthe spine-fan collapse, and overshoots, causing a collapse of the spine and fan magnetic fieldlinestructures in the opposite direction. The deflection current concentrations thus eventuallycoalesce to form a new current sheet by around t ≈
2. The process then repeats, restoring thesheet 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 in-vestigate its nature by tracking the changes to fieldline connectivity within selected magnetic igure 5: The reconnection jet at t = 0 . (a) The system of flow in the y = 0 plane, where the colour map illustratesdensity, blue curves are magnetic fieldlines in the outflow regions, and red curves are fieldlines in the inflow regions(with the inflow and outflow regions separated by the spine and fan). Black velocity arrows indicate flow directionand magnitude. The dashed white contour encloses the region of super-magnetosonic flow ( M A > (b) Thesystem of flow out of the y = 0 plane, where density enhancement and velocity arrows are instead shown in the planedefined by the current sheet’s lengthwise axis. The transparent yellow surface encloses the super-magnetosonic flow( M A > flux tubes during the simulation (Figure 4). Initially, we see that the two selected fieldlinebundles are carried in towards the null by the wave-generated inflow. As they enter the vicinityof the current sheet (diffusion region), the pairwise magnetic connectivity between the fluidelements at ends of the ropes is lost as fieldlines are continuously reconnected across both thespine line and fan plane. In this fashion, the reconnection resulting from 3D null collapse isdistinct from the 2D case. Unlike in 2D reconnection, where connectivity changes occur al-ways in a pair-wise sense at the null point only, here connectivity is changed instant-to-instantthroughout the whole of the non-ideal region [Priest et al., 2003]. The observed connectivitychange is indicative of the spine-fan mode of 3D reconnection [Pontin et al., 2005, Priest &Pontin, 2009], where magnetic flux is redistributed across both the spine and fan plane. Bythe point at which the first reversal has completed ( t ≈
2, cf. Figures 3 and 4) the selectedflux ropes have been fully redistributed into opposite lobes about the spine and fan. As allmagnetic connectivity is lost between the four groups of seed particles at the ends of the se-lected flux ropes during the first reconnection event (see Appendix E), the selected flux ropesare bifurcated into four distinct ropes of flux by the completion of the first reversal. As timefurther progresses during the period of negative current at the null point, the inner-most ofthe selected flux ropes are reconnected through the diffusion region again, but in the oppositedirection due to the reversal of current. The figure therefore demonstrates that the transfer offlux 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 alsoplasma, with the establishment of a system of flow through the current sheet (Figure 5). It isnot a simple ‘disc-shaped’ outflow, relieving the compression of the implosion symmetrically, ut rather plasma is directed predominantly along the current sheet axis as this is the dominantdirection of the Lorentz force associated with sharply curved reconnected fieldlines. Along theoutflow plasma is accelerated to super-magnetosonic speeds in two reconnection jets . The jetsare bounded by standing slow mode MHD shock fronts which separate the fieldlines in theinflow 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 ofthe reconnection jets leading to back-pressures that choke off the outflow. This is accompaniedby the relief of (plasma) compression in the current sheet itself via the expulsion of the excessplasma 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 continualdriving flows. The consequence of this plasma and flux redistribution is that the balance offorces acting across the spine and fan planes (the collapse of which is necessary to sustain thereconnection 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 jetheads exerting a force against the outflow. In conjunction with deleterious effects in the inflowregion (such as the diffusion/redistribution of the driving magnetic field and current sheetdecompression) this leads to a deceleration of the jets, which leads in turn to the formationof fast MHD shocks (termination shocks). These refract magnetic field towards the shock’stangent with increases in overall plasma pressure and field strength. This bending of the fieldaway from the normal gives rise to the opposite polarity deflection currents which were earlieridentified in Figure 3. These fast, nonlinear waves formed in the jet head travel in toward thenull, and act implosively with similar physics to the initial collapse. Thus, a nonlinear fastwave 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 secondarycollapse which pushes the spine and fan together in the opposite sense. This reverses the sign ofcurrent at the null, and thus reverses the sense of the reconnection. With the completion of thereversal and the secondary implosion, we see the establishment of new (weaker) reconnectionflows 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. How-ever, the global effects of the magnetic field restructuring are not artificially confined to thisarea by computational boundaries, but rather may escape the vicinity of the current layer asfreely propagating waves. We find that many escaping MHD waves are generated and highlighta selection here.Figure 6a,b shows time-distance diagrams which trace disturbances seen in | v | and ρ prop-agating along the positive z -axis. A number of different wavefronts are visible. The featureslabelled I and O respectively correspond to the incoming and outgoing waves generated byour initial perturbation, i.e., they are not generated by reconnection.We identify the features marked A as Alfv´en waves. They propagate at the linearly in-creasing background Alfv´en speed, which leads to the exponential profile along its fronts (bythe solution of dz/dt = | B z | = 2 z ). They are, as characteristic of Alfv´en waves, incompressibleand are as such not visible in the figure for ρ (Figure 6b). These waves are also manifest in thegeneration of oppositely signed vorticity tubes in the reconnection region which propagate upalong either side of the spine line (Figure 7). Vorticity propagation is a marker of Alfv´en wavesin 3D MHD, as it corresponds to a rotational disturbances of the magnetic field which allows igure 6: Time distance diagrams showing the evolution of velocity magnitude | v | (a) and density ρ (b) along the z -axis (the approximate position of the spine away from the collapse). The intensity of both quantities is scaled asthe square-root to enhance contrast between different features. igure 7: Vorticity tube formation about the spine and in the fan plane over time ( ω · ˆ B ) (see also the supplementarymovie). The vorticity is associated with propagating Alfv´en waves. igure 8: Propagating density fronts in the y = 0 plane at selected times corresponding, which correspond to slowmagnetoacoustic pulses. See also the corresponding movie. agnetic tension to act as the predominant restoring force giving rise to the key properties ofthe wave [Goossens et al., 2011, Thurgood & McLaughlin, 2012, Shelyag et al., 2013]. Thesewaves are generated by the swaying of the spine line and neighbouring fieldlines from side toside along the x -axis during each reversal event, which produces two counter rotational vorticesin the flow pattern surrounding the spine (somewhat reminiscent of the ‘ m = 1’ kink wave incoronal loops). In this way the vorticity produced by the magnetoacoustic modes driving theswaying of the spine to the Alfv´en mode in a process of mode conversion .The features marked S are identified as large-amplitude, anharmonic slow-mode waves –due to their magnetoacoustic propagation speed (always lower than the Alfv´en speed) – andhave signal present in both the velocity and density diagrams. They are excited in the outflowregions of the reconnection jets (i.e., both the slow waves and termination shock have thesame source). Owing to their strong nonlinearity these anharmonic waves not only propagateplasma compression but actually transport material from the reconnection outflow away fromthe null and out along the spine axis (and fan plane). The profile of material carried awayfrom 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 wavespropagate 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 ofthe material transport.
We have demonstrated unambiguously that magnetic reconnection at a 3D null point, driven bya 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 currentsheets by MHD waves in a process of 3D null collapse, (ii) flux and plasma redistribution vialocalised 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 ofenergy release in high magnetic Reynolds number plasmas, demonstrating that reconnectiontriggered in an aperiodic manner both produces periodic reconnection reversals and periodi-cally excited propagating MHD waves due to a self-generated oscillation. Both the inherentperiodicity of 3D reconnection events and the generation of waves may play a role in explainingthe observed quasi-periodicity in solar and stellar flare emission [Nakariakov & Melnikov, 2009,Nakariakov et al., 2016, Van Doorsselaere et al., 2016]. Understanding the physical mechanismunderpinning this quasi-periodicity is crucial for developing the next generation of flare mod-els and for exploiting these observations as a diagnostic tool (via magneto-seismology, see DeMoortel & Nakariakov 2012). At present, there are a number of competing physical mecha-nisms that may explain Quasi-Periodic Pulsations (QPPs) including Oscillatory Reconnection,the magnetic tuning fork [Takasao & Shibata, 2016], self-oscillatory process such as magnetic-dripping [Nakariakov et al., 2010], and dispersive wave trains [Nakariakov et al., 2004] (seealso the aforementioned review papers). Our work finds that Oscillatory Reconnection mayproceed at fully 3D nulls and therefore is one potential mechanism that may provide an ex-planation for QPPs, since a subset of flaring events involve truly 3D null point geometries asopposed 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 correspondphysically to an MHD pulse impinging on the null which is sufficiently energetic (relative to he 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 OscillatoryReconnection. Pulses that are not sufficiently energetic, e.g. lower amplitudes (for the samepressure 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 qualita-tively different fashion to the system presented here. As our experiment is free of global scalesits 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 observa-tional and experimental data. We expect that, in the solar atmosphere, given the abundanceof 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 tobe accessible and so Oscillatory Reconnection should proceed in a qualitative sense as pre-sented here. Determining the specific periods accessible by the mechanism in a solar scenariowill 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 [Tomczyket al., 2007, McIntosh et al., 2011, Morton et al., 2015]. Understanding these waves is essentialsince 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 mass-transporting slow modes due to reconnection may explain the origin of propagating intensitydisturbances observed in the open-field corona [DeForest & Gurman, 1998, Banerjee et al.,2011]. Additionally, models predict that the dissipation of torsional Alfv´en waves could providesignificant chromospheric and coronal heating [Antolin & Shibata, 2010]. Thus, the naturalexcitation of torsional Alfv´en waves by Oscillatory Reconnection is especially interesting giventhat twist and torsional motions have recently been shown to be highly-prevalent in the lowsolar atmosphere [De Pontieu et al., 2014].There are also significant implications for the restructuring of the global field in response tothe reconnection process. For example, reconnection at a 3D null point in the Sun’s atmospherecan facilitate the exchange of plasma between magnetic fieldlines that are closed (connectat both ends to the solar surface) and those that are open to interplanetary space [Pontinet al., 2013], driving material outflows [Bradshaw et al., 2011]. Our results thus provide apotential 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-likereconnection jets in a spine-fan system, replete with standing slow shocks and fast terminationshocks, with particular implications for non-thermal particle acceleration. We note that the steady-state Petschek reconnection in uniform-resistivity MHD is widely considered unattain-able, with the requirement of an anomalous resistivity, localised at the null, to act as anobstacle 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 connectedto the time-dependent nature of our problem.Finally, we note that previous studies involving (2D) null collapse were chiefly motivatedin terms of examining the possibility of the collapse to small scales and large current densitiesyielding favourable scaling laws of the reconnection rate with resistivity. Early results inthe zero- β limit suggested a ‘fast’ scaling, with decay rates of the perturbations scaling as ∼ | ln η | 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 haltthe collapse before it enters a regime of fast reconnection. There are a number of possible esolutions - in MHD this includes the possibility of strong nonlinear effects associated withhighly energetic perturbations mitigating the effects of pressure leading to secondary, fasterphases 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 microphysicalreconnection rates will be sufficiently fast. In the context of 2D null point collapse, the fastrates of collisionless reconnection have been demonstrated by particle-in-cell simulations in aseries 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 theefficiency of 3D null collapse as a reconnection and dissipation mechanism. In this paper, 3Dnull collapse features as the means of establishing a localised current sheet which precipitatesOscillatory Reconnection, which itself is our main focus. We do however note that we havefound that three-dimensional null collapse is at least qualitatively similar to the known 2Dbehaviour. 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. Acknowledgements
The authors acknowledges generous support from the Leverhulme Trust and this work wasfunded by a Leverhulme Trust Research Project Grant: RPG-2015-075. The authors acknowl-edge IDL support provided by STFC. The computational work for this paper was carriedout on HPC facilities provided by the Faculty of Engineering and Environment, NorthumbriaUniversity, UK.
Author Information
The authors declare no competing financial interests or other conflicts of interest. Correspon-dence should be addressed to [email protected]
References
Antolin, P., & Shibata, K. 2010, ApJ, 712, 494Arber, T. D., Brady, C. S., & Shelyag, S. 2016, ApJ, 817, 94Arber, T. D., Longbottom, A. W., Gerrard, C. L., & Milne, A. M. 2001, Journal of Computa-tional Physics, 171, 151Banerjee, D., Gupta, G. R., & Teriaca, L. 2011, Space Sci. Rev., 158, 267Biskamp, D., & Schwarz, E. 2001, Physics of Plasmas, 8, 4729Bradshaw, S. J., Aulanier, G., & Del Zanna, G. 2011, ApJ, 743, 66Caramana, E. J., Shashkov, M. J., & Whalen, P. P. 1998, Journal of Computational Physics,144, 70Chandrashekhar, K., Morton, R. J., Banerjee, D., & Gupta, G. R. 2014, A&A, 562, A98 hilds, H., Brugger, E., Whitlock, B., et al. 2012, in High Performance Visualization–EnablingExtreme-Scale Scientific Insight, 357–372Comisso, L., Lingam, M., Huang, Y.-M., & Bhattacharjee, A. 2016, Physics of Plasmas, 23,100702Craig, I. J., & Watson, P. G. 1992, ApJ, 393, 385Craig, I. J. D., & McClymont, A. N. 1991, ApJ, 371, L41—. 1993, ApJ, 405, 207De Moortel, I., & Nakariakov, V. M. 2012, Philosophical Transactions of the Royal Society ofLondon Series A, 370, 3193De Pontieu, B., Rouppe van der Voort, L., McIntosh, S. W., et al. 2014, Science, 346, 1255732DeForest, C. E., & Gurman, J. B. 1998, ApJ, 501, L217Forbes, T. G. 1982, Journal of Plasma Physics, 27, 491—. 1986, ApJ, 305, 553Goossens, M., Erd´elyi, R., & Ruderman, M. S. 2011, Space Sci. Rev., 158, 289Graf von der Pahlen, J., & Tsiklauri, D. 2014, Physics of Plasmas, 21, 012901—. 2016, A&A, 595, A84Gruszecki, M., Vasheghani Farahani, S., Nakariakov, V. M., & Arber, T. D. 2011, A&A, 531,A63Hassam, A. B. 1992, ApJ, 399, 159Leake, J. E., Linton, M. G., & Antiochos, S. K. 2014, ApJ, 787, 46Leake, J. E., Linton, M. G., & T¨or¨ok, T. 2013, ApJ, 778, 99Long, D. M., Bloomfield, D. S., Chen, P. F., et al. 2017, Sol. Phys., 292, 7Loureiro, N. F., Schekochihin, A. A., & Cowley, S. C. 2007, Physics of Plasmas, 14, 100703Masson, S., Pariat, E., Aulanier, G., & Schrijver, C. J. 2009, ApJ, 700, 559McClymont, A. N., & Craig, I. J. D. 1996, ApJ, 466, 487McIntosh, S. W., de Pontieu, B., Carlsson, M., et al. 2011, Nature, 477McLaughlin, J. A., De Moortel, I., Hood, A. W., & Brady, C. S. 2009, A&A, 493, 227McLaughlin, J. A., Ferguson, J. S. L., & Hood, A. W. 2008, Sol. Phys., 251, 563McLaughlin, J. A., & Hood, A. W. 2004, A&A, 420, 1129McLaughlin, J. A., Hood, A. W., & de Moortel, I. 2011, Space Sci. Rev., 158, 205 cLaughlin, J. A., Verth, G., Fedun, V., & Erd´elyi, R. 2012, ApJ, 749, 30Morton, R. J., Srivastava, A. K., & Erd´elyi, R. 2012, A&A, 542, A70Morton, R. J., Tomczyk, S., & Pinto, R. 2015, Nature Communications, 6, 7813Murray, M. J., van Driel-Gesztelyi, L., & Baker, D. 2009, A&A, 494, 329Nakariakov, V. M., Arber, T. D., Ault, C. E., et al. 2004, MNRAS, 349, 705Nakariakov, V. M., Inglis, A. R., Zimovets, I. V., et al. 2010, Plasma Physics and ControlledFusion, 52, 124009Nakariakov, V. M., & Melnikov, V. F. 2009, Space Sci. Rev., 149, 119Nakariakov, V. M., Pilipenko, V., Heilig, B., et al. 2016, Space Sci. Rev., 200, 75Ofman, L., Morrison, P. J., & Steinolfson, R. S. 1993, ApJ, 417, 748Parnell, C. E., Smith, J. M., Neukirch, T., & Priest, E. R. 1996, Physics of Plasmas, 3, 759Petschek, H. E. 1964, NASA Special Publication, 50, 425Pontin, D. I. 2012, Philosophical Transactions of the Royal Society of London A: Mathematical,Physical and Engineering Sciences, 370, 3169Pontin, D. I., Bhattacharjee, A., & Galsgaard, K. 2007, Phys. Plasmas, 14, 052106Pontin, D. I., & Craig, I. J. D. 2005, Phys. Plasmas, 12, 072112Pontin, D. I., Hornig, G., & Priest, E. R. 2005, Geophysical and Astrophysical Fluid Dynamics,99, 77Pontin, D. I., Priest, E. R., & Galsgaard, K. 2013, Astrophys. J., 774, 154Priest, E., & Forbes, T. 2000, Magnetic Reconnection, 612Priest, E. R., Hornig, G., & Pontin, D. I. 2003, J. Geophys. Res., 108, 1285Priest, E. R., & Pontin, D. I. 2009, Physics of Plasmas, 16, 122101Pucci, F., Onofri, M., & Malara, F. 2014, ApJ, 796, 43Roberts, G. O. 1971, in Lecture Notes in Physics, Berlin Springer Verlag, Vol. 8, NumericalMethods in Fluid DynamicsNumerical Methods in Fluid Dynamics, ed. M. Holt, 171–177Shelyag, S., Cally, P. S., Reid, A., & Mathioudakis, M. 2013, ApJ, 776, L4Takasao, S., & Shibata, K. 2016, ApJ, 823, 150Thurgood, J. O., & McLaughlin, J. A. 2012, A&A, 545, A9Tomczyk, S., McIntosh, S. W., Keil, S. L., et al. 2007, Science, 317, 1192Tsiklauri, D. 2008, Physics of Plasmas, 15, 112903 siklauri, D., & Haruki, T. 2007, Physics of Plasmas, 14, 112905—. 2008, Physics of Plasmas, 15, 102902Van Doorsselaere, T., Kupriyanova, E. G., & Yuan, D. 2016, Sol. Phys., arXiv:1609.02689Yamada, M., Kulsrud, R., & Ji, H. 2010, Reviews of Modern Physics, 82, 603Zhang, Q. M., Chen, P. F., Guo, Y., Fang, C., & Ding, M. D. 2012, ApJ, 746, 19Zhang, Q. M., & Ji, H. S. 2014, A&A, 567, A11Zhang, Q. M., Ji, H. S., & Su, Y. N. 2016, Sol. Phys., 291, 859Zuccarello, F., Romano, P., Farnik, F., et al. 2009, A&A, 493, 629 A Solver - Lare3d code
The simulation is the numerical solution of the nondimensional, resistive MHD equations:D ρ D t = − ρ ∇ · v (3)D v Dt = 1 ρ ( ∇ × B ) × B − ρ ∇ p + F shock (4)D B D t = ( B · ∇ ) v − B ( ∇ · v ) − ∇ × ( η ∇ × B ) (5)D (cid:15) D t = − pρ ∇ · v + ηρ j + H visc ρ (6) j = ∇ × B (7) E = − v × B + η j (8) p = (cid:15)ρ ( γ −
1) (9)which are solved using the
Lare3d code. All results presented are in non-dimensional units.Algorithmically, the code solves the ideal MHD equations explicitly using the Lagrangianremap 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-centred artificial viscosityapproach [Caramana et al., 1998], where shock viscosity is applied to the momentum equationthrough F shock and heats the system through H visc . Extended MHD options available withinthe code, such as the inclusion of Hall terms, were not used in these experiments. Full detailsof 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. Initial Conditions
A summary of the (nondimensional) initial conditions on the primitive variables, discussed inthe main document, is as follows: p = 0 .
005 (10) ρ = 1 (11) v = (12) B = B + ∇ × A (13) B = [ x, y, − z ] (14) A = ψ exp (cid:34) − (cid:0) x + y + z (cid:1) σ (cid:35) ˆ y (15) ψ = 0 .
05 (16) σ = 0 .
21 (17)with resistivity taken uniformly as η = 10 − throughout. The ratio of specific heats is takenas γ = 5 / C Boundary conditions and reflectivity testing
The boundary conditions are taken as no-slip ( v = ) with zero-gradient conditions on B , ρ and p . A simulation was run with no perturbation to check that the boundaries do not launchspurious 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 perturbationor those subsequently generated, do not reflect off the domain boundaries and return inwardsto influence our solution at the null. During our initial testing we found that open boundariesand artificial damping zones are imperfect strategies for removing outgoing waves in linearnull point geometries. As such the only robust method of avoiding reflections influencing thesolution is to simply place the boundaries sufficiently far away from the null point that theshortest possible signal travel time from the initialisation site to the boundary and back is inexcess of the period we wish to study. Computationally, this is not a trivial task. The shortestsignal travel time is governed by the Alfv´en speed which increases linearly with distance fromthe null point (at the center of the domain). Thus, the overall travel time to the boundaryincreases only logarithmically with increasing distances to the boundaries. Boundaries mustbe placed at vast distances whilst maintaining sufficient resolution in the dynamic and diffusiveregions close to the null point.For the results shown here, we placed the x , y and z boundaries at ±
40. The shortestpossible signal travel time permitted is that of a wave propagating at the Alfv´en speed directlyalong the spine (where B z increases as 2 z ), and so by separation of variables the travel timefrom the outer edges of the initial perturbation z ≈ ± .
25 out to the boundary is t out ≈ . / .
25) = 2 .
53. The fastest possible reflection will return to the diffusion dominatedregion about the null (estimated at the length of the initial current sheet, ≈ .
18) after afurther t in ≈ . / .
18) = 2 .
70, yielding the shortest possible theoretical reflection timeof t reflect ≈ .
23. This was also practically tested with a ‘boundary convergence test’ wheresimulations were run with boundaries which were closer to the null (hence having shorterreflection times). We found that results only differed at times in excess of the theoreticalreflection time calculated as above for each case. In these test cases it was often possible to track eflections along paths using time distance diagrams as per Figure 6a, which further confirmedthese theoretical times. The maximum boundary length used ( ±
40) was a compromise betweenmaximising the reflection time and maintaining the resolution available close to the null – notethat doubling the distance to the boundary will only increase the time by ln(2). Thus, bysimple wave theory, the fastest possible theoretical reflection time is t reflect ≈ .
23 and so thenon-interference of reflections during the initial implosion, the first reversal, and the initialstages of the second reversal is guaranteed. Therefore, the self-generated oscillations of thenull point are indeed due to reconnection itself. Further, we note in the later stages of thesecond reconnection reversal ( t > t reflect ) up to the end time of t = 6 that no incoming frontsof significant intensity are detected in Figure 6a, which is in fact taken along this shortestpossible signal travel path. D Grid geometry and convergence testing
To accommodate the vast distances to the boundary whilst sufficiently resolving the diffusionat the null point and the wave dynamics in the surrounding region a grid stretching scheme isemployed.The grid is stretched according to a variation on the scheme of Roberts [1971], namely thecell boundary positions x b along the x -direction are distributed according the transformation: x b = 40 (cid:26) λ ( ξ i − Γ)]sinh [ λ Γ] (cid:27) −
40 (18)Γ = 12 λ ln (cid:34) (cid:0) e λ − (cid:1) . − (1 − e − λ ) 0 . (cid:35) (19)where ξ i is a uniformly distributed computational coordinate ξ ∈ [0 ,
1] subdivided amongstthe number of cells used in the x direction. The degree of grid clustering at the origin iscontrolled by the stretching parameter λ . Likewise, the same form and parameters are usedfor the distribution of cells in y and z . This choice thus clusters cells near the origin (nullpoint).In our final simulations we take λ = 15 .
5, and use 512 computational cells (512 in eachdirection). This yields maximal resolution of ∆ x ≈ . × ) cells on auniform grid). The sufficiency of the resolution in the regions surrounding the null whichcontribute to the solution was tested in two ways. First, we ran quantitative convergence testson the initial implosive phase (up to approximately t = 0 .
6) on a version of the simulationusing 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 wereable to ensure that the collapse is adequately resolved, and the solution converges. Thesesimulations were also compared to the same stages of the final simulation using the stretchedgrid (showing good agreement). The resolution close to the null gives approximately 50 gridpoints across the width of the current sheet (diffusion region) itself, ensuring the current is notconcentrated at the grid scale, meaning that unphysical numerical diffusion is negligible. Oursecond approach was to perform lower resolution simulations using the grid stretching (128 and 256 cells), which produced similar results. We note that below 256 cells the solutionin the wave dominated region surrounding the null begins to deteriorate (e.g., the structureacross the jets). Additionally, the ability of the stretched geometry to maintain the self-similarnature of implosions and explosions was tested by comparison of the stretched grid to uniform rid problems in runs with Sedov and MHD blasts, and with the successful reconstruction ofthe 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 hyper-threaded 28 core 2 . element arrays fora simulation of 8 primitive variables (plus arrays for ancillary calculations) does not present asignificant RAM overhead for modern nodes, rather the large computation time is due to therequirement of a large number of simulation cycles due to the physical parameters and griddictating very small hydrodynamic and resistive timesteps. E Tracking connectivity change with tracer parti-cles
Figure 4 was created as follows: (initially) magnetically connected fluid elements at eitherside 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). Thismodule is called after each hydrodynamic step of the main code and simply updates eachparticle position with an Eulurian push based on the hydrodynamic timestep and tri-linearlyinterpolated velocity at the particle’s previous position ( x t+1 = x t + ∆t v ). The trackedparticle positions are used then as instantaneous seed points for fieldline tracing using theinbuilt features of VisIt [Childs et al., 2012] at each time frame shown in Figure 4 and itssupplementary movie. After the bifurcation, new seeds are traced from either end of theresultant four flux ropes which are corresponding to fluid elements which are initially connectedat t = 2 .
5. These seeds are used for tracing in the same fashion for t > .
5. The particle modulewas tested exactly in linear flow fields. The tracers were also found to behave in the expectedfashion for a number of test problems including MHD wave propagation, the Kelvin-Helmholtzinstability, and blast waves.5. The particle modulewas tested exactly in linear flow fields. The tracers were also found to behave in the expectedfashion for a number of test problems including MHD wave propagation, the Kelvin-Helmholtzinstability, and blast waves.