© ESO, 2011

1. Introduction

The solar photosphere is threaded by a complicated distribution of the magnetic field concentrations of opposite polarities. Above the photosphere in the equatorial region, the positive and negative flux concentrations typically connect to one another in a non-trivial way, creating complicated topological patterns. The simplest approach to investigating the topology of the magnetic field is to use magnetogram data as boundary conditions for potential or linear force-free magnetic field extrapolations. This has been done by a number of authors (Close et al. 2003; Schrijver & Title 2003; Longcope et al. 2003) and is found to lead to large numbers of separatrix surfaces dividing three-dimensional (3D) space into regions of different magnetic connectivity. These analyses show that 3D nulls are present in abundance at low altitudes, becoming increasingly rare high above the photosphere, and the number distribution depends on the distribution patterns of the photospheric magnetic field (Longcope & Parnell 2009). From observations there are also a few cases in which 3D nulls have been inferred from the structure of the emission in the corona (Filippov 1999), but using this indirect evidence for a 3D null is problematic without making attempts to derive the underlying magnetic field structure. This has recently been done by Masson et al. (2009), who verified their magnetic field extrapolations by comparing with the observed emission patterns. These independent approaches all show that 3D nulls are present in the solar atmosphere, and we must therefore investigate their role in the energy release process that takes place there.

Magnetic reconnection is thought to play a leading role in a large number of different types of magnetic energy release processes where fast conversion of magnetic energy into bulk motion, plasma heating and particle acceleration is observed. These processes have been observed in a wide variety of situations, for example in “explosive” events taking place in the outer solar atmosphere. To understand these phenomena in their various physical environments, a detailed understanding of the reconnection process is required. For many years two-dimensional (2D) reconnection models have formed the basis for this understanding. Analytical investigations have shown that various configurations with different characteristics are possible, and the driving profile controls these characteristic properties to a large extent. One key property that can be estimated is the scaling behaviour of the 2D reconnection, which informs us about how fast reconnection can take place in a given steady state situation.

One problem that has been encountered when starting to make numerical attempts to investigate the magnetic reconnection scenario is that the imposed magnetic resistivity (η) has significant implications for the development of the reconnection region. In the case where constant resistivity is imposed, the magnetic structure has a tendency to collapse into a current sheet, but it has a difficulty in reaching a steady state situation where there is a perfect balance between the inflow and the reconnection process with a constant size current sheet. Instead, the current sheet tends to continue growing towards a solution that resembles a Sweet-Parker current sheet, with close to “infinite” length and close to zero reconnection rate. The way to avoid the infinite current sheet is to change the η profile. As soon as η is made space, time or field dependent, the situation changes and the current sheet reaches a finite extent instead.

Over the last 10–15 years the focus of investigations has changed from 2D reconnection to the more challenging and realistic 3D reconnection processes. It is becoming clear that in 3D the possible reconnection scenarios are much more diverse than in 2D, and provide a series of new challenges to resolve. The process of 3D reconnection is not restricted to regions where the magnetic field vanishes as in 2D, but may also take place in regions with non-vanishing magnetic field. In this paper we will investigate in some detail the reconnection process that can take place in the vicinity of a single 3D null point. Null points constitute critical points of the magnetic field, with the only possible class being a hyperbolic critical point due to the solenoidal nature of B (Green 1989; Lau & Finn 1990; Parnell et al. 1996). A single pair of field lines approach (recede from) the null from opposite directions, and identify the so-called “spine” of the null. A set of field lines radiate away from (approach) the null in a surface known as the “fan” of the null (see Fig. 1).

Fig. 1 The structure of a single null. The e3 direction represents the spine axis, while the plane defined by e1 and e2 represents the fan plane.

Earlier investigations have shown that reconnection at 3D nulls can take place in a number of characteristic manners. First, the topological properties of the magnetic flux evolution have been shown to be crucially dependent on the orientation of the electric current at the null (Pontin et al. 2004, 2005). Second, depending on the motions that drive the formation of the current layer at the null, the reconnection may occur either in a tube aligned to the spine, a planar layer in the fan, or a layer focused at the null (Rickard & Titov 1996; Pontin & Galsgaard 2007). The various regimes have recently been categorized by Priest & Pontin (2009). Due to the increased complexity of the magnetic structure, analytical solutions to the problem of single null reconnection can only be achieved using various simplifying assumptions, and no-one has at present shown how steady-state reconnection takes place at a single null in a compressible plasma, or if indeed a steady-state solution is possible.

In the present investigation a 3D null is stressed by imposing a systematic driving velocity on the two boundaries through which the spine axis of the null penetrates. The driving is of a shear type, such that the angle between the spine and fan would be expected to change. Under the simplistic assumptions made by Priest & Titov (1996) this was expected to lead to current accumulation in the entire fan plane. Although it appears that such a situation would indeed occur in an incompressible plasma (Craig et al. 1995; Craig & Fabling 1998; Pontin et al. 2007b), when plasma compressibility is included, a local collapse of the null occurs, destroying the planar nature of the fan surface. A localised current sheet forms which is focused around the null itself, with the result that the geometry of the magnetic field around the null is significantly different from the initial linear profile. The resulting mode of reconnection has been dubbed “spine-fan reconnection” by Priest & Pontin (2009) (having the appearance of a hybrid of early spine and fan models), and involves magnetic field lines being transported through/around the spine and across the fan surface.

A series of numerical experiments have investigated the reconnection process at single null. These have used various forms of driving of the system, that generally have perturbed the null region for a finite period of time (Galsgaard et al. 2003; Pontin et al. 2007b,a; Pontin & Galsgaard 2007). This paper discusses the situation where a long time systematic driving is imposed on the system in an attempt to reach a steady state reconnection configuration. By continuously advecting the magnetic field in the vicinity of the spine axis, the magnetic field is compressed towards the null from two sides. Assuming that reconnection takes place, then due to the change in magnetic field strength away from the null, one would expect a constantly changing amount of magnetic flux to be advected towards the null per unit time. As a result, one might expect it to be difficult to obtain a true steady-state situation, while a continuously evolving quasi-steady-state situation may be more likely.

In this paper we will discuss, in Sect. 2, the generic set-up of the problem, the numerical approach and the imposed driving of the system. In Sect. 3 the experiments will be discussed and the results highlighted. This is followed, in Sect. 4, by a discussion of the consequences of the results for our general understanding of the null point reconnection process. Finally, in Sect. 5 we present our conclusions.

2. Numerical setup

The magnetic field configuration adopted here is the linear rotationally symmetric potential magnetic field defined by ${B}\mathrm{=}{{B}}_{\mathrm{0}}\mathrm{[}\mathrm{-}\mathrm{2}\mathit{x,y,z}\mathrm{]}\mathit{,}$(1)where B₀ is a scaling parameter, here of order unity. The null is located at (x,y,z) = (0,0,0) with the spine axis along the x-axis and the fan plane represented by the y − z plane, i.e. x = 0. This linear null point is simple to use in numerical experiments, but represents a non-physical situation where the magnetic field strength continues to increase with distance from the null point. In reality non-linear effects will become important at some distance and change the linear growth. To investigate the dynamical evolution of the null, the domain of interest has to be limited in size. For this investigation we have limited the domain to (x,y,z) =  ± (0.5,1.5,1.5), and have discretised it using a uniform grid with 200 × 300 × 300 points.

The plasma is initially assumed to be at rest and have a constant density and thermal energy. In the present experiment the MHD equations are non-dimensionalised and the density is set to unity, while the thermal energy is set to 0.025. This means that the surface on which the plasma beta equals unity is an ellipsoid with its minor axis along the x-axis of length 0.08 and major axis perpendicular to this with a length of 0.18. Inside this surface the plasma beta tends to infinity as the null point is approached.

The linear nature of the magnetic null point forces a limitation on the numerical domain size. As a Cartesian domain does not represent very well the initial rotational symmetry of the magnetic null structure, and because we are going to impose a symmetry breaking perturbation, it is a non-trivial problem to describe the boundary conditions of the imposed domain as the perturbation reaches the boundaries. In the present approach, we force the flow perpendicular to all of the boundaries to be zero, and restrict parallel velocities at the boundaries to be non-zero only in the two regions where dedicated stressing velocity flows are imposed.

Fig. 2 A surface plot of the imposed driving velocity profile.

In these experiments we limit the driving of the system to two areas on the x-boundaries, by advecting the magnetic field in the y-direction, with an additional spatial dependence on the z-direction. The imposed driver has the following expression, $\begin{array}{ccc}{\mathit{V}}_{\mathit{y}}\mathrm{\left(}\mathit{y,z}\mathrm{\right)}& \mathrm{=}& \mathrm{\pm }\mathrm{0.5}\hspace{0.17em}{\mathit{V}}_{\mathrm{d}}\hspace{0.17em}\mathrm{\left(}\tanh \mathrm{\right(}\mathrm{(}\mathit{y}\mathrm{-}{\mathit{y}}_{\mathrm{0}}\mathrm{)}\mathit{/}{\mathit{y}}_{\mathrm{h}}\mathrm{)}\mathrm{-}\tanh \mathrm{(}\mathrm{(}\mathit{y}\mathrm{+}{\mathit{y}}_{\mathrm{0}}\mathrm{)}\mathit{/}{\mathit{y}}_{\mathrm{h}}\mathrm{\left)}\mathrm{\right)}\\ & & \mathrm{\times } \mathrm{\left(}\tanh \mathrm{\right(}\mathrm{(}\mathit{z}\mathrm{-}{\mathit{z}}_{\mathrm{0}}\mathrm{)}\mathit{/}{\mathit{z}}_{\mathrm{h}}\mathrm{)}\mathrm{-}\tanh \mathrm{(}\mathrm{(}\mathit{z}\mathrm{+}{\mathit{z}}_{\mathrm{0}}\mathrm{)}\mathit{/}{\mathit{z}}_{\mathrm{h}}\mathrm{\left)}\mathrm{\right)}\mathit{,}\end{array}$(2)where V_d is the peak driving amplitude, and (y,z) = (± y₀, ± z₀) represent the locations (at the centre of the hyperbolic tangent function) where the driving velocity has half of its peak value. The variables y_h, z_h are the half widths of the hyperbolic tangent functions. This gives a near linear driving velocity between (y,z) = (± (y₀ − y_h), ± (z₀ − z_h)) while the velocity decreases outside this region, asymptotically approaching zero. For these experiments we use y₀ = 0.6, y_h = 0.2, z₀ = 0.3 and z_h = 0.2 (see Fig. 2). This profile is imposed on the two x-boundaries with opposite directions. The driving is switched on at the initiation of the experiment and reaches the peak velocity exponentially with a given time scale.

The experiments are conducted using the 3D non-ideal MHD code by Nordlund & Galsgaard (1997). This is a high order finite difference code using staggered grids to maintain conservation of physical quantities. The interpolation operators are fifth order in space while the derivative operators are sixth order. The solution is advanced in time using a third order explicit predictor-corrector method. Here we intend to derive scaling relations, so we adopt a constant η model.

3. Experiments

The aim of the experiments is to investigate the possibility of obtaining steady state reconnection at a 3D null point. We also seek to derive scaling relations to enable the extension of the results to higher Reynolds numbers. In what follows we investigate the effect of varying two different global parameters of the system, namely the magnetic resistivity (η) and the imposed driving velocity (V_d).

An overview of the various simulation runs is given in Table 1, together with three characteristic parameters for the experiments that will be discussed later. Adopting characteristic parameters for the experiment, we can determine the magnetic Reynolds numbers for these experiments Rm = V_aL/η, where V_a = 0.5 is the Alfvén velocity at the spine axis at the top/bottom boundary, L = 1 is the distance between the x-boundaries and η is given in Table 1. This leads to values for Rm between 50 and 5000.

The following subsections discuss different aspects of the experiments. We begin with a general description of the dynamical evolution, first discussing a 2D setup as a reference for the 3D experiment, since our knowledge of the behaviour in the 2D case is well established (at least within the MHD framework).

3.1. Evolution of a 2D system

Fig. 3 Shearing of a 2D X-point. Left: sheet length as a function of time for domain  [x,y]  ∈  [ ± 0.5, ± 1.5]  (dashed line) and  [x,y ]  ∈  [± 1, ± 3]  (solid line). Right: peak current as a function of time for the same two domain sizes.

We first briefly describe the evolution of a sheared 2D null point. The magnetic field is taken to be B =  [−x,y,0]  and the z-dependence of the driver is suppressed. As the driving is switched on two MHD waves propagate towards the null point from either of the driving boundaries. Close to the null they collide and a strong current layer forms in the region where the two separatrices collapse very close to one another. The length of the current sheet and the maximum current density (and therefore in 2D the reconnection rate) gradually build to peak values, before beginning a slow decay (see Fig. 3). During this period, the configuration can be described as being in a quasi-steady state.

There are a number of key points to note. First, the halt in the growth of the current layer (both length and modulus), is due to the finite size of our computational domain. If we repeat the simulations in a larger computational domain, with an extended driving region, the current layer continues to grow for substantially longer, before again being constrained by the boundaries (compare the full and dashed lines in Fig. 3). This is a well-known effect in 2D – that for uniform resistivity a “Sweet-Parker” type current layer develops, which continually grows as the forcing continues. The limitation on the current layer in these simulations is therefore equivalent to the limitation imposed when reconnection is forced in an initial Harris sheet with periodic boundaries (e.g. Birn et al. 2005). Another important point to note is that, were the driving really constant in our simulations, we would not expect such a noticeable decay of the current layer length and modulus after the peak is reached. The decay is due to a reduction in the quantity of flux being forced into the sheet per unit time as the simulations progress. This in turn is due to the evacuation of flux from the region upstream of the inflow, as flux from that region is constantly fed through the sheet. The above properties all have analogues in the 3D system, as discussed below.

3.2. Qualitative evolution of the 3D system

Fig. 4 Current amplitude at a number of characteristic instances in the dynamical evolution. The individual frames are scaled to their local dynamical range to enhance the visual appearance. The images represent the z = 0 plane for run “C3” at times t = (0.4,1.0,2.2,4.4,9.2,19.0).

Initiation of the driving on the two x-boundaries launches two wavefronts that propagate through the magnetised plasma. As the driving speed is slowly ramped up the wavefronts initially propagate towards the fan plane as simple plane Alfvén waves (with their x speed being independent of y and z). With the driving continuously stressing the magnetic field, the field lines behind the wavefront are systematically advected in the y-direction. This breaks the initial symmetry of the field, generating an asymmetry of the field on either side of the advected spine line. This causes the wavefront to tilt as its speed becomes dependent on y, due to an increase in wave speed in the region where the field lines are stretched, Fig. 4.

From the initial phase of the evolution, two points should be noted. First, the width of the current front decreases while its intensity increases on its way towards the fan plane. This is an effect of the decreasing Alfvén velocity in the x-direction, that for a pure Alfvén wave implies that it can never reach the fan plane of the null point in a strict mathematical sense. This results in a pileup producing a strong current front. As this occurs, the wavefront also spreads out increasingly fast along the y- and z-directions as it approaches the fan plane, due to the field geometry. One would expect this to act to reduce the current intensity in the front, were this effect the dominant one. However, due to the compression of the plasma, the wave is not a pure Alfvén wave, but couples to both slow and fast mode waves. The perturbation can therefore propagate across the magnetic field lines and advance towards the null point. This type of behaviour was seen in Galsgaard et al. (2003), where the null point was perturbed by two rotational twist waves centred on the spine axis. Here the fast mode wave wraps around the null and concentrates part of the energy of the perturbation there (see also McLaughlin et al. 2009, who studied the 2D case).

As the evolution proceeds, the two wavefronts initiated at opposite x-boundaries approach each other and eventually collide, with the two non-rotationally-symmetric current concentrations compressing into one central current distribution centred on the null. Within the resulting current layer, the spine and fan have collapsed towards one another, and the current spreads along the fan both along and across the driving direction. The tilt and magnitude of the arising current concentration depend on the imposed driving speed and imposed η as will be discussed later. Typically the tilt of the current sheet overshoots initially and the sheet experiences a few oscillations before settling to a stable angle. Due to the imposed constant η, reconnection takes place in the current sheet, and a near steady-state balance between the flux advected into the current sheet and the reconnected flux exiting the sheet through a reconnection jet is reached (see below). This picture is maintained with only minor changes until the termination of the experiments.

3.3. Steady state evolution

Fig. 5 The time evolution of the magnetic field in the inflow region (top left), the plasma inflow velocity (top right), the outflow velocity relative to the Alfvén velocity in the inflow region and the ratio between the inflow (lower left) and outflow velocities and sheet dimensions defined by L   vin/(l   vout) (lower right). For the run C3.

From the 2D experiment it is seen that the evolution reaches a quasi-steady-state. In order for it to be physically meaningful to search for scaling relations, we must first check if such a quasi-steady-state also exists in the 3D simulations. To determine this we examine some of the characteristic parameters addressed in the Sweet-Parker (Parker 1957; Sweet 1958) 2D reconnection model.

In order to define the characteristic size of the current layer we define the boundary of the layer to be the surface at which the current modulus  | J |  falls to 50% of its peak value (this peak occurring at the null). In 3D this provides three characteristic length scales that are here listed in descending order; the length (L) in the plane of the shear perturbation (here the xy-plane), the extent of the sheet in the direction perpendicular to this plane and parallel to the current at the null (here the z direction) which we denote here by the width (w), and finally the thickness, l, of the sheet (measured in the z = 0 plane).

First we estimate the amount of magnetic flux being transported into the sheet – for which one needs to know both the inflow velocity perpendicular to the current sheet and the magnetic vector aligned with the sheet. For the purposes of making such an estimate, we measure these values at the centre of the inflow region just outside the current sheet. We also consider the outflow jet velocity, measured toward the outflow edge of the current sheet. Figure 5 shows four combinations of these variables as functions of time from experiment “C3”. The top left frame represents the component of the magnetic field parallel to the current sheet. It shows an initial build-up followed by a near linear decline in magnitude. The latter is a consequence of the depletion of the magnetic field in the inflow region due to the limited length scale of the driving in the y-direction. The top right panel shows the inflow velocity, which is found to reach a stable level with only a slow decay with time. The lower left frame shows the ratio of the outflow velocity and the Alfvén speed in the inflow region (akin to the reconnection rate in the Sweet-Parker model), which reaches a nearly stable value at the same time as the magnetic field component reached its peak value. Finally the lower right frame shows the ratio of the length scale of the inflow region times the inflow velocity and the thickness of the outflow region times the outflow velocity. In the 2D Sweet Parker analysis this relation constitutes mass conservation, and here we can see that it reaches a stable value declining only slowly with time. It is interesting to note that this stable value is somewhat greater than 1 (the 2D value), being closer to a value around 1.5. This is likely due to the fact that the outflowing plasma may access the third dimension in 3D, rather than being constrained to flow out in the plane of the inflow.

In Fig. 6 the top left panel shows the time development of the maximum current density in the domain, which reaches a peak around t = 7 followed by a slow decay. The top right panel represents the time development of the jet velocity, which is slightly different from the current, taking longer to reach its maximum value. This is expected since the velocity is a collective effect of the Lorentz tension force and gas pressure in the newly reconnected field lines, acting on a “dense” plasma and therefore needing some additional time to accelerate the plasma to it cruising velocity. The lower left panel shows how the sheet length is very stable in time, while the lower right panel represents the tilt angle of the sheet that shows a slow decreases with time.

Fig. 6 The four panels show the time evolution of the peak value of the current in the sheet (top left), the peak velocity of the reconnection jet (top right), the length of the current sheet (lower left) and finally its angle relative to the y-axis (lower right), for the run C3.

From the above discussion it appears that the effect of the continual driving is to force the null structure into a near steady state, in which the reconnection region maintains its characteristics for a period of time much longer than the Alfvén travel time along the current sheet (Figs. 5 and 6). The qualitative behaviour discussed above is also observed in the other experiments listed in Table 1 above. We therefore go on in the following section to derive scaling relations between various parameters in the simulations.

3.4. Scaling relations

3.4.1. Peak current

It is well known that the peak current J_max reached in a numerical experiment depends on two key parameters, one being the numerical resolution. If the imposed value of η is too low then the numerical diffusion dominates and will limit J_max. The other key parameter is the chosen value of the resistivity η, which limits the current growth through the “correct” evolution of the induction equation. With sufficient numerical resolution one can resolve the structure of the current sheet and it becomes possible to make a meaningful analysis of the impact of changing η. In this investigation we have run simulations with a limited selection of η values.

One way to investigate whether the imposed resistivity dominates the numerical resistivity is to measure the width of the current sheet as a function of the imposed η value. In the incompressible case Heerikhuisen & Craig (2004) showed that this gives rise to a power law dependence with a power of 0.5 (although we note that the morphology of the current sheet is somewhat different in their incompressible models to what is observed here). Figure 7 shows the relation for the “C” experiments given in Table 1. The full line is a power law, but with a power of 0.35 for the higher η values. In the plot the diamonds represent the imposed values of η, with the point furthest to the left representing a case where η was set to zero. This indicates that the code introduces an efficient η when the width on the sheet becomes on the order of the numerical stencil. If one assumes that the power law dependence can be extended, then an asymptotic value of the numerical η can be estimated for the given resolution. Correcting the values of η correspondingly for the runs C5 and C6 (two left most points on the graph) brings those data points onto the line representing the scaling relation (diamonds in the plot). In the following we will investigate if this correction of the imposed resistivity will have any implication on the scaling relations.

Fig. 7 Scaling relation between the imposed η and the thickness of the current sheet (measured in units of the grid spacing in the x-direction), for the runs C1 − C6. The full line represents sheet thickness  ∝ η0.35.

Fig. 8 The scaling relation of the peak current with the imposed η value. The data points represent two series of experiments with different imposed driving velocities. The diamonds and the associated full line represent data with Vd = 0.05. The triangle represents the corrected η for the left most point as indicated in the previous graph, with the dot-dashed line being the correspondingly altered scaling. The crosses and the dashed line represent the data with Vd = 0.1 (with the star showing the corrected η value). The dashed line represents Eq. (3), the full line Eq. (4) and the dot-dashed line Eq. (5).

When data from the experiments listed in Table 1 are used, scaling relations of the peak current with η can be obtained. Figure 8 shows plots of peak current versus η for two different driving velocities. The triangle represents the “corrected” η value for the run “C4” and the star is the corrected value for “A3”. The scalings, corresponding to the lines plotted in the figure, are $\begin{array}{ccc}{\mathit{J}}_{\max }& \mathrm{=}& \mathrm{0.46}\hspace{0.17em}{\mathit{\eta }}^{-0.62} \mathrm{for}{\mathit{V}}_{\mathrm{d}}\mathrm{=}\mathrm{0.1}\mathit{,}\\ {\mathit{J}}_{\max }& \mathrm{=}& \mathrm{0.14}\hspace{0.17em}{\mathit{\eta }}^{-0.78} \mathrm{for}{\mathit{V}}_{\mathrm{d}}\mathrm{=}\mathrm{0.05}\mathrm{-}\mathrm{without}\mathit{\eta }\mathit{correction,}\\ {\mathit{J}}_{\max }& \mathrm{=}& \mathrm{0.1}\hspace{0.17em}{\mathit{\eta }}^{-0.85} \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for}{\mathit{V}}_{\mathrm{d}}\mathrm{=}\mathrm{0.05}\mathrm{-}\mathrm{with}\mathit{\eta }\mathit{correction}\mathit{.}\end{array}$Here one can see two different comparable fits for the slow (V_d = 0.05) driving experiments, with the scaling based on the corrected value leading to a slightly stronger dependence of J_max on η. However, the correction of the η value for the faster driver does not seem to lead to a linear scaling relation.

Fig. 9 The scaling relation between the imposed driving velocity and the peak current for the experiments with η = 10-3 listed in Table 1.

We now turn to the scaling of the peak current with the driving velocity, for a constant η value, which is investigated using the η = 10^-3 experiments listed in Table 1. The result is shown in Fig. 9, where it is seen that while the peak current clearly depends on the imposed driving velocity, there is no simple (e.g. linear or exponential) relation that can be obtained to fit the data.

3.4.2. Current layer dimensions and jet velocity

The dimensions and tilt angle of the current sheet have been measured for the same experiments. They are found to be independent of η for a constant driving velocity, to within our measurement accuracy. On the other hand, varying the driving velocity while maintaining a constant value for η, we find that there are significant changes of these dimensions, indicating that they are largely dependent on the amount of imposed stress. However, as discussed above, the present domain size and the extent of the driving region limit the current sheet from evolving freely, and thus with the present setup we are unable to provide scaling relations for these quantities.

Fig. 10 The top frame shows the time dependent evolution of the peak velocity, for the runs C1 (dot-dashed), C2 (dashed), C3 (triple-dot-dashed), C4 (long-dashed) and C5 (solid) (see Table 1). The middle frame shows the peak velocity as a function of the imposed driving velocity. The lower frame shows the peak velocity as a function of the imposed resistivity, η. The full line represents the scaling using the original values, Eq. (6), while the dot-dashed line included the corrected η values, Eq. (7).

The top frame of Fig. 10 shows the maximum outflow velocity as a function of time for different driving velocities V_d. This shows a significant variation both with time and imposed driving speed. The middle frame shows the peak outflow jet velocity (for all time) as a function of the imposed driving velocity. Here we can see that we have a non-linear dependence, which is also not fit well with an exponential scaling, but follows a similar pattern to the scaling of J_max with V_d. Finally, the lower frame of Fig. 10 shows the jet velocity as a function of the imposed resistivity, where we find approximately power law scalings, with $\begin{array}{ccc}{\mathit{V}}_{\mathrm{peak}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{8.1}}{\mathit{\eta }}^{-0.22}\mathrm{-}\mathrm{without}\mathit{\eta }\mathrm{correction,}& & \end{array}$(6)$\begin{array}{ccc}{\mathit{V}}_{\mathrm{peak}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{9.3}}{\mathit{\eta }}^{-0.25}\mathrm{-}\mathrm{with}\mathit{\eta }\mathrm{correction}\mathit{.}& & \end{array}$(7)The first scaling is for the case where we use the manually applied η values for the experiments, while in the second the corrected values are included.

3.5. Reconnection rate

From established theory (Schindler et al. 1988; Priest et al. 2003) we know that the correct measurement of the reconnection rate in 3D is given by the maximal value of Φ, the integrated electric field parallel to a magnetic field line (E_∥), where this maximum is taken over all field lines which thread the non-ideal region (assumed to be spatially localised). To determine the reconnection rate we therefore perform a systematic field line tracing for a collection of field lines that thread the current sheet, and seek the maximum in the integrated parallel electric field. In practice this occurs along a field line that is located very close to the z-axis (it would be expected to lie exactly along the z-axis based on the exactly symmetric analytical model of Pontin et al. 2005). In these experiments a constant η has been used, which implies that E_∥ = ηJ_∥, so that $\mathrm{\Phi }\mathrm{=}{\mathrm{\int }}_{\mathit{ϵ}}^{\mathit{L}}\mathrm{(}\mathit{\eta }{J}\mathrm{·}{B}\mathrm{)}\mathit{/}\mathrm{\left|}{B}\mathrm{\right|}\mathrm{d}\mathit{l,}$(8)where we integrate along the given field line starting close (ϵ) to the null out to a finite distance, L. To get the correct reconnection rate the integral has to be made in the z-direction in both senses, to  ± L. In previous simulations (e.g. Pontin et al. 2007a) this was a relatively simple process, as the diffusion region was limited to a finite size. However, the situation here is different. Due to the continual driving, the current not only has a strong component close to the null, but there also exists an extended (though weaker) concentration of current in the fan plane that reaches all the way to the z-boundaries.

Fig. 11 The top panel shows the parallel electric field (E∥ = E·B/|B | ) along the z-axis for the run “C3”. The bottom panel shows a 2D image of the parallel electric field in the x = 0 plane.

Fig. 12 The left panel shows the dependence of the reconnection rate on the imposed η with a constant driving velocity (Vd = 0.05), with the triangle showing the corrected value for the run “C5” as discussed above, Eq. (9). The right panel shows the dependence on the reconnection rate on the imposed driving velocity, for constant η (η = 10-3), Eq. (10).

This is demonstrated in the top image in Fig. 11, which shows a 1D plot of the parallel electric field along the z-axis (which is not a magnetic field line). In the lower frame a 2D (y,z) representation of the same quantity is presented, which again shows the extended parallel electric field. Plotting the same quantity in a plane perpendicular to this shows that the parallel electric field is by contrast limited to a very narrow layer in the x-direction (compare with Fig. 4). The actual reconnection process is therefore not limited to the local region around the strong current sheet, but extends all the way to the imposed boundary. This clearly implies that the value obtained using Eq. (8) depends on the choice of L. For the results below the value of L has been limited to L = 1.4, such that the domain boundary is not reached.

Performing the integration defined in Eq. (8) on the various data sets, we obtain the results shown in Fig. 12. The left frame shows the peak value of Φ as a function of η. If we correct the result obtained with the lowest η value (run C5, left most diamond in the plot) relative to the results obtained from Fig. 7, this point is shifted to the location marked by the triangle in the plot, that falls nicely on the straight line scaling relation that fits the other data points. This simple correction seems to work here, most likely because the current sheet in the diffusion region is all the way just being resolved with a minimum number of gridpoints.

In the right hand frame, an indication of a scaling relation between the driving velocity and reconnection rate can be seen, although with some fluctuations around the obtained line. The two power law relations are given by: $\begin{array}{ccc}\mathrm{\Phi }& \mathrm{=}& \int {\mathit{E}}_{\mathrm{\parallel }}\mathrm{d}\mathit{l}\mathrm{=}\mathrm{0.126} {\mathit{\eta }}^{\mathrm{0.25}} \mathrm{for} {\mathit{V}}_{\mathrm{d}}\mathrm{=}\mathrm{0.05}\\ \mathrm{\Phi }& \mathrm{=}& \mathrm{\int }{\mathit{E}}_{\mathrm{\parallel }}\mathrm{d}\mathit{l}\mathrm{=}\mathrm{0.15} {\mathit{V}}_{\mathrm{d}}^{\mathrm{0.65}} \mathrm{for} \mathit{\eta }\mathrm{=}{\mathrm{10}}^{-3}\mathit{.}\end{array}$

4. Discussion

The experiments described above were conducted to investigate whether a steady-state magnetic reconnection process could be set up in 3D by continually stressing a symmetric null point. The results presented in Sect. 3.3 show that a quasi-steady state exists where the inflow to the diffusion region is balanced by the ongoing magnetic reconnection (this is true for all values of the imposed driving velocity and η). Due to the magnetic geometry, the inflow magnetic field and the plasma inflow and outflow velocities (v_in and v_out) are not constant during this quasi-steady evolution. However, ratios between characteristic quantities familiar from 2D models do settle to approximately constant values. For example, the ratio of v_out to the inflow Alfvén speed is approximately fixed. So is the ratio Lv_in/lv_out, which characterises mass conservation in the sheet and settles to a constant value of around 1.5 (rather than 1 as in 2D, since the outflowing plasma may access the third dimension rather than being constrained to flow out in the plane of the inflow). It should be noted that in each of the simulations the quasi-steady configuration obtained is strongly influenced by the restricted dimensions of both the numerical domain and the region on which the boundary driving is imposed. If it was possible to extend the domain and driving region indefinitely, indications are that the current layer would extend indefinitely, until reaching a threshold for some instability such as the tearing instability.

We have investigated the dependence of different characteristic quantities in the quasi-steady state on both the driving velocity V_d and imposed resistivity η. Considering first the dependence on V_d, we found that the current maximum and outflow velocity increase as V_d increases, although not with a simple scaling relation. It is natural that these quantities increase when the rate at which flux is forced into the sheet is increased, and also that they behave in the same way since the outflow jet is accelerated by the J × B force.

For the simulation runs with different values of η we have been able to determine approximate scaling relationships. The peak current appears to scale like J_max ~ η^-0.6 − η^-0.8, which is consistent with the η ^− 3/4 scaling expected for planar incompressible fan current layers (Craig & Fabling 1998). The peak outflow velocity on the other hand follows an approximate scaling V_peak ~ η^-0.25, which mirrors the scaling of the reconnection rate, as we discuss below. It is also clear from Fig. 7 that the width of the current sheet depends on the imposed η value, and the graph indicates that as η approaches zero the sheet reaches a minimum thickness, controlled by numerical diffusion. Results for large η indicate a power law dependence for the thickness as a function of η, though with a smaller exponent than found in the incompressible models (Craig & Fabling 1998; Heerikhuisen & Craig 2004). One can use this to estimate the value of the numerical diffusivity, which takes effect only when structures collapse to the resolution limit, while its value decreases rapidly as the structures increase in size. Making use of this information to correct the scaling parameters in different plots is not simple, as indicated by the additional points included in a number of figures. In some cases it seems to provide a scaling relation, while in other cases, it makes no coherent correction to the result. The only way to provide better information for these values of η is through experiments with even higher numerical resolution. However, doubling the numerical resolution locally makes the experiments much more expensive. To make a serious difference in the scaling at least a factor of 10 variation in resolution is required – which gives a factor of 10⁴ in required computing time. It is therefore a significant computational challenge to improve these scaling relations significantly for lower η values.

From 2D reconnection theory the simple Sweet-Parker model gives a scaling of the reconnection rate that depends on η^0.5. With realistic parameters for astrophysical plasmas this reconnection rate is far too slow to make it an attractive mechanism for magnetic energy conversion, in for instance solar flares. It has therefore since then been a task to find reconnection setups that depend less strongly on η and therefore evolve much faster. Taking the scaling rate for the reconnection rate given in Eq. (9), the driven fan-spine reconnection is found to scale with a smaller power of η, namely 0.25. This scaling for the reconnection rate is significantly faster than the slow Sweet-Parker scaling, though slower than the fast Petcheck reconnection rate that scales with the logarithm of η. However, it is worth pointing out that over our range of η covering two orders of magnitude (η = 10^-2 − 10^-4) it is difficult to distinguish the power law and ln(η) scalings.

The obtained scaling relations should be compared with those obtained for the same setup, but with an impulsive (i.e. temporally localised) perturbation of the null point, described by Pontin et al. (2007a) and Priest & Pontin (2009). In these investigations, the boundary driving was applied for a fixed period of time and then reduced to zero. It is worth noting that this implies that the net displacement of the spine – under the assumption of an ideal evolution – increases as the driving speed increases. Their analysis showed that for fixed η, the peak current and reconnection rate increase linearly with the driving velocity, while the length and width of the sheet (L and w above) actually decrease linearly. This seems counter-intuitive, but is a result of the fact that the collapse of the field in the vicinity of the null is more efficient for stronger driving. Thus, the scaling of the current and reconnection rate follow similar patterns to those seen above in the continually driven case, however the sheet dimensions scale in the opposite way. Turning now to the scaling with resistivity η (for fixed driving speed) in the impulsively driven case, it was found that as η decreases, the peak current increases and the reconnection rate decreases, with a rate between power law and logarithmic, while here we may argue in favor for a power law dependence. Clearly there are some similarities, but also some differences, between the two different driving configurations.

5. Conclusions

The experiments have shown that it is possible for a 3D null point to reach a quasi-steady state, where the characteristic parameters of the current sheet are controlled by a combination of the imposed driving velocity and the effective constant resistivity. It is likely that other parameters, such as the dimensions of the plasma-β = 1 surface also play a role. It is found that scaling relations between characteristic parameters do not always follow simple exponential or power law expressions over the full range of η values. This is partly due to the limited numerical resolution in the simulations, which means that the code imposes its own diffusion when structures become too small. Only with substantially higher numerical resolution will it be possible to push the reliable range of η to significantly lower values. However, this preliminary study indicates that the reconnection rate scales relatively weakly with the plasma resistivity, as η^1/4.

The scaling relations found here are different from the ones derived for the impulsive perturbation of the single null. However, we find that even changing the amplitude of the fixed driving velocity changes the scaling with resistivity. It is therefore not easy to make any predictions regarding the global behaviour of 3D nulls in realistic situations, because too many parameters seem to influence the dynamical evolution and the scalings are linked in ways that appear to be non-linear, and that we therefore still need to unravel. Further to this, test experiments indicate that some of the obtained parameters are likely to change when the size of the imposed driver is changed.

It is clear from the results that more effort has to be invested in understanding null reconnection and its scaling behaviour before we are able to make any quantitative predictions regarding processes occurring at 3D nulls in the solar atmosphere or the Earth’s magnetosphere.

Acknowledgments

Thanks to the anonymous referee for making suggestions to view some of the results in a different light. Computing time was given by the DCSC at university of Copenhagen. Support by the European Commission through the Solaire Network (MTRNCT-2006-035484) is gratefully acknowledged. D.P. acknowledges support from the Royal Society.

References

All Tables

All Figures

	Fig. 1 The structure of a single null. The e3 direction represents the spine axis, while the plane defined by e1 and e2 represents the fan plane.
In the text

	Fig. 2 A surface plot of the imposed driving velocity profile.
In the text

	Fig. 3 Shearing of a 2D X-point. Left: sheet length as a function of time for domain  [x,y]  ∈  [ ± 0.5, ± 1.5]  (dashed line) and  [x,y ]  ∈  [± 1, ± 3]  (solid line). Right: peak current as a function of time for the same two domain sizes.
In the text

	Fig. 4 Current amplitude at a number of characteristic instances in the dynamical evolution. The individual frames are scaled to their local dynamical range to enhance the visual appearance. The images represent the z = 0 plane for run “C3” at times t = (0.4,1.0,2.2,4.4,9.2,19.0).
In the text

	Fig. 5 The time evolution of the magnetic field in the inflow region (top left), the plasma inflow velocity (top right), the outflow velocity relative to the Alfvén velocity in the inflow region and the ratio between the inflow (lower left) and outflow velocities and sheet dimensions defined by L   vin/(l   vout) (lower right). For the run C3.
In the text

	Fig. 6 The four panels show the time evolution of the peak value of the current in the sheet (top left), the peak velocity of the reconnection jet (top right), the length of the current sheet (lower left) and finally its angle relative to the y-axis (lower right), for the run C3.
In the text

	Fig. 7 Scaling relation between the imposed η and the thickness of the current sheet (measured in units of the grid spacing in the x-direction), for the runs C1 − C6. The full line represents sheet thickness  ∝ η0.35.
In the text

Fig. 8 The scaling relation of the peak current with the imposed η value. The data points represent two series of experiments with different imposed driving velocities. The diamonds and the associated full line represent data with Vd = 0.05. The triangle represents the corrected η for the left most point as indicated in the previous graph, with the dot-dashed line being the correspondingly altered scaling. The crosses and the dashed line represent the data with Vd = 0.1 (with the star showing the corrected η value). The dashed line represents Eq. (3), the full line Eq. (4) and the dot-dashed line Eq. (5).

In the text

	Fig. 9 The scaling relation between the imposed driving velocity and the peak current for the experiments with η = 10-3 listed in Table 1.
In the text

Fig. 10 The top frame shows the time dependent evolution of the peak velocity, for the runs C1 (dot-dashed), C2 (dashed), C3 (triple-dot-dashed), C4 (long-dashed) and C5 (solid) (see Table 1). The middle frame shows the peak velocity as a function of the imposed driving velocity. The lower frame shows the peak velocity as a function of the imposed resistivity, η. The full line represents the scaling using the original values, Eq. (6), while the dot-dashed line included the corrected η values, Eq. (7).

In the text

	Fig. 11 The top panel shows the parallel electric field (E∥ = E·B/|B | ) along the z-axis for the run “C3”. The bottom panel shows a 2D image of the parallel electric field in the x = 0 plane.
In the text

	Fig. 12 The left panel shows the dependence of the reconnection rate on the imposed η with a constant driving velocity (Vd = 0.05), with the triangle showing the corrected value for the run “C5” as discussed above, Eq. (9). The right panel shows the dependence on the reconnection rate on the imposed driving velocity, for constant η (η = 10-3), Eq. (10).
In the text