© ESO, 2012

1. Introduction

Traditionally, magnetic reconnection and magnetohydrodynamics (MHD) wave theory have been viewed as separate areas of solar physics (see, e.g., Priest & Forbes 2000; Roberts 2004; De Moortel 2005; Nakariakov & Verwichte 2005; De Moortel & Nakariakov 2012). However, this is a misconception: we know that (steady-state) reconnection models not only generate outflows/waves, but also require inflows/waves (e.g. Parker 1957; Sweet 1958; Petschek 1964). Several authors have already challenged this point-of-view (e.g. Craig & McClymont 1991; Longcope & Priest 2007; Murray et al. 2009; McLaughlin et al. 2009, 2012) and their investigations contribute to our understanding of dynamic or time-dependent models of magnetic reconnection. Of particular importance to this paper is the work of McLaughlin et al. (2009) which is the first demonstration of reconnection naturally driven by MHD wave propagation, via a process entitled oscillatory reconnection.

MHD wave propagation in inhomogeneous media is a fundamental plasma process and the study of MHD waves in the neighbourhood of magnetic null points directly contributes to this area (see review by McLaughlin et al. 2011). It is known that null points – weaknesses in the magnetic field where the field strength, and hence the Alfvén speed, is zero – and separatrices – topological features that separate regions of different magnetic flux connectivity – are an inevitable consequence of the distributed isolated magnetic flux sources at the photospheric surface, where the number of such null points will depend upon the magnetic complexity of the photospheric flux distribution (see, e.g., review by Longcope 2005; Close et al. 2004; Régnier et al. 2008; Longcope & Parnell 2009, for the statistics of coronal null points). It is also now known that MHD wave perturbations are omnipresent in the corona (e.g. Tomczyk et al. 2007). Thus, these two areas of scientific study; MHD waves and magnetic topology, will encounter each other in the corona, i.e. MHD waves will propagate into the neighbourhood of coronal null points (e.g. blast waves from a flare will at some point encounter a null point).

McLaughlin & Hood (2004, 2005, 2006a,b) investigated the behaviour of linear MHD waves (fast and slow magnetoacoustic waves and Alfvén waves) in the neighbourhood of a variety of 2D null points. It was found that the (linear) fast wave is focused towards the null point by a refraction effect and all the wave energy, and thus current density, accumulates close to the null, i.e. null points will be locations for preferential heating by (linear) fast waves. The Alfvén wave propagates along magnetic fieldlines and so accumulates along the separatrices (in 2D) or along the spine or fan-plane (in 3D). Waves in the neighbourhood of a single 2D null point have also been investigated using cylindrical models, in which the generated waves encircled the null point (e.g. Bulanov & Syrovatskii 1980; Craig & McClymont 1991, 1993; Craig & Watson 1992; Hassam 1992) and it was found that the wave propagation leads to an exponentially-large increase in the current density (see also Ofman 1992; Ofman et al. 1993; Steinolfson et al. 1995, and a comprehensive review by McLaughlin et al. 2011, for further details). 3D MHD wave activity about coronal null points has been investigated by various authors (e.g. Galsgaard et al. 2003; Pontin & Galsgaard 2007; Pontin et al. 2007; McLaughlin et al. 2008; Galsgaard & Pontin 2011a,b; Thurgood & McLaughlin 2012).

Reconnection can occur when strong currents cause the magnetic fieldlines to diffuse through the plasma and change their connectivity (Parker 1957; Sweet 1958; Petschek 1964). However, these papers did not include the effect of gas pressure, which would act to limit the growth of the current density. In considering the relaxation of a 2D X-type neutral point disturbed from equilibrium, Craig & McClymont (1991) found that free magnetic energy is dissipated by a physical mechanism which couples resistive diffusion at the null to global advection of the outer field, which they called oscillatory reconnection. An example of oscillatory reconnection generated by magnetic flux emerging into a coronal hole was reported by Murray et al. (2009) who found a series of “reconnection reversals” take places as the system searches for equilibrium, i.e. the system demonstrates oscillatory reconnection in a self-consistent manner. The physics behind oscillatory reconnection has been investigated by McLaughlin et al. (2009), Murray et al. (2009) and Threlfall et al. (2012).

McLaughlin et al. (2012) investigated the long-term evolution of an initially-buoyant magnetic flux tube emerging into a gravitationally-stratified coronal hole environment and reported on the resulting oscillations and outflows. They found that the physical mechanism of oscillatory reconnection naturally generates quasi-periodic vertical outflows with a transverse/swaying aspect. There is currently a great deal of interest in observations of transverse motions in the solar atmosphere (e.g. Tomczyk et al. 2007; De Pontieu et al. 2007, 2011; Cirtain et al. 2007; Erdélyi & Taroyan 2008; Nishizuka et al. 2008, 2011; He et al. 2009a,b; Liu et al. 2009, 2011; McIntosh et al. 2011; Okamoto & De Pontieu 2011; Morton et al. 2012; Yurchyshyn et al. 2012) and transverse/swaying motions have been observed over a range of wavelengths, speeds, temperatures and scales. However, the origin of these propagating, transverse oscillations remains a mystery, and these authors often note that the challenge remains to understand how and where these waves are generated in the solar atmosphere. Liu et al. (2011) summarises possible generation mechanisms for these transverse motions, including an oscillating wake from a coronal mass ejection or periodic reconnection (see, e.g., Chen & Priest 2006; Sych et al. 2009). The physical mechanism of oscillatory reconnection is another possible source of these transverse motions. As reported by McLaughlin et al. (2012), the transverse behaviour seen in the periodic jets originating from the reconnection region of the inverted Y-shaped structure is specifically due to the oscillatory reconnection mechanism, and would be absent for a single, steady-state reconnection jet. The physical mechanism also naturally generates periodic outputs even though no periodic driver is imposed on the system.

Thus, there is a clear interest in furthering our understanding of the periodic nature of oscillatory reconnection. In this paper, we investigate the periodic signal generated by the mechanism, with a specific interest in measuring periods and decay rates as well as the robustness of our results, i.e. how do the results vary with the strength of the driver.

1.1. Overview of McLaughlin et al. (2009)

This paper will closely follow the work of McLaughlin et al. (2009) as we investigate the periodic nature of oscillatory reconnection. These authors investigated the behaviour of nonlinear fast magnetoacoustic waves near a 2D X-type neutral point and found that the incoming wave deforms the null point into a cusp-like point which in turn collapses to a current sheet. The system then evolves periodically through a series of horizontal/vertical current sheets with associated changes in connectivity, i.e. the system demonstrates the mechanism of oscillatory reconnection.

More specifically, McLaughlin et al. (2009) found that the incoming (fast) wave propagates across the magnetic fieldlines and the initial profile, an annulus, contracts as the wave approaches the null point. This is the refraction behaviour that is typical of fast wave behaviour around magnetic null points (see, e.g., McLaughlin et al. 2011) and results from the spatially-varying (equilibrium) Alfvén-speed profile.

The incoming wave was observed to develop discontinuities (for a physical explanation, see Appendix B of McLaughlin et al. 2009 or, alternatively, Gruszecki et al. 2011) and these discontinuities form fast oblique magnetic shock waves, where the shock makes B refract away from the normal. Interestingly, the shock locally heats the initially β = 0 plasma, creating β ≠ 0 at these locations.

At a later time, the shocks overlap, forming a shock-cusp, which leads to the development of hot jets and in turn these jets substantially heat the local plasma and significantly deform the local magnetic field. By the time the shock waves reach the null, the (initially X-point) magnetic field has been deformed such that the separatrices now touch one another rather than intersecting at a non-zero angle (Priest & Cowley 1975, call this “cusp-like”). The osculating field structure continues to collapse and forms a horizontal current sheet. However, the separatrices continue to evolve: the jets at the ends of the (horizontal) current sheet continue to heat the local plasma, which in turn expands. This expansion squashes and shortens the current sheet, forcing the separatrices apart. The (squashed) current sheet thus returns to a “cusp-like” null point that, due to the continuing expansion of the heated plasma, in turn forms a vertical current sheet. In effect, the (net) restoring force acts to return the (deformed) null point to its equilibrium state, but overshoots the equilibrium. The phenomenon then repeats itself: jets heat the plasma at the ends of this newly-formed (vertical) current sheet, the local plasma expands, the (vertical) current sheet is shortened, the system attempts to return itself to equilibrium, overshoots and forms a (second) horizontal current sheet. The evolution proceeds through a series of horizontal and vertical current sheets, and the system clearly displays oscillatory behaviour. It is also interesting to note that the final state is non-potential, where this is because the plasma to the left and right of the null is (locally) hotter than that above and below. Consequently, a thermal-pressure gradient exists and causes the X-point to be slightly closed up in the vertical direction (i.e. generating a small, positive current). It is important to note that the non-potential final state is still in force balance and will eventually return to a potential state, but on a far greater timescale than our simulations (t_diffusion ~ R_m = 10⁴ Alfvén times, where R_m is the magnetic Reynolds number). We also note that there is nothing unique about the orientation of the first current sheet being horizontal followed by a vertical, this simply results from the particular choice of initial condition, and McLaughlin et al. (2012) use the more general terminology orientation 1 and orientation 2.

Fig. 1 Time evolution of jz(0,0,t) (measured in milliAmps) for 0 ≤ t ≤ 480   s. Red line indicates impulsive phase and black line indicates stationary phase. Insert shows the same time evolution over 158.9 ≤ t ≤ 480   s (same horizontal axis, different vertical axis). The black dotted lines indicate the formation times of all the horizontal current sheets and the green line indicates jfinal = 0.8165   j0 = 6.5 × 10-5   A. The blue dashed lines indicate an exponentially-damped envelope max(jz)|stationary × e−λt + jfinal and min(jz)|stationary × e−λt + jfinal where λ = −0.015   s-1.

McLaughlin et al. (2009) also present evidence of reconnection occurring in the system; reporting both a change in fieldline connectivity (qualitative evidence) and changes in the vector potential which directly showed a cyclic increase and decrease in magnetic flux on either side of the separatrices (see their Figs. 12 and 13). Hence, since the system displayed both oscillatory behaviour and reconnection, it was concluded that the system demonstrated the phenomenon of oscillatory reconnection.

Our paper has the following outline: the basic setup, equations and assumptions are described in Sect. 2, the periodic nature of oscillatory reconnection is detailed in Sect. 3 and the conclusions are given in Sect. 4.

2. Basic equations

We consider the nonlinear, compressible, resistive MHD equations: $\begin{array}{ccc}\mathit{\rho }\left[\frac{\mathit{\partial }{v}}{\mathit{\partial t}}\mathrm{+}\left({v}\mathrm{·}\mathrm{\nabla }\right){v}\right]& \mathrm{=}& \mathrm{-}\mathrm{\nabla }\mathit{p}\mathrm{+}\left(\frac{\mathrm{1}}{\mathit{\mu }}\mathrm{\nabla }\mathrm{\times }{B}\right)\mathrm{\times }{B}\mathit{,}\\ \frac{\mathit{\partial }{B}}{\mathit{\partial t}}& \mathrm{=}& \mathrm{\nabla }\mathrm{\times }\mathrm{(}{v}\mathrm{\times }{{B}}^{\mathrm{)}}\mathrm{+}\mathit{\eta }{\mathrm{\nabla }}^{\mathrm{2}}{B}\mathit{,}\\ \mathit{\rho }\left[\frac{\mathit{\partial ϵ}}{\mathit{\partial t}}\mathrm{+}\left({v}\mathrm{·}\mathrm{\nabla }\right)\mathit{ϵ}\right]& \mathrm{=}& \mathrm{-}\mathit{p}\mathrm{\nabla }\mathrm{·}{v}\mathrm{+}\frac{\mathrm{1}}{\mathit{\sigma }}{\left|{j}\right|}^{\mathrm{2}}\mathrm{+}{\mathit{Q}}_{\mathrm{shock}}\mathit{,}\\ \frac{\mathit{\partial \rho }}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\left(\mathit{\rho }{v}\right)& \mathrm{=}& \mathrm{0}\mathit{,}\end{array}$(1)where ρ is the mass density, v is the plasma velocity, B the magnetic induction (usually called the magnetic field), p is the plasma pressure, μ = 4π × 10^-7   Hm^-1 is the magnetic permeability, σ is the electrical conductivity, η = 1/μσ is the magnetic diffusivity, ϵ = p/ρ(γ − 1) is the specific internal energy density, γ = 5/3 is the ratio of specific heats and j = ∇ × B/μ is the electric current density.

We solve these governing equations numerically using a Lagrangian remap, shock-capturing code called LARE2D (Arber et al. 2001), which utilizes artificial shock viscosity to introduce dissipation at steep gradients. The details of this technique, often called Wilkins viscosity, can be found in Wilkins (1980). Thus, Q_shock represents the viscous heating at shocks.

We now introduce a change of scale to non-dimensionalise all variables. Letting v = v₀v^∗, B = BB^∗, x = Lx^∗, y = Ly^∗, z = Lz^∗, ρ = ρ₀ρ^∗, p = p₀p^∗, j = j₀   j^∗, ∇ = ∇^∗/L, t = t₀t^∗ and η = η₀, where * denotes a dimensionless quantity and v₀, B, L, ρ₀, p₀, j₀, t₀ and η₀ are constants with the dimensions of the variable they are scaling. We then set $\mathit{B}\mathit{/}\sqrt{\mathit{\mu }{\mathit{\rho }}_{\mathrm{0}}}\mathrm{=}{\mathit{v}}_{\mathrm{0}}$ and v₀ = L/t₀ (this sets v₀ as a constant background Alfvén speed). We also set j₀ = B/μL and ${\mathit{\eta }}_{\mathrm{0}}{\mathit{t}}_{\mathrm{0}}\mathit{/}{\mathit{L}}^{\mathrm{2}}\mathrm{=}{\mathit{R}}_{\mathrm{m}}^{-1}$, where R_m is the magnetic Reynolds number, and choose R_m = 10⁴. This process non-dimensionalises Eq. (1) and under these scalings, t^∗ = 1 (for example) refers to t = t₀ = L/v₀; i.e. the time taken to travel a distance L at the background Alfvén speed.

There is no fixed dimensional length scale to our X-point system (X-points are scale-free), and thus we have a great deal of freedom in choosing our dimensional constants. We choose L = 1   Mm and B = 1   G (for simplicity, and where these choices allow an intuitive understanding of Eq. (2) below) and we chose a coronal density of ρ₀ = 5 × 10^-13   kg/m³ and coronal temperature of T₀ = 10⁶   K. This sets ${\mathit{v}}_{\mathrm{0}}\mathrm{=}\mathit{B}\mathit{/}\sqrt{\mathit{\mu }{\mathit{\rho }}_{\mathrm{0}}}\mathrm{=}\mathrm{126.2}$km s^-1, t₀ = L/v₀ = 7.93   s and j₀ = B/μL = 8 × 10^-5   A. The values returned from Eq. (1) are made dimensional using these solar constants.

2.1. Basic equilibrium and numerical set-up

To set-up our system, we follow the numerical framework of McLaughlin et al. (2009). Thus, we consider a simple 2D X-type neutral point as our equilibrium magnetic field, where the initial field is taken as: $\begin{array}{ccc}{{B}}_{\mathrm{0}}\mathrm{=}\frac{\mathit{B}}{\mathit{L}}\left(\mathit{y,x,}\mathrm{0}\right)\mathit{,}& & \end{array}$(2)where B = 1   G is a characteristic field strength and L = 1   Mm is the length scale for magnetic field variations. This magnetic field can be seen in McLaughlin et al. (2009, their Fig. 1).

Fig. 2 a) Formation time (measured in seconds) and b) length (measured in km and evaluated at the formation time) of first horizontal current sheet versus initial velocity amplitude (measured in km s-1), where v0 = 25.2   km s-1 corresponds to C = 1 simulation.

Initially, we consider the equilibrium plasma to be cold: T = 0   K (i.e. β(t = 0) = 0) and, hence, ignore plasma pressure effects. However, McLaughlin et al. (2009) showed that magnetic shocks will heat the plasma and so the plasma will not remain cold (see e.g. Sect. 1.5 in Priest & Forbes 2000).

In order to excite a pure fast magnetoacoustic wave, we consider an initial condition that perturbs velocity purely across the equilibrium magnetic field. Using the terminology of McLaughlin et al. (2009), there are three distinguishing velocity components that can be considered in our system:

• ${\mathit{v}}_{\mathrm{\perp }}\mathrm{=}{v}\mathrm{\times }{{B}}_{\mathrm{0}}\mathrm{·}\hat{{z}}\mathrm{=}{\mathit{v}}_{\mathit{x}}{\mathit{B}}_{\mathit{y}}\mathrm{-}{\mathit{v}}_{\mathit{y}}{\mathit{B}}_{\mathit{x}}$, which corresponds to the velocity across the equilibrium field,and hence corresponds to the fast magnetoacoustic wave (theonly MHD wave that can cross fieldlines),

• v _∥  = v·B₀ = v_xB_x + v_yB_y, which corresponds to the velocity along the equilibrium field and corresponds to the propagation of the slow magnetoacoustic wave,

• ${\mathit{v}}_{\mathit{z}}\mathrm{=}{v}\mathrm{·}\hat{{z}}$, which corresponds to the velocity in the invariant direction and hence corresponds to the Alfvén wave.

Interestingly, these three velocity components – each isolating an individual MHD mode – are in good agreement with those reported by Thurgood & McLaughlin (2012) who used the equilibrium magnetic field and the flux function (which is parallel to the invariant direction) to define an orthogonal coordinate system to isolate and identify the propagation of each of the MHD modes. Using their convention, ${\mathit{v}}_{\mathrm{\perp }}\mathrm{=}{v}\mathrm{·}{{B}}_{\mathrm{0}}\mathrm{\times }\hat{{A}}$, where ${A}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{}^{\mathrm{(}}\mathit{y}^{\mathrm{2}}\mathrm{-}{\mathit{x}}^{\mathrm{2}}{}^{\mathrm{)}}\hat{{z}}$ is the flux function, and where perturbations in the ${{B}}_{\mathrm{0}}\mathrm{\times }\hat{{A}}$-direction were shown to correspond to those of the fast wave (for helicity-free systems, such as our 2D X-point).

To perturb our system, we consider an initial condition in velocity such that: $\begin{array}{ccc}{\mathit{v}}_{\mathrm{\perp }}\left(\mathit{x,y,t}\mathrm{=}\mathrm{0}\right)& \mathrm{=}& \mathrm{2}\mathit{C}\sin \left[\mathit{\pi }\left(\mathit{r}\mathrm{-}\mathrm{4.5}\right)\right]\mathrm{for}\mathrm{4.5}\mathrm{\le }\mathit{r}\mathrm{\le }\mathrm{5.5}\mathit{,}\\ {\mathit{v}}_{\mathrm{\parallel }}\left(\mathit{x,y,t}\mathrm{=}\mathrm{0}\right)& \mathrm{=}& \end{array}$(3)where r² = x² + y² and 2C is our initial amplitude and initial condition (3) describes a circular, sinusoidal pulse in v _⊥ . Thus, as argued above, we initially generate (only) a pure fast wave in our system. Note that the velocity profile prescribed by Eq. (3) appears as the symmetric m = 0 mode in v _⊥ , but corresponds to the asymmetric m = 2 mode in Cartesian components. This is why the first current sheet has horizontal orientation, as per Sect. 1.1.

This initial pulse will naturally split into two waves – an outgoing wave and an incoming wave – each of amplitude C. In this paper, we will focus on the incoming wave, i.e. the wave propagating towards the null point. In this paper, we will conduct a parameter study of initial wave amplitude C. Note that setting C = 1 recovers the results of McLaughlin et al. (2009) and choosing a small value for C, say C = 0.001, recovers the the linear results from McLaughlin & Hood (2004; i.e. see Appendix A of McLaughlin et al. 2009). Under our dimensionalisation, a choice of C = 1 corresponds to an incoming wave with maximum initial Cartesian velocity v _⊥ /r = v₀/5 = 25.2   km s^-1 at a distance of r = 5L = 5   Mm and an equilibrium magnetic field strength of 5   G.

The governing Eq. (1) with initial conditions (3) are solved computationally in a square domain x,y ∈  [− 20,   20]  with a numerical resolution of 6144 × 6144. Zero gradient boundary conditions are applied to the variables B, ρ, ϵ at the four boundaries, and v is set to zero on all boundaries, i.e. reflective boundaries. A numerical damping region exists for x² + y² ≥ 6 which gradually removes kinetic energy from the outgoing waves and so all oscillations that enter this region are slowly damped away, and hence they do not influence the behaviour about the null. The (equilibrium) Alfvén speed increases with distance from the null point and, hence, waves accelerate as they propagate outwards. Since we do not want reflected waves to influence our null point, implementation of such a damping region is essential.

3. Aperiodic driver leading to periodic behaviour

We set C = 1 in Eq. (3) and, as expected, we recover the results of McLaughlin et al. (2009) and readers are directed to that paper for full details (primarily their Figs. 2, 5 and 6) and also this paper’s Sect. 1.1. One of the key results from McLaughlin et al. (2009) was the production of a periodic response resulting from an aperiodic input, i.e. the physical mechanism of oscillatory reconnection naturally gave rise to periodic behaviour. This periodic response can be quantitatively measured by analysing the time evolution of the (electric) current density, specifically j = (0,0,j_z), at the null point itself. This can be seen in McLaughlin et al. (2009, their Fig. 10) and the analysis of this time signal is the primary focus of this paper. In recreating the simulations of McLaughlin et al. (2009), we recover this periodic time series, but at a higher numerical resolution, 6144 × 6144, and at a higher cadence, dt = 0.05   t₀. The time evolution of j_z(0,0,t) for 0 ≤ t ≤ 480   s can be seen in Fig. 1. Note that due to the symmetry of our system, the null is always located at the origin.

We identify two distinct regimes in Fig. 1: 0 ≤ t < 90   s (red line in Fig. 1) which we refer to as the impulsive phase and t ≥ 90   s (black line in Fig. 1) which we refer to as the stationary phase. The impulsive (or transient) phase is a spiky, irregular signal, with multiple local extrema indicating: the arrival time of the shock-cusp at the null point (at t = 20.6   s), the formation of the first horizontal current sheet (at t = 35.7   s, j_z < 0), the formation of the first vertical current sheet (at t = 59.1   s, j_z > 0) and the formation of the second horizontal current sheet (at t = 90.0   s). The (current-sheet) signal is further contaminated due to addition of small currents related to the propagation of shock waves across the null point. We define the end of the impulsive phase as the formation time of the second horizontal current sheet, which in our C = 1 simulation occurs at t = 90.0   s.

After t = 90.0   s, the evolution of j_z(0,0,t) is much cleaner and closer to (damped) sinusoidal. For t ≥ 90.0   s, the extrema exactly match the formation of the cyclic current sheets, with j_z < 0/>0 indicating horizontal/vertical current sheets, respectively. We define this “cleaner signal” regime as the stationary phase, i.e. the regime characterised by a (damped) sinusoidal signal after the transients of the impulsive phase have leaked away, and this phase starts at the formation time of the second horizontal current sheet. The black dotted lines in Fig. 1 indicate the formation times of all the horizontal current sheets (for both phases).

Note that in labelling these two regimes, i.e the impulsive phase and the stationary phase, we have adopted the terminology usually associated with the excitation and damping of trapped and leaky modes in coronal loop oscillations (see, e.g., Terradas et al. 2005, 2006; Luna et al. 2008; McLaughlin & Ofman 2008, and reference therein).

3.1. Impulsive phase

Let us consider the impulsive phase of the evolution of j_z(0,0,t), i.e. evolution over 0 ≤ t < 90 (indicated by red line in Fig. 1). We see that j_z(0,0,t) is zero (since the equilibrium is potential) until the arrival of the fast oblique magnetic shock (which has its own associated current density) at the null point, indicated by the first (negative) extrema at t = 20.6   s. This is followed by a second (local) minimum at t = 35.7   s, which is the formation time of the first horizontal current sheet (i.e. not at t = 20.6   s). The correspondence between extrema and current-sheet formation cannot be deduced from Fig. 1 alone, but can be determined by comparison with the evolution of the separatrices and contours of v. The value of j_z(0,0,t = 35.7) = −41.25j₀ = −3.3   mA is proportional to the length of the first horizontal current sheet, which is measured directly from the numerical simulation as 0.1162   L = 116.2   km (for the C = 1 simulation).

The formation time and length (at formation time) of the first horizontal current sheet is of key interest here, since this directly reflects the driver (akin to a forcing term) of the system. We investigate how the formation time and length of the first horizontal current sheet vary as functions of the initial driving amplitude (from Eq. (3)). Figure 2a shows the formation time (measured in seconds) of the first horizontal current sheet as a function of initial velocity amplitude (measured in km s^-1), where v₀ = 25.2    corresponds to C = 1 simulation. We see that the greater the driving amplitude, the earlier the first current sheet forms. This is intuitive as we expect fast oblique magnetic shocks with a larger amplitude to propagate faster and thus reach the null point more rapidly, compared to waves driven with a smaller amplitude.

Figure 2b shows how the current sheet length (measured in km) evaluated at the formation time of the first horizontal current sheet varies as a function of initial velocity amplitude. We see that the greater the driving amplitude, the longer the horizontal current sheet. The length of the current sheet is also directly proportional to the value of j_z(0,0), and so we also conclude that the greater the driving amplitude, the stronger the value of |j_z(0,0)| at the corresponding time. Again, this result is intuitive; it is the fast oblique magnetic shock that physically deforms the X-point into an osculating field structure (and ultimately into a horizontal current sheet) and thus we would expect stronger fast oblique magnetic shocks (i.e. with a larger amplitude) to deform, i.e. refract B away from the normal, and “squash” the magnetic field to a greater extent, and thus to form longer and stronger current sheets.

Finally, let us investigate the time taken to evolve from the first horizontal current sheet (at t = 35.7   s) to the formation time of the second horizontal current sheet (at t = 90.0   s for the C = 1 simulation), namely the time taken for one complete cycle (i.e. horizontal current sheet evolves to vertical, evolves back to horizontal). This time, which we refer to as the impulsive period, is calculated to be t = 54.3   s (for initial amplitude v₀ = 25.2   km s^-1 in the C = 1 system). Again, we now investigate how this impulsive period varies with the initial driving amplitude (∝ Cv₀) and this can be seen in Fig. 3. Here, we see that the greater the initial driving amplitude, the shorter the resulting impulsive period, i.e. the shorter the time taken to evolve from the first horizontal current sheet to the second. Recalling the dependency seen Fig. 2b, this means that for longer current sheets, the impulsive period is shorter. This means that the restoring force must be stronger and so we conclude that longer current sheets have a stronger restoring force. In this way, the system acts as a harmonic oscillator, i.e. the greater the displacement away from equilibrium, the stronger the restoring force.

3.2. Stationary phase

Let us now investigate the stationary phase of the oscillation seen in Fig. 1 (black line), i.e. the time evolution of j_z(0,0,t) for t₂ ≤ t ≤ 480   s, where t₂ is the formation time of the second horizontal current sheet (t₂ = 90.0   s for the C = 1 simulation). Figure 1 also has an insert showing the (same) time evolution of j_z(0,0,t) for 158.9 ≤ t ≤ 480   s, i.e. same horizontal time axis, different vertical axis. Note that the insert does not show the start of the stationary phase, only a later part of it (starting at the time of the third horizontal current sheet). The green line indicates j_final, i.e. the finite amount of current density left in the system at t = 480   s when the system has reached its final, non-potential state. For the C = 1 system, this final state occurs at t = 480   s (8 min) and j_final = 0.8165j₀ = 6.5 × 10^-5   A.

We see that there is clear oscillatory behaviour in the stationary phase and, moreover, the oscillation is exponentially decaying. This can be seen in Fig. 1 and the blue dashed lines indicate an exponentially-damped envelope max(j_z)|_stationary × e^−λt + j_final and min(j_z)|_stationary × e^−λt + j_final where λ = −0.015   s^-1 is determined experimentally.

Let us now investigate the period associated with the stationary phase, which we define as the time taken to evolve between horizontal current sheets. Specifically, we define the stationary period as the time taken to evolve from the second horizontal current sheet to the third, i.e. the first complete oscillation within the stationary phase. For the C = 1 simulation, these formation times are 90.0   s and 158.9   s resulting in a stationary period of 69.0   s (note we present results here correct to 1 decimal places, but calculate periods to a greater degree of accuracy). Similar results are obtained for alternative definitions of the stationary periods, e.g. time taken to evolve from one vertical current sheet to the next.

Fig. 3 Impulsive period, i.e. time taken to evolve from the first horizontal current sheet to the second horizontal current sheet (measured in seconds) versus initial velocity amplitude (measured in km s-1), where v0 = 25.2   km s-1 corresponds to C = 1 simulation.

Fig. 4 a) Stationary period, i.e. the time taken to evolve from the second horizontal current sheet to the third (measured in seconds) and b) jfinal, i.e. jz(0,0) evaluated at t = 480   s, (measured in milliAmps) versus initial velocity amplitude (measured in km s-1), where v0 = 25.2   km s-1 corresponds to C = 1 simulation.

We now investigate how the stationary period varies with the initial driving amplitude and this can be seen in Fig. 4a. Here, we see that the greater the initial driving amplitude, the shorter the resulting stationary period. Thus, these results are in agreement with those in Sect. 3.1, i.e. the stronger the initial driving amplitude, the longer the resulting current sheet, thus the stronger the restoring force, thus the shorter the resulting period. Coupled with the exponential decay, we see that in the stationary phase, the system acts akin to a damped harmonic oscillator.

We also measure all the proceeding oscillations in the stationary phase, i.e. time taken to evolve from the second/third horizontal current sheet to the third/fourth horizontal current sheet, and obtain similar periods of 67.8   s and 66.6   s respectively. Interestingly, this means that the stationary period appears to be slightly decreasing by roughly 1.8% per oscillation.

Finally, we investigate j_final, the finite amount of current density left in the system when the system has reached its final, non-potential state and this can be seen in Fig. 4b. For the C = 1 system, j_final = 0.8165j₀ = 6.5 × 10^-5   A which is measured at t = 480   s. In all our simulations, j_final > 0 indicating that the (final) X-point is very slightly closed up in the vertical direction, i.e. j_z > 0 is associated with vertical current sheets. This is because the (local) plasma to the left and right of the X-point is slightly hotter, since that is where the initial, strongest jet heating occurs. Thus, the existence of this thermal-pressure gradient coupled with force balance requires the final state to be non-potential.

From Fig. 4b, we see that the stronger the initial driving amplitude, the greater the value of j_final. Again, this is intuitive: fast oblique magnetic shocks with a larger amplitude will intersect to form stronger, hotter jets to the left and right of the X-point. Thus, this local plasma will be hotter at the end of the simulation, indicating a stronger thermal-pressure gradient and thus, in order to achieve force balance, a greater absolutely value of the Lorentz force, i.e. a greater value of j_final.

4. Conclusions

This paper describes an investigation into the periodicity of oscillatory reconnection, specifically oscillatory reconnection initiated by a nonlinear fast magnetoacoustic wave deforming a 2D magnetic X-point. We have solved the compressible, resistive, nonlinear MHD equations using a Lagrangian remap, shock-capturing code (LARE2D) and have followed the numerical set-up of McLaughlin et al. (2009). As in that paper, we find that the fast magnetoacoustic wave develops into a fast oblique magnetic shock wave which significantly deforms the local magnetic fieldlines, to the extent that the incoming wave deforms the null point into a cusp-like point which in turn collapses to a current sheet. The system then evolves periodically through a series of horizontal and vertical current sheets with associated changes in connectivity, i.e. the system demonstrates the mechanism of oscillatory reconnection.

The main focus of this paper is on the periodic nature of this oscillatory cycle of horizontal and vertical current sheets. For the first time, we identify two distinct phases in the oscillation: a transient, impulsive phase, encompassing the development of the first horizontal current sheet, the formation of the first vertical current sheet, and ending with the formation of the second horizontal current sheet. We define the stationary phase to begin at the formation of the second horizontal current sheet and thus this phase includes all the proceeding cyclic behaviour.

In the impulsive phase, we find that greater the driving amplitude (Cv₀) of the velocity initial condition (Eq. (3)), [a] the earlier the first horizontal current forms, [b] the longer its maximum length and [c] the greater its maximum current density. These results are intuitive since we would expect magnetic shocks with larger amplitudes to propagate faster and thus arrive at the null point more rapidly than those with smaller amplitudes. We would also expect magnetic shocks with larger amplitudes to deform the pre-existing magnetic field to a greater extent (specifically to refract B away from the normal to a greater extent) and thus, ultimately, to form longer and stronger current sheets.

We also investigate the time taken to evolve from the first horizontal current sheet to the second, which we labelled as the impulsive period. We find that the greater the initial velocity amplitude (Cv₀) the shorter the resultant impulsive period. Coupled with the results on current sheet length, this means that longer current sheets have shorter corresponding impulsive periods. In this way, the system is acts as a harmonic oscillator, i.e. the greater the displacement away from equilibrium, the stronger the restoring force, and thus the shorter the impulsive period.

For a driving amplitude of 25.2   km s^-1 (corresponding to C = 1 simulation) we measure an impulsive period of 54.3   s. We also investigate the resultant impulsive periods for driving amplitudes 6.3−126.2   km s^-1 and find associated impulsive periods in the range 46.9−65.4   s.

The stationary phase is found to be dominated by an exponentially-decaying oscillation, tending to a finite value, j_final. We also investigate the stationary period, namely the time taken to evolve from the second horizontal current sheet to the third, i.e. the first complete cycle of the stationary phase. As in the impulsive phase, we find that the greater the initial velocity amplitude (Cv₀) the shorter the resultant stationary period. This is explained just as before, the greater the initial amplitude, the longer and stronger the current sheets at each stage, and thus the greater restoring force, leading to shorter periods (compared to smaller initial amplitude, shorter resultant current sheets, weaker restoring force and thus longer periods). Hence, again the system acts as a harmonic oscillator and, coupled with the exponential decay, we relate the oscillatory reconnection mechanism to that of a damped harmonic oscillator during the stationary phase.

For a driving amplitude of 25.2   km s^-1 (corresponding to C = 1 simulation) we measure a stationary period of 69.0   s. We also investigate the resultant stationary periods for driving amplitudes 6.3−126.2   km s^-1 and find associated stationary periods in the range 56.3−78.9   s, i.e. these are high frequency (0.0127−0.0178   Hz) oscillations.

It is also prudent at this stage to ask what determines this stationary period and what determines the exponentially-decaying timescale. To this end, we can consider the work of Craig & McClymont (1991) who investigated the relaxation of a 2D X-point disturbed from equilibrium. By neglecting both nonlinear and thermal pressure effects, Craig & McClymont (1991) derived an analytical prediction for two timescales: $\begin{array}{ccc}{\mathit{t}}_{\mathrm{oscillation}}\mathrm{\approx }\mathrm{2}\ln {\mathit{R}}_{\mathrm{m}}\mathit{,} {\mathit{t}}_{\mathrm{decay}}\mathrm{\approx }{\mathit{t}}_{\mathrm{oscillation}}^{\mathrm{2}}\mathit{/}\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}& & \end{array}$where we identify t_oscillation as our stationary period and t_decay as our decay time; 1/λ. For a driving amplitude of 25.2   km s^-1, these correspond to t_oscillation ≈ 109.6   s, compared to our measured stationary period of 69.0   s, and t_decay ≈ 76.7   s compared to our measured decay time of 1/λ = 1/0.015 = 66.7   s, given that in our investigations R_m = 10⁴ and time is made dimensional using t₀ = 7.93   s. Thus, given the (relative) simplicity of the Craig & McClymont (1991) system, these estimates are in fair agreement with our results. Note, however, that these simple analytical formulae cannot predict the variation in period versus amplitude of the initial velocity driver. This suggests that nonlinear effects and thermal-pressure gradients play a crucial role, which seems reasonable given that the restoring force of oscillatory reconnection has been shown to be a dynamic competition between the thermal-pressure gradients and the Lorentz force (see Sect. 3.2 of Murray et al. 2009; Sect. 3.3 of McLaughlin et al. 2009; Fig. 7 of Threlfall et al. 2012).

It is also important to note that there is a significant difference between the impulsive period (e.g. 54.3   s for C = 1 system) and the stationary period (e.g. 69.0   s for the C = 1 system), and that in every numerical experiment we find that the stationary period is longer than the impulsive period. This indicates that different physical processes dominate in each phase (i.e. the deformation of the X-point by the shock and the importance of jet heating in the impulsive phase, and the elastic motion of the magnetic field trying to get back to equilibrium in the stationary phase) and validates our approach of dividing the whole evolution into two distinct phases.

This difference in impulsive and stationary periodicities actually has an intriguing caveat: it is important to note that when oscillatory reconnection is seen in, say, a numerical simulation, one must be careful to interpret which phase one is actually observing and to observe several oscillations, i.e. if only two periods are seen, say the impulsive period followed by a single stationary period, then one would conclude that the period was actually increasing between oscillations. A similar result would pertain in solar observations of oscillatory reconnection, i.e. the first period measured would be shorter than the proceding periods (assuming the first, i.e. impulsive, period is also observed). This is a clear prediction for the oscillatory reconnection mechanism.

It is also important to note that, as shown by McLaughlin et al. (2012), the mechanism periodically generates v_x and v_y but that these are generated exponentially damped. Thus, if such signals are detected, then they may be decaying not due to a particular damping mechanism, but due to the generation mechanism itself.

In addition to the stationary period, we also measured all the proceeding periods in the stationary phase, e.g. time taken to evolve from third horizontal current sheet to fourth, etc. Interestingly, it was found that the period very slightly decreases by roughly 1.8% per oscillation. The exact reason for this decrease in the period is uncertain and may be a numerical effect. This will be investigated further in future work.

As in McLaughlin et al. (2009), it was found that the final state (i.e. velocity zero, oscillatory behaviour ceased) is in force balance but is non-potential and a small, finite amount of current density exists in the system, j_final. The (final) X-point is very slightly closed up in the vertical direction, i.e. j_z > 0 is associated with vertical current sheets. This is because the (local) plasma to the left and right of the X-point is slightly hotter, since that is where the initial, strongest, jet heating occurs. Thus, the existence of this thermal-pressure gradient in force balance requires the final state to be non-potential. We find that the greater the initial velocity amplitude (Cv₀) the larger the value of j_final. Again, this is intuitive: fast oblique magnetic shocks with a greater amplitude will overlap to form stronger, hotter jets to the left and right of the (equilibrium) X-point. Thus, this local plasma will be hotter at the end of the simulation, indicating a stronger thermal-pressure gradient and thus, in order to achieve force balance, a greater value of the Lorentz force, i.e. a larger value of j_final.

We have presented an investigation into the periodic nature of oscillatory reconnection and have found that an aperiodic driver can naturally generate a period signal via the physical mechanism of oscillatory reconnection. We have found that the system behaves akin to a damped harmonic oscillator. Again, this is not surprising: effectively our velocity initial condition can be thought of as injecting a finite amount of energy into the oscillatory reconnection mechanism and so intuitively the resultant periodic behaviour must be finite in duration, i.e. this is a dynamic reconnection phenomena as opposed to the classical steady-state, time-independent reconnection models.

Oscillatory Reconnection may also play a role in generating quasi-periodic pulsations (see, e.g., reviews by Aschwanden 2003; Nakariakov & Melnikov 2009). Oscillatory behavior has been reported in a number of solar and stellar flare observations (e.g. Mathioudakis et al. 2003, 2006; McAteer et al. 2005; Inglis et al. 2008; Inglis & Nakariakov 2009; Nakariakov et al. 2010; Nakariakov & Zimovets 2011; Inglis & Dennis 2012; Shen & Liu 2012) but the generation mechanism responsible remains an open question.

We believe the physical mechanism of oscillatory reconnection described in this paper is a robust, general phenomenon that will be observed in other systems that demonstrate finite-duration/non-steady-state reconnection (although we have only presented a specific example of oscillatory reconnection in this paper). For example, evidence of oscillatory reconnection in 3D flux emergence simulations has been reported by Archontis et al. (2010).

Acknowledgments

The authors acknowledge IDL support provided by STFC. JAM and JOT acknowledge financial assistance from the Royal Astronomical Society. The computational work for this paper was carried out on the joint STFC and SFC (SRIF) funded cluster at the University of St Andrews (Scotland, UK).

References

All Figures

Fig. 1 Time evolution of jz(0,0,t) (measured in milliAmps) for 0 ≤ t ≤ 480   s. Red line indicates impulsive phase and black line indicates stationary phase. Insert shows the same time evolution over 158.9 ≤ t ≤ 480   s (same horizontal axis, different vertical axis). The black dotted lines indicate the formation times of all the horizontal current sheets and the green line indicates jfinal = 0.8165   j0 = 6.5 × 10-5   A. The blue dashed lines indicate an exponentially-damped envelope max(jz)|stationary × e−λt + jfinal and min(jz)|stationary × e−λt + jfinal where λ = −0.015   s-1.

In the text

	Fig. 2 a) Formation time (measured in seconds) and b) length (measured in km and evaluated at the formation time) of first horizontal current sheet versus initial velocity amplitude (measured in km s-1), where v0 = 25.2   km s-1 corresponds to C = 1 simulation.
In the text

	Fig. 3 Impulsive period, i.e. time taken to evolve from the first horizontal current sheet to the second horizontal current sheet (measured in seconds) versus initial velocity amplitude (measured in km s-1), where v0 = 25.2   km s-1 corresponds to C = 1 simulation.
In the text

	Fig. 4 a) Stationary period, i.e. the time taken to evolve from the second horizontal current sheet to the third (measured in seconds) and b) jfinal, i.e. jz(0,0) evaluated at t = 480   s, (measured in milliAmps) versus initial velocity amplitude (measured in km s-1), where v0 = 25.2   km s-1 corresponds to C = 1 simulation.
In the text