© J. Quinn et al. 2021

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 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 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 current-vortex 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 spine-fan reconnection and collapse.

The KHI has been well-studied 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 current-vortex 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 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 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 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 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 high-resolution 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 high-resolution 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 non-dimensionalised visco-resistive MHD equations,

(1)

(2)

(3)

(4)

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,

(5)

where the specific heat capacity ratio is given by γ = 5/3. The terms Q_ν = σ : ∇u and Q_η = η|ȷ|² are the viscous and ohmic heating contributions, respectively. The material derivative is D/Dt = ∂/∂t + (u ⋅ ∇).

The non-dimensionalisation scheme is identical to that used in the code Lare3d (Arber et al. 2001), where a typical magnetic field strength B₀, density ρ₀, and length scale L₀ are chosen and the other variables non-dimensionalised appropriately. Velocity and time are non-dimensionalised using the Alfvén speed and Alfvén crossing time t_A = L₀/u_A, respectively. Temperature is non-dimensionalised 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 non-dimensional variables by their respective reference value (e.g. B_dim = B₀B). The reference values used here are B₀ = 5 × 10⁻³ T, L₀ = 1 Mm, and ρ₀ = 1.67 × 10⁻¹² kg m⁻³, giving reference values for the Alfvén speed of u_A = 3.45 Mm s⁻¹, Alfvén time t_A = 0.29 s, and temperature T₀ = 1.73 × 10⁹ 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,

(6)

where ν is the viscous parameter (termed the viscosity throughout this paper) and W is the rate-of-strain tensor,

(7)

In the anisotropic case, we use the switching viscosity model of MacTaggart et al. (2017),

(8)

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,

(9)

where

(10)

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

(11)

where n is the ion number density. For typical active region conditions, T = 2 × 10⁶ K and n = 3 × 10³ m⁻³ giving τ = 0.773 s and α ≈ 10⁸. 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 sub-grid 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⁸. 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 z-axis is determined by an imposed initial magnetic field given in non-dimensional units by

(12)

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⁹ 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 x-boundaries,

(13)

and similarly, the y- and z-derivatives are zero on their respective boundaries. On the x- and y-boundaries, the velocity and its (normal) derivative are zero, and on the z-boundaries 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

(14)

Here u_r(r) describes the radial profile of the twisting motion in terms of the radius r² = x² + y²,

(15)

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

(16)

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₀. The parameters used in all simulations are u₀ = 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 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¹. The fast mode Mach number is given by

(17)

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

(18)

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

(19)

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 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 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 high-resolution simulations is presented, both performed using resistivity of η = 10⁻⁴ and viscosity ν = 10⁻⁴, 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.

4.1. Evolution of a typical case

We detail the evolution of the high-resolution, typical cases in stages, exploring the formation, stability, and breakup of the current-vortex 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 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 (Fig. 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 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 current-vortex 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 current-vortex 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 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 (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 out-of-plane velocity uz 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 current-vortex sheet

As expected from the linear stability measures shown in Figs. 2b and 3 shows the development of the out-of-plane velocity from t = 2 to 10 and reveals that the current-vortex 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 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 (Fig. 4a) and the resultant rolls create strong, small-scale current sheets, enhancing the local reconnection rate.

Fig. 4. Breakup of the current-vortex 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 current-vortex sheet remains stable in the isotropic case while that in the switching case has been fragmented by the KHI. The resultant small-scale 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 current-vortex 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 small-scale current sheets created by the KHI, and the enhanced reconnection far from the null can be seen as ripple-like structures at the fringes of the plot.

4.1.3. 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 that generates a magnetic shear, permitting spine-fan 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 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 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 spine-fan 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: |ux(x)| − |ux(−x)|, the difference in ux 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 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 (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 spine-fan 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 KHI-unstable current-vortex sheet as presented in detail previously (Fig. 7a). The initial injection of the Alfvén waves and formation of the current-vortex 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 small-scale 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 current-enhancing 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 current-vortex 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⁻⁵, 10⁻⁴, or 10⁻³ and η takes values of 10⁻⁴ or 10⁻³. 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 high-resolution 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⁻³ 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 current-vortex 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⁻⁵, 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⁻³, 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 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; in these the KHI is directly observed in the out-of-plane 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⁻³, ν = 10⁻⁴ 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.

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 KHI-unstable 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⁻⁴ 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⁻³ and ν = 10⁻⁵ (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⁻⁴ and ν = 10⁻⁵, does show some growth but it is notably less than the fully unstable case (Fig. 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 that acts to stabilise the sheet. This conclusion can similarly be drawn for the outlying switching case where η = 10⁻⁴ and ν = 10⁻³ (Fig. 9f). The remaining case of interest is where η = 10⁻³ and ν = 10⁻⁴, 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⁻⁴. 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 upwards-facing 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.

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⁻³ 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 current-vortex 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 non-trivial and results in the observed difference in trends.

Figure 11 reveals how the contributions from the current-vortex 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.

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 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.

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 Fig. 10b. When η = 10⁻⁴, 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⁻³ (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 current-vortex 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 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 current-vortex 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 current-vortex 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 spine-fan 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 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.

We also explored the effect of the form of viscosity on the collapse of the null point in the two high-resolution 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 spine-aligned 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⁻⁵ and η = 10⁻³) 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 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.

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 current-vortex 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 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 (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 small-scale 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 small-scale 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 kink-unstable 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 build-up 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.

Acknowledgments

Results were obtained using the ARCHIE-WeSt High Performance Computer (www.archie-west.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

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², 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³.

All simulations were performed on a single, multi-core 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,

(A.1)

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 second-order Runge-Kutta scheme to iteratively calculate the discrete positions x_i along a field line passing through some seed position x₀,

(A.2)

(A.3)

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,

(A.4)

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, second-order Runge-Kutta, 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 current-vortex 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 out-of-plane velocity uz 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 current-vortex 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 current-vortex sheet remains stable in the isotropic case while that in the switching case has been fragmented by the KHI. The resultant small-scale 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 spine-fan 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: |ux(x)| − |ux(−x)|, the difference in ux 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 upwards-facing 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.

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 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.

In the text