Kelvin-Helmholtz instability and collapse of a twisted magnetic null point with anisotropic viscosity

Context: Magnetic null points are associated with high-energy 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 Kelvin-Helmholtz instability (KHI) within its fan plane and, under continued twisting, can instigate spine-fan reconnection and an associated collapse of the null point. Aim: 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: Simulations were performed using the 3D magnetohydrodynamics code Lare3d with a custom anisotropic viscosity module. A pair of high resolution simulations was performed, one using isotropic viscosity and another using anisotropic viscosity, keeping all other factors identical, and the results analysed 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 Kelvin-Helmholtz instability 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, leading to a greater overall heating rate when using anisotropic viscosity. The collapse of the null point occurs significantly sooner when anisotropic viscosity is employed.


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 Kelvin-Helmholtz 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 to also 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 current-vortex sheet which forms in the fan plane and can be unstable to the KHI, given appropriate parameter choices. All experiments are carried out using both isotropic and anisotropic viscosity over a range of parameter choices. Anisotropic viscosity is modeled following MacTaggart et al. (2017). Continued driving after the moment at which the KHI occurs causes the null point to spontaneously undergo spine-fan reconnection and collapse.
In general, the effect of a magnetic field is stabilising; when the wavevector 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 current-vortex sheet, where a velocity shear coincides with a magnetic shear, the balance of shear layer strength and thickness dictates if the KHI, the tearing instability, or some mixture of them, 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 current-vortex 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 sug-Article number, page 1 of 12 arXiv:2102.00761v1 [astro-ph.SR] 1 Feb 2021 A&A proofs: manuscript no. paper gest 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 coinciding 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 high-energy 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 (Moreno-Insertis & 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 much-studied form of reconnection in 3D null points is spine-fan reconnection, where a strong current sheet forms in the vicinity of the null point and enables efficient reconnection between the magnetic field making 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. The numerical setup of the simulations is presented in section 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 section 3. In the first part of section 4 the results of a high-resolution pair of simulations for a single choice of viscosity and resistivity parameters are presented and the effects of the viscosity model used are compared. In the second part, the results of a parameter study are described, generalising the high-resolution results. The chapter concludes with a discussion of findings in section 5 and conclusions in section 6.

Governing equations
We consider the non-dimensionalised visco-resistive 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. D/Dt = ∂/∂t + (u · ∇) is the material derivative. The non-dimensionalisation 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 non-dimensionalised appropriately. Velocity and time are non-dimensionalised using the Alfvén speed u A = B 0 / √ ρ 0 µ 0 and Alfvén crossing time t A = L 0 /u A , respectively. Temperature is non-dimensionalised via T 0 = u 2 Am /k B , where k B is the Boltzmann constant andm is the average mass of ions, here taken to bem = 1.2m p (a mass typical for the solar corona) where m p is the proton mass. Dimensional quantities can be recovered by multiplying the non-dimensional 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 u A = 3.45 Mm s −1 , Alfvén time t A = 0.29 s and temperature T 0 = 1.73 × 10 9 K.

Models of viscosity
Both isotropic and anisotropic models of viscosity are used and compared here. 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 rate-of-strain tensor, In the anisotropic case, the switching viscosity model of Mac-Taggart et al. (2017) is used, where b = B/|B| is the unit vector in the direction of the magnetic field and s(α|B|) is the switching function, 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 to better model viscosity in the vicinity of magnetic null points (Hollweg 1986 where f (a) = 6 5 a 2 + 2.23 a 4 + 4.03a 2 + 2.23 , 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 ion-ion 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 Fig. 1: Field configuration after 4 Alfvén times. The driver speed is also shown as a slice, where white indicates peak driving speed.
and α ≈ 10 8 . Were this value to be used in (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 sub-grid scale at the resolutions used here. To properly resolve the region of isotropic viscosity, we choose α = 12.

Null point model
The magnetic structure of the null point with the spine parallel to the z-axis is determined by an imposed initial magnetic field given in non-dimensional 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. 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. 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 figure 1. The imposed acceleration of the twisting motion is described by u t (t) with 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 footpoints of the spine of the null point, dragging the field and introducing twist throughout the entire null point (figure 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 which acts as a perturbation. As in Wyper & Pontin (2013), there is no prescribed perturbation which 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 9 parameter choices. To limit the time required to complete the study, a relatively low resolution of 320 grid points in each direction is used for these runs. A single, higher-resolution pair of simulations were run, one for each viscosity model, at the 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.

Stability measures
Following Wyper & Pontin (2013), two quantities are used in understanding the stability of the current-vortex 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 current-vortex 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 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 and 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 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 current-vortex 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 current-vortex 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 are 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 .

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, the reconnection rate local to a given magnetic field line is calculated 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).

Results
In the first subsection of the results section, the evolution of a pair of high-resolution 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 current-vortex 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 high-resolution 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.

Evolution of a typical case
The evolution of the high-resolution, typical cases is detailed in stages, exploring the formation, stability and breakup of the current-vortex sheet before investigating the collapse of the null point. Then, the evolution is summarised through an analysis of the energy budget and the reconnection rate in time.

Formation of the current-vortex 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 (figure 2a). These shear layers are jointly called the current-vortex sheet. Without any diffusion in the system the waves would travel far along the fan plane be- Rings of vorticity and current density and associated linear stability criteria. Subfigure (a) plots the vorticity and current density for both viscosity models at t = 3 and z = 0. Subfigure (b) plots the 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 current-vortex sheet which 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). fore 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 current-vortex sheet forms. The hole in the sheet is due to magnetic tension forces opposing the twisting motion, as illustrated in figure 3 of Wyper & Pontin (2013). This also gives rise to counterflows (not shown) similar to those seen in Wyper & Pontin (2013) and Galsgaard (2003).
In the switching case, the reduced effective viscosity produces a vortex ring which 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 figure 3b). The measure Λ confirms that the current-vortex 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 (figure 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.

Disruption of the current-vortex sheet
As expected from the linear stability measures shown in figure 2b, figure 3 shows the development of the out-of-plane velocity from t = 2 to 10 and reveals that the current-vortex sheet in only the switching case is unstable to the KHI. Both cases exhibit a similar morphology of u z , despite a growing difference in strength, until t = 6 when the KHI appears only in the switching case, initially along the diagonals (figure 3b) before spreading azimuthally (figure 3c). There is no evidence of the KHI in the isotropic case.
In both cases the current-vortex 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 current-vortex sheet (figure 4a) and the resultant rolls create strong, small-scale current sheets, enhancing the local reconnection rate. 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. The reconnection rate is greatest close to the origin, 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 which spiral outwards from the origin. The switching case shows a greater peak reconnection rate due to the small scale current sheets created by the KHI, and the enhanced reconnection far from the null can be seen as ripple-like structures in the fringes of the plot.

Spine-fan 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 which generates a magnetic shear, permitting spine-fan reconnection. The results of the isotropic case are presented first, in detail, then the effect of the KHI is explored in the switching case.
In typical studies of spine-fan 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 which 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 which 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.
The twist in the spines creates a current which 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 (figure 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 ( figure 6b).
The velocity shear around the null point shears the magnetic field accordingly, creating a current sheet through the null point ( figure 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 (figure 7a). The development of the current sheet and the resultant spine-fan reconnection is similar to that of Pontin et al. (2007) with the exception that the twist in the field unravels as the reconnection proceeds. 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.

Energy budget and reconnection rate
The kinetic energy in the switching case is a measure of the main evolution of a KHI-unstable current-vortex sheet as presented in detail previously (figure 7a). The initial injection of the Alfvén waves and formation of the current-vortex sheet can be seen at 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 small-scale current sheets form and Ohmic heating begins to drain energy from the instability (figure 7e). This is also reflected in the reconnection rate (figure 7f) 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 (figure 7a). The ohmic heating is similarly damped without the influence of the KHI (figure 7e). This results in the switching model extracting more energy from the field (figure 7b) and heating the plasma much more effectively (figure 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 (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 (figure 7f). One Article number, page 5 of 12 A&A proofs: manuscript no. paper (c) t = 10 Fig. 3: Development of the KHI in the out-of-plane velocity u z at t = 2, 6 and 10 for both viscosity models. Note the isotropic results have been multiplied by as much as 1000 in order to compare 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.
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 current-enhancing effect of the switching model or 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 has broken up the current-vortex sheet, introduced inhomogeneities throughout the fan plane and generated small current sheets. This results in a collapse which 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.

Study of parameter dependencies
The results shown in section 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. The simulations were performed 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 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 more unstable KHI. 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 current-vortex sheet. The condition on Λ for instability is Λ > 1. All layers with the switching model are linearly unstable to the KHI (figures 8d and 8b) 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. The two middle cases, when ν = 10 −4 show stability when η = 10 −4 and instability when η = 10 −3 . These findings can be attributed to the way in which diffusion affects the magnitude and thickness of the current-vortex sheet.

Shear layer stability
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, resulting 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 current-vortex sheet to the KHI in each case is determined via inspection of the out-of-plane 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, that is the KHI is directly 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. Figure 9 shows the kinetic energy as a function of time for all parameter choices and for both viscosity models. The strongly KHI-unstable cases show a similar kinetic energy profile to that of the unstable typical case (figure 7a). In the switching cases, the peak kinetic energy is larger when η = 10 −4 in two cases (figures 9d and 9e). This is a result of the reduced diffusion of the magnetic field resulting in a stronger vorticity layer. The isotropic cases in figure 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 (figure 9a) is the only isotropic case where the KHI is significantly unstable. In this case the kinetic energy profile shares a similar   Table 1: Stability in the isotopic and switching cases for different choices of ν and η. Both linear stability ("lin", as predicted by Λ > 1 in figure 8) and observed stability ("obs") are shown. Entries marked as unstable* show growth of the KHI but the growth rate of the perturbation is close to zero. The isotropic model mostly results in stability while the switching model mostly results in instability.

Kinetic energy profiles
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. Table 1 reveals three cases of interest which are explored through the 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 (figure 9d). Given that the current-vortex 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 which acts to stabilise the sheet. This conclusion can similarly be drawn for the outlying switching case where η = 10 −4 and ν = 10 −3 (figure 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 (figure 9b). In this case, the kinetic energy plateaus as the viscosity dissipates kinetic energy as it is generated by the instability. 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 (figure 10a), Q ν generally decreases with decreasing viscosity ν, as one may expect, with the exception of the isotropic cases where the resistvity η = 10 −4 . Instead, in these cases the viscous heating shows little dependence on ν. This reveals the complex, nonlin-ear relationship between viscous heating, the value of ν and the flows generated.

Heating profiles
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 (figure 10b. The stark difference in trends can be explained by considering the spatial distribution of ohmic heating which 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 one here). These are the main sources of ohmic heating and the balance of contributions from each source, for different values of η and ν, is non-trivial 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 xy-plane (representing the heating within the current-vortex sheet) and that in the yz-plane (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.
For any value of η, the spine heating increases and the current-vortex 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 current-vortex 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 figure 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 (figure 11b). The opposite is true when η = 10 −3 (figure 11a).

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. It has been found that, when using a model of viscosity appropriate in the solar corona, i.e. the switching model, 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 spine-fan reconnection specifically has been 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 current-vortex sheets studied here.
The simulations detailed here have been performed with a model of anisotropic viscosity which 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 current-vortex sheet could be accelerated by perpendicular vis-cosity, 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 spine-fan reconnection and null point collapse) is primarily driven by a shearing motion at A&A proofs: manuscript no. paper 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 downwards-facing 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. 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 which 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 null-directed 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.
The effect of the form of viscosity on the collapse of the null point is also explored in the two high-resolution simulations (of section 4.1) and it is 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 spine-aligned flows, so one may conjecture that, in the switching case, the KHI introduces its own asymmetries which 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 back-pressure which halts and reverses the spine-fan reconnection process. The simulations performed here do not run for a long enough time to investigate the generation of back-pressure. It may be that a collapsing twisted null point is unable to produce the required back-pressure if it unravels during its initial collapse. Running the high-resolution simulations reported here for a longer time would reveal if the particular setup studied here can generate oscillatory spine-fan reconnection. Alternatively, using a pre-twisted 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.  Fig. 11: Ohmic heating contributions from separate current structures in the spine and fan. Shown are mean ohmic heating contributions from the spines (dotted) and current-vortex 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.

Conclusions
In this paper two models of viscosity have been applied to a magnetic null point which has been dynamically twisted at its footpoints in such a way that a current-vortex sheet forms in the fan plane. This sheet has the potential to become unstable to the KHI. It was 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 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 spine-directed, pressure-driven flows. This was found to occur sooner when the KHI is present. The general development of the collapse and associated spine-fan 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 current-vortex 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 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 which would be absent if only isotropic viscosity were considered.