© ESO, 2011

1. Introduction

There has been significant progress in recent years towards understanding the properties of magnetic reconnection in three dimensions (3D). In particular, it is now appreciated that the intense current layers necessary for reconnection in astrophysical plasmas may form at a number of different characteristic structures of the magnetic field. One such structure is a 3D magnetic null point. Such null points have been demonstrated to be present in abundance in the solar corona (e.g. Longcope & Parnell 2009), and the same is highly likely to be true in other astrophysical environments such as the coronae of other stars and of accretion disks. Recent results from observations and simulations suggest that reconnection at these nulls may play an important role in the coronal dynamics (e.g. Luoni et al. 2007; Lynch et al. 2008; Pariat et al. 2009; Masson et al. 2009). Note also that the importance of reconnection at 3D nulls is not restircted to astrophysical plasmas, but plays a role both in the Earth’s magnetosphere (e.g. Xiao et al. 2006) and some laboratory plasmas (e.g. Gray et al. 2010).

Recent studies have revealed a number of characteristic modes of reconnection that may occur at such 3D nulls. These have recently been categorised by Priest & Pontin (2009) into “torsional spine reconnection”, “torsional fan reconnection” and “spine-fan reconnection”. Torsional spine and torsional fan reconnection occur when an equilibrium magnetic null point field is disturbed by a rotational perturbation (the rotation being in a plane perpendicular to the spine). Spine-fan reconnection occurs when a shear perturbation is applied that disturbs the locations of the spine and fan – this leads to a localised current layer forming around the null and flux transport across the fan separatrix surface.

Close to the null the magnetic field may be written (1)and the eigenvalues and eigenvectors of the matrix ∇B determine the locations of the spine and fan of the null (see e.g. Parnell et al. 1996). Previous studies of the generation of current layers at 3D nulls due to rotational motions have considered only the case where the background equilibrium null point has a rotational symmetry, i.e. in which the two eigenvalues associated with the fan are equal. During the resulting evolution in which torsional spine and torsional fan reconnection takes place, the respective current layers, plasma flows, etc., have exhibited azimuthal symmetry due to the azimuthally symmetric background magnetic fields and perturbations.We have previously investigated the effect of varying the symmetry of the background null point field on the spine-fan reconnection mode, and shown that while the topological properties of the reconnection are unchanged, the dimensions and intensity of the current layer, and the reconnection rate, are strongly dependent on the degree of asymmetry (Al-Hachami & Pontin 2010), as well as the relative angle between the shear driving and the strong/weak field directions (Galsgaard & Pontin 2011). In this paper we generalise existing models for torsional spine and torsional fan null point reconnection as follows. We begin in each case by introducing a new kinematic analytical solution for the corresponding reconnection mode in which a localised current layer is present at the null. We then go on to consider the effect of varying the symmetry of the background magnetic field, by varying the ratio of the fan eigenvalues of the null. Lastly we perform numerical simulations of the full system of MHD equations in which null point configurations with varying degrees of symmetry are subjected to rotational disturbances, to complement the analytical models.

The paper is organised as follows. In Sect. 2 we review previous modelling of torsional spine and torsional fan reconnection. In Sect. 3 we introduce our new analytical and numerical models for torsional spine reconnection, and discuss their implications. We go on in Sect. 4 to do the same for torsional fan reconnection. Finally in Sect. 5 we present our conclusions.

2. Existing models for torsional spine and fan reconnection

Both torsional spine and torsional fan reconnection involve the formation of a current layer in which the current vector is directed parallel to the spine line at the null. It was first shown by Pontin et al. (2004) that the corresponding reconnection takes the form of a rotational slippage of magnetic flux threading the non-ideal region. (This is in contrast to the case where the current vector is parallel to the fan, in which case magnetic flux is reconnected across the spine and fan (Pontin et al. 2005).) The magnetic flux undergoes this rotational slippage in response to rotational flows in the ideal region in which the rate of flux transport in the azimuthal direction is different for field lines entering the non-ideal region than it is for field lines exiting the non-ideal region. The original model of Pontin et al. (2004) is based on the magnetic field in cylindrical polar coordinates. This results in a spatially uniform current parallel to the spine (z-axis). In order to obtain a 3D-localised non-ideal region (relevant in astrophysical plasmas which are approximately ideal almost everywhere), it was therefore necessary to impose a localised resistivity profile.

In the following sections we introduce two new analytical kinematic solutions. In these solutions we have been able, for the first time, to include a current density which is fully spatially localised with its peak intensity focussed at the spine or fan. This current localisation allows a spatially uniform resistivity to be used, adding a degree of physical plausibility to the models, since in practice in an astrophysical plasma a localised non-ideal region is associated with a localised current layer. In each of our solutions, the structure chosen for the magnetic field and the resulting current layer is based on the results of resistive MHD simulations in which the dynamic formation of these current layers has been observed. It is worth noting that additional kinematic solutions with a spatially varying current density have been presented by Wyper & Jain (2010). In their solutions the current is spatially localised, but in contrast to our new solutions is focussed away from the null point.

Resistive MHD simulations have demonstrated that the form of the current layer is different depending on whether the rotational perturbation primarily disturbs the fan field lines or field lines around the spine. The perturbation behaves essentially as an Alfvén wave, travelling along the magnetic field lines, which due to the hyperbolic structure of the field around the null leads to the perturbation accumulating either in the vicinity of the spine or the fan. When the fan field lines are subjected to a rotational disturbance, torsional spine reconnection takes place in a tubular current structure that forms at the spine (Rickard & Titov 1996; Pontin & Galsgaard 2007). Within this tube, the magnetic field lines spiral around the spine line. The radius of the tube decreases, and the current intensifies, until the twisting of the field being driven by the perturbation is balanced by rotational slippage facilitated by resistive diffusion. The reconnection rate, defined in 3D as the maximal value of (2)along any field line threading the non-ideal region, then quantifies this rotational slippage.

When field lines in the vicinity of the spine line are disturbed, a current layer forms on the fan surface leading to torsional fan reconnection (Rickard & Titov 1996; Galsgaard et al. 2003). Again field lines spiral within the current layer, whose magnitude intensifies as the twisting of the field is concentrated in an increasingly narrow sheet over the fan surface. Once the sheet becomes sufficiently thin resistive diffusion dissipates the twist leading once again to a rotational slippage of field lines.

3. Torsional spine reconnection

3.1. Torsional spine reconnection: kinematic model

In this section we describe a kinematic model for torsional spine reconnection within a current tube that is localised to the spine of the null point. The form of the magnetic field and resulting current structure is chosen to match behaviour observed in the numerical simulations described by Pontin & Galsgaard (2007). The method of solution is the same as that described by Hornig & Priest (2003) and Pontin et al. (2004). Specifically, we solve the steady-state, kinematic, resistive MHD equations (3)We specify the form of the magnetic field B and solve for the electric field (E) and plasma flow perpendicular to B (v _⊥ ) via (4)(the component of v parallel to B being arbitrary). Note that s is a parameter along magnetic field lines, i.e. the integration of Φ is performed along these field lines. The resistivity η = η₀ is taken to be uniform.

We consider first the rotationally symmetric case. We set (5)where in cylindrical polar coordinates (r,θ,z) where j, a, b are positive real numbers and n, m, μ, ν are positive integers. B_P defines the potential “background” null point component, and B_J defines the non-potential component of B associated with the tubular current structure, which extends to radius r = a and to z =  ± b. In accordance with previous simulation results we assume an extended tube of current aligned to the spine so that b ≫ a. For m = n = μ = ν = 1, B and J are continuous and differentiable. However, in order that all physical quantities in the final solution are continuous and differentiable we find that it is necessary to choose m = 3, μ = n = 2 and ν = 1. The current density is represented by the plot in Fig. 1.

Fig. 1 Left: magnetic field lines for the torsional spine model defined by Eqs. (5), (8), (9) for a = 1, b = 5, j = 3, p = 1.5, q = 1. The shaded surface shows a current isosurface. Right: current vectors in the y = 0 plane for a = 1, b = 2, q = 1.

Fig. 2 Current density (left) and plasma flow vectors (right) for the torsional spine kinematic solution. Plotted in the plane z = 2 for B0 = L0 = j = a = 1, b = 4 and p = q = 1 (above), p = q = 2 (middle) p = q = 5 (below). The grey curve (red online) in the right-hand plots marks the boundary of the non-ideal region.

Parametric equations for the magnetic field lines associated with Eqs. (5)–(7) can be found in a straightforward way by solving ∂X/∂s = B(X(s)). Since the integrand in the first equation in (4) is independent of θ we require only r = r₀exp(B₀s/L₀) and z = z₀exp( − 2B₀s/L₀) (it is straightforward to also obtain an expression for θ(r₀,θ₀,z₀,s)). Solving Eqs. (4) with s = 0 at z =  ± z₀ we find that the qualitative properties of the electric field and flow field are similar to those discovered by Pontin et al. (2004). In particular, rotational plasma flows are still present around the spine axis of the null. Using the freedom of arbitrary flow parallel to B in the model we can for illustrative purposes choose to add a component v_∥ such that v_z = 0. We then see a purely azimuthal flow. Due to the return current region close to r = a (see Fig. 1) the sense of rotation of the plasma changes at an intermediate radius, as shown in the top-right frame of Fig. 2.

We now investigate how the properties of the solution vary when the rotational symmetry of the above system is broken. When the rotational symmetry is lost it is no longer possible to find closed-form expressions for the field line mapping. We therefore numerically integrate B to find field lines and solve Eqs. (4) on a rectangular grid. We may break the symmetry either in the potential component B_P defining the magnetic null or in the component B_J defining the current tube. Our new potential component of the magnetic field is given by (8)in Cartesian coordinates where p > 0 is a parameter. As p varies the magnetic field along the spine direction is fixed while the ratio between the fan eigenvalues (associated with the eigenvectors along the  and  directions) varies (see Fig. 1). We choose to break the symmetry in B_J by converting to Cartesian coordinates and setting (9)where R² = x² + qy² (note that this reduces to expression (7) when q = 1). This has the effect of distorting the current into a cylinder with elliptical cross-section, with major and minor axes along the x- and y-axes, extending to x =  ± a, (see the images on the left in Fig. 2). Pre-empting the results of the following section, we present here results for p = q, such that as p increases the current tube narrows along the direction associated with the large fan eigenvalue, i.e. the strong field direction in the fan. We set B₀ = L₀ = η₀ = j = 1, a = 1, b = 4, and solve Eqs. (4) on a rectangular grid with 81 gridpoints in each direction covering the volume  −2 < x,y < 2,0 < z < 4 with the solution being symmetric about z = 0. We restrict our attention to the range p ≥ 1, which simply selects the  direction as the strong field direction in the fan.

The results of the above analysis are presented in Figs. 2–4. As p is increased, the current tube shrinks in the y-direction, with the dominant current component J_z intensifying in the part of the tube close to the y-axis (i.e. the direction of the short axis of the ellipse) – see Fig. 2. The stronger current in this region results in an enhanced plasma flow speed. The direction of the flow is also distorted from the circular pattern at p = 1, but continues to flow on closed elliptical paths around the spine (z-axis). As the fan plane is approached the radii of the elliptical shells of positive and negative azimuthal flow increase, owing to the hyperbolic nature of the field structure – see Fig. 3.

Fig. 3 Intensity map of vx in the x = 0 plane for the kinematic torsional spine solution with p = q = 1.

In order to determine the reconnection rate we calculate Ψ as defined in Eq. (2). Due to the breaking of the symmetry it is no longer clear that the maximal value of Ψ should occur along the spine field line, as was found in previous studies (note that the current modulus has maximum value away from the spine for large p). However, it turns out that indeed the maximum occurs along field lines asymptotically close to the spine for all p. Figure 4 displays the peak value of the current density (which we impose) and the reconnection rate as a function of the degree of asymmetry. It is clear that the peak current scales linearly with p( = q) and that correspondingly the reconnection rate scales linearly with p. Note however that all of the above solutions were obtained with a fixed value of the parameter j which also contributes to controlling the peak current density, and that the velocity and reconnection rate will increase proportional to this parameter. As we shall see below, however, we do not have a linear increase of the peak current with p in our simulations, and it may therefore be more realistic to take j = j(p).

Fig. 4 Dependence on the anisotropy parameter q = p of the maximum values of |J| and the reconnection rate Ψmax for the kinematic torsional spine model.

3.2. Torsional spine reconnection: simulations

To complement the above kinematic model we now describe the results of numerical simulations of the full system of resistive MHD equations. The code has been extensively used for other similar simulations, the details of the scheme being described by Nordlund & Galsgaard (1997). The MHD equations are non-dimensionalised by setting the magnetic permeability μ₀ = 1 and the gas constant equal to the mean molecular weight. The result is that for a volume in which |ρ| = |B| = 1, time units are such that an Alfvén wave would travel one space unit in one unit of time. The equations are solved on a numerical grid of 256³ gridpoints over  [x,y,z] ∈  [ ± 0.75, ± 0.75, ± 3] with uniform η = 10^-4. The grid is stretched to give higher resolution around the null, where Δx = Δy = 0.0034, Δz = 0.014. All boundaries are closed and line-tied. We repeat the simulations described by Pontin & Galsgaard (2007), in which a localised rotational perturbation of the magnetic field is imposed on a background equilibrium null point. Specifically, we begin with a potential magnetic field given by Eq. (8), plasma velocity zero everywhere, and initialise the plasma density and thermal energy to be spatially uniform with values ρ = 1 and e = 0.025, respectively. As a result there exist regions of both high and low plasma-β in the domain. β approaches infinity as one approaches the null, and takes a value of 1 on an ellipsoidal surface which cuts the coordinate axes at x =  ± 0.09(p + 1), y =  ± 0.09(p + 1)/p, z =  ± 0.09.

In addition to the above, we impose at t = 0 a magnetic field perturbation composed of a ring of magnetic flux centred on the null point and lying in the fan plane: (10)with b₀ = 0.05, r₀ = 0.18, α = 0.08, ζ = 0.06 (see Figs. 1 and 2 of Pontin & Galsgaard 2007). (We have also performed simulations where this initial perturbation is elliptical rather than circular, but found no change to the qualitative results – thus here we confine our discussion to the circular perturbation.)

For t > 0 the perturbation splits, with wavefronts travelling both toward and away from the null. We focus on the behaviour of the ingoing pulse, which for p = 1 gradually stretches out to form a cylindrical tube of intense current around the spine. When p ≠ 1, the azimuthal symmetry of the perturbation wavefront is broken as soon as the evolution begins. The Alfvén speed in the radial direction now depends on x and y, and so the wavefront travels toward the spine faster along the y-axis where |B| is stronger. Propagation in the z-direction is essentially unaffected, and the current distribution associated with the perturbation forms into a cylinder with elliptical cross-section, whose length and eccentricity both increase as the pulse steepens towards the spine and the peak current correspondingly intensifies (see Figs. 5, 6).

Fig. 5 Time sequence showing contours of Jz in the z = 0 (left) and x = 0 (right) planes, for the torsional spine simulation with p = 2.

Fig. 6 Above: current isosurface at 50% of maximum value, for the torsional spine simulations with p = 1 (left), p = 2 (middle) and p = 5 (right). Taken in each case at the time when the positive value of Jz reaches a maximum. Below: contours of Jz in the z = 0 (fan) plane over  [x,y] ∈  [ ± 0.2, ± 0.2] for the same runs at the same times.

Fig. 7 Evolution of the peak value of |J| for the different torsional spine simulations: p = 1 (solid line), p = 2 (dotted), p = 3 (dashed), p = 5 (triple-dot-dashed) and p = 10 (dot-dashed).

As shown in Fig. 7, there is a strong (at least 100%) increase in the peak current density as the disturbance reaches the spine. It is also clear from the plot that this occurs earlier for larger p due to the increased speed of propagation along the y-axis. In general, the peak current is higher for larger p. However, this is complicated by a competing effect – namely that there are two distinct maxima in |J|_max during the simulations, one corresponding to the localisation of the leading edge of the pulse and the other to the trailing edge of the pulse. As can be readily seen in Fig. 5, they correspond to opposite signs of J_z. Thus, while the maximum positive value of J_z strictly increases as p increases, it is found that for p = 1,2 the trailing edge of the pulse (corresponding to J_z < 0) dominates. The dimensions of the current layer are shown in the final three columns of Table 1. They are measured at a time corresponding to the localisation of the leading edge of the pulse, i.e. the time when J_z > 0 reaches its temporal maximum (which occurs at a time marked by the peak on each curve in Fig. 7 between t = 1 and t = 2, except for the run with p = 1 in which J_z reaches its maximum at t = 2.1). The dimensions of the current tube in the xy-plane demonstrate its elliptical nature (centred on the spine) with eccentricity increasing with p as shown in Fig. 6. The length in the z-direction (L_z) decreases with increasing p. This is because the z-component of B_P (and thus the propagation speed in z) is independent of p, and because the pulse localises in the xy-plane at an earlier time for larger p as explained above. It is also worth noting from Fig. 6 that for p ≠ 1 the maximum current is attained not exactly on the spine, but in two locations displaced symmetrically from the spine along the x-axis.

Fig. 8 Evolution of the reconnection rate for the different torsional spine simulations: p = 1 (solid line), p = 2 (dotted), p = 3 (dashed), p = 5 (triple-dot-dashed) and p = 10 (dot-dashed).

We turn now to consider the reconnection rate, calculated as described in Eq. (2). The maximal value of Ψ is found over all field lines that thread the current layer, these being field lines that pass close to the null and its spine and fan. For p = 1, due to symmetry, the maximum can be found on any field line which runs asymptotically close to the null. However, for p > 1 the current is strongest in the weak field regions around the x-axis, and so field lines that pass through these regions attain the highest values of Ψ. The evolution of the reconnection rate is plotted for the different simulations in Fig. 8. The maximum value attained does not depend strongly on p (see the fourth column of Table 1), except that it is significantly lower for p = 1. However, this is likely down to the fact that for p = 1 the disturbance extends all the way to the z-boundaries, and so we “miss” some of the length of the current layer. Note that the reconnection rate calculated here is a net effect of integrating through regions of both positive and negative E·B – we have different senses of reconnection (rotational slippage) occurring on the leading and trailing edges of the pulse, and here we measure the net effect.

3.3. Torsional spine reconnection: discussion

The numerical simulations discussed above demonstrate the localisation of a rotational perturbation towards the spine of a non-symmetric linear 3D null. The resulting current intensification is associated with torsional spine magnetic reconnection. The current tube that forms around the spine has a structure that is closely matched by the kinematic steady-state model described in Sect. 3.1. In particular, the current is dominated by the component parallel to the spine (J_z), and is localised within a tube of elliptical cross-section. The short axis of the ellipse is aligned with the strong field direction in the fan plane, along which the current is most intense. The eccentricity of the ellipse increases as the magnetic field asymmetry increases. This dynamic evolution should be compared with the stationary analytical solution at the stage where the current layer has become fully localised on the spine, as in the analytical solution – i.e. when the maximum (or minimum) of J_z is achieved. If it were feasible to drive the perturbation from the boundaries in a steady fashion, we would expect this spine-localised current tube to persist in a quasi-steady state for a period controlled by the period of steady driving, as in Sect. 4 below – see Fig. 13 (driving a rotational motion in a smooth way around the fan boundaries in our Cartesian domain is not practical).

The plasma flow also has a similar qualitative structure in the kinematic model and simulations. This structure is that of a rotation along elliptical paths around the spine, with these elliptical paths closely following the current density contours. One difference between the model and simulations is that both signs of rotational flow are only seen during a certain period in the simulations. Specifically, as the twist associated with the perturbation propagates towards the null it drives flow predominantly in the positive rotational sense, and when the reconnection process is completed (t ≳ 3.5) the field then “untwists” leading to a large-scale flow in the opposite direction. It is only approximately between the times that J_z reaches its maximum positive and negative values (approximately 1.5 < t < 3, see Fig. 7) that the Lorentz force accelerates the plasma in opposite rotational senses at different distances from the null. By contrast, in the steady model where momentum balance is neglected, rotational plasma flows of both senses are required to maintain the steady state.

In the numerical simulations it is found that the peak reconnection rate is approximately independent of p. We can understand why this should be the case by noting that, while J_z increases with increasing p, the length of the current layer along z decreases with increasing p. Thus for larger p we have a larger integrand in Eq. (2), but it is non-zero over a shorter distance. Note also that the reconnection rate is found to be given by the integral of E _∥  along field lines lying asymptotically close to the spine. Since J is dominated by J_z and the current layer has minimal extent in the xy-plane, then it is natural that the reconnection rate should not depend strongly on the B_xy field components away from the null. In order to match these simulation results for the p-dependence of the reconnection rate, the magnetic field B_J in our kinematic model could be normalised by a factor proportional the current modulus.

It is worth considering here the interpretation of the reconnection rate as measured above. Within the framework of general magnetic reconnection (Schindler et al. 1988), the flux reconnected is measured by calculating the differential transport of flux on opposite sides of the non-ideal region, for a flux surface passing through the peak of the pseudo-potential Ψ (see also Hesse et al. 2005). This theory is based on the assumption that no nulls are present within the non-ideal region, while here of course we consider a non-ideal region containing a null point. In the rotationally symmetric case the maximal value of |Ψ| is found on any field line passing asymptotically close to the null. This measures the differential rate of flux transport through any flux surface bounded by the spine and fan in the half-space (above or below the fan), by symmetry (Pontin et al. 2004; Pontin 2011). That is, it measures the net magnetic flux reconnected at the given time. However, in the present simulations for p ≠ 1 the value of Ψ no longer exhibits azimuthal symmetry. Therefore the maximal value of |Ψ| obtained above measures the maximum differential rate of flux transport for any flux surface bounded by the spine and fan, though one should bear in mind that this now varies around the spine. Of course in a steady-state situation, this differential rate of flux transport must be the same for all such flux surfaces (by continuity), and indeed it is found to be in the steady-state kinematic solution.

The behaviour of the disturbance – spreading out along the spine in all of the simulations – suggests that the dominant wave mode is an Alfvén wave. However, it is highly likely that other wave modes are present. Owing to the differing waves speeds approaching the null for different directions, one may speculate that if the disturbance were sufficiently strong as to be considered non-linear, then some of the effects discussed by McLaughlin et al. (2009) would be present. Those authors considered the propagation of a non-linear fast-mode wave towards a 2D X-point, and observed the formation of cusp-shaped structures as the wavefront collapsed, and the formation of both fast and slow mode shock waves. However, since our simulation is three-dimensional, it is not feasible at present to use the numerical resolution required to properly resolve such features. Note further that the perturbation in our simulations is initialised in the vicinity of the β = 1 surface. It is known that wave-mode conversion may take place in the vicinity of this surface. Thus the propagation of the disturbance towards the null, as well as the final form of the current layer, may be different if we were to repeat the series of simulations with different values of the plasma-β.

Fig. 9 Above: magnetic field lines for the torsional fan model defined by Eqs. (5), (8), (12) for a = 5, b = 1, j = 50, p = 1.5, q = 1. The shaded surface shows a current isosurface. Below: current vectors in the y = 0 plane for a = 3, b = 1, q = 1.

4. Torsional fan reconnection

4.1. Torsional fan reconnection: kinematic model

We now turn our attention to modelling of the torsional fan reconnection mode, which involves rotational slippage of field lines in a current layer localised around the fan surface. We proceed to solve Eqs. (3) in the same way as described in Sect. 3.1. Again, we first analyse a model for the cylindrically symmetric case, in which for the first time a localised current layer is included in the fan plane. The structure of the magnetic field is chosen by comparing with the numerical simulations of Galsgaard et al. (2003); Pontin & Galsgaard (2007). We again construct our magnetic field as the sum of a potential part (B_P) and non-potential part (B_J) as in Eq. (5), with B_P =  [r,0, − 2z], and this time (11)see Fig. 9. Note that B_θ is now odd in z, and since we are modelling a current layer focused on the fan plane we assume that b ≪ a. We again choose the integers m,n,μ,ν in such a way that all physical quantities in our solution are continuous and differentiable, specifically m = 3,μ = n = 2,ν = 6. The field line mapping for B is again given by r = r₀exp(B₀s/L₀) and z = z₀exp( − 2B₀s/L₀), along with a lengthy expression for θ(r₀,θ₀,z₀,s) which is not required to obtain the solution. This time we solve Eqs. (4) setting s = 0 at r = r₀ = a. We choose to do this because it leads to non-zero flow for |z| > |b|, which is consistent with the observed result from the simulations that this reconnection mode is set up by rotational driving flows in the regions around the spine footpoints.

Fig. 10 Current density (left) and plasma flow vectors (right) in the kinematic torsional fan model. Plotted for B0 = L0 = j = b = 1, a = 4 and p = q = 1 (above), p = q = 2 (middle) p = q = 5 (below). The grey curve (red online) in the right-hand plots marks the boundary of the non-ideal region. The current density is plotted at z = 0, while the plasma flow is plotted in the plane z = 0.5.

The new torsional fan solution is represented in Figs. 10, 11. The plasma flow has a rotational structure, as in the solution of Pontin et al. (2004). That is, when we subtract a component of v parallel to B such that v_z = 0, then the remaining flow is non-zero only in the azimuthal direction (see the top-right frame of Fig. 10). Thus, field lines traced from comoving footpoints in the ideal region at z > z₀ (or z <  − z₀) rotate around the spine at a fixed radius, while field lines traced from the ideal region at r > r₀ remain fixed (v = 0 there), and we have a change of connectivity in the form of a counter-rotational slippage. Owing to the fact that J_r changes sign for different levels of z, the rotational flow within the current layer has regions of both clockwise and anti-clockwise rotation, in a similar way to the torsional spine model discussed above.

As before, we now investigate the dependence of the properties of our solution on the symmetry of the magnetic field. We proceed as in Sect. 3.1 to solve Eqs. (4) using the semi-analytical method described there. 81³ gridpoints are used over the domain  − 4 ≤ x,y ≤ 4, 0 ≤ z ≤ 2, and we use parameters B₀ = L₀ = j = η₀ = 1, b = 1,a = 4.

We take the potential component of our magnetic field (B_P) to be given by Eq. (8), with the non-potential component taken to be (12)where R² = qx² + y² (which reduces to Eq. (11) when q = 1). The current layer now has the shape of an elliptical disc, with major and minor axes along the x- and y-axes, extending to , y =  ± a. Pre-empting the results of the following section, we present here results for p = q, such that as p increases the current disc shrinks along the direction associated with the small fan eigenvalue, i.e. the weak field direction in the fan. As shown in Fig. 10, the current density is enhanced in the regions around the short axis of the ellipse. As in the torsional spine solution, if we set the parallel flow in such a way as to eliminate v_z, then the plasma flow in the xy-plane follows closed elliptical paths, being strongest in magnitude where the current is enhanced.

Examining the dependence of the peak current and reconnection rate on the degree of asymmetry, we find that both increase with increasing q, as shown in Fig. 11.

Fig. 11 Dependence on the anisotropy parameter q = p of the maximum values of |J| and the reconnection rate Ψmax for the kinematic torsional fan model.

4.2. Torsional fan reconnection: simulations

We now perform numerical simulations similar to those described in Sect. 3.2, designed to investigate the effect of the background field asymmetry on the torsional fan reconnection mode (see Table 2). The setup of the simulations is the same as before, except that this time we take 256³ gridpoints over  [x,y,z] ∈  [±2.5, ± 2.5, ± 0.5]. This time we perturb the initial equilibrium (Eq. (8)) by applying a rotational driving velocity on the boundaries around the spine footpoints both above and below the null, as in Galsgaard et al. (2003). The sense of rotation is opposite above and below the fan plane. Specifically, we apply the following azimuthal velocity in the z =  ± 0.5 planes; (13)where , C = 10, τ = 1.6 and A =  ∓ 0.1 at z =  ± 0.5. We perturb the system using a different method to Sect. 3.2 for the following reasons. First, the rotation of the fan plane (Sect. 3.2) is done with an internal perturbation within the domain since driving from the boundary in this way in our Cartesian simulation domain would involve driving around the corners of the domain which is difficult to do in a smooth and stable way. Second, we choose to drive from the boundary in this section since applying an internal perturbation it is difficult to develop a significant current increase at the fan. This is because owing to the spatially uniform resistivity (which we cannot significantly reduce due to the limitations of numerical resolution) an internal perturbation (which must be well localised initially) suffers significant diffusion before reaching the fan, as discussed by Pontin & Galsgaard (2007). Note that this is not such an issue for the case where the fan field lines are rotated (Sect. 3.2) since in that case the incoming pulse spreads only in one direction (z) owing to the hyperbolic nature of the field while it contracts in two (x,y) – however, for the case of a spine rotation the in-coming pulse spreads in two directions (x,y) and contracts only in one (z).

Fig. 12 Contours of Jy (left) and Jz (right) in the x = 0 plane over  [y,z] ∈  [±2.5, ± 0.5] for t = 0.8 (top), t = 1.8 (middle) and t = 4.2 (bottom), for the torsional fan simulation with p = 1.

For the symmetric case p = 1, the evolution of the system following the initiation of the driving velocity is described in detail by Galsgaard et al. (2003). A disturbance that is dominated by a torsional Alfvén wave is launched from the driving z boundaries towards the null. The wave front spreads along the x- and y-directions travelling along the field lines, with its velocity in the z-direction being independent of x and y. The wavefront steepens as it approaches the fan surface due to the hyperbolic structure of the field, eventually forming a planar current layer in the fan. This process is demonstrated in Fig. 12 (see also the left-hand image in Fig. 14), from which one can see that the current associated with the current layer is oriented parallel to the fan surface and flows radially inwards toward the null (note that in this x = 0 plane, J_x is zero by symmetry). It is worth noting that there also remains a large distributed current near the boundaries where the twisting was applied – observed as strong concentrations of J_z – with the modulus of the current in this region being approximately equal to that of the current in the layer at the fan. We have repeated the simulation with different values of η, and found that for lower η the peak current density at the fan is stronger relative to the concentration near the boundary. Hence, if we were able to run the simulation with a realistic value of η for an astrophysical plasma the current layer at the fan would dominate.

As a result of the dominance of the J_z component near the driving boundaries (as discussed above), it is most clear to observe the formation and evolution of the fan current layer by plotting the evolution of the peak value of |J_xy|, as in Fig. 13 (the dominant current component in the current layer is parallel to the fan surface for all p). We can see that for p = 1 this quantity rises steadily as the pulse approaches the fan, after which there is a period of around t = 1.5 when the value remains steady, after which it decreases. Note that this steady period is consistent with the period during which the driving velocity remains steady – see the crosses in the figure. It is clear that the overall qualitative behaviour is similar between the simulations. However, the overall peak current attained is largest when the initial field is most asymmetric (p = 10). Note that for the simulations with the largest values of p the peak current does not remain steady for such a long period as for small p. This may be influenced by the fact that the disturbance interacts with the y-boundaries at later time for large p.

Fig. 13 Evolution of the maximum of |Jxy| for the different torsional fan simulation runs: p = 1 (solid line), p = 2 (dashed), p = 3 (dotted), p = 5 (dot-dashed) and p = 10 (triple-dot-dashed). For comparison, the time variation of the amplitude of the driving velocity is also represented, with the crosses.

Fig. 14 Contours of |J| for the torsional fan simulations, plotted in the z = 0 (fan) plane over  [x,y] ∈  [ ± 2.5, ± 2.5] for p = 1 (left), p = 2 (middle) and p = 5 (right). Taken in each case at the time when the value of |Jxy| reaches a maximum.

We now examine the spatial structure of the current distribution for simulation runs with different p. When the boundary driving is initiated the azimuthal symmetry of the wavefront that propagates into the domain is broken. The wavefront remains approximately planar (z constant), spreading into an elliptical shape with long axis along the y-direction (along which the Alfvén speed is greater) and short axis along x. However, as the wavefront gets closer to the null, it becomes steadily more inhomogeneous, with the current density focussing in the weak field regions near the x-axis. That is, although the current layer that forms eventually at the fan is more extended along y (as a simple consequence of mapping a circular driving region on the boundary along B), this current is most intense along the short axis of the ellipse, as shown in Fig. 14. As p is increased, the peak current increases as described above, and this increasing current maximum is localised in a gradually narrower “channel” around the x-axis.

We turn now to consider the rate of reconnection in the different simulations. As discussed above, at the limited magnetic Reynolds number we are able to use the fan current layer does not dominate over the distributed current near the driving boundaries – though indications are that it would for more realistic astrophysical parameters. Therefore, in order to examine only the effect of reconnection in the thin current layer, we calculate the reconnection rate by integrating E _∥  along the magnetic field line in the fan plane that passes through the peak of the current density – which in practice lies very close to the x-axis by symmetry. When the current layer has fully formed at the fan this genuinely measures the reconnection rate associated with the dynamically forming fan current layer. Clearly at earlier stages when the perturbation is yet to reach the fan there will still be a modest amount of reconnection, which is not measured by this approach. Therefore, when analysing the plots of reconnection rate versus time shown in Fig. 15, one should bear in mind that the curves do not accurately portray the reconnection rate prior to t ≈ 2.5. We see that the peak reconnection rate associated with the fan current layer occurs in each of the simulations at t ≈ 4, just before the peak current density starts its decline. A clear pattern emerges that the reconnection rate is greatest for the symmetric case with p = 1, and steadily decreases for simulations with larger p.

Fig. 15 Evolution of the reconnection rate for the different torsional fan simulation runs: p = 1 (solid line), p = 2 (dashed), p = 3 (dotted), p = 5 (dot-dashed) and p = 10 (triple-dot-dashed).

4.3. Torsional fan reconnection: discussion

The numerical simulations discussed above demonstrate the localisation of a rotational perturbation towards the fan of a non-symmetric linear 3D null, resulting in torsional fan magnetic reconnection. The planar current structure that forms around the fan qualitatively resembles that found in the kinematic steady-state model described in Sect. 3.1. In particular, the current is dominated by a component parallel to the fan where the current vector is directed (at its maximum intensity when p ≠ 1) radially towards the null. Current isosurfaces have an elliptical shape, with the short axis of the ellipse aligned with the weak field direction in the fan plane, along which the current is most intense. The eccentricity of the ellipse increases as the magnetic field asymmetry increases. The plasma flow also has a similar qualitative structure in the kinematic model and simulations. This structure is that of a rotation along elliptical paths around the spine, with these elliptical paths closely following the current density contours.

In the torsional fan reconnection simulations – in contrast to the torsional spine case – the peak reconnection rate depends strongly on the magnetic field asymmetry (p parameter). This is because in the torsional fan case, the contribution to the integrand in the reconnection rate definition (2) comes largely from field lines closely aligned to the xy-plane, which is the plane in which the magnetic field varies as p is varied. The interpretation of the reconnection rate when p ≠ 1 should be the same as that described above in Sect. 3.3. Perhaps counter-intuitively, the scaling of the reconnection rate with p is opposite to the dependence of the peak current on p (as p increases the peak current increases while the reconnection rate decreases – see Figs. 13, 15). This is because, unlike in 2D, the reconnection rate is not a local quantity defined at a point, but rather is defined as an integral along a field line. So if the peak value of the integrand in Eq. (2) is larger, but the integrand is non-zero over a shorter section of the field line, this may result in a lower net value for the integral, or vice versa. Note that this behaviour was also observed in the study of spine-fan reconnection carried out by Al-Hachami & Pontin (2010). Note also that the effect was not captured in the kinematic solution – however, it is worthwhile noting that the scaling of the reconnection rate with p is much weaker than that of the peak current. Specifically, while the current increases by a factor  ≈ 3 between p = 1 and p = 10 in the kinematic solution, the reconnection rate increases only by a factor  ≈ 1.2.

In the simulations after the current layer that formed dynamically at the fan has started to dissipate, the distributed current near the driving boundaries dominates to a steadily greater extent, since the large-scale current concentrations dissipate only very slowly. This is consistent with the results of Pontin & Craig (2005), who found that only when shear perturbations of the spine/fan are made does the lowest energy state of the system involve a current sheet at the null – when rotational perturbations are performed around the spine then the lowest energy state is achieved when the twist of the field is distributed along the field lines, in concentrations that extend outwards from the driving boundaries.

Finally, note that here we have chosen to model an instance of torsional fan reconnection in which the vorticity of the driving flow around the spine has opposite sign on the opposite boundaries. As a consequence, a strong current – directed predominantly in the radial direction – is present at the fan surface. However, in our kinematic model we could equally have chosen B_θ to be even in z to model the situation where the driving flows have the same sign of vorticity, which leads to cancellation at the fan plane of the currents generated (Galsgaard et al. 2003). It should also be noted that the return currents present close to z =  ± b are not present in the simulations as we drive only in one sense, but could be induced by reversing the sign of the driving velocity at some intermediate time in the simulations.

5. Conclusion

Here we have presented analytical and numerical models for torsional spine and torsional fan magnetic reconnection at 3D null points. The analytical models included for the first time fully localised current layers – focussed at the spine or fan – that determine the boundary of the non-ideal region, thus alleviating the requirement in previous models to have an artificially localised (“anomalous”) resistivity. We also for the first time investigated the generic case where the null point is not radially symmetric, i.e. where the fan eigenvalues are not equal.

3D null points have been demonstrated to be present in abundance in the solar corona, and the same is likely to be true in other astrophysical environments. Recent results from solar observations and from simulations suggest that reconnection at such 3D nulls may play an important role in the coronal dynamics. The fan separatrix surfaces of these 3D nulls divide the coronal magnetic field into distinct topological domains, i.e. distinct regions of magnetic connectivity between the photospheric flux concentrations. The torsional spine and torsional fan reconnection modes do not act to transfer magnetic flux between these topological domains, but rather permit a rotational slippage within the domains of the magnetic flux lying close to the nulls and therefore the domain boundaries. These null point reconnection modes are therefore unlikely to be involved in energetic events that involve a large-scale restructuring of the magnetic flux between topological domains as the coronal field seeks a lower energy state (e.g. during solar flares). Rather, they are a mechanism to dissipate energy and reduce stress associated with the dynamic perturbation of the coronal field by the turbulent boundary driving from the photosphere.

We have shown that the geometry of the current layers within which torsional spine and torsional fan reconnection occur is strongly dependent on the symmetry of the magnetic field defining the null point. Torsional spine reconnection still occurs in a narrow tube around the spine, but with elliptical cross-section when the fan eigenvalues are different. The eccentricity of the ellipse increases as the degree of asymmetry increases, with the short axis of the ellipse being along the strong field direction. Furthermore, the current profile is not azimuthally symmetric around the spine, but is peaked in these strong field regions. The numerical simulations suggest that the spatiotemporal peak current, and the peak reconnection rate attained, do not depend strongly on the degree of asymmetry. For torsional fan reconnection, the reconnection occurs in a planar disk in the fan surface, which is again elliptical when the symmetry of the magnetic field is broken. The short axis of the ellipse is along the weak field direction, with the current being peaked in these weak field regions. The peak current and peak reconnection rate in this case are clearly dependent on the asymmetry, with the peak current increasing but the reconnection rate decreasing as the degree of asymmetry is increased.

While we have relaxed the rotational symmetry of the magnetic field in these studies, the field structure of a linear null is still relatively simple. In the future it will be important to understand how these reconnection modes are modified – and how they release the energy associated with imposed stresses – when the null point is embedded in a more realistic coronal geometry. One recent study suggests that other features present in the magnetic field may attract the current preferentially over the nulls and therefore may in some cases inhibit the formation of the torsional spine and torsional fan current layers (Santos et al. 2011). Other studies have focussed on explaining explosive events and thus have been motivated by identifying reconnection via flux transfer between topological domains, facilitated by the spine-fan reconnection mode (e.g. Masson et al. 2009; Pariat et al. 2009; Török et al. 2009). In practice one would expect a combination of the pure reconnection modes to appear, and significant further study is required to determine the manifestations of, interactions between, and relative importance of the spine-fan mode and the two torsional modes discussed here. The importance of a number of other parameters such as the magnitude of the perturbation, the plasma-β and the resistivity η is also yet to be explored.

Acknowledgments

The authors would like to acknowledge fruitful discussions with G. Hornig, A. Wilmot-Smith and A. Yeates. D.I.P. gratefully acknowledges support from the Royal Society. A.K. Al-H. was supported by a grant from

the Iraqi Government. Support by the European Commission through the Solaire Network (MTRNCT-2006-035484) is gratefully acknowledged. Computational simulations were run on the UKMHD Computing Consortium’s Beowulf cluster.

References

All Tables

All Figures

	Fig. 1 Left: magnetic field lines for the torsional spine model defined by Eqs. (5), (8), (9) for a = 1, b = 5, j = 3, p = 1.5, q = 1. The shaded surface shows a current isosurface. Right: current vectors in the y = 0 plane for a = 1, b = 2, q = 1.
In the text

	Fig. 2 Current density (left) and plasma flow vectors (right) for the torsional spine kinematic solution. Plotted in the plane z = 2 for B0 = L0 = j = a = 1, b = 4 and p = q = 1 (above), p = q = 2 (middle) p = q = 5 (below). The grey curve (red online) in the right-hand plots marks the boundary of the non-ideal region.
In the text

	Fig. 3 Intensity map of vx in the x = 0 plane for the kinematic torsional spine solution with p = q = 1.
In the text

	Fig. 4 Dependence on the anisotropy parameter q = p of the maximum values of |J| and the reconnection rate Ψmax for the kinematic torsional spine model.
In the text

	Fig. 5 Time sequence showing contours of Jz in the z = 0 (left) and x = 0 (right) planes, for the torsional spine simulation with p = 2.
In the text

	Fig. 6 Above: current isosurface at 50% of maximum value, for the torsional spine simulations with p = 1 (left), p = 2 (middle) and p = 5 (right). Taken in each case at the time when the positive value of Jz reaches a maximum. Below: contours of Jz in the z = 0 (fan) plane over  [x,y] ∈  [ ± 0.2, ± 0.2] for the same runs at the same times.
In the text

	Fig. 7 Evolution of the peak value of |J| for the different torsional spine simulations: p = 1 (solid line), p = 2 (dotted), p = 3 (dashed), p = 5 (triple-dot-dashed) and p = 10 (dot-dashed).
In the text

	Fig. 8 Evolution of the reconnection rate for the different torsional spine simulations: p = 1 (solid line), p = 2 (dotted), p = 3 (dashed), p = 5 (triple-dot-dashed) and p = 10 (dot-dashed).
In the text

	Fig. 9 Above: magnetic field lines for the torsional fan model defined by Eqs. (5), (8), (12) for a = 5, b = 1, j = 50, p = 1.5, q = 1. The shaded surface shows a current isosurface. Below: current vectors in the y = 0 plane for a = 3, b = 1, q = 1.
In the text

	Fig. 10 Current density (left) and plasma flow vectors (right) in the kinematic torsional fan model. Plotted for B0 = L0 = j = b = 1, a = 4 and p = q = 1 (above), p = q = 2 (middle) p = q = 5 (below). The grey curve (red online) in the right-hand plots marks the boundary of the non-ideal region. The current density is plotted at z = 0, while the plasma flow is plotted in the plane z = 0.5.
In the text

	Fig. 11 Dependence on the anisotropy parameter q = p of the maximum values of |J| and the reconnection rate Ψmax for the kinematic torsional fan model.
In the text

	Fig. 12 Contours of Jy (left) and Jz (right) in the x = 0 plane over  [y,z] ∈  [±2.5, ± 0.5] for t = 0.8 (top), t = 1.8 (middle) and t = 4.2 (bottom), for the torsional fan simulation with p = 1.
In the text

	Fig. 13 Evolution of the maximum of |Jxy| for the different torsional fan simulation runs: p = 1 (solid line), p = 2 (dashed), p = 3 (dotted), p = 5 (dot-dashed) and p = 10 (triple-dot-dashed). For comparison, the time variation of the amplitude of the driving velocity is also represented, with the crosses.
In the text

	Fig. 14 Contours of |J| for the torsional fan simulations, plotted in the z = 0 (fan) plane over  [x,y] ∈  [ ± 2.5, ± 2.5] for p = 1 (left), p = 2 (middle) and p = 5 (right). Taken in each case at the time when the value of |Jxy| reaches a maximum.
In the text

	Fig. 15 Evolution of the reconnection rate for the different torsional fan simulation runs: p = 1 (solid line), p = 2 (dashed), p = 3 (dotted), p = 5 (dot-dashed) and p = 10 (triple-dot-dashed).
In the text