A&A 493, 227-240 (2009)
DOI: 10.1051/0004-6361:200810465

## Nonlinear fast magnetoacoustic wave propagation in the neighbourhood of a 2D magnetic X-point: oscillatory reconnection

J. A. McLaughlin1 - I. De Moortel1 - A. W. Hood1 - C. S. Brady2

1 - School of Mathematics and Statistics, University of St Andrews, KY16 9SS, UK
2 - Physics Department, University of Warwick, CV4 7AL, Coventry, UK

Received 25 June 2008 / Accepted 15 October 2008

Abstract
Context. This paper extends the models of Craig & McClymont (1991, ApJ, 371, L41) and McLaughlin & Hood (2004, A&A, 420, 1129) to include finite  and nonlinear effects.
Aims. We investigate the nature of nonlinear fast magnetoacoustic waves about a 2D magnetic X-point.
Methods. We solve the compressible and resistive MHD equations using a Lagrangian remap, shock capturing code (Arber et al. 2001, J. Comp. Phys., 171, 151) and consider an initial condition in (a natural variable of the system).
Results. We observe the formation of both fast and slow oblique magnetic shocks. The nonlinear wave deforms the X-point into a cusp-like'' point which in turn collapses to a current sheet. The system then evolves through a series of horizontal and vertical current sheets, with associated changes in connectivity, i.e. the system exhibits oscillatory reconnection. Our final state is non-potential (but in force balance) due to asymmetric heating from the shocks. Larger amplitudes in our initial condition correspond to larger values of the final current density left in the system.
Conclusions. The inclusion of nonlinear terms introduces several new features to the system that were absent from the linear regime.

Key words: magnetohydrodynamics (MHD) - waves - shock waves - Sun: corona - Sun: magnetic fields - Sun: oscillations

### 1 Introduction

It is now known that MHD wave motions (e.g. Roberts 2004; De Moortel 2005; Nakariakov & Verwichte 2005) are omnipresent throughout the solar corona (Tomczyk et al. 2007). Many solar instruments have observed various MHD wave motions in the solar atmosphere: slow magnetoacoustic (MA) waves have been seen in SOHO data (Berghmans & Clette 1999; Kliem et al. 2002; Wang et al. 2002) and TRACE data (De Moortel et al. 2000). Fast MA waves have been seen with TRACE (Aschwanden et al. 1999, 2002; Nakariakov et al. 1999; Wang & Solanki 2004) and Hinode (Ofman & Wang 2008). Non-thermal line narrowing/broadening due to Alfvén waves has been reported by Harrison et al. (2002)/Erdélyi et al. (1998) and O'Shea et al. (2003). More recently, Alfvén waves have been observed in the corona (Okamoto et al. 2007; Tomczyk et al. 2007) and chromosphere (De Pontieu et al. 2007), although these claims are currently subject to intense discussion (Erdélyi & Fedun 2007; Van Doorsselaere et al. 2008).

It is clear that the coronal magnetic field plays a fundamental role in the propagation and properties of MHD waves, and to begin to understand this inhomogeneous, magnetised environment, it is useful to look at the topology (structure) of the magnetic field itself. Potential-field extrapolations of the coronal magnetic field can be made from photospheric magnetograms, and such extrapolations show the existence of important features of the topology: null points - locations in the field where the magnetic field, and hence the Alfvén speed, is zero, and separatrices - topological features that separate regions of different magnetic flux connectivity. Detailed investigations of the coronal magnetic field, using such potential field calculations, can be found in e.g. Brown & Priest (2001), Beveridge et al. (2002), Régnier et al. (2008) or a more comprehensive review by Longcope (2005).

The propagation of fast magnetoacoustic waves in an inhomogeneous coronal plasma has been investigated by Nakariakov & Roberts (1995), who showed that the waves are refracted into regions of low Alfvén speed. In the case of null points, the Alfvén speed actually drops to zero.

McLaughlin & Hood (2004, hereafter referred to as Paper I) solved the linearised, MHD equations using a two-step Lax-Wendroff numerical scheme. They found that in the neighbourhood of a single 2D X-point, the fast MA wave refracted around and accumulated at the null point. These key results have been found to carry over from the simple 2D single magnetic null point to two null points (McLaughlin & Hood 2005) and to a more realistic magnetic configuration of a null point created by two dipoles (McLaughlin & Hood 2006a). It should be noted that the behaviour of the fast wave is entirely dominated by the Alfvén-speed profile, and since the magnetic field drops to zero at the null point, the wave will never reach the actual null for a  plasma. McLaughlin & Hood (2006b) extended the model of Paper I to include plasma pressure effects. This naturally led to the inclusion of the slow MA wave and the introduction of a layer around the null point; representing a high  environment inside and low outside. Coupling and mode conversion is observed at locations where the sound speed and Alfvén speed become comparable in magnitude (e.g. Zhugzhda & Dzhalilov 1982; Cally 2001; Bloomfield et al. 2006; McDougall & Hood 2007).

Waves in the neighbourhood of a single 2D null point have been investigated by various authors. Bulanov & Syrovatskii (1980) provided a detailed discussion of the propagation of harmonic fast and Alfvén waves using cylindrical symmetry. Craig & Watson (1992) mainly considered the radial propagation of the m=0 mode (where m is the azimuthal wavenumber) using a mixture of analytical and numerical solutions. They showed that the propagation of the m=0 wave towards the null point generates an exponentially large increase in the current density and that magnetic resistivity dissipates this current in a time related to . Craig & McClymont (1991,1993), Hassam (1992) and Ofman et al. (1993) investigated the normal mode solutions for both m=0 and modes with resistivity included. Again, they emphasise that the current builds up as the inverse square of the radial distance from the null point. All these investigations were carried out using cylindrical models in which the generated waves encircled the null point.

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). In 2D, reconnection can only occur at null points (Priest & Forbes 2000). Dungey (1953) reported that a perturbed X-point can collapse if the footpoints of the field are free to move, Mellor et al. (2002) looked at the linear collapse of a 2D null point, and Imshennik & Syrovatsky (1967) described the collapse with an exact, non-linear solution of the ideal MHD equations. 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 oscillatory reconnection, which couples resistive diffusion at the null to global advection of the outer field. An example of oscillatory reconnection generated by flux emergence within a coronal hole was recently detailed by Murray et al. (2008). Finally, Longcope & Priest (2007) investigated the diffusion of a 2D current sheet subject to suddenly enhanced resistivity. They found that the diffusion couples to a fast MA mode which propagates the current away at the local Alfvén speed.

The aim of this paper is to extend the model used in Paper I and Craig & McClymont (1991) to include finite  and nonlinear effects. To realise this, we will be solving the compressible and resistive MHD equations using a Lagrangian remap, shock capturing code: LARE2D (Arber et al. 2001). The key results from Paper I, i.e. for the linear fast wave, are that it demonstrates refraction, that the wave energy accumulates at the null, and that current density builds up exponentially at that point. We believe it is important to extend this work to include nonlinear effects since Paper I indicated a preferential topological location for (ohmic) heating and we need to see if this observational prediction persists when we consider larger (and hence nonlinear) wave amplitudes. This paper will address three main questions that naturally arise when extending Paper I into the nonlinear regime:

(1)
Does the fast wave now steepen to form shocks, and can these propagate across or escape the null?
(2)
Can the refraction effect drag enough magnetic field into the null to initiate X-point collapse or reconnection?
(3)
Has the rate of current density accumulation changed, and is the null still the preferential location of wave heating?
The paper has the following outline: the basic setup, equations and assumptions are described in Sects. 2, 3 details the nonlinear fast MA wave behaviour and the resultant oscillatory reconnection is discussed in Sect. 4. We consider different amplitudes for our initial condition in Sect. 5 and the conclusions are given in Sect. 6.

### 2 Basic equations

We consider the 2D compressible and resistive MHD equations appropriate to the solar corona:

 = = = 0 , = (1)

where is the mass density, is the plasma velocity, the magnetic induction (usually called the magnetic field), p is the plasma pressure, is the magnetic permeability, is the electrical conductivity, is the magnetic diffusivity, is the specific internal energy density, where is the ratio of specific heats and is the electric current density.

The LARE2D numerical code utilises artificial shock viscosity to introduce dissipation at steep gradients. The details of this technique, often called Wilkins viscosity, can be found in Wilkins (1980) and Arber et al. (2001).

#### 2.1 Basic equilibrium and non-dimensionalisation

The equilibrium magnetic field structure is taken as a simple 2D X-type neutral point. The aim of studying waves in a 2D configuration is one of simplicity. There are many complicated effects including mode conversion and coupling, and a 2D geometry allows us to understand and explain these behaviours better, before the extension to 3D. The initial magnetic field is taken as

 (2)

where B is a characteristic field strength and L is the length scale for magnetic field variations. This magnetic field can be seen in Fig. 1. Equation (2) is slightly different to that used in Paper I, where the authors considered . This simply represents a  rotation of our magnetic field and the key results are still valid. Note that this magnetic configuration is no longer valid far from the null point, as the field strength tends to infinity. However, McLaughlin & Hood (2006a) looked at a magnetic field that decays far from the null and they found that the key results from Paper I remain true close to the null.

We will also find it useful to consider , the vector potential, where . In 2D and for the coordinate system used in this paper, and our equilibrium vector potential is given by:

 (3)

We now consider a change of scale to non-dimensionalise all variables. Let , , x = L x*, y=L y*, , p = p0 p*, , t=t0 t*, and , where we let * denote a dimensionless quantity and v0, B, L, , p0, t0 and  are constants with the dimensions of the variable they are scaling. We then set and v0 = L / t0 (this sets v0 as a constant background Alfvén speed). We also set , where is the magnetic Reynolds number, and set , where is the plasma- at a radius L from the origin. This process non-dimensionalises Eqs. (1) and under these scalings, t*=1 (for example) refers to t=t0= L / v0; i.e. the time taken to travel a distance L at the background Alfvén speed. For the rest of this paper, we drop the star indices; the fact that all variables are now non-dimensionalised is understood.

We take the equilibrium density to be uniform, i.e. ; a spatial variation in  can cause phase mixing (Heyvaerts & Priest 1983; De Moortel et al. 1999; Hood et al. 2002). We set .

Finally, we consider the equilibrium plasma to be cold: T=0 K (i.e. ) and, hence, ignore plasma pressure effects (as in Paper I). However, as we will see below, magnetic shocks heat the plasma and so the plasma will not remain cold (see e.g. Sect. 1.5 in Priest & Forbes 2000).

 Figure 1: The equilibrium magnetic field. The red lines denote the separatrices and arrows indicate the direction of the magnetic field. The null point is located at the origin, where the separatrices intersect. Open with DEXTER

#### 2.2 Initial and boundary conditions

Equations (1) are solved numerically using a Lagrangian remap, shock-capturing code called LARE2D (Arber et al. 2001). The equations are solved computationally in a square domain with a numerical resolution of . Zero gradient boundary conditions are applied to the variables , , at the four boundaries, and  is set to zero on all boundaries, i.e. reflective boundaries. A damping region exists for and so all oscillations that enter this region are slowly damped away. 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.

As seen in Paper I, there are two natural variables to consider in our system: and . Here, and are related to the perpendicular and parallel velocity, respectively and, as seen in Paper I, their implementation naturally simplifies the governing equations, aids in MHD mode interpretation and (for  ) led to an analytical solution to the linear, cold plasma equations.

In cartesian coordinates: and . In polar coordinates: , , where and we take to be initially zero. We note that with our choice of magnetic null point (Eq. (2)), if we drive any of the velocity variables, the system will naturally develop a  dependence.

Paper I clearly demonstrates that the Alfvén speed ( ) plays a vital role. Hence, it is natural to consider either a polar coordinate system or a circular pulse. In addition, as commented by McClements et al. (2004), a disturbance initially consisting of a plane wave is refracted as it approaches the null in such a way that it becomes more azimuthally symmetric. Thus, it is appropriate to consider the evolution of azimuthally symmetric perturbations. We expect the nonlinear behaviour to be more complicated than the linear equivalent and so, in order to clearly demonstrate the differences between the two systems, in this paper we consider an initial condition in velocity, such that:

 = (4) = 0

where 2C is our initial amplitude. Initial condition (4) describes a circular, sinusoidal pulse in  . When the simulation begins, this initial pulse will naturally split into two waves, each of amplitude C, travelling in different directions: an outgoing wave and an incoming wave. In this paper we will focus on the incoming wave, i.e. the wave travelling towards the null point. The damping regions will remove kinetic energy from the outgoing waves, and so they do not influence the null. The initial condition produces a propagating disturbance that crosses magnetic fieldlines. Hence, we identify these waves as (initially) fast MA waves.

By choosing a small value for C in Eq. (4), we can recover the linear results from Paper I. This can be seen in Appendix A where we set C=0.001. For the nonlinear work described in Sects. 3 and 4, we set C=1. Different values of C will be considered in Sect. 5.

 Figure 2: Contours of for a fast wave pulse initially located at a radius r=5, and its resultant propagation at (Alfvén) times . The black lines denote the separatrices and the null point is located at their intersection (the origin). The full evolution, , is available as an mpeg animation in the online edition of the Journal. Open with DEXTER

### 3 Nonlinear fast MA behaviour

#### 3.1 Development of shocks (0  t  1)

The evolution of can be seen in Fig. 2 ( ), Fig. 5 ( ) and Fig. 6 ( ), and the full evolution, , is available as an mpeg animation in the online edition of the Journal.

From Fig. 2, we see that the initial pulse has split into two oppositely travelling wave pulses, where both waves can be seen at t=0.1. In order to clearly show the wave behaviour, we have plotted , i.e. a square subsection of our total computational box. Hence, only the incoming wave can be seen after time t=0.1. This figure can be compared with Fig. A.1, which corresponds to the linear regime.

There are two key features to note from Fig. 2. Firstly, we see that the incoming wave propagates across the magnetic fieldlines and keeps its initial pulse profile, i.e. an annulus, and the maximum amplitude remains at C=1. The annulus contracts as the wave approaches the null point and this is the same refraction effect described in Paper I. This refraction effect occurs since the (equilibrium) Alfvén speed is spatially varying.

Secondly, we note that the incoming wave pulse is developing an asymmetry: the wave peaks are propagating faster (relative to the footpoints) in the y-direction than the x-direction. Thus, the wave pulse is developing discontinuities, where in the y-direction the wave peak is catching up with the leading footpoint, and in the x-direction the trailing footpoint is catching up with the wave peak. These discontinuities can be clearly seen in Fig. 3a (trailing edge in vx[x,0]) and Fig. 3b (leading edge in vy[0,y]). We have also indicated the grid resolution (using stars) on the right-hand side of each subfigure which shows that the developing discontinuity is well resolved.

 Figure 3: a) Plot of vx(x,0) for at t=1. b) Plot of vy(0,y) for at t=1. The red line indicates the location of the null point. For , we have plotted stars to indicate the grid resolution. Open with DEXTER

 Figure 4: Our choice of initial condition for a) vx and b) vy at t=0. This choice of initial condition prescribes a background velocity profile. c) Cartoon representation of the vx (red arrows) and vy (blue) contributions to this background velocity profile, where the black lines denote the separatrices and the null point is located at their intersection. Open with DEXTER

This development of discontinuities was not reported in Paper I as it is an entirely nonlinear effect, which arises because of our choice of a velocity initial condition (Eq. (4)). In the nonlinear regime, specifying an initial velocity condition also prescribes a background velocity profile, and this profile can be seen in Fig. 4. It is this background velocity profile that leads to the development of discontinuities on the leading edges in the y-direction and on the trailing edges in the x-direction. The vx (red arrows) and vy (blue arrows) contributions to this background velocity profile can be seen in Fig. 4c. This phenomenon is most easily understood by means of a simple 1D example, details of which can be found in Appendix B. Note that the background velocity profile prescribed by Eq. (4) appears as the m=0 mode in but corresponds to the m=2 mode in cartesian components.

#### 3.2 X-point collapse (1.4  t  4)

The wave evolution between can be seen in Fig. 5. From the first row, , we see that the asymmetry seen in Fig. 2 has now lead to the formation of shock waves. From the second and third rows ( ) of Fig. 5, we see that the shocks above and below y=0 have started to overlap (starting at ). This overlap leads to the development of hot jets.

The white box in Fig. 5 at t=2 is analysed in Fig. 7a. By considering the physical quantities along a line perpendicular to the shock front, we see that there is an abrupt increase in density, temperature and, consequently, pressure (Fig. 7b). Bx increases in magnitude, whereas By decreases, and both preserve their original directions (Fig. 7c). Hence, the shock makes refract away from the normal and so we identify it as a fast oblique magnetic shock. In the idealised limit , this would be called a switch-on shock (see e.g. Sect. 1.5.2 in Priest & Forbes 2000). It is interesting to note that the shock has heated the plasma and so in these locations.

The white box in Fig. 5 at t=2.4 is analysed in Fig. 8. Figure 8a shows a contour of for , , and we see that the overlap of the shock waves forms a triangular cusp'', which we call the shock-cusp. Figure 8b shows that these hot jets heat the plasma, substantially more than the previous shock heating (compare magnitudes of T from Figs. 7b and 8b). These hot jets also significantly bend the magnetic fieldlines.

The hot jets set up new shock waves emanating from the shock-cusp, and by considering the changes of physical quantities perpendicular to the shock front (similar to before), we see that there is an abrupt increase in Bx and decrease in By, although both preserve their original directions (Fig. 8c). Thus, the shock makes refract towards the normal, and we identify this as a slow oblique magnetic shock. In the idealised limit , this would be called a switch-off shock.

The structure of the hot jet is in good agreement with that described by Forbes (1988, see his Fig. 15). In addition to the slow shocks downstream of the tip of the jet, there is evidence of slow shocks along the sides of the jet upstream of the tip (see Fig. 8b). There are also kinks in the fieldlines at the tip of the jet, indicative of a fast oblique magnetic shock. Thus, the jet heating itself is accomplished by a combination of slow and fast shocks. The abrupt jump in temperature and magnetic field at the tip of the jet is also consistent with the jump predicted by Soward and Priest (1982).

The fourth row of Fig. 5 shows that the shocks, both the fast shocks and their overlap (the tails of the jets) have reached the null point (at t=2.6). The shocks have deformed the magnetic field such that the separatrices now touch one another rather than intersecting at a non-zero angle (called cusp-like'' by Priest & Cowley 1975). That the perturbation can reach the neutral point is a phenomenon not seen in Paper I (where the perturbation reached the null after an infinite amount of time).

After time t=2.6, we see that the shock wave can now pass through the null point: entirely different behaviour to that seen in Paper I. This is more clearly seen in Fig. 9a, which shows a plot of vy (0,y) at t=2.6 (just before null point is crossed) and t=3.0 (after null has been crossed).

 Figure 5: Contours of at various (Alfvén) times between t=1.4 and t=2.8. Note that the amplitude varies substantially throughout the evolution, and hence each row is assigned its own colour bar. The black lines denote the separatrices. The white boxes at t=2 and 2.4 are analysed in Figs. 7 and 8, respectively. Open with DEXTER

#### 3.3 Oscillatory Behaviour v (4.1  t  60)

The evolution of for can be seen in Fig. 6. We can see that the evolution proceeds in two separate ways. Firstly, some of the wave now escapes the system: this propagation can be seen above and below the location of the null point. Secondly, the (deformed) neutral point itself continues to collapse and forms a horizontal current sheet. This can be clearly seen at t=4, and contours of jz are shown in Fig. 9b. Here, current density can be seen at the four slow oblique magnetic shocks, along (approximately) (i.e. a horizontal current sheet), at the locations of our two shock-cusps and due to the wave propagating away. Thus, we see that current density forms at several locations. The formation of a current sheet was not seen in Paper I (which reported the formation of a current density line at the null point).

Returning to Fig. 6, we see that this horizontal current sheet shortens in length (t=4.1). This is because the jets to the left and right of the origin heat the plasma, which begins to expand. This expansion squashes the horizontal current sheet, forcing the separatrices apart. The (squashed) horizontal current sheet then returns to a cusp-like'' null point which, due to the continuing expansion from the heated plasma, in turn forms a vertical current sheet (t=7). The evolution now proceeds through a series of horizontal and vertical current sheets and displays oscillatory behaviour (as seen by Craig & McClymont 1991). At t=60, we see that the majority of the initial wave pulse has propagated away from the X-point, and that the associated velocities are negligible. It is also interesting to note that the final state (t=60) is non-potential; the separatrices for t=0 are overplotted in red and there is clearly a small offset. This is due to the finite amount of current left in our system (see below).

The oscillatory nature of the system can be clearly seen from the time evolution of jz(0,0) shown in Fig. 10. The red/blue lines indicate maxima/minima in the system and the green line shows jz(0,0)=0.8615, which is the limiting value of the oscillation. We can see that there is no current density at x=y=0 (i.e. the location of the potential null point) before t=2.6. This is as expected since our initial wave pulse does not reach the null before this time. At t=2.6 we see a strong spike in the current density: this is the result of our nonlinear wave reaching the null point for the first time (see Fig. 5 at t=2.6). At t=4.4, the current density associated with the triangular shock-cusps (from the slow oblique magnetic shocks) reaches the origin; giving rise to a second sudden spike in current density.

At t=6.4, there is a third large spike in current density (opposite sign to previous two peaks). This indicates the formation of the first vertical current sheet. The peak shortly following this (at t=7.4) is due to the current density associated with the triangular shock-cusps from this vertical current sheet reaching the origin.

The next peak in current density occurs at t=11.8. This represents the formation of a second horizontal current sheet (see Fig. 6 at t=12). However, there is no secondary peak associated with this horizontal current sheet, unlike the previous two current sheets.

From Fig. 10 we can see that at later times, the cycle of maxima and minima continues, with each maxima (red lines) associated with the formation of vertical current sheet and each minima (blue) associated with a horizontal one. This oscillation has a decreasing period, and the time taken to go from a horizontal to vertical current sheet is shorter than vice-versa. This is because the system finds it easier to form vertical current sheets due to the asymmetric heating around the null point (see below).

We can also see that the oscillation in current density is tending to a constant value of jz(0,0)=0.8615. This is in agreement with our statement above that the final state is non-potential. This is because at the end of the simulation, the plasma to the left and right of the neutral point is hotter than above and below, due to the hot jets that formed and heated the (initially cold) plasma after t=2. Consequently, the plasma pressure in greater to the left and right of the null and hence the system finds it easier to form vertical current sheets (due to the asymmetric heating). The final state (Fig. 6 at t=60) shows that the X-point is very slightly closed up in the vertical direction, in agreement with jz tending to a small positive value.

If the final state is in force-balance, we would expect the pressure to be a function of Az, i.e. the plasma pressure should lie along contours of Az if P=P(Az). Figure 11a shows that a plot of Az(x,0) and Az(0,y) against P(x,0) (red) and P(x,0) (blue) yields a single curve, and from Fig. 11b, we can see that the agreement between the pressure and contours of Az(t=60) is excellent. This implies that although the final state is non-potential, it is still in force-balance.

 Figure 6: Contours of at various (Alfvén) times between t=2.9 and t=60. Note that the amplitude varies substantially throughout the evolution, and hence each row is assigned its own colour bar. The black lines denote the separatrices. The red lines at t=60 denote the potential separatrices. Open with DEXTER

 Figure 7: a) Contour of at t=2 for . The black line denotes the separatrix, the dotted, y= 4.382x-0.938, and dashed, y=-0.228x+0.341, white lines denote the lines perpendicular and parallel to the shock front, respectively. b) Plots of Temperature, T (red), Pressure, P (blue), and density, (green), and c) plots of Bx (red) and By (blue) perpendicular to the shock front. Open with DEXTER

 Figure 8: a) Contour of at t=2.4 for , . The black line denotes the separatrices, the dotted, y= -1.111x+1.427, and dashed, y=0.9x-0.685, white lines denote the lines perpendicular and parallel to the shock front, respectively. b) Contour of temperature at t=2.4 for , . The white lines denote magnetic fieldlines (thick white line denotes the separatrices). c) Plots of Bx (red) and By (blue) perpendicular to the shock front. Open with DEXTER

 Figure 9: a) Plot of vy(0,y) for at t=2.6 (red) and t=3.0 (blue). The perturbation has crossed the null point. The dashed line indicates vy(0,y)=0. b) Contour of jz at t=4. The black lines denote the separatrices. Open with DEXTER

 Figure 10: Plot of time evolution of jz(0,0) for . Insert shows plot of time evolution of jz(0,0) for (same x-axis, different y-axis). The dashed lines indicate maxima (red) and minima (blue) in the system and the green line shows jz(0,0)=0.8615. Open with DEXTER

 Figure 11: a) Plot of P[x,0] (red) against Az[x,0] and plot of P[0,y] (blue) against Az[0,y]. b) Contour of plasma pressure at t=60. Contours of Az(t=60) are overplotted in white and the separatrices are in black. Open with DEXTER

### 4 Reconnection

From Sect. 3.3, it is clear that we have oscillatory behaviour in our system. However, do we have reconnection? We demonstrate that reconnection is occurring via the following two pieces of evidence: firstly, in Fig. 12 we have plotted a selection of fieldlines in our system at four different time slices, where red fieldlines originate from and blue fieldlines originate from , all with equal separation of 0.01. The separatrices (which change in time) are plotted in black. In order to clearly show the results, we have plotted a subsection of our computational box: .

At t=0, all the red fieldlines from the top/bottom left corner connect to the top/bottom right corner, and all the blue fieldlines connect from the top left/right corner to the bottom left/right corner. The (potential) separatrices pass through and the origin, as expected from Eq. (2), and separate the red and blue fieldlines. Thus, if we have a change in connectivity, we should see a red/blue fieldline cross the (evolving) separatrices. Since we have reflecting boundaries and the value of the vector potential, Az, on the boundaries does not change with time, we can be confident that we are always plotting the same fieldlines in each subfigure.

 Figure 12: Selection of fieldlines in our system at times . Red fieldlines originate from and blue fieldlines originate from , all with equal separation of 0.01. The separatrices is plotted in black. The same fieldines are shown in each subfigure. Open with DEXTER

 Figure 13: a) Plot of the time evolution of Az(0,0), where we have overplotted a straight line fit: . b) Plot of the rate of change of Az at x=y=0 with straight line trend removed: , where are shown in blue. c) (Dashed line) Plot of the time evolution of the reconnection rate, i.e. , (dotted line) overplot of evolution of jz(0,0). The green line shows . Open with DEXTER

Figure 12 (t=4) shows our choice of magnetic fieldlines shortly after the formation of the first horizontal current sheet (compare to Fig. 5 at t=4). Here, we can see that one of our red fieldlines has crossed the separatrices; indicating a change in connectivity. Figure 12 (t=8) shows our choice of magnetic fieldlines shortly after the formation of the first vertical current sheet, and now we see that some of the blue lines have crossed the separatrices and the red fieldlines are once again all on the same side of the separatrices. At the end of our simulation (t=60), we see that more blue fieldlines have crossed the separatrices. Thus, throughout the evolution of our system, we have several changes in connectivity, giving us qualitative evidence for reconnection.

Quantitative evidence for reconnection in our system can be seen in Fig. 13. Here, we see a plot of the time evolution of Az(0,0), where changes in the vector potential at the origin indicate changes in connectivity. Figure 13a shows a plot of the time evolution of Az(0,0). Before t=2.6, Az(0,0)=0 as expected (no disturbance to the null point). After this time, we see that Az(0,0) oscillates but also displays a clear trend that is tending towards a straight line: . This straight line is associated with the final current density remaining in our system, i.e.:

 (5)

Hence, if our final state contains constant current density, say jC=1, then we expect Az to change linearly in time. Assuming and integrating Eq. (5) and comparing this to our straight line, we see that jC=1 = 0.8615, which is in excellent agreement with our estimate from Fig. 10.

Figure 13b shows a plot of the time evolution of , i.e. we have removed the straight line trend. We have actually plotted -A1 to aid comparison between Figs. 13a and 13b. We can clearly see the oscillatory nature of A1(0,0) and that it tends towards a constant value of -1 (which gives confidence that our straight line fit is appropriate). We can also see that A1(0,0) undergoes significant changes at (overplotted in blue), where these times correspond to the nonlinear wave reaching the X-point, the first triangular shock-cusps reaching the null and the formation of the first vertical current sheet, respectively.

Figure 13b shows us that A1(0,0) is continuously changing as it tends to -1, and thus reconnection is continuously occurring until this time ( ). Thus, whenever -A1(0,0) oscillates through -1, we are increasing or decreasing flux on one side or the other of our final state.

Finally, we can calculate the reconnection rate in our system, defined as , and we have plotted in Fig. 13c. We can clearly see the reconnection rate varies throughout the simulation and tends towards a constant value. However, from Eq. (5) we see that the reconnection rate is the same as and thus is expected to tend towards . Indeed we have overplotted the current density evolution (dotted line) and the agreement is excellent.

### 5 Amplitudes C = 0.5 and C = 2

In this section, we examine the effect of changing the amplitude, C, of our initial condition in (Eq. (4)). Here, we consider C=0.5 and C=2, and the results can be seen in Fig. 14. Comparing to Fig. 13, we can see that the overall behaviour is very similar that of C=1, where we observe oscillatory reconnection and a cycle of horizontal and vertical current sheets. Figure 14a shows a plot of the time evolution of AC=0.5(0,0), i.e. the flux function at the origin for C=0.5. We see that AC=0.5(0,0) oscillates but displays a clear downwards trend that is tending to a straight line: . As explained in Sect. 4, this straight line is associated with the final current density remaining in our system. Figure 14b shows a plot of the time evolution of , i.e. we have detrended AC=0.5(0,0). We can clearly see the oscillatory nature of AC=0.5(0,0). As before, the repeated changes in AC=0.5(0,0) implies that reconnection is continuously occurring in our system.

 Figure 14: a) Plot of the time evolution of AC=0.5(0,0), where we have overplotted a straight line fit: . b) Plot of therate of change of detrended AC=0.5 at x=y=0, where , where t=2.6 is shown in blue. c) Plot ofthe time evolution of the reconnection rate, i.e. , where the green line indicates . d) Plot of the time evolution of AC=2(0,0), where we have overplotted a straight line fit: . e) Plot of the rate of changeof detrended AC=2 at x=y=0, where , where t=2.6 is shown in blue. f) Plot of the time evolution of the reconnection rate, i.e. , where the green line indicates . Open with DEXTER

We can also calculate the reconnection rate in the C=0.5 system, defined as and we have plotted in Fig. 14c. We see that the reconnection rate varies throughout the simulation and tends to a constant value. From Eq. (5), we see that the reconnection rate is also an excellent measure of the current density in the system, tending to a final current density of jC=0.5= 0.6672. Comparing to Fig. 13c, we see that the current evolution is comparable in nature but slightly smaller in magnitude. In addition, the final current left in the system is smaller than jC=1=0.8615.

Figure 14d shows a plot of the time evolution of AC=2(0,0), i.e. the flux function at the origin for C=2. As before, we see that AC=2(0,0) oscillates and displays a clear downwards trend that is tending to a straight line: . Figure 14e shows a plot of the time evolution of . Again, we can clearly see the oscillatory nature of the system. Figure 14d shows the reconnection rate, and therefore the current density evolution in the C=2 system. Comparing to Fig. 13c, we see the current evolution is comparable in nature but slightly larger in magnitude. The final current left in the system, jC=2=1.1943, is larger than jC=1=0.8615.

We have also indicated t=2.6 in Figs. 14b and 14e, i.e. the time taken for the nonlinear wave to first reach the X-point in the C=1 system. We see that in the C=0.5/C=2 system, the wave reaches the X-point after/before this time, since the shock forms earlier and travels faster the larger the value of C.

### 6 Conclusions

This paper describes an investigation into the nature of nonlinear fast magnetoacoustic waves in the neighbourhood of a 2D magnetic X-point. We have solved the compressible and resistive MHD equations using a Lagrangian remap, shock capturing code (LARE2D). We consider a circular, sinusoidal pulse in as our initial condition in velocity (Eq. (4)), which naturally splits into two waves and we focus on the wave travelling towards the null point. Implemented damping regions remove the kinetic energy from the outgoing wave. Initially, we consider an incoming wave of amplitude C=1.

Between , we find that the incoming wave propagates across the magnetic fieldlines and keeps its initial pulse profile (an annulus). The annulus contracts as the wave approaches the null point, and this is the same refraction behaviour reported in Paper I, due to the spatial variation of the (equilibrium) Alfvén speed. We also note that the incoming wave pulse is developing an asymmetry: the wave peaks are propagating faster (relative to the footpoints) in the y-direction than in the x-direction. The incoming wave develops discontinuities, where in the y-/x- direction the wave peak/trailing footpoint is catching up with the leading footpoint/wave peak. The development of shocks around t=1 was not reported in Paper I, as it is a nonlinear effect, and arises because of our choice of a velocity initial condition: in the nonlinear regime, specifying an initial condition in velocity also prescribes a background velocity profile (Fig. 4). This phenomenon is demonstrated for a simple system in Appendix B. In addition, our initial condition in appears to excite the m=0 mode, but this corresponds to the m=2 mode in cartesian components.

For , we follow the asymmetry reported above in the development of fast oblique magnetic shock waves. By considering the physical quantities along a line perpendicular to the shock front, we found an abrupt increase in density, temperature and (consequently) pressure. It is interesting to note that the shock has heated the initially  plasma, thus creating at these locations.

At , we observed that the shocks above and below y=0 began to overlap (starting at ), forming a triangular cusp'': the shock-cusp. This led to the development of hot jets, which substantially heated the local plasma and significantly bent the local magnetic fieldlines. The hot jets set up slow oblique magnetic shock waves emanating from the shock-cusp. In addition, there is evidence of slow shocks along the sides of the jet upstream of the tip and we see kinks in the fieldlines at the tip of the jet, indicative of a fast shock. Thus, the jet heating itself is accomplished by a combination of slow and fast shocks. It is interesting to note that the jet has a bimodal structure consisting of a hot, narrow jet incased within a broader, lower temperature jet, which is a feature that is not predicted by steady-state reconnection theory.

The nonlinear wave, both the fast shocks and their overlap (the tails of the jets), reach the null point at t=2.6. The shocks have deformed the magnetic field such that the separtrices now touch one another rather than intersecting at a non-zero angle (called cusp-like'' by Priest & Cowley 1975). However, as can be seen at later times, the separatrices continue to evolve and so this osculating field structure is not sustained for any length of time. That the perturbation can reach the null point is a phenomenon not seen in Paper I. The nonlinear wave then passes through the null point, again entirely different behaviour to that seen in Paper I.

After t=2.6, The evolution subsequently proceeds in two separate ways. Firstly, some of the wave now escapes the system, corresponding to the wave that has passed through the null point. Secondly, the (deformed) null point itself continues to collapse and forms a horizontal current sheet. Current density exists in this horizontal current sheet, and also at the four slow oblique magnetic shocks, at the location of the shock-cusps and in the wave propagating away (as opposed to in the linear system, where the current density accumulated at the null point).

After t=4, the evolution proceeds as follows: the jets to the left and right of the origin continue to heat the plasma, which in turn expands. This expansion squashes and shortens the horizontal current sheet, forcing the separatrices apart. The (squashed) horizontal current sheet then returns to a cusp-like'' null point which, due to the continuing expansion from the heated plasma, in turn forms a vertical current sheet. The evolution then proceeds through a series of horizontal and vertical current sheets and displays oscillatory behaviour (as reported by Craig & McClymont 1991).

At the end of our simulation (t=60), we see that the majority of the wave pulse has propagated away from the null point, and that the associated velocities are negligible. It is also interesting to note that our final state is non-potential: there is a finite amount of current density left in our system ( jC=1=0.8615). This non-potential state occurs because the plasma to the left and right of the null point is hotter than that above and below, due to the hot jets that formed and heated the initially  plasma after t=2. Consequently, the plasma pressure is greater to the left and right and hence the final state shows that the X-point is slightly closed up in the vertical direction. This also explains why the time taken to go from a horizontal to vertical current sheet is shorter than the reverse: the system finds it easier to form vertical current sheets due to this asymmetric heating around the null. We also found that the pressure is approximately a function of the vector potential, i.e. P=P(Az), in the final state, indicating that although the final state is non-potential, it is still approximately in force balance. Of course, the system will eventually return to a potential state due to diffusion, but this will occur on a far greater timescale than that of our simulation ( ).

We provide two pieces of evidence for reconnection in our system. Qualitatively, we observed changes in fieldline connectivity (Fig. 12) and quantitatively we looked at the evolution of the vector potential at the origin. We found that Az(0,0) oscillates and displays a clear trend that was tending towards a straight line: . This straight line is associated with the final current density left in the system. Thus, since we have both oscillatory behaviour in our system and evidence for reconnection, we conclude that the system displays oscillatory reconnection (detailed by Craig & McClymont 1991, and reported more recently by Murray et al. 2008).

We then extended our study to look at the effect of changing the amplitude, C, of our initial condition (Eq. (4)). Considering both C=0.5 and C=2, we found that the overall behaviour was very similar to the C=1 study, i.e. oscillatory reconnection and a cycle of horizontal and vertical current sheets. Looking at the evolution of Az(0,0) for the two systems, labelled AC=0.5 and AC=2, we saw that both oscillated and displayed a clear downward tend. Again, this was associated with the final current left in the system, where jC=0.5=0.6672 was smaller than jC=1=0.8615, which were both less than jC=2=1.1943. Thus, we conclude that a larger initial amplitude results in a larger amount of current being left in the system at the end of the simulation, i.e. in the non-potential final state. Note that in the linear limit, jC=0.001=0, since the wave cannot reach the null point in a finite time.

Thus, we can now fully answer our three original questions (posed in Sect. 1):

(1)
Does the fast wave now steepen to form shocks, and can these propagate across or escape the null?

The nonlinear behaviour is completely different to the linear regime. We observe the formation of both fast and slow oblique magnetic shocks, and the nonlinear wave can now cross, and thus escape, the null point. We have also seen that the shocks (asymmetrically) heat the plasma such that .

(2)
Can the refraction effect drag enough magnetic field into the null to initiate X-point collapse or reconnection?

The nonlinear wave deforms the X-point into a cusp-like'' point, which in turn collapses into a horizontal current sheet. Expanding plasma to the left and right squashes the horizontal current sheet and the system evolves through a series of horizontal and vertical current sheets. Changes in the value of the vector potential at the origin as well as changes in connectivity demonstrate we have reconnection in our system. Our final state is non-potential, but in force balance. Larger amplitudes in our initial condition correspond to larger values of the final current left in the system.

(3)
Has the rate of current density accumulation changed, and is the null still the preferential location of wave heating?

The current density now accumulates at many locations, such as along horizontal or vertical current sheets, along slow oblique magnetic shocks and at the location of shock-cusps. Current density can now also leave the system since the nonlinear wave can escape the null point. This was not the case in the linear regime, where all the current density accumulated at the null point exponentially in time.

The aim of this paper is to contribute to the understanding of how nonlinear MHD waves behave in inhomogeneous, magnetised environments. Future work in this area will consider the consequences of different initial conditions, such as utilising a localised increase in pressure or internal energy in a plasma. Such an initial condition is expected to introduce slow waves and create a region around the null point where the sound speed and Alfvén speed become comparable in magnitude, and hence a layer where mode conversion could occur (e.g. Cally 2001; McLaughlin & Hood 2006b; McDougall & Hood 2007). Finally, this work will also be extended to nonlinear wave behaviour in the neighbourhood of a 3D null point, and compared to linear work by, for example, Galsgaard et al. (2003), Pontin & Galsgaard (2007), Pontin et al. (2007) and McLaughlin et al. (2008).

Acknowledgements
The authors would like to thank the referee, Terry Forbes, for his valuable comments that helped to improve this paper. J.A.M. and I.D.M. acknowledge financial assistance from the Leverhulme Trust and Royal Society, respectively. J.A.M. also wishes to thank Michelle Murray, Valery Nakariakov and Neil Symington for helpful and insightful discussions. 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).

### Appendix A: Linear regime

In this Appendix, we numerically solve Eqs. (1) subject to initial condition (4) and set C=0.001. This recovers the linear results for the fast MA wave seen in Paper I. This appendix has been included because (i) it will be useful to directly compare the linear and nonlinear systems subject to the same initial condition; and (ii) Paper I described a wave pulse driven in from one of the boundaries, which is somewhat different from the setup studied in the present paper. The results described below (simulations using the LARE2D numerical code) are identical to those in Paper I (simulations using a two-step Lax-Wendroff numerical scheme), and the results can be seen in Fig. A.1. We find that the initial condition, as expected, splits into an outgoing wave and an incoming wave. We concentrate on the incoming wave. This wave propagates across the magnetic fieldlines and travels at the Alfvén speed, and we identify it as a linear fast MA wave, in a cold plasma (). The linear fast wave keeps its initial pulse profile, i.e. an annulus, and the maximum amplitude remains at C, and the annulus contracts as the wave approaches the null point. Since the Alfvén speed is spatially varying, a refraction effect focuses the wave into the null point.

 Figure A.1: Contours of for a fast wave pulse initially located at a radius r=5, and its resultant propagation at (Alfvén) times t=0, 0.1,1 and 3. At t=0.1, we can see the initial pulse has split into two, oppositely travelling wave pulses. The black lines denote the separatrices and the null point is located at the origin. Open with DEXTER

As the length scales decrease, there is a build up of current density, growing exponentially in time, while the velocity remains finite in magnitude. Ohmic heating occurs preferentially at the null point. The refraction effect and the preferential heating at the null point are the key results for linear, plasma fast wave propagation (Paper I). Note that since the magnetic field, and hence the Alfvén speed, is zero at the null point, the fast wave cannot cross the null point and never actually reaches it.

### Appendix B: Isothermal 1D HD equations

 Figure B.1: Evolution of Eq. (B.1) at times t=0 (black lines), t=1 (blue) and t=5 (red), for three choices of initial amplitude a) D=0.001 (linear), b) D=0.4 and c) D=-0.4. Note we have plotted - v in c) to aid comparison. Open with DEXTER

Consider the isothermal, one-dimensional, nonlinear hydrodynamic equations:

which is equivalent to Eq. (1) with , . Now let us non-dimensionalise, i.e. v=v0 v*, , p= p0 p*, x = x0 x*, t=t0 t*, where we let * denote a dimensionless quantity and v0, , p0, x0 and t0 are constants with the dimensions of the variable they are scaling. We set v0 = x0 / t0 and , i.e. non-dimensionalise with respect to the sound speed. For the rest of this appendix, we drop the star indices. Letting gives:

This can be solved in general by the Method of Characteristics:

where F and G are arbitrary functions, specified by the initial conditions.

Now let us consider an initial condition in velocity:

 (B.1)

where D is our initial amplitude. This has a general solution

 (B.2)

where F(x)=v(x,0). Thus, we can see that the general solution consists of two waves, moving with speeds v +1 and v -1, which can be thought of as and due to our choice of non-dimensionalisation. We can see the evolution of initial condition (B.1) in Fig. B.1, at times t=0 (black lines), t=1 (blue) and t=5 (red). In Fig. B.1a we see that the velocity pulse (initial amplitude = D=0.001) naturally splits into two waves, as expected from the solution, each propagating in opposite directions with amplitude D / 2. which, since v is small, can be thought of as corresponds to a disturbance travelling in the increasing x-direction (to the right), whereas corresponds to a disturbance travelling in the decreasing x-direction (to the left). This choice of D represents the linear regime.

In Fig. B.1b we see the evolution of a velocity pulse of initial amplitude = D=0.4. Again, the initial pulse splits into two oppositely travelling waves: one propagating to the right (located x=5) and one to the left (located x=-5). However, we can now clearly see the nonlinear effects: the two wave pulses are beginning to tip over'' and shock, as expected for nonlinear waves. However, non-intuitively, the waves are both developing discontinuities, i.e. shock fronts, on the same faces. The pulse propagating to the right is developing a discontinuity at , where the peak is catching up with the the leading footpoint. Conversely, the discontinuity developing in the pulse propagating to the left is located at , i.e. where the trailing footpoint is catching up with the peak.

This asymmetry in the development of discontinuities can be explained by Eq. (B.2) and can be thought of as follows: when we set our initial condition in velocity, we are effectively prescribing a background velocity profile, and it is this profile that leads to asymmetry in our system (as opposed to the symmetry present in Fig. B.1a). A choice of D>0 corresponds to a positive background velocity profile that explains the development of the discontinuities on the rightmost faces of both the left and right propagating waves.

This explanation is confirmed by the evolution seen in Fig. B.1c. Here, we set D=-0.4 and have plotted -v to aid comparison with (a) and (b). We see that again the initial pulse splits into two oppositely propagating disturbances (each with amplitude D /2) and that both these pulses are developing discontinuities on the leftmost faces. This is in agreement with the above interpretation, since a choice of D<0 corresponds to a negative background velocity profile.

Interestingly, the footpoints in all three subfigures are located at the same locations. This is because the nonlinear effects are only apparent at large amplitudes and so the footpoints propagate at a speed of unity, i.e. .