Implosive collapse about magnetic null points: A quantitative comparison between 2D and 3D nulls

Null collapse is an implosive process whereby MHD waves focus their energy in the vicinity of a null point, forming a current sheet and initiating magnetic reconnection. We consider, for the first time, the case of collapsing 3D magnetic null points in nonlinear, resistive MHD using numerical simulation, exploring key physical aspects of the system as well as performing a detailed parameter study. We find that within a particular plane containing the 3D null, the plasma and current density enhancements resulting from the collapse are quantitatively and qualitatively as per the 2D case in both the linear and nonlinear collapse regimes. However, the scaling with resistivity of the 3D reconnection rate - which is a global quantity - is found to be less favourable when the magnetic null point is more rotationally symmetric, due to the action of increased magnetic back-pressure. Furthermore, we find that with increasing ambient plasma pressure the collapse can be throttled, as is the case for 2D nulls. We discuss this pressure-limiting in the context of fast reconnection in the solar atmosphere and suggest mechanisms by which it may be overcome. We also discuss the implications of the results in the context of null collapse as a trigger mechanism of Oscillatory Reconnection, a time-dependent reconnection mechanism, and also within the wider subject of wave-null point interactions. We conclude that, in general, increasingly rotationally-asymmetric nulls will be more favourable in terms of magnetic energy release via null collapse than their more symmetric counterparts.


Introduction
Magnetic reconnection is an important process for energy conversion throughout astrophysical plasmas, in the Sun and planetary magnetospheres as well as further afield -for example in γ-ray bursts [e.g. Zweibel andYamada, 2009, Pontin, 2012]. One particular location in which reconnection can preferentially take place is in the vicinity of a magnetic null point (locations at which the magnetic field strength is zero). Recent studies show that such null points exist in abundance in the solar atmosphere [Régnier et al., 2008, Longcope andParnell, 2009], and there is growing evidence that reconnection at nulls plays an important role in solar flares, CMEs, jets and bright points [e.g. Barnes, 2007, Sun et al., 2013, Moreno-Insertis and Galsgaard, 2013, Zhang et al., 2016, Wyper et al., 2017, Chitta et al., 2017. This is further supported by observations of flare ribbons that appear to result from particle acceleration during null point reconnection [Zuccarello et al., 2009, Liu et al., 2011.
One particular mechanism involving reconnection at nulls that may be important for energy conversion is 'null point collapse': an implosive process whereby MHD waves form large current densities by concentrating flux at small scales. If the ambient plasma pressure is sufficiently low, the collapse is thought to be halted at a scale determined by resistive diffusion, yielding favourable scalings for reconnection rates with decreasing plasma resistivity Watson, 1992, McClymont andCraig, 1996]. The basic idea that perturbations tend to collect at Xpoints (or X-type neutral lines) leading to a growth in current densities was first realised by Dungey [1953]. Null point collapse in 2D and 3D has subsequently been considered from various perspectives. Dynamic collapse studies in the close vicinity of the null (within which the magnetic field and flow can be approximated as linear) have been performed in both the linear and nonlinear MHD regimes by Imshennik and Syrovatskii [1967], Bulanov and Olshanetsky [1984], Klapper et al. [1996], Mellor et al. [2003]. These studies tend to indicate unbounded growth of the current in the absence of dissipation. However, since they explicitly exclude the surrounding field, it is unclear whether sufficient energy could accumulate at the null in the full system to sustain this current blowup.
In this paper we take a different approach, which is to study the process computationally including the full nonlinear field and plasma flow geometries. This approach received significant attention during the 1990s in advocating null collapse as a possible mechanism for obtaining fast reconnection rates in two dimensions [Craig and McClymont, 1991, Hassam, 1992, Craig and Watson, 1992, Craig and McClymont, 1993, McClymont and Craig, 1996, Priest and Forbes, 2000, but is yet to be considered in 3D. In the 2D case, it was eventually realised that for the solar corona the plasma pressure would be likely sufficient to restrict the process and so questions were raised over its viability as a fast reconnection mechanism, at least from a solar physics perspective [e.g. Craig, 1996, Priest andForbes, 2000]. A possible resolution is hinted at in the work of McClymont and Craig [1996] in that nonlinear effects may overcome this limitation, permitting a secondary stage of fast reconnection after the initial halting, although to our knowledge this possibility has not since been investigated. More recently, PIC simulations of 2D null collapse find fast rates occur due to collisionless effects Haruki, 2007, 2008], and that null collapse seems to be able to efficiently accelerate particles which has been proposed as a source of γ-ray flares in the Crab Nebula [Lyutikov et al., 2016]. However, since these simulations begin with a magnetic field in which the null is 'pre-collapsed' at a kinetic scale, there remains the question of whether an external perturbation would initiate a collapse that forms a thin-enough current sheet to promote collisionless reconnection before being pressure-limited at an MHD-scale.
An alternative perspective on null collapse is obtained by making arguments based on lowest energy states in a domain in which the magnetic field lines are line-tied at the boundaries.
Starting from an equilibrium, one perturbs the magnetic field in such a way as to displace some of the separatrix (2D) or spine and fan (3D) field lines at the boundaries, and then considers the properties of the new equilibrium that is accessible through an ideal dynamics with these boundaries held fixed. Studies using different analytical and computational approaches indicate that the lowest energy state contains a singular current layer at the null in both 2D and 3D [Syrovatskii, 1971, Craig and Litvinenko, 2005, Pontin and Craig, 2005, Fuentes-Fernández and Parnell, 2012, Fuentes-Fernández and Parnell, 2013. Thus, in any real plasma with finite magnetic Reynolds number, resistive effects must eventually become important during the collapse, making 2D and 3D nulls favourable sites for reconnection. Notably, pressure forces are shown to weaken the scaling of the singularity, but cannot remove it Litvinenko, 2005, Pontin and.
In a recent study [Thurgood et al., 2017], we used MHD pulses to initialise a localised and finite-duration collapse of 3D null points, to form a current sheet embedded within a largerscale (global) field and discovered that a phenomenon known as Oscillatory Reconnection can occur at 3D nulls which are undergoing spine-fan reconnection . We noted that a number of qualitative aspects of dynamic null collapse as described in the papers above appeared to carry-over to the case of 3D null collapse, but were at the time unable to quantitatively investigate the collapse phenomena by way of parameter study due to restrictive computational requirements associated with those particular simulations. Thus, in this paper it is our primary aim to present the results of such a study and investigate the extension of known null collapse scalings with resistivity, and associated magnetic reconnection efficiency, to 3D for the first time. We also aim to present a fresh perspective on both 2D and 3D null collapse which suggests that it may in-fact be a phenomenon of some importance even in the case of relatively 'high' coronal back-pressures. The paper is structured as follows; first, we out outline the set-up of our numerical experiments (section 2). Then, in Section 3 we briefly discuss physical aspects of null collapse based on past 2D results and the qualitative similarity and differences in extension to 3D. We then in Section 4 present the scaling of key measurements with variable resistivity in 2D and 3D collapses in low-β plasmas as, calculated from our simulations, in order to investigate empirically the applicability of previous analysis to the 3D case. We then consider the possibility of null collapse as a fast reconnection mechanism in astrophysics, and the effects of more appreciable ambient plasma pressure in Section 5. We summarise our results and draw conclusions in (Section 6).

Simulation Setup
The simulations involve the numerical solution of the 2.5D and 3D single-fluid, resistive MHD equations using the LareXd code [Arber et al., 2001]. Here we outline the simulation setup (initial conditions), with full technical details deferred to the appendix. All variables in this paper are nondimensionalised (see appendix A), unless units are explicitly stated.
We consider the collapse of magnetic fields containing null points, of the 2D Cartesian form and the 3D form each of which is a potential field, free from electric currents, and so constitutes a minimum energy state. These fields are therefore force free (j × B 0 = 0). These prescriptions, often referred to as 'linear null points', arise from the first-order Taylor series expansion near a null point embedded within a generic field, and so represent an approximation to realistic fields sufficiently close to the null point itself [see, e.g., Parnell et al., 1996, Priest andForbes, 2000, for full details]. We present results for four different geometries in this paper, namely the 2D null (equation 1) in conjunction with enforcing ∂/∂y = 0 throughout the solution (reducing the equations to so-called 2.5D) and the 3D null (equation 2) with the fully 3D MHD equations for k = 0.25, k = 0.5, and k = 1. Note, that setting k = 0 in equation (2) and rotating (the transformation x→z, z→ − x) recovers the same 2D null point as equation (1). In practice, we use the form of equation (1) preferentially in the 2D simulations for computational reasons, although the underlying physical problem is essentially the same. The fieldline structure of these fields is illustrated in Figure 1. Topologically, the fields consist of the null point itself (located at the origin), and (in 2D) a set of separatrices, or (in 3D) 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, connectively separated by the separatrices (or spine and fan in 3D), have a hyperbolic structure. The parameter k is an eccentricity parameter, introducing an azimuthal asymmetry to these fieldlines about the z-axis, taking on a preferential direction within the fan plane (the x-direction), and an associated rescaling of magnetic field strength in the z-direction. The computational domain is the cube |x, y, z| ≤ 1 (or equivalent square in 2D), with the boundaries taken to be closed with a fixed magnetic through-flux. The ratio of specific heats is γ = 5/3 throughout. Plasma resistivity η is taken as uniform and is a variable of the parameter study, and under our normalisation its value corresponds to an inverse Lundquist number as defined by the normalisation constants (see Appendix A for further details). Given that our outer boundaries are located at a nondimensional distance of 1 from the null, where the local Alfvén speed is of the order unity under the normalisation, this Lundquist quantifies the relative strength of the diffusivity of the whole domain. We take the plasma to be initially at rest (v = 0), of uniform density (ρ = 1) and a uniform gas pressure, chosen such that a fixed plasma-β defined at radius r = 1 may be set as a variable, denoted β 0 . Thus, the background state of the simulations is given by the force-free fields (equations 1, 2D, and 2, 3D) and a uniform gas pressure (zero pressure gradient), and is therefore an exact equilibrium (it has furthermore been verified to be a numerical equilibrium, as described in the Appendix B).
This background state is then subject to finite amplitude perturbations to the total field B = B 0 + B of the Cartesian form and in 3D, The perturbation, given by equations 3 and 4, corresponds to a superimposed, uniformly distributed current density in the y-component j y of initial magnitude j 0 . This current means that the separatrices (or the equivalent spine fieldline and fan separatrix plane in 3D) are no longer perpendicular. This perturbation to the initial background therefore disrupts force balance and, immediately after initialisation, the perturbation begins to focus towards the null, establishing a plasma flow that drives the collapse, as described qualitatively in the following section. We note that, due to the nature of our boundary conditions (Appendix B), that this corresponds to a perturbation which although uniform in the domain extends to the boundaries only. The boundary conditions are v = 0 and B line-tied (normal component fixed). As such there is no inflow of either plasma or magnetic energy. This is a crucial difference between our study and the linear collapse studies of e.g. Imshennik and Syrovatskii [1967], Bulanov and Olshanetsky [1984], Klapper et al. [1996]: in our case there is no energy inflow through the boundaries to drive unlimited collapse. This perturbation is thus of fixed spatial extent (extending to the edge of the domain) and finite total energy. The effects of boundary conditions on collapse, which is a subtle issue, is discussed further in Chapter 7.1 of Priest and Forbes [2000] and also by, e.g., Forbes and Speiser [1979], Klapper [1998]. As such, we stress the overall setup in this paper is that of a finite perturbation to a stable system. Presupposing this localised initial disturbance to the magnetic flux is motivated by the well-established result that MHD waves are generically attracted to null points, and so we expect that external perturbations to the larger-scale magnetic field will preferentially collect at nulls [see McLaughlin et al., 2011, for a review]. Specific examples of this in application can be seen in the work by Santamaria et al. [2015] and Tarr et al. [2017] where photospheric motions lead to current accumulation at nulls in realistic model solar atmospheres. In our paper, where we focus on the details of the current sheet formation (i.e. dynamics very close to the null rather than those in the external field), we thus presuppose this disturbance as both a matter of computational feasibility and as a modelling simplification (in line with previous null collapse studies, e.g. Craig 1996, Tsiklauri andHaruki 2007).
3 A brief overview of the physical and qualitative aspects of 2D and 3D null collapse in low-β plasma We begin by outlining the qualitative and physical aspects of our various simulations, before going on to quantify these results in Section 4. The perturbations considered in this paper (Section 2) correspond to an initial enhancement of free magnetic energy centred around the null. The magnetic field is not in force balance, and after t = 0 the evolution can be understood in terms of the propagation of the perturbation towards the null as a MHD wave. Due to the arrangement of the Lorentz force, the incoming wave immediately drives a fluid flow typical of reconnection, with the null itself coinciding with a stagnation point separating symmetric and anti-symmetric regions of inflow and outflow, delineated by the separatrices, or the spine and fan in 3D (i.e. streamlines of this flow field resemble rectangular hyperbolae). It is the focusing of this incoming wave (namely, its excess magnetic flux, and the flow it forces through its Lorentz force, both of which will increase in magnitude in time) which is at the heart of null collapse.
In the β < 1 case here, the incoming perturbation propagates as a fast (magnetic) MHD wave. Since first discussed by Dungey [1953], it has become well-established that such waves are, in-general, attracted to null points due to a refraction effect close to nulls, with fast waves propagating across surfaces of constant Alfvén speed, or down the potential well of the background field [e.g. McLaughlin et al., 2011, for a review]. Accordingly, the energy of the perturbation is propagated in towards the null point, transporting its flux/magnetic energy and associated mass-flow. As such, null collapse is a class of MHD implosive process with a null point being the center of converging magnetic flux, and also the aforementioned convergingdiverging flow (in this sense, it differs somewhat to most implosions which involve only a convergence of mass at a symmetry point, line or plane, dependent on dimensionality). In linear MHD such an implosion evolves in a relatively simple manner. Close to a generic null the field grows linearly away from the null: external disturbances accumulate at the null according to the linear Alfvén speed profile, and therefore the volume in which the energy is contained decreases exponentially (viz. gradients across this volume will increase exponentially). Therefore, as perturbations focus their energy in the vicinity of the null an exponential increase in current Figure 1: Representative fieldlines about the 3D equilibrium magnetic fields B 0 , for variable eccentricity parameter k (Equation 2) . In the k = 1 case (a 'proper' null point) topologically consists of a spine fieldline (red line), running along the z-axis into the null (at the origin), and the fan plane (transparent grey), sitting in the z = 0-plane, in which fieldlines run radially outward from the null. Elsewhere, fieldlines run hyperbolically, parallel to the spine or fan when far from the null. As k is decreased, the fieldlines become increasingly eccentric and take on a prefered direction in the fan plane. At k = 0, the null reduces to a translationally invariant 2D null structure, which is a 45 • rotation of Equation 1, consisting of two sets of separatricies (red fieldlines). In the k = 0 cases, the light blue fieldlines illustrate planar projections of the selected fieldlines. density at the null point occurs, a result that has been demonstrated both in the null collapse literature [e.g. Watson, 1992, McClymont andCraig, 1996] and in papers which consider current build up due to more generalised, externally originating MHD waves [e.g. McLaughlin and Hood, 2004, McLaughlin et al., 2011. This indicates that, at least in the low-β limit, the process is not too dependent on the rather symmetric initial conditions often employed in the collapse studies. As the waves cannot reach the null by propagation at the background Alfvén speed (v A → 0 as r → 0) this focusing is expected to continue until resistive diffusion and/or plasma backpressure becomes sufficiently large to allow the perturbations to propagate/diffuse through to the null, possibly reflecting the wave. In the related case of the collapse of a 2.5D X-line, the process is essentially the same although the guide-field also provides a magnetic backpressure which may oppose the collapse, in a manner analogous to the plasma back-pressure [McClymont and Craig, 1996].
In the more physically realistic case of a finite-amplitude (nonlinear) perturbation, the increasing amplitude of the focussing wave means that the excess flux carried by the wave may overwhelm the background field, so that the wave undergoes nonlinear evolution. Indeed, this will inevitably occur unless the collapse is first halted by resistive diffusion and backpressure as described above. As such, for our nonlinear MHD simulations, we expect there to be two distinct regimes in which the implosion evolves, depending on the initial wave energy. If the perturbation energy is sufficiently large that nonlinear evolution occurs, the implosion becomes characterised by quasi-1D behaviour (and as we will see in Section 4, different scaling laws as a consequence). This is because the magnetic field of the perturbation reinforces the background field in two quadrants -where the wavefront undergoes nonlinear evolution -and partially cancels it in the other two quadrants, where the wavefront stalls [Craig and Watson, 1992, McClymont and Craig, 1996, Gruszecki et al., 2011. This self-reinforcing process means that nonlinear collapse naturally creates sheet-like current distributions, as we will soon demonstrate. The quasi-1D phase of the collapse is closely related to the case of imploding 1D Harris-like current sheets (i.e., anti-parallel field about a null line) such as those considered by Forbes [1982] and Takeshige et al. [2015], where even, in certain limits, previous 2D null collapse and 1D collapse solutions have been shown to be equivalent [see e.g., the Appendix of Forbes, 1982].
We noted in our discussions of Thurgood et al. [2017], where null collapse was used to trigger Oscillatory Reconnection, that collapse for a k = 1 3D null appeared to proceed in a qualitatively similar way to the 2D collapse, within a particular plane containing the null point due to the same nonlinear evolution process of 2D collapse. We also showed (for a specific case) that the resulting reconnection was of the spine-fan type [Thurgood et al., 2017, Figure 4], which is consistent with related results of 3D reconnection excited by disturbances emanating at the computational boundary [e.g. . In the following section we aim to more broadly quantify this 3D collapse and the associated reconnection, for a variety of plasma parameters and null geometries. We note that during the process, the spine and fan of the null point locally collapse towards one another (see Figure 2). The plane in which this collapse occurs (the plane that contains the deformed spine line) is determined both by the perturbation that drives the collapse and by the null-point structure. In the remainder of this paper, we refer to this plane as the plane of collapse. Figure 2 shows some example current sheets resulting from this parameter study. It shows integrated values of the current density through the xz-planes ('side-view') which illustrates the aforementioned similarities to 2D collapse in planes of fixed-y, and also integrated values through xy-planes ('top-view') which shows the out-of-plane distribution of current. Importantly, as k is increased from k = 0 (the 2D null with no guide-field) to k = 1 (the rotationally 1 −1 |j(x, y, z)| dy illustrating the current distribution in the x − z halfplanes (a 'side-on' view of the null), where the top row corresponds the nonlinear regime (j 0 = 10 −1 ) and the bottom the linear (j 0 = 10 −1 ). Colours are fixed to a constant, saturated, linear scale for each row. We find that, in planes of fixed-y sufficiently close to the null that the collapse is qualitatively as per the 2D case. Bottom: Corresponding distributions of 1 −1 |j(x, y, z)| dz illustrating the current distribution in the x − y halfplanes (a 'top-down' view of the null). We find that the enhanced magnetic pressure for increasingly rotationally symmetric nulls (k → 1) inhibits collapse away from the null by providing a magnetic back-pressure via the growing out of plane field component B y = ky. symmetric limit), the out-of-plane extent of the current cylinder/ring (linear regime) or sheet (nonlinear) is reduced by a factor commensurate with the change in k -consistent with the behaviour found in boundary-driven simulations by Al-Hachami and Pontin [2010], Galsgaard and Pontin [2011]. This is due to the influence of the magnetic field component, B z = ky, which provides a magnetic back-pressure to resist the collapse (away from y = 0) in a manner analogous to the guide-field in the 2.5D X-line case [e.g. McClymont and Craig, 1996]. We will see that this has an important impact on the overall reconnection rate. These current concentrations can also be further contextualised in terms of the collapse field by comparison of Figure 2 to Thurgood et al. [2017, Figure 2b].

Quantitative scaling in low-β, compressible plasma
We now present the results of a parameter study of collapsing null points of the form given by equations (1) and (2) in a low pressure plasma such that β 0 = 10 −8 , for a variety of uniform plasma resistivity values in the range 10 −4 ≤ η ≤10 −2 (i.e. global Lundquist numbers in the range 10 2 − 10 4 ), and (in 3D) for field line eccentricities k = 0.25, 0.5, 1.0. The data shown in this section is for such simulations where the nulls are subjected perturbations (equations 3 and 4) which correspond to (initially) uniform current densities of magnitude j 0 in the j y component. We consider two different fixed-energy perturbation amplitudes, the first (j 0 = 10 −5 ) is expected to reside in the linear regime for the range of η considered, and the second (j 0 = 10 −1 ) is expected to behave nonlinearly, according to the arguments of Craig and McClymont [1993]. In the following subsections we first consider the critical times (the halting times) of each implosion, then present measurements of current sheet geometries along the width-wise, length-wise, and out-of plane (3D) axes, then finally the peak reconnection rates and evaluate also the efficiency of overall flux transfer by the critical time. Due to the extremely low resistivities of astrophysical plasmas, which are so small that they cannot be directly simulated, we pay particular attention to the scaling of these quantities across the computationally accessible range of resistivity. This mirrors the standard approach of previous (2D) collapse studies, such as McClymont and Craig [1996]. These scalings, by extrapolation, allow for an estimate of collapse behaviour as resistivity is further reduced. Figure 3 shows the time of peak current density (i.e. critical time, or the time of halting) as a function of plasma resistivity for different 2D and 3D runs in the linear and nonlinear collapse regimes. As pointed out by Craig and Watson [1992], a simple calculation suggests for the 2D null that the peak time in the linear case will simply be the time for the components of the perturbation at the boundary to travel to the critical radius r c of the diffusion region, which is r c = a √ η where a ∼ 1. Thus, for the 2D null, the shortest such time is for a radially propagating mode, giving t peak = 0.5 ln |η|, which is indicated by the straight line in the figure. The 2D simulations in the linear regime are in good agreement with this prediction, and we find that the principal effect of considering 3D nulls of decreasing eccentricity from the translationally invariant k = 0 case to the rotationally symmetric k = 1 case is a decrease in the peak time. This is naturally explained by the increasing background Alfvén speed as k → 1 (where B 2 0 = x 2 + k 2 y 2 + [k + 1] 2 z 2 ) decreasing the signal travel times. These critical times are found to be independent of the perturbation energy, so long as it is sufficiently small that the perturbation remains always in the linear regime (i.e., they are unchanged for other small values of j 0 , such as j 0 = 0.1η, which we do not show here). This amplitude-independent peak time confirms the linearity of the dynamics involved in these particular collapses. The perturbations which are expected to behave nonlinearly are indicated by unfilled symbols in Figure 3. Our simulations confirm that this is also the case for the 3D null collapses, and again, we find that the primary effect of decreasing 3D null point eccentricity is the uniform reduction in overall peak time across the range of η considered. This is, like the linear case, simply a consequence of faster background Alfvén speeds at the 3D nulls (which still influences the collapse and the focusing of the pulse before it enters its nonlinear stage). We note that we also considered a set of runs with perturbation j 0 = 10 −2 , which was be expected to transition between the linear and nonlinear collapses for the range of η considered. We indeed found that such runs begin to depart from the straight line when j 0 becomes of the order η for both 2D and 3D collapses in a manner reminiscent of McClymont and Craig [1996, their Figure 1], which we have not shown here to avoid excess clutter in the figure.

Current sheet geometry
Next, we consider geometrical scalings of the current sheet formed at the critical/peak times identified. Figure 4 shows measurements of the current sheet width w and length l for the different linear and nonlinear runs. Here, with regards to 3D nulls, width and length are defined  Figure 5: Left: Measured current sheet 'depth' d, for eccentricities k = 0.25, 0.5, 1.0 with corresponding symbols; blue squares, green circles and red diamonds respectively. Between different values of k, the measured d scales as approximately k −1 . Right: representative profiles a the peak time along the y-axis of the product ηj y , normalised to j 0 , for resistivity η = 10 −4 , 3 · 10 −4 , and 10 −3 for dashed-dot, dashed and dotted lines, respectively, with black (red) lines corresponding to the linear (nonlinear) case.
in the y = 0 plane (the plane of collapse) by analogy to the 2D null, and we will consider the out-of-plane distribution shortly (current sheet 'depth' d). As per Section 3 for runs in the linear regime the current distribution is not a 'true current sheet' but rather forms a more uniform ring-like distribution in the plane of collapse. For the 2D linear runs, we find that our simulations conform to the expected η 0.5 radial scaling discussed by McClymont and Craig [1996] with w = l. For the 3D linear runs, w was (arbitrarily) chosen to lie along the x-axis and l along the z-axis. The measurement conforms closely to the η 0.5 line with no discernible difference with eccentricity in the case of measurements along the x-axis (w), whereas the length l along the z-axis uniformly decreased for decreasing eccentricity due to the B z =B x imbalance. Thus, linear collapse at decreasingly eccentric 3D nulls produces increasingly ellipsoidal current distributions within the plane of collapse , departing from the cylindrical symmetry of 2D due to the nonuniformity of the background Alfvén speed (this feature is visible in the second row of Figure 2). For the nonlinear runs, w and l are defined more meaningfully as the short and long current sheet axis (in the y = 0 plane, for the 3D case). We find that for both our 2D and 3D runs, the width scales as η 0.89 (unfilled symbols, Figure 4a). McClymont and Craig [1996, their section 2.2 and 2.3] proposed that this scaling should be η 1 in the absence of density inhomogeneities in the current layer, or η 0.89 if such inhomogeneities play a role in halting the collapse. Their argument was based on comparison with 1D analytical results by Forbes [1982] describing the ideal implosion of a planar current sheet (taking into account the formation of density inhomogeneities) to the scale at which the diffusion speed should become comparable to wavespeeds in the inflow region. As such, this indicates that even in these initially low-β simulations growing current sheet inhomogeneity plays a role in halting the collapse and setting the smallest width w it can access -although we stress that it is ultimately halted by diffusive effects, as opposed to, say, a back-pressure associated with adiabatic heating in the compressed current sheet which would be η-independent. Figure 4 also indicates that null point eccentricity has little influence on current sheet width. We also find that for the nonlinear runs, the measured current sheet length l shows only a very weak (if any) dependence on η. This is consistent with the notion of Craig and Watson [1992] that the length is set by the point at which nonlinear evolution occurs in the width-wise direction (and so controls length-wise stalling). This is ultimately determined by the energy of the imploding perturbation, which is fixed across the different nonlinear runs. There is some slight variation with eccentricity, but unlike for other parameters there is no clear pattern. As such, we do not attribute any specific physics to this variability.
Turning now to the current sheet 'depth' d along the y-axis in the 3D simulations, the geometry of the out-of-plane current distribution is qualitatively as discussed in Section 3 and shown in Figure 2. To investigate the scaling, we plot the extent (full width at tenth maximum -FWTM -the maximum always being located at the null) of the j y distribution along the y-axis (the component contributing to 3D reconnection along that fieldline) in Figure 5. This measurement does not appear to follow a power-law scaling, perhaps due to the number of competing effects. We can understand the decrease in the FWTM as η decreases as follows. At large |y| there is a large magnetic pressure, that acts to halt the collapse. As η is reduced, the current sheet -in the absence of this magnetic pressure -should become progressively thinner across its width-wise axis (see Figure 4). However, for a given perturbation energy, the minimum current sheet width allowed by the magnetic pressure is fixed (though increases proportionally to |y|). Therefore as η is reduced, the point at which the magnetic backpressure halts the collapse moves towards progressively smaller |y|, such that the collapse is only resistively halted in a small region around the null (in y), being halted by magnetic backpressure for larger |y|. This is seen in the line plots to the right in Figure 5. All of the preceding implies that the current becomes increasingly peaked around y = 0 as η is decreased, leading  Figure 6: Left: Measures of the peak 3D reconnection rate, for the 3D runs with k = 0.25,0.5 and 1.0 (blue square, green circle and red diamond, respectively). A comparative 3D reconnection rate for the 2D runs (black triangles) is obtained by multiplying the 2D reconnection rate by the equivalent out-of-plane length L y = 2. Right: The product ηj at the null itself, which is similarly enhanced by collapse in the plane of collapse regardless of k.
to the observed scaling of the FWTM.

3D Reconnection rate, ηj at the null, and net flux transfer
We next consider the peak reconnection rate attained during the collapse. For the 2D runs, the reconnection rate is defined simply by ηj y at the null point. In three dimensions, the reconnection rate is the maximal value of the integration being performed along magnetic fieldlines [Schindler et al., 1988]. Following , by symmetry this coincides in our simulations with x=z=0 E y dy (verified by calculating the integral (5) along ∼ O(10 4 ) fieldlines traced from seed points clustered near the null). To compare the 2D and 3D rates, we simply multiply the recorded 2D rate by a factor of L y = 2 (the domain size in y for the 3D runs). These quantities are shown in Figure 6. For the 2D runs, we find in the linear case that the reconnection rate becomes independent of η in this range, which is an expected scaling from McClymont [1991, 1993] and was also found in the numerical experiments of McClymont and Craig [1996]. However, contrary to the simulations of McClymont and Craig [1996] our 2D nonlinear peak reconnection does seem to exhibit some (albeit, weak) scaling, decreasing with η. We hypothesise that this is due to the presence of Ohmic heating in our simulations (absent on those of McClymont and Craig 1996), which leads to an η-dependent increase in temperature, and therefore pressure, in the collapsing current layer.
In the 3D simulations, we find that for both linear and nonlinear collapse the reconnection rate depends on both η and the degree of magnetic field asymmetry, k. We note that the reconnection rate has a much stronger scaling with η that the value of ηj at the null (compare the left and right frames of Figure 6). Thus, the strong scaling of the reconnection rate in 3D with η is due primarily to the increasingly peaked current distribution along y discussed in the previous section. That is, the depression of the 3D rate is mainly due to the strong scaling of j y with η at large y, since it is essentially independent of η, E = ηj ∼ η 1 . This directly determines affects the 3D reconnection rate, it being obtained by integrating E along y. As such, we find that the rotationally symmetric 3D null produces a smaller peak reconnection rates than the asymmetric nulls (with smaller k).

Accessing fast reconnection through collapse in astrophysical plasmas
The main historical motivation of previous 2D collapse studies was assessing the viability of collapse as a fast reconnection mechanism, with a view towards, say, modelling flare energy release. We have found in the previous section that the basic physics of the implosion for maintaining a favourable (fast) scaling of the product ηj is unchanged in 3D; however, the introduction of a spatially-dependent magnetic back-pressure in 3D ultimately limits the spatial extent of the non-ideal, high-ηj region, which adversely impacts the reconnection rate. It is currently unclear if this is a significant limiting factor of reconnection efficiency in real astrophysical applications involving 3D nulls. This is in part due to the fact that the linear null employed here is a simplified model -in any real application the magnetic pressure would not grow monotonically away from the null. Therefore the impact of this stalling of the collapse will depend on the relative scales involved in the given problem. Nevertheless, we  , ηj at the null in the k = 1 3D simulations (center) and the corresponding true 3D reconnection rate (right) for variable β = 10 −8 , 10 −2 , 10 −1 , represented by +, × and * symbols, respectively. The linear (j 0 = 10 −5 ) and nonlinear (j 0 = 10 −1 ) runs are simply distinguished as being the lower and upper cluster of points along the vertical axes, respectively, as the plotted quantities contain no normalisation by j 0 . We find that in both cases a raised plasma-β curtails current build up a the null (as the current sheet does not narrow as much), decreasing the reconnection rate in the 2D case. However, we find that in the 3D case, the null current makes only a small contribution to the overall reconnection rate. For k = 1, the magnetic backpressure exerted by the B y component is sufficiently dominant mechanism curtailing collapse out of the y = 0 plane for the β considered. Figure 9: Demonstration of nonlinear tearing in an unstable, high-aspect ratio current sheet produced by nonlinear null collapse. This particular setup is as per the 2D parameter study, but with larger amplitude and more appreciable ambient pressure (β 0 = 10 −2 , η = 10 −5 and j 0 = 2). The instability develops after the halting time of the initial collapse.
can conclude from our results that more (rotationally) asymmetric nulls (that are 'closer to 2D') would generally be more favourable sites for magnetic energy release due to collapse than rotationally symmetric nulls.
Regardless of the fast instantaneous reconnection rates reported, it has been pointed out by Priest and Forbes [2000] (Chapter 7) it is also important to consider the total reconnected flux during the collapse in assessing the efficacy of magnetic energy release. Figure 7 shows the total flux reconnected up to the critical time for the various runs in our low-β parameter study. We find that only in the linear 2D case is the amount of reconnected flux resistivity-independent. In this case, curves of ηj(t) exhibit a simple translation, with the area underneath the curve remaining fixed. This translation is determined by the time required to advect the pulse at the background Alfvén speed to the diffusion-dominated scale (which is longer for smaller η). In the linear 3D case, this η-independence is no longer observed due to the action of the non-zero out of plane field component limiting the collapse at large |y|. The reconnected flux thus inherits the (relatively-weak) scaling with resistivity and the propensity to be lower as k → 1. For both the 2D and 3D nonlinear collapses, however, the total flux reconnected by peak time is strongly dependent on resistivity due to the fact that the implosion only produces significantly enhanced ηj (the peak value of which is η-dependent) for a short duration in a highly-impulsive fashion, which can also be seen in the curves in Figure 7.
An important further possible limitation of null collapse from the perspective of fast reconnection in astrophysical plasmas is the adverse effect of finite plasma pressure, first reported by McClymont and Craig [1996], Priest and Forbes [2000]. The observation is that pressure increases in the imploding current sheet until a sufficiently large outward pressure grows, halting the collapse. Given the nature of our results, there is no reason to expect this problem not to persist in the 3D case. To verify this, we repeated the 2D simulations and the k = 1 3D simulations in the linear and nonlinear perturbation regimes for β 0 = 10 −2 and β 0 = 10 −1 . Figure 8 shows the (2D) reconnection rate for the 2D simulations, the 3D reconnection rate for the 3D simulations, and the equivalent value of ηj at the null in the 3D simulations. We see that, for the 2D case, the overall dependence of peak reconnection rate departs from (near)-η independence in the low-β case (horizontal line) and begins to scale more strongly, although the curves are always less steep than a simple η 1 dependence which may be expected to be the limiting case for a failed collapse where reconnection occurs at the static rate. In the 3D, k = 1 case, the current at the null itself behaves similarly however we note that there is not much impact on the overall 3D reconnection rate (as the current at the null itself only makes a small contribution to the reconnection rate). This is because the collapse at large |y|, as discussed in Section 4, is strongly limited by a magnetic back-pressure. This indicates that, for these nulls, this magnetic back-pressure, which here is 3D null equivalent of the 'guide-field' back-pressure discussed in the 2D literature, remains the dominant throttling mechanism. Based upon the arguments given in Priest and Forbes [2000, Chapter 7.1], for (linear 2D) collapse to not enter this ambient pressure-limited regime requires β 0 < η 0.56 , and so for solar resistivity of the order of η ∼ 10 −10 this requires β 0 ≤10 −6 on the boundary. Although there are inherent difficulties in ascribing real length scales to the boundaries of problems involving linear nulls, this value is smaller than typical β inferred throughout various regions of the solar atmosphere via observation (however, other plasmas can have higher resistivities and therefore less prohibitive β 0 restriction). Thus, there remain doubts as to whether simple MHD null collapse is an efficient mechanism for magnetic energy release in the solar atmosphere, and we have found that considering the MHD collapse of a 3D null as opposed to 2D does not by itself remedy this problem.
The unfavourable scalings with resistivity in the nonlinear regime reported above thus suggest that the initial implosion may not provide a viable mechanism for fast energy release for coronal parameters. However, this initial collapse does set up a current sheet geometry that could lead to rapid energy release following some secondary process. Three candidates that we discuss below are secondary current sheet thinning as an MHD process, nonlinear tearing, or a collapse to collisionless scales. First, a 'secondary thinning' has been proposed by McClymont and Craig [1996]. The essential idea is that once the current sheet becomes highly pressurised, the pressure gradient may drive strong outflows that relieve the pressure enhancement, permitting further collapse. Accounting for thermal conductivity (excluded from our simulations) may also allow for a reduction in current sheet temperatures, also relieving internal current sheet pressure which sustains it against the inwardly directed Lorentz force of the implosion. To our knowledge, secondary thinning processes have not yet been further considered, and we note that although it is tempting to suppose that this may occur more readily in 3D (due to the additional dimension for outflow), we actually have previously found that for these 3D nulls the plasma is predominantly ejected in a collimated jet near the plane of collapse due to the curvature of fieldlines during spine-fan reconnection [Thurgood et al., 2017, Figure 5]. It is also the case that, if future attempts are made to investigate secondary thinning simulation, care should be taken that appropriate boundary conditions are employed (i.e. those allowing the passage of outgoing waves), otherwise secondary thinning could in fact be caused by the reflection and return of the outgoing fast shocks rather than being a self-consistent process within the locality of the current sheet. Alternatively, the collapse only needs to reach a scale which at which the effective local resistivity grows anomalously (i.e. a current-dependent resistivity), or a scale at which reconnection becomes collisionless [e.g. Haruki, 2007, 2008]. Otherwise, firmly within the realm of single-fluid, uniform resistivity MHD, it may be possible to promote fast reconnection via collapse regardless of β 0 by having a sufficiently energetic collapse that a current sheet forms with an aspect ratio that is susceptible to nonlinear tearing. This is now a well established route to fast reconnection [e.g. Huang et al., 2017, and references therein], and has been demonstrated by occur in a modified form in 3D by Daughton et al. [2011], Wyper and Pontin [2014]. At 3D nulls, it is expected that a nonlinear tearing of the current sheet occurs for Lundquist numbers greater than 2 × 10 4 and current sheet aspect ratios greater than 100 [Wyper and Pontin, 2014]. Since the instability takes some time to set in, and the current sheet formed during our collapse simulations gradually broadens following the initial implosion, it may be that the current layer at the critical time should have a larger aspect than this for nonlinear tearing to set in; however, this requires more careful study in a simulation with 'open' boundaries. What is clear is that -examining the nonlinear scalings for w and l obtained from Figure 4, for a given value of η, a sufficiently energetic perturbation should yield a current sheet beyond the critical aspect ratio. We have ourselves been able to observe current sheet tearing in 2D collapses after the initial phase of collapse ( Figure 9). In this simulation, we consider the same 2D setup as in the parameter study but with β 0 = 10 −2 , η = 10 −5 , and j 0 = 2 (larger than perturbations previously considered). The same nonlinear implosion process proceeds, rapidly producing a high-aspect ratio current sheet by t peak ≈ 0.6 which subsequently undergoes nonlinear tearing. This particular simulation however suffers from questions relating to the applicability of the closed boundary, namely that the current sheet has been impacted upon by the reflected fast waves before the instability develops have and also that, once it does, ejected plasmoids are artificially confined near the boundary due to the no-flow through boundary conditions. At present, it is unclear whether these effects of non-physical confinement significantly effect the evolution of the instability (and most crucially, whether it occurs at all). As such, we caution that this is a preliminary result intended primarily as a conceptual demonstration, and hope to study this further in the future (ideally, with an open system).

Conclusions and Discussion
In this paper we have considered a detailed parameter study of collapsing 3D magnetic null points of variable eccentricity (k = 0.25, 0.5, 1), alongside 2D nulls (equivalent to k = 0) with variable resistivity and variable initial perturbation amplitude, for both low and moderate plasma-β. The Key Findings are as follows: 1. The collapse of 3D nulls, across variable-k, is found to be both qualitatively and quantitatively as per the 2D case within the plane of collapse (here, the y = 0 plane). In both 2D and 3D, there exist two regimes of collapse: linear and nonlinear. These are characterised by self-similar and quasi-1D evolution, respectively, where the nature of the implosion is dependent on the relative energy of the perturbation and ability of the plasma to diffuse or resist the perturbation. For both regimes we find that the implosion proceeds to increasingly small length-scales as resistivity is decreased (here, seen and measured as the achieved current sheet widths, Figure 4). This length scale, which determines the current density via ∇ × B, is found to scale with resistivity as η 0.5 (linear) and η 0.89 (nonlinear). This leads to an independent (linear) or weak dependence (nonlinear) of the product ηj at the null on the resistivity, which is the reconnection rate in 2D rate and also contributes to the 3D rate.
2. The crucial difference between 2D and 3D null collapse occurs out of the y = 0 plane. For k > 0, the magnetic field component (B y = ky) increases away from the xz-plane in which the collapse proceeds. This field component acts to provide a magnetic backpressure that grows with distance from the null, opposing the collapse, and is analogous to the effect of guide-field back-pressure considered in 2D by McClymont and Craig [1996]. Thus, in the third dimension, the growth of ηj for y = 0 is increasingly limited for more rotationally symmetric 3D null points (k → 1). As the 3D reconnection rate (equation 5) is determined by the parallel electric field throughout the non-ideal volume (as opposed to ηj at the null itself), this magnetic back-pressure limits the overall reconnection rate for 3D nulls and so we find that the overall reconnection rate appreciably decreases as η is reduced, despite the favourable scalings for ηj at the null itself discussed in Key Finding 1 (cf. the left and right panels of Figure 6).
3. Increasing the plasma-β to realistic values for the solar corona (β 0 ∼ 10 −2 −10 −1 ) inhibits the collapse, leading to reduced reconnection rates. This effect is less severe in 3D since the magnetic back-pressure tends to be the dominant throttling mechanism across the majority of the current layer volume.
There are a number of astrophysical implications of these findings. Firstly, from the perspective of low-β null collapse as a fast reconnection mechanism, although we find that the basic physics of the implosion for maintaining favourable scalings of the product ηj is unchanged in 3D (Key Finding 1), it appears that a further complication of 3D collapse is that the introduction of a spatially-dependent magnetic back-pressure (Key Finding 2) ultimately limits the spatial extent of the non-ideal, high-ηj region, which negatively impacts upon the reconnection. Thus, we can conclude from our results that more eccentric, rotationally asymmetric nulls that would generally be more favourable in terms of magnetic energy release than a more rotationally symmetric counterpart. Regardless, we show directly that even if the peak reconnection rate obtains a favourable scaling with resistivity, that the total reconnected flux by the critical time is still limited by decreasing resistivity (Figure 7). It may be the case that higher rates of reconnection are maintained after the stalling of the collapse, but we cannot directly consider this from our simulations for the parameter study. As β 0 (the value of β at the boundary) is raised to values that are thought to be representative of the solar atmosphere, we have found that in the case of 3D null collapse the increasing initial plasma pressure interferes with the favourable scaling of the reconnection rate as resistivity is lowered (Key Finding 3, Section 5). Thus, there remain old and new doubts as to whether simple MHD null collapse is an efficient mechanism for magnetic energy release in the solar atmosphere, and we have found that considering the MHD collapse of a 3D null as opposed to 2D does not by itself provide any remedy to this problem. However, in Section 5 we have discussed several ways in which null collapse could bring about efficient reconnection and magnetic energy release as a secondary process even if the scaling of the initial reconnection rate itself is limited (either by the 3D out-of-plane magnetic pressure, or by too-high ambient plasma pressures). These processes include secondary thinning, accessing sufficient scales for current-dependent resistivity, microphysics and collisionless reconnection, and finally, by creating a current sheet that is susceptible to tearing after the implosion. We have been able to demonstrate tearing as the result of a 2D collapse using the simulation setup considered in these papers for relatively high β 0 = 10 −2 for the case of a closed system.
Outside of the fast reconnection perspective, we were also motivated to investigate 3D null collapse as we expect the results to be useful in future theoretical studies of Oscillatory Reconnection (OR). OR [McLaughlin et al., 2009, Thurgood et al., 2017 is a time-dependent and oscillatory magnetic reconnection system which is currently considered as a candidate for explaining quasi-periodic pulsations (QPPs) in solar and stellar flares [Nakariakov and Melnikov, 2009, Van Doorsselaere et al., 2016, McLaughlin et al., 2018, where a crucial question in assessing its applicability is; what periods can it produce for solar flare parameters (and are they compatible with QPPs)? It has been shown that the period of OR is dependent upon the initial disturbance to the null point field, behaving akin to a damped harmonic oscillator , and therefore a deeper understanding of null collapse can enhance our understanding of OR. This may be achieved by either simple assumptions regarding the proportion of available wave energy reaching close to the null, which would be equated to the null collapse system considered here as the initial perturbation energy, or with complementary numerical modelling efforts to estimate how much wave energy makes it to the immediate vicinity of the null whilst accounting for effects such as mode-conversion and atmospheric stratification [e.g., Tarr et al., 2017]. Indeed, our results could be used in tandem with globalscale simulations of these processes, which by necessity cannot resolve the details of the current sheet and reconnection dynamics.
Finally, we note that it may be possible to investigate the processes of implosive current sheet formation as described in this paper (and its predecessors) in the laboratory with devices such as CS-3D [Frank, 1999, Frank and Bogdanov, 2001, Frank and Kyrie, 2017, which can investigate current sheet formation and subsequent reconnection in a variety of null-containing fields including null lines (2D null points), null lines with guide fields, and truly 3D nulls. In these experiments a over-dense current sheet forms after some time delay which correlates with times expected for radial propagation of converging magnetoacousic waves from the edge to the centre of the chamber, which is then proceeded by 'metastable' stage, which is then followed by an eventual 'explosive' release of magnetic energy. We suggest that the theory of null collapse may describe the physics of the initial current sheet formation within such devices, explaining the compression ratios and sheet thickness achieved. Indeed, some of the reported results are suggestive of behaviour as predicted by null collapse theory, for example, it has been observed that electrical current is less effectively concentrated at small scales (i.e. measured current sheets are thicker) and the plasma within is less effectively compressed in response to growing guide fields in the case of the 2D null line with guide field [Frank et al., 2005].

A Nondimensionalisation and the Solver (LareXd code)
Following the details in the LareXd user manual, the normalisation is through the choice of three basic normalising constants, specifically: where quantities with and without a hat symbol are dimensional and nondimensional, respectively. These are then used to define the normalisation of quantities with derived units , t = t 0t and P = P 0P etc. Applying this normalisation to the ideal MHD equations simply removes the vacuum permeability µ 0 . In resistive MHD, this scheme leads naturally to a resistivity normalisation: Since v 0 is the normalised Alfvén speed this means thatη = 1/S where S is the Lundquist number as defined by the basic normalisation constants.
The simulation is the numerical solution of the nondimensional, resistive MHD equations: (NB: we drop the hat from this point onwards in the appendix, and throughout the main paper all quantities are nondimensional) which are solved on a Cartesian grid using the Lare2d and Lare3d codes. All results presented are in non-dimensional units. Algorithmically, the code solves the ideal MHD equations explicitly using a Lagrangian remap approach and includes the resistive terms using explicit subcycling [Arber et al., 2001[Arber et al., , 2016. The solution is fully nonlinear and captures shocks via an edge-centred artificial viscosity approach [Caramana et al., 1998], where shock viscosity is applied to the momentum equation through F shock and heats the system through H visc . Extended MHD options available within the code, such as the inclusion of Hall terms, were not used in these simulations. Full details of the code can be found in the original paper [Arber et al., 2001] and the users manual.

B Boundary conditions
The calculations presented in this paper represent the solution for the case of perturbed nulls contained within the cube |x, y, z| ≤ 1. The faces are subject to boundary conditions that permit no flow through or along the boundary (v = 0) with zero-gradient conditions taken ρ and , and also on magnetic field components which are tangential to a given face. The normal component of the field is held fixed (line-tied) through the boundary. The suitability of these boundary conditions, and overall stability of the setup, was checked by runs with no perturbation where we found that there was no undesirable behaviour such as the launching of spurious waves or erroneous current formation at the boundary. In practice, we apply the aforementioned conditions only on 'external' computational boundaries, and exploit appropriate symmetry/antisymmetry conditions on the 'internal' computational boundaries. Specifically, for the 3D setup we solve only for the half-cube y ≥ 0, and in 2D for the quarter-plane 0 ≤ (x, z) ≤ 1 (y = 0, arbitrarily). The more favourable reduction in 2D is facilitated by the fact that the form of equation (1) results in the current sheets length-wise and width-wise axes aligning with a computational boundary, hence why we have rotated the field in the 2D case. The implementation and accuracy of the symmetry conditions were checked simply by re-running some simulations in the whole domain, and we find perfect agreement. We note that the symmetry conditions are not used for the tearing mode simulation (Figure 9), but rather we simulate the full domain order to permit the symmetry breaking expected to occur during the instability.

C Grid geometry, resolution and testing
To adequately resolve the small-scale features produced by the collapse, especially in the lower resistivity cases, grid stretching is employed to concentrate resolution in the vicinity of the current sheets. 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: Γ = 1 2λ ln 1 + e λ − 1 0.5 1 − (1 − e −λ ) 0.5 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 (the null point) is controlled by the stretching parameter λ. Likewise, the same form and parameters are used for the distribution of cells in y and z. In our final 3D simulations of the parameter study, presented here we chose λ = 8 then performed simulations with increasing numbers of cells up to a maximum of nx = nz = 1024, ny = 512 (effectively, 1024 3 given the symmetry). Generally, we found that provided the resolution is sufficient to stop the current sheet collapsing to the grid-scale (i.e., capture the physics of the pressurisation and resistive heating of the current sheet which facilitates the halting process) the solution as measured by the maximal values of current density, density and other variables at the null itself demonstrates convergent behaviour as the numerical resolution is increased. In practice, only the simulations for the smallest resistivity can be ran within a reasonable time at 1024 3 , due to the unfavourable impact of smaller cell sizes upon the resistive timestep (∆t η ∝ ∆ x /η). Conveniently, however, higher values of η correspond to much wider current sheets at the critical time which therefore do not require such a fine grid (see Figure 4). The final resolution as used for the data presented in the parameter study is as follows (all k, linear and nonlinear); 10 −4 ≤ η ≤ 10 −3 use 1024 3 cells yielding ∆x min ≈ 0.00029, η = 3 × 10 −3 uses 512 3 cells yielding ∆x min ≈ 0.00058, and η = 10 −2 uses 256 3 cells, yielding ∆x min ≈ 0.00115. Each of these final simulations is in good agreement with a simulation at half the stated resolution (half of the cells in each dimension), in a qualitative sense during the evolution of the implosion and in the sense of producing the same scaling laws (which are in agreement with analytical results). They are also in an acceptable level of quantitative agreement with lower resolution simulations, with difference in measured current at the null being less than one percent when compared to the half-resolution case, except in the case of η = 10 −4 (the most challenging to resolve) which has the largest difference (∼ 3%). Overall, given the excellent agreement with analytic results for collapse scaling, which provide an independent means of verification where applicable, we are confident our simulations have faithfully captured the key aspects of the collapse up to the critical time.
In 2D, equivalent stretching schemes are utilised in x and z, however in test runs we also accessed much higher resolutions than possible for 3D for the sake of further testing (similar tests are also performed regarding variable stretching factors λ, to test the stretching).