Issue 
A&A
Volume 650, June 2021



Article Number  A143  
Number of page(s)  12  
Section  The Sun and the Heliosphere  
DOI  https://doi.org/10.1051/00046361/202140460  
Published online  22 June 2021 
KelvinHelmholtz instability and collapse of a twisted magnetic null point with anisotropic viscosity
^{1}
Research Software Development Group, University College London, Gower Street, London WC1E 6BT, UK
email: jamiejquinn@jamiejquinn.com
^{2}
School of Mathematics and Statistics, University of Glasgow, Glasgow G12 8QQ, UK
Received:
31
January
2021
Accepted:
3
April
2021
Context. Magnetic null points are associated with highenergy coronal phenomena such as solar flares and are often sites of reconnection and particle acceleration. Dynamic twisting of a magnetic null point can generate a KelvinHelmholtz instability (KHI) within its fan plane and can instigate spinefan reconnection and an associated collapse of the null point under continued twisting.
Aims. This article aims to compare the effects of isotropic and anisotropic viscosity in simulations of the KHI and collapse in a dynamically twisted magnetic null point.
Methods. We performed simulations using the 3D magnetohydrodynamics code Lare3d with a custom anisotropic viscosity module. A pair of highresolution simulations were performed, one using isotropic viscosity and another using anisotropic viscosity, keeping all other factors identical. We analysed the results in detail. A further parameter study was performed over a range of values for viscosity and resistivity.
Results. Both viscosity models permit the growth of the KHI and the eventual collapse of the null point. Over all studied parameters, anisotropic viscosity allows a faster growing instability, while isotropic viscosity damps the instability to the extent of stabilisation in some cases. Although the viscous heating associated with anisotropic viscosity is generally smaller, the ohmic heating dominates and is enhanced by the current sheets generated by the instability. This leads to a greater overall heating rate when using anisotropic viscosity. The collapse of the null point occurs significantly sooner when anisotropic viscosity is employed.
Key words: instabilities / plasmas / Sun: corona / magnetohydrodynamics (MHD)
© J. Quinn et al. 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
This paper presents the results of a series of numerical experiments intended to develop an understanding of the effect of anisotropic viscosity on the KelvinHelmholtz instability (KHI) in the fan plane of a magnetic null point, reproducing and extending the work of Wyper & Pontin (2013). We continue to stress the null point beyond the time investigated in Wyper & Pontin (2013), which allows us also to study the effect of anisotropic viscosity on the spontaneous collapse of the null point. The experiments take the form of a dynamic twisting of an initially static magnetic null point at the footpoints of its spine, resulting in a currentvortex sheet that forms in the fan plane and can be unstable to the KHI, given appropriate parameter choices. All experiments were carried out using both isotropic and anisotropic viscosity over a range of parameter choices. Anisotropic viscosity was modelled following MacTaggart et al. (2017). Continued driving after the moment at which the KHI occurs causes the null point to spontaneously undergo spinefan reconnection and collapse.
The KHI has been wellstudied in the magnetohydrodynamic (MHD) context and can be found in a number of coronal situations, both in numerical simulations (Howson et al. 2017; Wyper & Pontin 2013) and in observations (Foullon et al. 2011; Yang et al. 2018). Faganello & Califano (2017) offer a recent review of the KHI in MHD and Chandrasekhar (1961) provide a classical treatment.
In general, the effect of a magnetic field is stabilising; when the wave vector of a perturbation in a velocity shear layer is parallel or at an oblique angle to a magnetic field, magnetic tension acts to stabilise the KHI (Chandrasekhar 1961; Ryu et al. 2000). Otherwise, the KHI acts as an interchange instability and the magnetic field does not affect its linear stability (Chandrasekhar 1961).
In a currentvortex sheet, where a velocity shear coincides with a magnetic shear, the balance of shear layer strength and thickness dictates whether the KHI, the tearing instability, or a combination of the two is excited. Generally, when the magnetic shear is strong compared to the velocity shear, the KHI is suppressed and the tearing instability grows (Einaudi & Rubini 1986). This can be somewhat modified by the inclusion of viscosity (Einaudi & Rubini 1989). The nonlinear development of the KHI is known to enhance reconnection by the local distortion of magnetic field lines, the generation of current sheets (Min et al. 1997) and by generating local turbulence in conjunction with the tearing instability (Kowal et al. 2020).
The effect of (anisotropic) viscosity on the stability of a currentvortex sheet is to suppress the growth of the KHI, although viscosity is found to enhance the linear growth of the tearing instability when the KHI is stabilised by a strong magnetic field (Einaudi & Rubini 1989). A number of studies suggest isotropic viscosity can also slow and even suppress the KHI (Howson et al. 2017; Roediger et al. 2013; Wyper & Pontin 2013).
Magnetic null points are an abundant feature in the topologically complex coronal magnetic field (Edwards & Parnell 2015). Given they are sites that coincide with changes in topology, they are strongly associated with reconnection processes (Yang et al. 2020; Sun et al. 2013). Additionally they are inferred to participate in a number of highenergy phenomena, such as in the generation of flare ribbons in compact solar flares (Masson et al. 2009; Pontin et al. 2016) and the production of jets (MorenoInsertis & Galsgaard 2013) and of coronal mass ejections (Barnes 2007; Zou et al. 2020), particularly through their key involvement in the breakout model of solar eruptions (Antiochos et al. 1999; Maclean et al. 2005).
A muchstudied form of reconnection in 3D null points is spinefan reconnection, where a strong current sheet forms in the vicinity of the null point and enables efficient reconnection between the magnetic field that makes up the spine and fan plane, collapsing the field around the null in the process (Thurgood et al. 2018). The collapse of a null point has the potential to develop into a form of oscillatory reconnection (McLaughlin et al. 2012; Thurgood et al. 2017).
The layout of this paper is as follows. We present the numerical setup of the simulations in Sect. 2, including a description of the model of the linear null point and footpoint driver. The methods of calculating stability measures, shear layer properties, and reconnection rate are described in Sect. 3. In the first part of Sect. 4 we present the results of a highresolution pair of simulations for a single choice of viscosity and resistivity parameters and we compare the effects of the viscosity model used. In the second part, the results of a parameter study are described, generalising the highresolution results. The article concludes with a discussion of findings in Sect. 5 and conclusions in Sect. 6.
2. Model and numerical setup
2.1. Governing equations
We consider the nondimensionalised viscoresistive MHD equations,
posed in a rectangular domain Ω, where ρ is the mass density, u is the plasma velocity, p is the thermal pressure, ȷ is the current density, B is the magnetic field, σ is the viscous stress tensor, η is the resistivity, equivalent to the inverse Lundquist number, and ε is the specific energy density, given by the equation of state for an ideal gas,
where the specific heat capacity ratio is given by γ = 5/3. The terms Q_{ν} = σ : ∇u and Q_{η} = ηȷ^{2} are the viscous and ohmic heating contributions, respectively. The material derivative is D/Dt = ∂/∂t + (u ⋅ ∇).
The nondimensionalisation scheme is identical to that used in the code Lare3d (Arber et al. 2001), where a typical magnetic field strength B_{0}, density ρ_{0}, and length scale L_{0} are chosen and the other variables nondimensionalised appropriately. Velocity and time are nondimensionalised using the Alfvén speed and Alfvén crossing time t_{A} = L_{0}/u_{A}, respectively. Temperature is nondimensionalised via , where k_{B} is the Boltzmann constant and is the average mass of ions, here taken to be (a mass typical for the solar corona) where m_{p} is the proton mass. Dimensional quantities can be recovered by multiplying the nondimensional variables by their respective reference value (e.g. B_{dim} = B_{0}B). The reference values used here are B_{0} = 5 × 10^{−3} T, L_{0} = 1 Mm, and ρ_{0} = 1.67 × 10^{−12} kg m^{−3}, giving reference values for the Alfvén speed of u_{A} = 3.45 Mm s^{−1}, Alfvén time t_{A} = 0.29 s, and temperature T_{0} = 1.73 × 10^{9} K. These reference values are defined away from the null point, where the Alfvén speed is zero and the Alfvén time is not defined.
2.2. Models of viscosity
Here we use and compare both isotropic and anisotropic models of viscosity. In the isotropic case, the Newtonian viscosity model is used,
where ν is the viscous parameter (termed the viscosity throughout this paper) and W is the rateofstrain tensor,
In the anisotropic case, we use the switching viscosity model of MacTaggart et al. (2017),
where b = B/B is the unit vector in the direction of the magnetic field and s(αB) is the switching function; this is an interpolation function which controls the degree of anisotropy in the tensor based on the local magnetic field strength. This model focuses on the most important components (for the solar corona) of the full Braginskii tensor (Braginskii 1965), the parallel and isotropic components, in order better to model viscosity in the vicinity of magnetic null points (Hollweg 1986). While previous work (MacTaggart et al. 2017) has used a phenomenological form for the switching function s, here we use a form based on the coefficients of the Braginskii tensor,
where
and a = αB, where α is a parameter used to control the dependence of s on B.
Physically, α = eτ/m, where e is the electron charge, m is the ion mass, and τ is the ionion collision time
where n is the ion number density. For typical active region conditions, T = 2 × 10^{6} K and n = 3 × 10^{3} m^{−3} giving τ = 0.773 s and α ≈ 10^{8}. Were this value to be used in Eq. (9), s would change so rapidly with B that the region near the null point where s ≈ 0, corresponding to where the viscosity is isotropic, would be of subgrid scale at the resolutions used here. To properly resolve the region of isotropic viscosity, we choose α = 12. It should be noted that s(a) with α > O(10) has a very similar profile to that with α = 10^{8}. The main results that we will present will depend primarily on the region of anisotropic viscosity, outside the region of isotropic viscosity. Therefore, it is important for our results that they are not influenced by any potential numerical errors, at the null, which are outside our control. Our choice of α allows us to avoid such errors whilst still representing an improved model of plasma viscosity.
2.3. Null point model
The magnetic structure of the null point with the spine parallel to the zaxis is determined by an imposed initial magnetic field given in nondimensional units by
The domain Ω is a box of dimension [ − 3.5, 3.5]×[−3.5, 3.5]×[−0.25, 0.25] in the x, y, and z directions, respectively. We performed test simulations in both larger and smaller domains at resolutions of 320 grid points per dimension and the results are qualitatively similar. The initial velocity is uniformly zero, the initial density is uniformly ρ = 1, and the internal energy is uniformly ε = 5/4, corresponding to a temperature of 1.44 × 10^{9} K and a plasma beta of β ≈ 0.017. While this temperature is at the high end for the solar corona, the value was chosen to align with parameters used by Wyper & Pontin (2013). Additionally, test simulations run with lower temperatures resulted in prohibitively long running times.
On the domain boundaries, the fluxes of the magnetic field, the density, and the specific energy density are set to zero. That is, on the xboundaries,
and similarly, the y and zderivatives are zero on their respective boundaries. On the x and yboundaries, the velocity and its (normal) derivative are zero, and on the zboundaries the velocity is given by a driver, described below.
The driver takes the form of a slowly accelerating, rotating flow at the upper face of the box z = 0.25 given by
Here u_{r}(r) describes the radial profile of the twisting motion in terms of the radius r^{2} = x^{2} + y^{2},
where r_{d} controls the radial extent of the driver and u_{r0} is a normalising factor. The extent of the driving region can be seen in Fig. 1. The imposed acceleration of the twisting motion is described by u_{t}(t) with
Fig. 1. Field configuration after four Alfvén times. The driver speed is also shown as a slice, where white indicates peak driving speed. 
where the parameter t_{r} controls the time taken to reach the final driver velocity u_{0}. The parameters used in all simulations are u_{0} = 0.09, u_{r0} = 5.56, t_{r} = 0.25, and r_{d} = 36. At the lower boundary z = −0.25, the flow is in the opposite direction.
This driver twists the magnetic field around the spine of the null point, introducing twist throughout the entire structure (Fig. 1). The form of driver allows the system to be accelerated slowly enough that the production of disruptive shocks and fast waves is minimal. It is unavoidable that some waves are produced during the boundary acceleration, however these provide a useful source of noise that acts as a perturbation. As in Wyper & Pontin (2013), there is no prescribed perturbation that results in either the KHI or the null collapse; all perturbations are dynamically generated due to noise in the system. This entire setup is similar to that of Wyper & Pontin (2013).
The main parameter study required 18 simulations to be run in total; one per viscosity model for each of the nine parameter choices. To limit the time required to complete the study, a relatively low resolution of 320 grid points in each direction was used for these runs. A single, higher resolution pair of simulations were run, one for each viscosity model, at a resolution of 640 grid points for a single parameter choice. As well as allowing a detailed analysis of this case, these higher resolution simulations provide evidence that the lower resolution simulations have suitably converged.
3. Methods of analysis
3.1. Stability measures
Following Wyper & Pontin (2013), two quantities are used to understand the stability of the currentvortex sheet: the fast mode Mach number M_{f}, associated with the velocity shear, and a parameter Λ describing the balance of stability between the tearing mode and the KHI in a currentvortex sheet^{1}. The fast mode Mach number is given by
where Δu is the velocity change across the shear layer and c_{s} and v_{A} are the local sound and Alfvén speeds, respectively. The parameter Λ measures the relative strength of the velocity shear to magnetic shear and is given by
where L_{b} is the width of the magnetic shear layer, L_{u} is the width of the velocity shear layer, and M_{A} is the projected Alfvén Mach number
Since the shear layer occurs in the presence of a guide field (that of the initial magnetic null point), which is not included in the linear stability study of the KHI, the difference in magnetic field ΔB across the shear layer is used in the Alfvén Mach number as opposed to the full magnetic field strength B. In this way the Alfvén Mach number can be considered to be projected on to the shear layer.
Plotting the radial dependence of these quantities over the shear layers gives an indication of the local linear stability based on the stability analysis performed by Einaudi & Rubini (1986). The analysis predicts that a currentvortex sheet is linearly unstable to the KHI where M_{f} < 2 and Λ > 1. Where Λ < 1, the analysis predicts that the sheet is unstable to the tearing instability instead. It should be noted that the analysis of Einaudi & Rubini (1986) is 2D so can only be approximately used in the study of the KHI here, where there is an additional guide field in the system.
To calculate the stability measures, the peak vorticity and current density within the currentvortex sheets are measured, along with the radii at which the peaks occur. These radii are then used as the locations at which the absolute difference in azimuthal velocity Δu and magnetic field ΔB across the shear layers is measured, calculated as the difference between the maximum and minimum values of velocity or magnetic field either side of the shear layer. The distance between the maximum and minimum points gives a measurement of the thickness of the shear layers, L_{u} and L_{B}.
3.2. Reconnection rate
The reconnection rate is calculated using the same method employed in previous work by the same authors (Quinn et al. 2020). In summary, we calculate the reconnection rate local to a given magnetic field line as the local parallel electric field (that is, parallel to the magnetic field) integrated along the field line. By choosing a grid of starting points and integrating along each associated field line, an image is constructed of reconnection rates projected onto the grid of field line seed points. This is used to explore the spatial distribution of reconnection. The maximum value across all seed points gives the conventionally accepted measure of reconnection rate, the maximum integrated value (Galsgaard & Pontin 2011; Priest et al. 2003; Schindler et al. 1988).
4. Results
In the first subsection of the results section, the evolution of a pair of highresolution simulations is presented, both performed using resistivity of η = 10^{−4} and viscosity ν = 10^{−4}, but with different viscosity models. The purpose of these simulations is to capture the main features of the null point dynamics in response to the driver: the formation of a currentvortex sheet in the fan plane, the appearance of counterflows, the (potential) growth of a KHI, and the eventual collapse of the null point. This pair of highresolution simulations also highlights, in detail, the differences between the isotropic and the switching viscosity models, in particular the suppression of the KHI in the isotropic case, and the quicker collapse of the null point in the switching case. These results are then extended to other parameter choices in the second subsection.
4.1. Evolution of a typical case
We detail the evolution of the highresolution, typical cases in stages, exploring the formation, stability, and breakup of the currentvortex sheet before investigating the collapse of the null point. Then, we summarise the evolution through an analysis of the energy budget and the reconnection rate in time.
4.1.1. Formation of the currentvortex sheet
Initially, the torsional Alfvén waves injected by the driver trace out the field surrounding the null, moving first along the spine then out across the fan plane. This occurs from above and below. The upper and lower waves diffuse via viscosity and resistivity and eventually meet, creating shear layers in the velocity and magnetic field in the form of rings of vorticity and current density centred around the null point (Fig. 2a). These shear layers are jointly called the currentvortex sheet. Without any diffusion in the system, the waves would travel far along the fan plane before meeting. The presence of both viscosity and resistivity diffuses the waves as they travel along the field, allowing the upper and lower waves to meet around r = 1, where the currentvortex sheet forms. The hole in the sheet is due to magnetic tension forces opposing the twisting motion, as illustrated in Fig. 3 of Wyper & Pontin (2013). This also gives rise to counterflows (not shown) similar to those seen in Wyper & Pontin (2013) and Galsgaard (2003).
Fig. 2. Rings of vorticity and current density and associated linear stability criteria. Panel a: vorticity and current density for both viscosity models at t = 3 and z = 0. Panel b: linear stability measures as functions of radius at t = 6. The switching model permits rings of greater radial extent and notably stronger vorticity, resulting in a currentvortex sheet that is linearly unstable to the KHI. The thin dashed grey line Λ = 1 is the boundary between the KHI and the tearing instability from the 2D linear analysis of Einaudi & Rubini (1986). 
In the switching case, the reduced effective viscosity produces a vortex ring that is larger in radius and stronger in magnitude. The current density ring is somewhat larger in the switching case, but of equivalent peak magnitude to that in the isotropic case. Since the viscosity diffuses velocity directly and affects the magnetic field only indirectly, the vorticity is naturally affected by the change in viscosity model more than the current density. This difference in vorticity but not current density affects the relative size of the stability measures.
Figure 2b shows the relevant stability measures as functions of radius across the fan plane at t = 6, a time when the fan plane has become unstable to the KHI in the switching case but remains stable in the isotropic case (see Fig. 3b). The measure Λ confirms that the currentvortex sheet is linearly stable to the KHI in the isotropic case and unstable in the switching case for r > 0.6. This linear prediction matches the location where the KHI is observed to develop. In the switching case, the peak of M_{f} aligns with the observed region of initial growth of the instability. In the switching case, Λ and M_{f} are significantly larger due to the greater vorticity (Fig. 2a). In the isotropic case the more efficient dissipation of velocity results in a generally weaker vorticity ring and thus lower values of the stability measures.
Fig. 3. Development of the KHI in the outofplane velocity u_{z} at t = 2, 6, and 10 for both viscosity models. The isotropic results have been multiplied by as much as 1000 in order to compare them to the switching results. In the switching case, the KHI appears initially along the diagonals before extending azimuthally. In the isotropic case, there is no evidence of the instability. 
4.1.2. Disruption of the currentvortex sheet
As expected from the linear stability measures shown in Figs. 2b and 3 shows the development of the outofplane velocity from t = 2 to 10 and reveals that the currentvortex sheet is only unstable to the KHI in the switching case. Both cases exhibit a similar morphology of u_{z}, despite a growing difference in strength, until t = 6 when the KHI only appears in the switching case, initially along the diagonals (Fig. 3b) before spreading azimuthally (Fig. 3c). There is no evidence of the KHI in the isotropic case.
In both cases the currentvortex sheet grows in radius and magnitude with time, more in the switching case than in the isotropic. The shearing action of the counterflows produces a secondary ring of strong current density closer to the spine, which is greater in magnitude in the isotropic case. By t = 10 the KHI has disrupted the currentvortex sheet (Fig. 4a) and the resultant rolls create strong, smallscale current sheets, enhancing the local reconnection rate.
Fig. 4. Breakup of the currentvortex sheet and associated reconnection at t = 10. Panel a: current and vorticity density at z = 0. Panel b: spatial distribution of reconnection rate measured as the parallel electric field integrated along field lines traced from locations in the plane z = 0.23. Both panels show both viscosity models with isotropic shown on the left half of each figure and switching on the right. The currentvortex sheet remains stable in the isotropic case while that in the switching case has been fragmented by the KHI. The resultant smallscale reconnection in the rolls produces localised pockets of strong vorticity and current density. 
Figure 4b shows the spatial distribution of the reconnection rate for both viscosity models. Each pixel in the image represents one field line passing through that pixel along which the parallel electric field has been integrated. The colour of the pixel is given by the value of the integration and the starting locations of the field lines are a grid of points in the plane z = 0.23. The reconnection rate is greatest close to the centre of the plot, corresponding to regions of slippage reconnection due to the strong currents in the spine and currentvortex sheet. The effects of the boundary can be seen as long dark lines that spiral outwards from the origin. The switching case shows a greater peak reconnection rate due to the smallscale current sheets created by the KHI, and the enhanced reconnection far from the null can be seen as ripplelike structures at the fringes of the plot.
4.1.3. Spinefan reconnection
This section presents the results of driving the magnetic null point to the moment at which it undergoes spontaneous collapse. The collapse is instigated by a velocity shear across the null that generates a magnetic shear, permitting spinefan reconnection. We present the results of the isotropic case in detail first, then explore the effect of the KHI in the switching case.
In typical studies of spinefan reconnection (such as Pontin et al. 2007), the spines of a null point are dragged in opposite directions at the boundaries. This motion pulls the field above and below the null point in opposite directions and creates a current sheet that acts to reconnect field lines between the spine and fan. Here, the field near the null is shifted not because of motions at the footpoints of the spine, but due to imbalances in the velocity above and below the null point. These imbalances arise due to small pressure differences generated during the course of the initial driving. Figure 5 presents the magnetic field lines before and during the reconnection.
Fig. 5. Collapse of the null point visualised with field lines in the isotropic case. Field lines are plotted from a circle of radius 0.05 around the upper and lower spine footpoints. Contours of j = 60 are also plotted and reveal the strong current within the spine as well as the formation of the central sheet associated with the spinefan reconnection. At t = 18.5 the bulk of the field lines making up the core of the spine have reconnected. 
The twist in the spines creates a current that heats the contained plasma via ohmic heating and generates a small pressure force directed towards the null point, driving two oppositely directed streams of plasma along the spine (Fig. 6a). Where these streams meet (at the null point), they form a stagnation point flow, compressing the plasma in the vicinity of the null point and flowing out along the fan plane. Due to small asymmetries in the pressure that accrue during the simulation, an imbalance in the velocity appears above and below the null point (Fig. 6b).
Fig. 6. Velocity imbalance above and below the null point. Panel a: slice of the pressure through y = 0 overlaid with fluid velocity, where the longest arrows correspond to a fluid velocity of approximately 0.1. Panel b: u_{x}(x) − u_{x}(−x), the difference in u_{x} between the left and right sides of the plane x = 0. This gives a measure of the asymmetry in the velocity around the null point. 
The velocity shear around the null point shears the magnetic field accordingly, creating a current sheet through the null point (Fig. 5b). This current sheet enables reconnection between the spine and fan, which further extends, thins, and strengthens the sheet, continuing the reconnection process until the field around the null point collapses. The collapse itself can be seen in the kinetic energy as a dramatic increase starting at t ≈ 18 (Fig. 7a). The development of the current sheet and the resultant spinefan reconnection is similar to that of Pontin et al. (2007) with the exception that the twist in the field unravels as the reconnection proceeds (see also the similarity in behaviour compared to a model for solar polar jets in Pariat et al. 2009 which is driven by twisting motions). In the switching case, the development of the spinefan reconnection and associated collapse is qualitatively similar to that in the isotropic case with the exception that the reconnection occurs notably earlier and evolves over a shorter timescale.
Fig. 7. Energy components and reconnection rate as functions of time. (a) Kinetic energy. (b) Magnetic energy. (c) Internal energy. (d) Viscous heating. (e) Ohmic heating. (f) Reconnection rate. 
4.1.4. Energy budget and reconnection rate
The kinetic energy in the switching case is a measure of the main evolution of a KHIunstable currentvortex sheet as presented in detail previously (Fig. 7a). The initial injection of the Alfvén waves and formation of the currentvortex sheet can be seen at t ≈ 3. As the null point continues to be driven, the sheet becomes unstable to the KHI and the kinetic energy grows accordingly from t ≈ 3 to 8. At t ≈ 8 the KHI saturates as smallscale current sheets form and Ohmic heating begins to drain energy from the instability (Fig. 7f). This is also reflected in the reconnection rate (Fig. 7g) where the small current sheets in the rolls of the KHI in the fan plane enhance reconnection locally. Around t = 14, a transient increase in the kinetic energy reveals the start of the null point collapse.
In the isotropic case, the increased kinetic energy and enhanced reconnection rate associated with the KHI are absent, however the collapse of the null point produces significantly more kinetic energy at t ≈ 17 than in the switching case (Fig. 7a). The ohmic heating is similarly damped without the influence of the KHI (Fig. 7f). This results in the switching model extracting more energy from the field (Fig. 7b) and heating the plasma much more effectively (Fig. 7c). One significant finding is that the velocity shears created by the KHI allow anisotropic viscous heating of comparable levels to that of isotropic viscosity (in contrast to much larger differences observed in other situations e.g. the kink instability of a twisted flux tube presented in Quinn et al. 2020).
The reconnection rate reveals some interesting features about the nature of reconnection within the system and how the presence of the KHI affects the null point collapse (Fig. 7g). One interesting observation is that the switching case shows a greater reconnection rate than that of the isotropic case even before the onset of the KHI (i.e. for t < 6), suggesting the switching model allows for enhanced reconnection. As in Quinn et al. (2020), here too this is due to the switching model permitting greater velocities, greater compression, and thinner, stronger current sheets. It is then unclear whether the generally enhanced reconnection rate in the switching case for times t = 5 to 10 can be attributed to the currentenhancing effect of the switching model or is an effect of the KHI. Certainly, the spiky nature of the reconnection rate from t = 8 to 15 can be attributed to the small, strong current sheets produced in the rolls of the KHI, which do not appear in the stable isotropic case. The collapse of the null point is observed in the reconnection rate in the switching case around t = 15 and in the isotropic case around t = 17, however it differs significantly between the two cases. In the isotropic case, the reconnection rate increases during the collapse, while in the switching case, it decreases. This is due to the difference in the state of the nulls in each case as the collapse occurs.
In the isotropic case, where the KHI has not been excited, the flows and magnetic field are relatively simple and smooth such that the collapse is able to form large current sheets and reconnect many field lines at once. In contrast, in the switching case the KHI broke up the currentvortex sheet, introduced inhomogeneities throughout the fan plane, and generated small current sheets. This results in a collapse that struggles to reconnect with the same efficiency as in the smoother, simpler isotropic case. Additionally, there is simply more free magnetic energy in the system where the KHI remains stable, allowing current sheets to form more effectively during the collapse. In essence, the KHI places the null point in a more complex state where the collapse is less efficient at reconnecting field lines.
4.2. Study of parameter dependencies
The results shown in Sect. 4.1 change dramatically when the resistivity ν and the viscosity η are varied. This section presents results of simulations where ν takes values of 10^{−5}, 10^{−4}, or 10^{−3} and η takes values of 10^{−4} or 10^{−3}. This results in six pairs of simulations, each choice being run with switching viscosity or isotropic viscosity. We performed the simulations at a resolution of 320 grid points per dimension, half that of the highresolution cases described in detail above, and run to t = 15 instead of to t = 20 as in the latter. The null point collapse occurs sooner than at higher resolution and shows behaviour more typical of fast reconnection, which is indicative of inadequate resolution (Miyama et al. 2012). For this reason, focus is placed on the development of the KHI rather than on the null point collapse, leaving a parameter study of the null point collapse itself open as an avenue of future research. Generally, increasing the resistivity η to 10^{−3} produces a null point that is more unstable to the KHI (even in isotropic cases). Increasing the viscosity ν damps the KHI but does not totally suppress it, while decreasing ν leads to a faster growing KHI.
4.2.1. Shear layer stability
Figure 8 shows the stability measures as functions of radius for every studied parameter choice and for both viscosity models at t = 8. In every case M_{f} < 2, a necessary condition for instability of the currentvortex sheet. The condition on Λ for instability is Λ > 1. All layers with the switching model are linearly unstable to the KHI (Figs. 8b and d) while the isotropic cases show a mix of linear stability. When ν = 10^{−5}, the viscosity is weaker and the linear stability analysis predicts that the layers should be unstable for either value of η. The opposite is true for ν = 10^{−3}, when the isotropic viscosity is at its most dissipative.
Fig. 8. Stability measures as functions of radius r for all parameter choices at t = 8. Measures are plotted for ν = 10^{−5} (solid), ν = 10^{−4} (dashed), and ν = 10^{−3} (dotted) for each value of η and each viscosity model. The cases where Λ > 1 are unstable with the exception of the isotropic cases where ν = 10^{−4}, η = 10^{−3} (dashed line in Fig. 8c). (a) Isotropic; η = 10^{−4}. (b) Switching; η = 10^{−4}. (c) Isotropic; η = 10^{−3}. (d) Switching; η = 10^{−3}. 
In general, increased diffusion leads to a thicker, weaker ring, due to the Alfvén waves diffusing more before meeting in the fan plane. The switching model, being generally less diffusive than the isotropic model, permits velocity shear layers with much greater peak vorticity. Due to the coupling between the magnetic field and the velocity in an Alfvén wave, the isotropic model appears to provide some diffusion to the magnetic field during the formation of the magnetic shear layer. This results in a layer with weaker peak current, however the switching model affects the magnetic layer very little. Lower resistivity results in a stronger vorticity layer (in the switching case) but also a stronger current layer, and vice versa. for larger resistivity.
The observed stability of the currentvortex sheet to the KHI in each case is determined via inspection of the outofplane velocity for each parameter choice and is summarised for each parameter choice in Table 1. Some entries are marked as unstable*, referring to their being marginal cases; in these the KHI is directly observed in the outofplane velocity but the growth rate is close to zero and the perturbation remains negligibly small even at the final time of t = 15. The observed stability is well matched by the theoretical conditions of instability Λ > 1 and M_{f} < 2 in all but one case. This exception at η = 10^{−3}, ν = 10^{−4} is close to marginal and the result is most likely affected by the choice of resolution. This indicates that, despite the difference in geometry, the stability analysis of Einaudi & Rubini (1986) is of practical use in predicting the stability of the KHI in magnetic null points. This condition even accurately predicts the stability of the marginal cases.
Stability in the isotopic and switching cases for different choices of ν and η.
4.2.2. Kinetic energy profiles
Figure 9 shows the kinetic energy as a function of time for all parameter choices and for both viscosity models. The strongly KHIunstable cases show a similar kinetic energy profile to that of the unstable typical case (Fig. 7a). In the switching cases, the peak kinetic energy is larger when η = 10^{−4} in two cases (Figs. 9d and e). This is a result of the reduced diffusion of the magnetic field resulting in a stronger vorticity layer. The isotropic cases in Fig. 9 show an interesting trend in that the kinetic energy becomes stronger for smaller values of the viscosity ν (as expected) or for larger values of the resistivity η. In particular, the isotropic case where η = 10^{−3} and ν = 10^{−5} (Fig. 9a) is the only isotropic case where the KHI is significantly unstable. In this case, the kinetic energy profile shares a similar shape to the associated switching case, but the enhanced dissipation prevents the KHI from generating similar levels of kinetic energy. Instead, the profile is flatter and saturates at a later time.
Fig. 9. Kinetic energy as function of time for each parameter choice and viscosity model. An increase in ν damps the energy released during the KHI in the switching cases and totally suppresses the KHI in most isotropic cases. (a) ν = 10^{−5}; η = 10^{−3}. (b) ν = 10^{−4}; η = 10^{−3}. (c) ν = 10^{−3}; η = 10^{−3}. (d) ν = 10^{−5}; η = 10^{−4}. (e) ν = 10^{−4}; η = 10^{−4}. (f) ν = 10^{−3}; η = 10^{−4}. 
Table 1 reveals three cases of interest that are explored through their kinetic energy profiles. The marginally unstable isotropic case, where η = 10^{−4} and ν = 10^{−5}, does show some growth but it is notably less than the fully unstable case (Fig. 9d). Given that the currentvortex sheets in both these cases share similar strengths of vorticity, it is the combination of viscous dissipation of perturbations and enhanced magnetic shear in the lower η case that acts to stabilise the sheet. This conclusion can similarly be drawn for the outlying switching case where η = 10^{−4} and ν = 10^{−3} (Fig. 9f). The remaining case of interest is where η = 10^{−3} and ν = 10^{−4}, the single case where the linear prediction disagrees with the observed stability (Fig. 9b). In this case, the kinetic energy plateaus as the viscosity dissipates kinetic energy as it is generated by the instability.
4.2.3. Heating profiles
Figure 10 presents the total heat generated by viscous and ohmic dissipation and the total internal energy at t = 13, prior to any null collapse. Looking first at the viscous heating (Fig. 10a), Q_{ν} generally decreases with decreasing viscosity ν, as one may expect, with the exception of the isotropic cases where the resistivity is η = 10^{−4}. Instead, in these cases the viscous heating shows little dependence on ν. This reveals the complex, nonlinear relationship between viscous heating, the value of ν, and the flows generated.
Fig. 10. Internal energy, total viscous, and total ohmic heating as functions of viscosity ν at t = 12.5, before the onset of any null collapse, for all parameter choices. Isotropic (blue, solid) and switching viscosity (orange, dashed) are both shown. An upwardsfacing triangle denotes the higher value of η = 10^{−3} and a downwardsfacing triangle, η = 10^{−4}. The anisotropic viscous heating can become significant for larger values of ν yet, when ν is smaller, the lack of viscous heating is compensated by enhanced ohmic heating. 
In the switching cases, generally an increase in ν increases viscous heating and decreases ohmic heating. The decrease in ohmic heating is due to two complementary effects. Firstly, viscosity generally slows flows and limits the compression of current sheets, consequently limiting ohmic heating; thus a larger ν produces less ohmic heating. Secondly, the nonlinear phase of the KHI enhances ohmic heating in the fan plane and, since the instability is more unstable for smaller ν, ohmic heating increases with decreasing ν. The overall effect is a decrease in internal energy with increasing ν. This is also true for the η = 10^{−3} isotropic cases.
The ohmic heating profile similarly reveals complex behaviour in the isotropic cases (Fig. 10b). The stark difference in trends can be explained by considering the spatial distribution of ohmic heating that mirrors that of the current density. The two main current structures in a twisted null point are the currentvortex sheet and the structure associated with the twisted spines. Although the spine currents are two separate regions of current density, they contribute equally to the ohmic heating so are considered as one here. These are the main sources of ohmic heating and the balance of contributions from each source, for different values of η and ν, is nontrivial and results in the observed difference in trends.
Figure 11 reveals how the contributions from the currentvortex sheet and spines change with ν, and the effect that has on the total ohmic heating. These measures are calculated as the mean of the ohmic heating in the xyplane (representing the heating within the currentvortex sheet) and that in the yzplane (representing the heating in the spines). These are not true measurements of the ohmic heating within each current structure, however they provide a useful proxy.
Fig. 11. Ohmic heating contributions from separate current structures in the spine and fan. We show the mean ohmic heating contributions from the spines (dotted) and currentvortex sheet (dashed) and their sum (solid) for η = 10^{−3} (a) and η = 10^{−4} (b) in only the isotropic cases at t = 10. The value of η dictates how rapidly the balance of contributions shifts from fan to spine with ν, resulting in different trends in total ohmic heating. 
For any value of η, the spine heating increases and the currentvortex heating decreases as ν increases. This is due to greater viscosity dissipating the initial Alfvén waves more effectively and reducing the magnetic shear in the currentvortex sheet while retaining magnetic shear in the spines. The difference in how rapidly the relative contributions change with ν gives rise to the difference in total ohmic heating trends found in Fig. 10b. When η = 10^{−4}, the heating in the sheet decreases faster than the spine heating increases with ν, resulting in a drop in total ohmic heating (Fig. 11b). The opposite is true when η = 10^{−3} (Fig. 11a).
5. Discussion
While the values of the resistivity η used in the simulations performed here are orders of magnitude greater than typical coronal estimates, the values of viscosity ν are certainly within realistic bounds. When using a model of viscosity appropriate in the solar corona (the switching model), we found that the KHI is unstable, regardless of parameter choice. This strongly suggests that in the real corona the KHI can be excited in currentvortex sheets similar to those studied here.
This paper details the investigation of the KHI and null point collapse around an axisymmetric, linear null point, an idealised model of a real null point (such as that observed in Masson et al. 2009). The impact of different null point configurations such as those with asymmetry (e.g. those investigated by Thurgood et al. 2018; Pontin et al. 2016) is unclear. Similarly, the simplicity of the driver used here is unlikely to reflect the true nature of drivers in the real solar corona. The impact of driver complexity on spinefan reconnection has been specifically investigated by Wyper et al. (2012) who, however, focus on sheared drivers as opposed to the torsional drivers employed here. It would be of interest to understand how different magnetic field configurations and forms of driver affect the formation and stability of the kind of currentvortex sheets studied here.
The simulations detailed here have been performed with a model of anisotropic viscosity that only captures the parallel component of viscosity. As discussed by Einaudi & Rubini (1989), perpendicular components can become significant in strong velocity shears (such as those found in the fan plane of a twisted null point) despite the small size of the associated transport coefficient. A similar set of experiments exploring the effect of perpendicular viscosity could provide useful insight, particularly in ascertaining if the growth of the tearing instability in the currentvortex sheet could be accelerated by perpendicular viscosity, as is found in the linear analysis performed by Einaudi & Rubini (1989). That being said, we expect that the inclusion of perpendicular effects will lead to differences that are quantitative (e.g. affecting the numerical size of the KHI growth rate) rather than qualitative. This is because the perpendicular components of viscosity act in localised regions (e.g. MacTaggart et al. 2017) and not everywhere in the domain like isotropic viscosity. Thus, we expect the effects of perpendicular viscosity to play a role more similar to parallel viscosity than isotropic viscosity.
An important finding of this investigation is the spontaneous collapse of the null point without shearing drivers (as in Pontin et al. 2007) or prescribed current density perturbations (as in Thurgood et al. 2018). The formation of the current sheet at the null point (which facilitates spinefan reconnection and null point collapse) is primarily driven by a shearing motion at the null point, itself a result of oppositely directed streams of plasma flowing along the spine towards the null point. These flows rely on pressure gradients that are generated by ohmic heating in the twisted spine. It may be that the pressure generated under ideal conditions (or using realistically small coronal resistivity) is not enough to drive the nulldirected flows and, thus not enough to collapse the null point. Further investigation of this form of collapse within a twisted null point is required to ascertain if such a collapse is possible in the real corona.
We also explored the effect of the form of viscosity on the collapse of the null point in the two highresolution simulations (of Sect. 4.1) and we found that in the switching case, where the KHI is unstable, the null point collapses notably earlier than in the isotropic case, where the KHI is stable. It is unclear if the early null collapse is a consequence of the KHI or the use of the switching model. From the results of the isotropic case, the null point collapse appears to be ultimately caused by slight asymmetries in the spinealigned flows, so one may conjecture that, in the switching case, the KHI introduces its own asymmetries that cause the early collapse of the null point. A higher resolution version of the unstable, isotropic case (where ν = 10^{−5} and η = 10^{−3}) would provide clarity.
It is unclear how the unravelling of the null point as it collapses affects the ability of the null point to undergo the kind of oscillatory reconnection found by Thurgood et al. (2017). One observed phase in the oscillatory process is the generation of backpressure, which halts and reverses the spinefan reconnection process. The simulations performed here do not run for a long enough time to investigate the generation of backpressure. It may be that a collapsing twisted null point is unable to produce the required backpressure if it unravels during its initial collapse. Running the highresolution simulations reported here for a longer time would reveal if the particular setup studied can generate oscillatory spinefan reconnection. Alternatively, using a pretwisted null point as an initial condition with the perturbation used to collapse the null point found in Thurgood et al. (2017) would provide a similar experiment.
6. Conclusions
In this paper two models of viscosity have been applied to a magnetic null point that has been dynamically twisted at its footpoints in such a way that a currentvortex sheet forms in the fan plane. This sheet has the potential to become unstable to the KHI. We found that increased viscous dissipation, particularly in the form of isotropic viscosity, has a stabilising effect on the sheet, to the point of complete suppression of the instability. This is primarily due to viscosity thickening the sheet and thus increasing its stability. The presence of the instability enhances reconnection and viscous heating within the sheet.
After some time, the null point spontaneously collapses due to an imbalance in spinedirected, pressuredriven flows. This was found to occur sooner when the KHI is present. The general development of the collapse and associated spinefan reconnection is similar to that of previous work with the exception that the twist in the spine unravels during the collapse.
The investigation of the stability of the currentvortex sheet was extended with a parameter study over an order of magnitude difference in resistivity and two orders of magnitude in viscosity. The results show that the KHI is mostly unstable when using anisotropic viscosity and mostly stable when using isotropic viscosity.
The general qualitative behaviour that we have found in this study matches well with our previous study on the effect of anisotropic viscosity on the kink instability in a flux tube (Quinn et al. 2020). The more localised influence of anisotropic viscosity compared to isotropic viscosity allows for the creation of much smaller length scales, both in the magnetic and velocity fields. The weaker effect of direct damping, by anisotropic viscosity compared to isotropic viscosity, means that the shorter length scales also coincide, in general, with higher magnitudes, thus enhancing the possibility of instability. Even with different dynamical drivers (the null point is twisted in this study whereas the flux tube in Quinn et al. 2020 is initially unstable and allowed to relax), the cases with anisotropic viscosity exhibit the fast development of instabilities and reconnection in the way described above. Although anisotropic viscosity does not affect the magnetic field directly, it does have a significant indirect influence through nonlinear interaction. Therefore, in the study of nonlinear dynamics in the solar corona, taking account of anisotropic viscosity is important as it can result in the development of behaviour that would be absent if only isotropic viscosity were considered.
Finally, our results have important implications for understanding energy release in the corona. The more efficient reconnection described above, which is an indirect consequence of anisotropic viscosity, allows for two important effects in the corona. First, the amount of ohmic heating produced is greater with anisotropic viscosity compared to isotropic viscosity. This is simply due to the production of many smallscale current sheets, as described above. Second, magnetic instabilities develop more quickly with anisotropic viscosity. By this, we mean that the build up or development (Quinn et al. 2020) of instabilities occurs faster. This is because more efficient reconnection, with anisotropic viscosity, allows for the easier release of energy (see Fig. 7a for example). With isotropic viscosity, the smallscale current sheets allowed by anisotropic viscosity do not form, and so it is more difficult for the magnetic field to release its energy. Quinn et al. (2020) showed that in an evolving kinkunstable flux tube, magnetic relaxation occurs through secondary instabilities when only isotropic viscosity is present. That is, the field has to wait until an ideal instability occurs to release its energy as reconnection is not efficient enough to release it sooner. In the present study, the instabilities in the magnetic field develop faster for the anisotropic case, whilst there is a slower build up for the isotropic case, leading to a later phase of instability (see Fig. 7a again for a clear picture of this). There are many models of coronal dynamics that are based on the slow buildup of magnetic energy through photospheric deformations leading to a sudden release via an ideal instability (an example that is pertinent to our study can be found in Pariat et al. 2009). Our results suggest that instabilities of this type may be more efficient due to the indirect influence of anisotropic viscosity on reconnection. It will be an important avenue of research to model the effects we have found in this work on an active region scale.
Wyper & Pontin (2013) additionally use the projected Alfvén Mach number alongside the two measures used here, however it is our opinion that Λ captures the same information.
https://github.com/jamiejquinn/khi_null_point_code and have been archived at https://doi.org/10.5281/zenodo.4483221
Acknowledgments
Results were obtained using the ARCHIEWeSt High Performance Computer (www.archiewest.ac.uk) based at the University of Strathclyde. JQ was supported by the Engineering and Physical Sciences Research Council [grant number EP/N509668/1].
References
 Antiochos, S. K., DeVore, C. R., & Klimchuk, J. A. 1999, ApJ, 510, 485 [Google Scholar]
 Arber, T., Longbottom, A., Gerrard, C., & Milne, A. 2001, J. Comput. Phys., 171, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Arber, T., Bennett, K., Quinn, J., & Brady, C. 2020, https://doi.org/10.5281/zenodo.4155646 [Google Scholar]
 Barnes, G. 2007, ApJ, 670, L53 [NASA ADS] [CrossRef] [Google Scholar]
 Bennett, K., Arber, T., Quinn, J., & Brady, C. 2020, https://doi.org/10.5281/zenodo.4155546 [Google Scholar]
 Braginskii, S. I. 1965, Rev. Plasma Phys., 1, 205 [Google Scholar]
 Chandrasekhar, S. 1961, Hydrodynamic and Hydromagnetic Stability (Oxford: Clarendon Press) [Google Scholar]
 Edwards, S. J., & Parnell, C. E. 2015, Sol. Phys., 290, 2055 [Google Scholar]
 Einaudi, G., & Rubini, F. 1986, Phys. Fluids, 29, 2563 [Google Scholar]
 Einaudi, G., & Rubini, F. 1989, Phys. Fluids B: Plasma Phys., 1, 2224 [Google Scholar]
 Faganello, M., & Califano, F. 2017, J. Plasma Phys., 83, 535830601 [Google Scholar]
 Foullon, C., Verwichte, E., Nakariakov, V. M., Nykyri, K., & Farrugia, C. J. 2011, ApJ, 729, L8 [NASA ADS] [CrossRef] [Google Scholar]
 Galsgaard, K. 2003, J. Geophys. Res., 108, 1042 [Google Scholar]
 Galsgaard, K., & Pontin, D. I. 2011, A&A, 529, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hollweg, J. V. 1986, ApJ, 306, 730 [NASA ADS] [CrossRef] [Google Scholar]
 Howson, T. A., Moortel, I. D., & Antolin, P. 2017, A&A, 602, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kowal, G., FalcetaGonçalves, D. A., Lazarian, A., & Vishniac, E. T. 2020, ApJ, 892, 50 [Google Scholar]
 Maclean, R., Beveridge, C., Longcope, D., Brown, D., & Priest, E. 2005, Proc. R. Soc. A: Math. Phys. Eng. Sci., 461, 2099 [Google Scholar]
 MacTaggart, D., Vergori, L., & Quinn, J. 2017, J. Fluid Mech., 826, 615 [NASA ADS] [CrossRef] [Google Scholar]
 Masson, S., Pariat, E., Aulanier, G., & Schrijver, C. J. 2009, ApJ, 700, 559 [NASA ADS] [CrossRef] [Google Scholar]
 McLaughlin, J. A., Thurgood, J. O., & MacTaggart, D. 2012, A&A, 548, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Min, K. W., Kim, T., & Lee, H. 1997, Planet. Space Sci., 45, 495 [Google Scholar]
 Miyama, S. M., Tomisaka, K., & Hanawa, T. 2012, Numerical Astrophysics: Proceedings of the International Conference on Numerical Astrophysics 1998 (NAP98), Held at the National Olympic Memorial Youth Center, Tokyo, Japan, March 10–13, 1998 (Springer Science& Business Media) [Google Scholar]
 MorenoInsertis, F., & Galsgaard, K. 2013, ApJ, 771, 20 [Google Scholar]
 Pariat, E., Antiochos, S. K., & DeVore, C. R. 2009, ApJ, 691, 61 [Google Scholar]
 Pontin, D. I., Bhattacharjee, A., & Galsgaard, K. 2007, Phys. Plasmas, 14, 052106 [NASA ADS] [CrossRef] [Google Scholar]
 Pontin, D., Galsgaard, K., & Démoulin, P. 2016, Sol. Phys., 291, 1739 [Google Scholar]
 Priest, E. R., Hornig, G., & Pontin, D. I. 2003, J. Geophys. Res., 108, 1285 [Google Scholar]
 Quinn, J., MacTaggart, D., & Simitev, R. D. 2020, Commun. Nonlinear Sci. Numer. Simul., 83, 105131 [Google Scholar]
 Roediger, E., Kraft, R. P., Nulsen, P., et al. 2013, MNRAS, 436, 1721 [NASA ADS] [CrossRef] [Google Scholar]
 Ryu, D., Jones, T. W., & Frank, A. 2000, ApJ, 545, 475 [NASA ADS] [CrossRef] [Google Scholar]
 Schindler, K., Hesse, M., & Birn, J. 1988, J. Geophys. Res., 93, 5547 [NASA ADS] [CrossRef] [Google Scholar]
 Sun, X., Hoeksema, J. T., Liu, Y., et al. 2013, ApJ, 778, 139 [NASA ADS] [CrossRef] [Google Scholar]
 Thurgood, J. O., Pontin, D. I., & McLaughlin, J. A. 2017, ApJ, 844, 2 [Google Scholar]
 Thurgood, J. O., Pontin, D. I., & McLaughlin, J. A. 2018, ApJ, 855, 50 [NASA ADS] [CrossRef] [Google Scholar]
 Wyper, P. F., & Pontin, D. I. 2013, Phys. Plasmas, 20, 032117 [NASA ADS] [CrossRef] [Google Scholar]
 Wyper, P. F., Jain, R., & Pontin, D. I. 2012, A&A, 545, A78 [EDP Sciences] [Google Scholar]
 Yang, H., Xu, Z., Lim, E.K., et al. 2018, ApJ, 857, 115 [Google Scholar]
 Yang, S., Zhang, Q., Xu, Z., et al. 2020, ApJ, 898, 101 [CrossRef] [Google Scholar]
 Zou, P., Jiang, C., Wei, F., et al. 2020, ApJ, 890, 10 [CrossRef] [Google Scholar]
Appendix A: Associated software
A.1. Anisotropic viscosity module
A custom version of Lare3d (Arber et al. 2001) has been developed where a new module for anisotropic viscosity has been included. The new version can be found online^{2}, and also archived at Bennett et al. (2020), and should be simple to merge into another version of Lare3d for future research. The version of Lare3d used in the production of the results presented here, including initial conditions, boundary conditions, control parameters, and the anisotropic viscosity module, can be found at Arber et al. (2020). The data analysis and instructions for reproducing all results found in this report may be also found online^{3}.
All simulations were performed on a single, multicore machine with 40 cores and 192 GB of RAM, although this amount of RAM is much higher than was required; a conservative estimate of the memory used in the largest simulations is around 64 GB. Most simulations completed in under two days, although the longest running simulations (the highest resolution cases shown here) completed in around two weeks.
A.2. Field line integrator
As described in Sect. 3.2, the reconnection rate local to a single field line is given by the electric field parallel to the magnetic field, integrated along the field line. The global reconnection rate for a given region of magnetic diffusion is the maximum value of the local reconnection rate over all field lines threading the region. A field line integrator was developed specifically for this calculation and is detailed here.
Magnetic field lines lie tangential to the local magnetic field at every point x(s) along the line,
where s is a variable that tracks along a single field line and b is the unit vector in the direction of B. This equation is discretised using a secondorder RungeKutta scheme to iteratively calculate the discrete positions x_{i} along a field line passing through some seed position x_{0},
where h is a small step size. Since b is discretised, the value at an arbitrary location x_{i} is calculated using a linear approximation. The integration of a scalar variable y is carried out along a field line given by a sequence of N locations x_{i} using the midpoint rule,
where Y is the result of the integration. In practice, N is not specified and the discretised field line contains the required number of points to thread from its seed location to the boundary of the domain.
While the linear interpolation, secondorder RungeKutta, and midpoint rule are all low order methods, testing higher order methods showed little change in results but dramatically increased the runtime of the analysis. The lower order methods used offer an acceptable compromise between speed and accuracy. The above algorithm is implemented in Python and can be found in the code directory linked in the Associated software section above, in shared/field_line_integrator.py with examples of use in main/field_line_integrator.Rmd. The integration of multiple field lines is an embarrassingly parallel problem and is parallelised in a straightforward manner using a pool of threads supplied by the Pool feature of the Python library multiprocessing. Although the integrator is used solely to integrate the parallel electric field along magnetic field lines in this study, the tool can be easily applied to arbitrary vector and scalar fields.
All Tables
All Figures
Fig. 1. Field configuration after four Alfvén times. The driver speed is also shown as a slice, where white indicates peak driving speed. 

In the text 
Fig. 2. Rings of vorticity and current density and associated linear stability criteria. Panel a: vorticity and current density for both viscosity models at t = 3 and z = 0. Panel b: linear stability measures as functions of radius at t = 6. The switching model permits rings of greater radial extent and notably stronger vorticity, resulting in a currentvortex sheet that is linearly unstable to the KHI. The thin dashed grey line Λ = 1 is the boundary between the KHI and the tearing instability from the 2D linear analysis of Einaudi & Rubini (1986). 

In the text 
Fig. 3. Development of the KHI in the outofplane velocity u_{z} at t = 2, 6, and 10 for both viscosity models. The isotropic results have been multiplied by as much as 1000 in order to compare them to the switching results. In the switching case, the KHI appears initially along the diagonals before extending azimuthally. In the isotropic case, there is no evidence of the instability. 

In the text 
Fig. 4. Breakup of the currentvortex sheet and associated reconnection at t = 10. Panel a: current and vorticity density at z = 0. Panel b: spatial distribution of reconnection rate measured as the parallel electric field integrated along field lines traced from locations in the plane z = 0.23. Both panels show both viscosity models with isotropic shown on the left half of each figure and switching on the right. The currentvortex sheet remains stable in the isotropic case while that in the switching case has been fragmented by the KHI. The resultant smallscale reconnection in the rolls produces localised pockets of strong vorticity and current density. 

In the text 
Fig. 5. Collapse of the null point visualised with field lines in the isotropic case. Field lines are plotted from a circle of radius 0.05 around the upper and lower spine footpoints. Contours of j = 60 are also plotted and reveal the strong current within the spine as well as the formation of the central sheet associated with the spinefan reconnection. At t = 18.5 the bulk of the field lines making up the core of the spine have reconnected. 

In the text 
Fig. 6. Velocity imbalance above and below the null point. Panel a: slice of the pressure through y = 0 overlaid with fluid velocity, where the longest arrows correspond to a fluid velocity of approximately 0.1. Panel b: u_{x}(x) − u_{x}(−x), the difference in u_{x} between the left and right sides of the plane x = 0. This gives a measure of the asymmetry in the velocity around the null point. 

In the text 
Fig. 7. Energy components and reconnection rate as functions of time. (a) Kinetic energy. (b) Magnetic energy. (c) Internal energy. (d) Viscous heating. (e) Ohmic heating. (f) Reconnection rate. 

In the text 
Fig. 8. Stability measures as functions of radius r for all parameter choices at t = 8. Measures are plotted for ν = 10^{−5} (solid), ν = 10^{−4} (dashed), and ν = 10^{−3} (dotted) for each value of η and each viscosity model. The cases where Λ > 1 are unstable with the exception of the isotropic cases where ν = 10^{−4}, η = 10^{−3} (dashed line in Fig. 8c). (a) Isotropic; η = 10^{−4}. (b) Switching; η = 10^{−4}. (c) Isotropic; η = 10^{−3}. (d) Switching; η = 10^{−3}. 

In the text 
Fig. 9. Kinetic energy as function of time for each parameter choice and viscosity model. An increase in ν damps the energy released during the KHI in the switching cases and totally suppresses the KHI in most isotropic cases. (a) ν = 10^{−5}; η = 10^{−3}. (b) ν = 10^{−4}; η = 10^{−3}. (c) ν = 10^{−3}; η = 10^{−3}. (d) ν = 10^{−5}; η = 10^{−4}. (e) ν = 10^{−4}; η = 10^{−4}. (f) ν = 10^{−3}; η = 10^{−4}. 

In the text 
Fig. 10. Internal energy, total viscous, and total ohmic heating as functions of viscosity ν at t = 12.5, before the onset of any null collapse, for all parameter choices. Isotropic (blue, solid) and switching viscosity (orange, dashed) are both shown. An upwardsfacing triangle denotes the higher value of η = 10^{−3} and a downwardsfacing triangle, η = 10^{−4}. The anisotropic viscous heating can become significant for larger values of ν yet, when ν is smaller, the lack of viscous heating is compensated by enhanced ohmic heating. 

In the text 
Fig. 11. Ohmic heating contributions from separate current structures in the spine and fan. We show the mean ohmic heating contributions from the spines (dotted) and currentvortex sheet (dashed) and their sum (solid) for η = 10^{−3} (a) and η = 10^{−4} (b) in only the isotropic cases at t = 10. The value of η dictates how rapidly the balance of contributions shifts from fan to spine with ν, resulting in different trends in total ohmic heating. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.