© ESO, 2018

1. Introduction

Magnetic reconnection is a topological restructuring of a magnetic field causing a change in connectivity of its fieldlines (Priest & Forbes 2000). In this process, the magnetic energy is converted into kinetic energy and thermal energy of the plasma. Reconnection plays an important role in many dynamical phenomena in the solar atmosphere including the photosphere, for example Ellerman bombs (e.g. Ellerman 1917; Reid et al. 2016), chromosphere, for example penumbral microjets (e.g. Katsukawa et al. 2007) and calcium jets (e.g. Shibata et al. 2007; Morita et al. 2010), and the corona, for example X-ray jets (e.g. Shibata et al. 1992; Savcheva et al. 2007) and flares (e.g. Moore et al. 2001). Comprehensive reviews of magnetic reconnection can be found in Zweibel & Yamada (2009), Yamada et al. (2010) and Pontin (2012).

There are several models for steady-state reconnection which can be separated broadly into two categories: slow and fast. The model of Sweet-Parker reconnection (Parker 1957; Sweet 1958) has a long diffusion region between oppositely directed magnetic fieldlines. The rate at which the magnetic fields diffuse predicts an Alfvén Mach number of the inflow region, M_A, equal to the inverse root of the Lundquist number, that is M_A = S^{− 1/2}. In astrophysical plasmas, the Lundquist number can be many orders of magnitude larger than unity, resulting in a very small Alfvén Mach number, hence the name slow reconnection. Slow reconnection is not sufficient to explain observed phenomena, e.g. the energy release timescales necessary for solar and stellar flares (Aschwanden 2005, p. 410). Amendments have been made on the Sweet-Parker model to create fast reconnection models, for example Petschek (1964) and flux pile-up models (Priest & Forbes 1986). The Petschek model achieves fast reconnection by including two pairs of slow-mode shocks which act to carry flow away from the diffusion region. This results in a diffusion region that shortens as the inflow rate increases. Flux pile-up models have a rise in magnetic field energy entering the diffusion region and a slow-mode expansion. In the flux pile-up models, the diffusion region becomes long and thin, and the reconnection rate can exceed the maximum Petschek reconnection rate (Priest & Forbes 1986). These models all relate to steady-state reconnection, whereas in the solar atmosphere reconnection is time dependent.

The time evolution of a reconnection event has been studied via the incompressible Taylor problem, where a small boundary perturbation is applied to a stable slab plasma, as reviewed in Bhattacharjee (2004). This is typically referred to as forced or driven reconnection. The linear and non-linear phases of this process can be studied separately to analyse the evolution of a reconnection event and the formation of current sheets and magnetic islands (Wang & Bhattacharjee 1992). Numerical simulations of the Taylor problem use an external perturbation to develop a localised current (Fitzpatrick 2003) yielding results consistent with analytical solutions (Hahm & Kulsrud 1985; Dewar et al. 2013). The non-linear phase is only weakly dependant on the resistivity (Wang et al. 1996). The external velocity perturbation reflects in the enclosed domain producing narrow spikes in the current (Fitzpatrick et al. 2003).

When investigating reconnection in the solar atmosphere it is common to study the initialisation and development of magnetic reconnection through an internal trigger, that is various perturbations specified as initial conditions that are centred on the current sheet. These include invoking locally-enhanced resistivity or small velocity perturbations inside the current sheet (e.g. Ugai & Tsuda 1977; Arber & Haynes 2006; Malyshkin & Kulsrud 2010). However, such an initialisation may not always be applicable for reconnection in the highly dynamic solar atmosphere.

Stable current sheets can form in a number of ways throughout the solar atmosphere (Priest & Forbes 2000, Chap. 2). For example, current sheets can form around rational surfaces of a force-free magnetic field, through collapse of magnetic structures, and by compression of the plasma. Analytically, the induction equation can be solved to find a stable current sheet with velocity inflow (Priest & Forbes 2000, p. 94).

It is known that waves and flows are ubiquitous in the solar atmosphere (e.g. Edwin & Roberts 1983; Nakariakov & Verwichte 2005; De Moortel 2005; Nakariakov 2007; Tomczyk et al. 2007). Sub-Alfvénic flows are observed at coronal (e.g. Nakariakov & Verwichte 2005), chromospheric (e.g. Jess et al. 2015) and photospheric (e.g. Priest 1984, p.19) atmospheric heights, and flux cancellation has been observed at photospheric levels (e.g. Nelson et al. 2016). Thus, these waves and flows can, at some time, encounter a current sheet and this can initiate reconnection. However, the process by which reconnection can then develop as a result of this external driver is still to be studied in detail.

In the solar chromosphere, the temperature and density allow both neutral and ionised particles to exist, hence the plasma is partially ionised. This results in Cowling resistivity that acts perpendicular to magnetic field, in addition to Spitzer resistivity. Previous work investigating the implications of partial ionisation includes studies of flux emergence (Leake & Arber 2006; Arber et al. 2007), wave dissipation (Leake et al. 2005), tearing mode instabilities (Leake et al. 2012), chromospheric current sheet collapse (Arber et al. 2009), and the evolution of slow mode shocks (Hillier et al. 2016). For partially-ionised coalescing loops it has been shown that the ionisation increases the amount of reconnected magnetic flux, but the reconnection rate remains unchanged (Smith & Sakai 2008). The effects of partial ionisation can also be used to explain the existence of penumbral microjets (Sakai & Smith 2008). For strong chromospheric magnetic field strengths, the length scale of a tearing mode instability can become comparable to the kinetic length scales that have been hypothesised to be necessary for fast reconnection (Singh et al. 2015). For chromospheric reconnection the Hall term, J × B, in Ohm’s law can generally be neglected (Malyshkin & Zweibel 2011). A review of the effects of partial ionisation can be found in Zweibel et al. (2011).

This paper describes an investigation into the onset of 2D time-dependent magnetic reconnection that is triggered using an (external or non-local) sub-Alfvénic velocity driver specified perpendicular to an equilibrium Harris current sheet, removed from the centre of the domain. The physical motivation of this model is to study how reconnection events can be triggered by converging flows of the dynamic solar atmosphere itself. Reconnection start-up and development is investigated in coronal (fully-ionised), chromospheric (partially-ionised) and photospheric (weakly-ionised) atmospheric conditions. In order to understand the differences between coronal, chromospheric and photospheric reconnection, the coronal model is first presented and analysed as a reference model, and then a parameter study is performed. This allows the effects of each parameter to be studied independently. We note that this paper is studying the time-dependent onset of reconnection, and should not be confused with steady-state reconnection, which has been well studied from both an analytical and numerical description, for example (Priest & Forbes 2000).

The structure of this paper is as follows. First the computational model is described (Sect. 2). Then a reference coronal case is investigated (Sect. 3). A parameter study is then formed on this reference model, investigating the velocity dependence (Sect. 4) and plasma-β (Sect. 5). The role of ambipolar diffusion is investigated in Sect. 6 for a chromospheric case. Finally, comparisons are made regarding the onset of reconnection in the photosphere, chromosphere and corona (Sect. 7).

2. Computational set-up and diagnostics

Simulations have been performed using Lare3D (Arber et al. 2001) in the Cartesian yz-plane, with an invariant x-direction. Lare3D solves the compressible, resistive magnetohydrodynamic (MHD) equations given by for density ρ, time t, velocity v, current density J, magnetic field B, internal energy ϵ, Spitzer resistivity η, electric field E, pressure P, temperature T, neutral fraction ξ_n and reduced mass μ_m. J_⊥ is the current density perpendicular to the magnetic field. η_⊥ is an additional resistivity due to partial ionisation that acts on the perpendicular current density only. The universal constants are Boltzmann’s constant k_B = 1.381 × 10^-23 m² kg s^-2 K^-1, ionisation energy of hydrogen χ_i = 2.179 × 10^-18 J and permeability of free space μ₀ = 4π × 10^-7 m kg s^-2 A^-2. Spitzer resistivity η is uniform across the domain. Lare3D normalises parameters based on a normalisation density ρ₀, length l₀ and magnetic field strength B₀. In this paper, l₀ = 10⁴ m, B₀ = 0.002 T and ρ₀ varies in the corona, photosphere and chromosphere. The numerical resistivity has been tested and is of the order 10^-8Ωm (see Appendix A).

The neutral fraction ξ_n is calculated using the modified Saha equation (Athay & Thomas 1961). The neutral fraction is a function of density and temperature, that is ξ_n = ξ_n(ρ,T), as implemented in Lare3D by Leake et al. (2005).

Our paper investigates time-dependent magnetic reconnection in a similar set-up as the Taylor problem, however we use careful treatment of the boundaries to prevent reflections of our initial perturbation (whereas the Taylor problem has reflecting boundaries). We apply this methodology using atmospheric conditions at various heights of the solar atmosphere allowing us to consider a naturalistic onset of reconnection occurring throughout the Sun.

A sketch of the computational domain is shown in Fig. 1, where the damping region is indicated by the shaded region. The initial magnetic field is specified using a Harris current sheet defined as (11)where B_y is the magnetic field component in the Cartesian yz-plane. B_z = 0 initially. The maximum current is in the centre of the domain indicated by the dashed line at z = 0 in Fig. 1.

For a uniform temperature the pressure balance equation is solved analytically to find the equilibrium density as (12)The physical extent of the domain is 320 km × 160 km, with 512 × 512 cells and a grid size is Δy = 0.625 km by Δz = 0.3125 km. This domain size is sufficiently large in the y-direction that any shocks produced in the simulation do not reach the outflow boundaries. It also allows us to consider an isolated parameter regime in the z-direction, neglecting variations with altitude, i.e., the box size is smaller than the pressure scale height. The normalised grid dimensions are −16 ≤ y ≤ 16 and −8 ≤ z ≤ 8. A grid convergence test was performed using 256, 512 and 1024 cells and there are no significant quantitative or qualitative changes in the output for increased resolution, therefore 512 × 512 is sufficient.

2.1. Velocity driver

The reconnection is triggered by sub-Alfvénic velocity drivers specified far away from the centre of the current sheet at 4 ≤ z ≤ 8 and −8 ≤ z ≤ −4. Both velocity drivers are of the form (13)This velocity propagates across the domain dragging magnetic field energy towards the centre of the domain and triggers the reconnection.

When the velocity fronts meet, a fast-mode reflection occurs. It propagates towards the driven boundary and if left untreated bounces off the upper and lower boundaries and back into the reconnection region causing unwanted phenomena. To limit the influence of the reflection on the simulation, a damping region is specified between −8 ≤ z ≤ −4 and 4 ≤ z ≤ 8. One can damp out the unwanted perturbation by acknowledging the velocity in this region has two components, the driver v_d and the perturbation v_p, i.e. v = v_d + v_p. Kinetic energy damping is applied to v_p only and the driven velocity is applied uniformly between −8 ≤ z ≤ −4 and 4 ≤ z ≤ 8. This has the effect of damping out the unwanted perturbations whilst maintaining the imposed driven velocity.

2.2. Reconnection region

A reconnection region of width 2δ and length 2L is fitted to the domain, see Fig. 1. δ is estimated from the half width at half maximum (HWHM) of the current density J. L is estimated by the point along the centre line (z = 0) at which the outflow velocity is maximum. This defines L as the point at which the plasma ceases to be accelerated away from the reconnection region.

2.3. Electric field

The electric field E is calculated using Ohm’s law (Eq. (5)). In our computational domain there is only one non-zero component of E that is perpendicular to the yz-plane, that is E = (E_x,0,0). The electric field can be separated into a diffusion term ηJ + η_⊥J_⊥ and an advection term v × B. We note that the only non-zero component of J is J_x. This is equivalent to J_⊥ for our 2D simulation.

In 2D systems, reconnection can only occur at a null point. Therefore, in steady-state reconnection, the reconnection rate is given by the resistive electric field ηJ at the null point. This also provides a good measure of reconnection rate during time-dependent simulations of reconnection, for example tearing-mode studies. In this paper, we are interested in the onset of reconnection, namely the transition from ideal to resistive dynamics. At early times of the simulation, the system is not reconnecting despite current being present. As the simulation progresses, there is a continuous change moving from the ideal phase to the resistive phase, and the flux inflow (v × B in the inflow region) does not equate to the flux transfer (ηJ at the X-point). Therefore, using ηJ is potentially misleading since there is an imbalance of the fluxes. To capture the dynamics of the reconnection region, including the onset of reconnection, the maximum electric field in the diffusion region is used. This takes into account the motional electric field v × B in addition to the resistive electric field ηJ. By looking at the maximum electric field in the diffusion region, we gain insight into the dynamics in and around the reconnection site. Note that as the reconnection evolves, the motional electric field v × B stagnates and the resistive electric field ηJ increases. This results in the fluxes tending towards an equilibrium and means that towards the end of the simulation, the reconnection rate ηJ becomes an appropriate metric since the system is tending towards a steady-state.

3. Coronal (fully-ionised) environment

We first consider a coronal-type atmosphere driven with a sub-Alfvénic velocity. The key parameters for this study are in Table 1.

The onset process can be separated into three phases: reflection (Fig. 2), development (Fig. 3) and advection-dominated (Fig. 4). In this section each phase is discussed in order to understand the evolution of the system. More quantitative analysis of the coronal (reference) model is provided in subsequent sections.

3.1. Reflection phase

The driven velocity propagates from z = ± 4 at the local fast-mode speed towards the current sheet located at z = 0. The equilibrium Alfvén speed decreases towards the centre of the domain due to the increase in gas pressure and the decrease in magnetic field strength. The initial sound speed is constant across the domain because the temperature is constant. The fast-mode speed is a function of the equilibrium Alfvén speed and the sound speed and decreases towards the centre of the domain. Thus, the driven wave steepens as it approaches the centre of the Harris current sheet. The velocity amplitude of the driver is chosen to prevent this wave steepening into a shock. (We note that Lare3D has shock-capturing capabilities, Arber et al. 2001). When these wavefronts hit the centre of the Harris current sheet they trigger a pair of fast-mode reflections that propagate backwards towards the drivers (Fig. 2). The majority of the kinetic energy of the fast-mode reflection is damped out in the damping region and has no significant impact on the dynamics in the current sheet.

3.2. Development phase

In the wake of this fast-mode reflection, the current at the centre of the domain increases and the reconnection region begins to form (Fig. 3). There is a small amount of outflow at the start of this phase, however there is no electric field increase associated with this. Therefore this is only a fluid process as a result of the driver, as opposed to a magnetic process. The current rises as a result of the inflow pushing magnetic field together. Towards the end of this stage, magnetic reconnection begins to accelerate the plasma and there is an associated rise in electric field. Note that there is uniform resistivity in the model so diffusion is always present. The magnetic field pile-up in the inflow region creates a larger current and allows plasma to be significantly accelerated by reconnection.

3.3. Advection-dominated phase

For late times, there is clear acceleration of the plasma exiting the reconnection region (Fig. 4). The reconnection reaches a linear phase in terms of reconnection region size, with L increasing linearly with time and δ decreasing linearly with time. The reconnection event has qualities of both Sweet-Parker and Petschek models. The reconnection region is elongated and narrow like in Sweet-Parker reconnection, however there are weak slow-mode shocks that form on the interface between the inflow and outflow flow regions similar to Petschek. The inflow velocity is constant at late times, however the magnetic flux entering the reconnection region is increasing. This creates a reconnection process that behaves like the flux pile-up model of Priest & Forbes (1986). This flux pile-up on the inflow region is also present in numerical simulations of the Taylor problem and acts to enhance the current sheet (Wang et al. 1996). A colour map of the perturbation current J_x−J_x(t = 0) and streamlines of magnetic field at time t = 600Δt are shown in Fig. 5. During this phase, the physics is dominated by the advection term in Ohm’s law, as discussed in Sect. 4.2.

4. Velocity dependence

The coronal model in the previous section acts as our reference case and provides a qualitative analysis of the onset process. We now consider variation from this reference model to determine the role of various parameters on the onset of reconnection. First the dependence on the magnitude of the velocity driver will be considered. Three different velocity magnitudes are presented here. The velocities have been chosen to prevent shocks occurring on the inflow. The driver Alfvén speed is shown is Table 2. The high driver considered here is the same as the reference coronal case in Sect. 3.

4.1. Evolution of the reconnection region

The half-width δ and half-length L of the reconnection region are shown in Fig. 6. After the initial development phase, the half-length L of the reconnection region increases and the half-width δ decreases. All tested driver velocities influence L similarly until approximately t = 400Δt. There are two linear post-development phases in the diffusion half-length L. The first between 250Δt and 350Δt during which all velocities yield the same growth rate for L. However after t = 400Δt the high driver demonstrates a steeper growth rate for L. The diffusion half-width δ decreases linearly after t = 250Δt. The low driver velocity does not appear to affect δ.

The velocity driver acts to bring magnetic flux into the diffusion region, resulting in an increased current density in the inflow region. For the low velocity driver, the current density rises by only a small amount and hence the half-width δ (calculated as the HWHM of current density) remains constant. For the other velocity drivers, the current density increases sufficiently in the diffusion region to produce a smaller δ.

4.2. Electric field magnitude and evolution

The maximum electric field occurs at the centre of the domain when diffusion dominates. When advection dominates, the maximum electric field moves along the line z = 0 with the outflow. Note however the electric field variation along this line is small inside the diffusion region (significantly less than an order of magnitude).

The maximum electric field in the reconnection region shows a significant difference between the three different velocity drivers in terms of magnitudes. However, all have reasonably similar behaviour during the reflection and reconnection phases, see Fig. 7. There is a sudden change in magnitude of the electric field that is dependent on the magnitude of the velocity driver. The electric field increases by approximately four, two and one orders magnitude for the high, medium and low drivers respectively. By looking at the components of the electric field at the location of the maximum electric field in Fig. 7, one can see that the steep gradient change occurs in the v × B term. This implies the advection process dominates at late times. η is small in this case so the diffusion term is also small. The v × B term shows the three distinct phases: reflection phase, development phase and advection-dominated phase.

During the reflection phase, the wavefront from the driver hits the centreline (at time t = 80Δt) and causes a fast-mode reflection. In the wake of this reflection the current sheet fluctuates in width towards a new equilibrium. This same process occurs for all three velocities tested. These fluctuations are present in the v × B component of Ohm’s law, shown in Fig. 7 by the damped oscillations between t = 80Δt and 100Δt.

One can also identify the end of the development phase in Fig. 7 as the point where the v × B term becomes larger than the ηJ term. This occurs at a different time in all cases. Here some magnetic reconnection is occurring, however the process appears to be different in each case. To understand this, one must consider the components of v × B in the x-direction, that is v_yB_z and v_zB_y. The magnetic field in the y-direction is far stronger than the z-direction and hence the effects of v_z are magnified when looking at v × B. The high driver shows a plateau in the development region, at approximately times t = 120Δt to 200Δt in Fig. 7. The pressure of the incoming velocity is sufficient to result in a velocity directed predominantly in the y-direction at the location of the maximum electric field. The low velocity driver increases in v × B during the development phase. This is because the pressure from the incoming velocity is insufficient to direct the flow in the y-direction. Instead the velocity flows predominantly along the magnetic fieldlines. The medium driver is a combination of these two effects.

During the advection-dominated phase, the v × B term dominates the diffusion term and slow-mode shocks begin to form. The difference between the magnitudes of the advection and diffusion terms is dependent on the driver velocity. For the high driver the advection term is approximately four orders larger than the diffusion term. A larger inflow velocity produces a larger reconnection signature since more magnetic flux is transported into the reconnection region.

4.3. Resistivity dependence

A diffusion test was performed on the simulation grid determining that the numerical resistivity is of order η_n ≈ 10^-8Ωm (see Appendix A), that is far smaller than the Spitzer resistivity for the coronal case η = 10^-6Ωm (see Table 1). However, performing the simulation with a Spitzer resistivity of η = 10^-3Ωm and η = 10^-4.5Ωm did not significantly change the structure of the current density. The reconnected magnetic field along the line z = 0 does increase for the higher resistivity however remains several orders of magnitude smaller than the inflow magnetic field. This shows that changing the resistivity does have an effect on the amount of reconnected magnetic field, however the reconnected flux remains too small to have an effect on other parameters. The velocity driver carries magnetic field towards the diffusion region and hence, the rate at which magnetic flux enters the diffusion region is determined by the velocity amplitude. This acts as a limiting factor on the reconnection rate during the onset of magnetic reconnection. Towards the end of the simulation, the advection term stabilises but the current term is still increasing. One would expect that the diffusion term ηJ would eventually dominate the physics of the reconnection process as time increases beyond our simulation. It was not possible to extend the simulation time due to numerical stability issues. A weak dependence of the resistivity is also present in results from the Taylor problem (e.g. Wang et al. 1996).

5. Plasma-β dependence

To investigate the effect of the plasma-β, the high driver (Alfvén Mach number 0.003) from the previous section is used (note that this is also used in the reference coronal model in Sect. 3). Three plasma-β values have been tested: 0.002,0.1 and 1. The change in plasma-β changes the density distribution (see Eq. (12)) and hence the propagation time of the driven velocity. This changes the time at which the wavefronts collide. Note that the change in plasma-β affects the initial equilibrium sound speed, and thus the fast-mode propagation speed, via the density and wave propagation time.

5.1. Reconnection region

As before, a reconnection region of width 2δ and length 2L can be fitted to the data. The half-width δ and half-length L are shown in Fig. 8.

The half-width δ shows nearly identical values for late times. The main differences appear to be the magnitude of δ during the reflection phase, when the wavefront hits the current sheet. However all three cases show the same qualitative behaviour. The absolute magnitude of the velocity driver is the same across the three cases and hence the rate at which magnetic flux is entering the diffusion region is the same. This results in similar increases in current density in the inflow region and hence similar behaviour in δ.

The half-length L shows differences between the three tested plasma-β values. All three cases show a general increasing trend however the rates of this increase vary. In the case where plasma-β = 0.002, L increases in an approximately linear fashion. For the high plasma-β case, β = 1, the half-length L is relatively constant after the initial spike due to the waves colliding. L starts to increase for this case towards the end of the simulation. As the plasma-β increases, plasma pressure becomes more important and the outflow velocity is determined more by fluid motion, resulting in a slower increase in L for higher plasma-β values.

5.2. Electric field: variation with plasma-β

The maximum electric field in the reconnection region is shown in Fig. 9. There are significant differences between the electric field signature for different plasma-β values.

The low plasma-β case (β = 0.002) was analysed in Sect. 3 as the reference coronal model and Sect. 4 as the high velocity case. A peak in electric field occurs during the reflection phase, when the two wavefronts collide at approximately t = 80Δt. After this, the electric field is relatively constant throughout the development phase. Finally there is a rise of four orders of magnitude in electric field as the reconnection phase starts up at approximately time t = 200Δt in Fig. 9. This is due to the v × B term in Ohm’s law becoming far larger than the ηJ term as was seen in Fig. 7 and discussed in Sect. 4.2.

The other plasma-β values behave quite differently. As the plasma-β value increases, the plasma pressure becomes more important. This results in a larger electric field spike when the two wavefronts collide. Following this, the behaviour is different than in the reference model (black line in Fig. 9). The development phase occurs at a later time and has multiple peaks as oppose to being fairly constant. The advection-dominated phase is also very different from the reference coronal model (Sect. 3); there is no large, smooth increase, instead the electric field is fairly constant and noisy. For the plasma-β = 1 case the advection term is far larger than the diffusion term. For the plasma-β = 0.1 case, the advection and diffusion terms are of similar orders of magnitude.

Increasing the plasma-β appears to inhibit the reconnection. The electric field only demonstrates a sharp rise, indicative of reconnection, for high plasma-β values when the inflow velocity approaches the Alfvén speed.

6. Ambipolar diffusion

In the solar chromosphere, the plasma is partially ionised. This introduces Cowling resistivity η_⊥ that acts perpendicular to the magnetic field and is calculated by Eq. (9). η_⊥ appears in the internal energy equation (Eq. (4)) and Ohm’s law (Eq. (5)). Ambipolar diffusion refers to the η_⊥J_⊥ term in Ohm’s law (Eq. (5)).

To investigate the effects of Cowling resistivity, simulations were performed using chromospheric conditions (see Table 3). Three initial temperature values have been chosen at 7200, 9200 and 10 300 K such that the respective neutral fraction ξ_n is 0.9, 0.5 and 0.1. The driver velocity is chosen to avoid shocks in the inflow region at all three temperatures. We note that the change in temperature changes the propagation speed and hence the time at which the wavefronts collide.

6.1. Reconnection region

The development and evolution of the reconnection region (half-length L and half-width δ) are shown in Fig. 10. At late times, all three cases show a linear increase in L and a linear decrease in δ. There is only minimal difference between the considered cases in terms of the reconnection region half-width and half-length. The atmospheric conditions are fairly similar in the three cases so the behaviour of the start-up of reconnection is similar across the three cases, that is a narrowing and elongating reconnection region with comparable gradients across the three cases.

Previous results indicate a correlation between the current layer thickness and neutral fraction (e.g. Yamada et al. 2006). Our results in Fig. 10b show similar rates of change in the width δ during the onset of reconnection. This difference in these results originates from the stage of reconnection considered. Ambipolar diffusion is usually considered in the context of a tearing mode instability, where reconnection is achieved directly using an interior trigger (e.g. Zweibel 1989). In this paper we consider the early onset of reconnection, using an exterior trigger in the form of a velocity driver. It was established in Sect. 4.3 that the resistivity plays a minimal role in the onset of reconnection in our configuration and since the ambipolar diffusion manifests in the equations as an additional resistivity term one would expect that the ionisation fraction would have little role in the onset of reconnection. As time advances, resistivity begins to play a more dominant role and the effects of partial ionisation would become more pronounced. Note that outside the reconnection region, ambipolar diffusion shapes the inflowing current density (Sect. 6.3).

6.2. Electric field

Figure 11 shows the evolution of the maximum electric field in the reconnection region. Again there is very little difference between the cases. All three cases show near identical behaviour, both quantitatively and qualitatively. The time off-set of the results is due to the difference in temperature altering the propagation speed of the driven wavefronts, and thus the time at which the wavefronts collide. The large increase in electric field occurs when the v × B term becomes larger than the ηJ term in Ohm’s law.

6.3. Current distribution

For a partially-ionised plasma, the distribution of the current should follow a power law of −2/3 (Brandenburg & Zweibel 1994). This power law region allows sharper features to form in partially-ionised plasmas, for example the solar chromosphere. Figure 12 shows the distribution of the current J_x in the z-direction at time t = 500Δt along the line y = 0, that is the distribution of the current along the inflow direction. At this time the diffusion region is located at z ≈ 0.2, 0.3 and 0.4 for ξ_n = 0.1, 0.5 and 0.9 respectively (Fig. 10b). All three cases demonstrate a power law region with roughly the expected −2/3 gradient. This indicates that the inclusion of partial ionisation alters the distribution of the current along the inflow direction. In contrast, there is no power-law region in the distribution of the current for the coronal case. Figure 12 also shows that there are differences between the distribution of current for different neutral fractions. As time advances, one would expect these to become more pronounced as the resistivity begins to play a larger role in the reconnection.

6.4. Heating

There is a slight difference in the maximum percentage temperature increase when changing the neutral fraction. For a neutral fraction of ξ_n = 0.5 and ξ_n = 0.1 the temperature in the centre of the domain increases by approximately 280 K. For ξ_n = 0.9 there is temperature increase of approximately 410 K. The perpendicular resistiviy depends on ξ² (Eq. (9)) and hence the additional diffusion due to the η_⊥J_⊥ term is significantly higher when ξ_n = 0.9 than when ξ_n = 0.5 or ξ_n = 0.1, producing more of a temperature increase. This implies that current sheets lower in the solar atmosphere (i.e. photosphere or lower chromosphere) produce more local heating. The distribution of the heated plasma is roughly identical in all cases and confined to the reconnection region.

7. Onset of reconnection: photospheric vs. chromospheric vs. coronal environments

We can compare the onset of a reconnection event occurring at different atmospheric heights now that the different parameters have been investigated separately. The driver magnitude in all cases has been chosen such that the Alfvén Mach number on the inflow approaches but never exceeds unity. The key parameters in the different cases are in Table 4.

The output data for the simulations have been scaled such that the time for a fast-mode wave to propagate across the domain is approximately 100Δt. This means the wavefronts collide at time t ≈ 50Δt. The values of Δt are 100, 38.5 and 24 s for photosphere, chromosphere and corona respectively.

7.1. Size and extent of the reconnection region

The half-length L and half-width δ of the reconnection region are shown in Fig. 13. In the figure time is scaled by Δt and the photospheric case appears to narrow at the fastest rate. However, by looking at the rates of the linear change in δ and using Δt from Table 4, one can calculate the rate Δδ/ Δt, that is the speed at which the reconnection region is narrowing. These values are approximately 0.213, 0.4132 and 0.208 km s^-1 for the photosphere, chromosphere and corona respectively. Therefore the fastest rate of change is in the chromosphere. The perpendicular diffusion in the chromospheric case produced a sharp current structure and hence the half-width δ (calculated as the HWHM of the current density) narrows at a faster rate than the other atmospheres.

7.2. Electric field evolution

Figure 14 shows the electric field in each case normalised by its electric field at time t = 0. The three cases demonstrate very different behaviour during the advection-dominated phase. There are orders of magnitude difference between the maximum electric field in each case. In the photosphere, the electric field exhibits exponential behaviour (Fig. 14) but has a very low gradient. We observe two separate exponential phases with different gradients: one between 180Δt and 280Δt corresponding to the tail end of the development phase, and one from 280Δt and 400Δt with a slightly steeper gradient that corresponds to the advection-dominated phase. The chromosphere has a very flat electric field until t = 270Δt when the advection-dominated phase starts and the v × B term dominates. The corona has the highest electric field gradient, again corresponding to the v × B term dominating the electric field.

The electric field signature also depends on the driver velocity. This was analysed for the coronal case in Sect. 4. The electric field for different driver velocities in the photosphere and chromosphere are shown in Figs. 15a and 15b respectively. The different driver velocities are shown in Table 5 for the photospheric case and Table 6 for the chromospheric case. In the photosphere, only the fastest driver yields a steep gradient in electric field at late times. This occurs when the advection and diffusion terms are of the same order. In the chromosphere, one can produce behaviour in the electric field similar to both photospheric and coronal cases depending on the magnitude of the velocity driver. The fastest driver has comparable behaviour to the coronal case: a steep gradient in electric field and advection dominates at late times. However for slow drivers there is very little change in the electric field and the diffusion dominates, similar to the slow photospheric case. Therefore reconnection occurring at chromospheric levels is highly dependent on the magnitude of the velocity driver.

7.3. Heating: change in internal energy

The heating can be calculated in kelvin. In all cases the maximum heating occurs at the centre of the domain. The original values of the temperature are 6000, 9200 and 10⁶ K for the photosphere, chromosphere and corona (Table 4). The temperature change at the end of the simulation is approximately 186, 267 and 39 456 K in the three cases. This corresponds to a 3.1%, 2.9% and 3.9% rise in temperature over the simulation for the photosphere, chromosphere and corona.

The shape of the heating region also varies in the different cases, see Fig. 16. In all cases, the bulk of the heating is created in the diffusion region. For the photospheric case, the heating is a combination of two effects: Ohmic heating generated in the diffusion region is carried out along the double Y-shaped slow-mode shocks that exist on the interface between inflow and outflow regions, and some heat is created directly by the shock. This double Y-shape is similar to the shape of the outflow jet and termination shocks in Forbes (1988). The coronal case produces a fairly evenly distributed heating region with very little heat being carried out along the magnetic field lines. The chromospheric case has a very localised heating region of approximately the same size as the reconnection region. This is due to the power-law region in the current density creating a highly localised heating region.

Reconnection events generated in this setup deposit heat very differently in the surrounding plasma. In the coronal case, the heat is spread uniformly in a large area around the reconnection region. Chromospheric reconnection creates a heating signature that is entirely localised to the diffusion region. In the photosphere, the heat is distributed by shocks and flows around the reconnection site.

8. Conclusions

In this paper the onset of 2D magnetic reconnection has been investigated at different atmospheric heights of the Sun, namely the solar corona, chromosphere and photosphere, by using an external sub-Alfvénic velocity driver specified perpendicular to a Harris current sheet. This allows us to investigate the early behaviour of magnetic reconnection in a naturalistic manner. As waves and flows are ubiquitous in the solar atmosphere, these waves and flows can, at some time, encounter a current sheet and this can initiate reconnection as described in this paper. This is the physical interpretation of our velocity driver. Furthermore, by choosing different physical parameters, we can investigate the signatures of the reconnection onset as a result of this driver at different solar atmospheric layers, namely the photosphere, chromosphere and corona.

The process has been separated into three phases: reflection, development and advection-dominated. When the velocity wavefront hits the centre of the equilibrium Harris current sheet there is a fast-mode reflection. This causes the current sheet to fluctuate in width. During the development phase, the current density starts to rise due to the inflow generated by the driver and magnetic field lines begin to reconnect. Finally during the advection-dominated phase there is clear acceleration of the plasma exiting the reconnection region. With reconnection initiated with this type of velocity driver, the reconnection appears to be fairly independent of the resistivity; the current density J_x and outflow velocity v_y remain identical when the resistivity is changed. Note that the specified resistivity is resolved in our grid (Appendix A). The fundamental reconnection structure is formed from the collision of the wavefronts. Acceleration of the plasma occurs, however the advection term in Ohm’s law is dominant over the diffusion term. The rate at which magnetic flux can enter the diffusion region is determined by the velocity amplitude. This implies that for this configuration in the solar atmosphere, a large increase in electric field can only occur for large velocity amplitudes.

For the photospheric case, a large increase in electric field was only obtained as the velocity flowing into the reconnection region approached the Alfvén speed. The high gas pressure results in a very narrow current peak at the centre of the domain and hence a narrow reconnection region. This high gas pressure also means that more energy is required to accelerate the plasma, producing a weak electric field. The heat is carried out along slow-mode shocks that exist on the interface between the inflow and outflow regimes. The dense plasma limits the rate at which magnetic flux can enter the reconnection region during the onset of magnetic reconnection. The implication for the solar photosphere is that reconnection is less readily achieved, compared to other atmospheric layers, and requires a very large inflow amplitude.

In the corona, the electric field shows a far larger increase, compared to the other atmospheric layers. The advection term becomes four orders of magnitude larger than the diffusion term in Ohm’s law at late times. The low plasma-β means that plasma can easily be accelerated away from the reconnection region, resulting in a large velocity and hence a large advection term.

The chromospheric case has atmospheric conditions somewhere between the corona and the photosphere. The components of Ohm’s law behave much like the coronal case where the advection term is far larger than the diffusion term at late times. However the difference is approximately two orders of magnitude in the chromosphere, compared to four orders of magnitude in the corona. The inclusion of partial ionisation results in a power-law region in the distribution of the current on the inflow that was not present in the coronal or photospheric cases. However, the partial ionisation does not appear to have a significant effect on the outflow velocity or electric field signatures when using this driver. The initial onset of reconnection in the solar chromosphere behaves similarly to the fully-ionised coronal case, however the effects of partial ionisation become more important as the reconnection evolves.

There is a fundamental difference in the heating region in the three cases. The corona has a large regular-shaped heating region around the reconnection region. The low density in the corona allows the heat to spread fairly uniformly around the reconnection region. The photosphere appears to have heating that is guided along the double Y-shaped slow-mode shocks, surrounding the reconnection region. The chromospheric case has a very localised heating region almost the same size as the diffusion region. This is a result of the partial ionisation; the structure of the current density is a peak in the centre, surrounded by a power law region. This narrow peak in J_x results in localised Ohmic heating and hence a heating region that is localised to the diffusion region. The heat generated by reconnection in the corona is distributed over a large area, whereas reconnection at chromospheric and photospheric levels produced localised heating as a signature of reconnection.

The reconnection event behaves Sweet-Parker-like in some ways, and Petschek-like in others. The reconnection region is narrow and elongated, as in Sweet-Parker. However weak slow-mode shocks form on the interface between the inflow and outflow and the outflow velocity is closer to the Petschek speed, where . The inflow magnetic field is stratified, thus the general behaviour is similar to the flux pile-up models (Priest & Forbes 1986).

In this paper, we consider the onset of 2D magnetic reconnection. A test was also performed in 2.5D using a constant magnetic field in the invariant x-direction. The main result of this was that the Alfvén speed at the centre of the domain does not reduce to zero, resulting in less flux pile-up on the inflow region. Qualitatively the results were the same between the 2D and 2.5D cases. Expanding to 3D opens up the potential for more complex configurations, however we expect that qualitatively the results would be comparable to this 2D study if they occur is a similar setup.

Reconnection is ubiquitous in the solar atmosphere and plays a key role in many phenomena occurring at photospheric, chromospheric and coronal levels. In this paper, numerical simulations were performed investigating the transient onset of magnetic reconnection at different atmospheric heights. This was achieved by using an external driver to trigger reconnection from an equilibrium Harris current sheet. The amount of flux that can enter the diffusion region is determined by the velocity of the inflow. A lower plasma-β means a low gas pressure and hence magnetic field can be pushed together using lower velocities. This creates a larger change in the electric field signature of the reconnection event. This sharp rise in electric field is obtained for all tested velocity drivers. However the magnitude of this rise is dependent on the amplitude of the velocity driver; a higher amplitude driver produces a larger rise in electric field. The high plasma pressure (relative to the magnetic pressure) in the photospheric case leads to a narrower current sheet, compared to coronal and chromospheric cases. For a high-β plasma, a sharp rise in electric field is obtained only as the inflow into the reconnection region approaches the Alfvén speed. The rate at which magnetic flux enters the diffusion region is the key parameter determining the onset of magnetic reconnection.

Acknowledgments

The authors acknowledge IDL support provided by STFC. A. Hillier is supported by his STFC Ernest Rutherford Fellowship grant number ST/L00397X/2. J.A. McLaughlin acknowledges generous support from the Leverhulme Trust and this work was funded by a Leverhulme Trust Research Project Grant: RPG-2015-075.

References

Appendix A: Resistivity dependence

Numerical resistivity:

the numerical resistivity on our grid was calculated by specifying a Harris current sheet as the initial condition and running Lare3D specifying η = 0. The current sheet will collapse due to the numerical resistivity η_n. The numerical resistivity can then be calculated by solving the diffusion equation: (A.1)Simulations are performed with two invariant directions and a grid density 512 cells over 32 Mm, corresponding to the spatial resolution used in the rest of the paper. This simulation was performed for a physical time of 20 000 s (corresponding to approximately 15 000 numerical iterations). The initial and final states are used to calculate η_n from Eq. (A.1). The numerical resistivity calculated is of the order η_n = 10^-8Ωm.

Order of magnitude argument:

the calculated value of numerical resistivity is very small (). An order of magnitude analysis can be used to compare this to previous work. The resistivity (using normalisation values) is effectively: (A.2)Coronal simulations in the literature that use Lare3d typically use a length scale of 10⁶ m and have a resistivity of order −3 to −4Ωm. Here we are interested in localised reconnection so our length scale is smaller, l₀ = 10⁴ m. This effectively means that we can use a resistivity that is of order −5 to −6Ωm without the output becoming dominated by numerical resistivity.

Physical resistivity:

the coronal simulation was repeated using a specified resistivity of η = 10^-3Ωm. (The coronal resistivity used elsewhere in this paper is 10^-6Ωm.) The only variable this changes is B_z on the outflow plane: it becomes larger as the simulation advances. This means that there is a slight change in behaviour when the resistivity is increased, further validating our calculation of numerical resistivity. The variable B_z on the outflow region is significantly smaller, by several orders of magnitude, than the inflow B_y. As a result, the change in B_z, whilst validating our choice of resistivity, does not correspond to changes one may expect from changing the resistivity. The reconnection would have to become significantly more developed in order to notice changes in other parameters from changing the resistivity. This further suggests that the onset of reconnection is dominated by advection effects.

B_z scaling with resistivity:

Fig. A.1 gives the reconnected magnetic field (B_z) along the z = 0 line at the time t = 600Δt for three different resistivities (η = 10^-3, 10^-4.5 and 10^-6Ωm). It is clear from this figure that the amount of reconnected flux is different for each resistivity value. It is possible to create a simple model to estimate the change in the reconnected flux. Firstly we assume that the current is dominated by the B_y field component, so we set J_x ≈ −∂B_y/∂z. As the early stages of reconnection have weak flows, then we can estimate the temporal evolution of B_z is dominated by diffusion, i.e. ∂B_z/∂t ~ −η∂J_x/∂y. If J_x remains relatively constant (and determined by the inflow magnitude and structure) once the reconnection has started, then B_z(t) ~ −ηt∂J_x/∂y. Using this we can estimate the change in magnitude of the reconnected flux (ϵ) for a given change in resistivity as: (A.3)From this we would expect that if we calculate the value of ϵ for η₁ = 10^-3Ωm and η₂ = 10^-4.5Ωm, or η₁ = 10^-4.5Ωm and η₂ = 10^-6Ωm we would get a difference of approximately 10^1.5 ~ 30 between them. For the simulations shown in Fig. A.1 with η₁ = 10^-3Ωm and η₂ = 10^-4.5Ωm we get ϵ = 6.5 × 10^-5, and with η₁ = 10^-4.5Ωm and η₂ = 10^-6Ωm we get ϵ = 1.9 × 10^-6 which differ by a factor of ~ 30 as we predict from our simple model. This test does imply that we are able to see differences in the reconnection behaviour down to very small resistivity values, highlighting that we are not in a regime dominated by numerical resistivity for the simulations shown in this paper. We would like to note that we are not claiming that this would hold into the fully developed stage of reconnection, because in the later stages of reconnection when the current sheet dynamically thins the smaller the resistivity value the sooner the system will become dominated by numerical resistivity.

All Tables

All Figures