FAST MAGNETIC RECONNECTION IN THE SOLAR CHROMOSPHERE MEDIATED BY THE PLASMOID INSTABILITY

It is well known that magnetic reconnection plays an important role in the dynamic solar corona, where it is a key process in eruptions (flares, filament eruptions, coronal mass ejections), jets, and other phenomena. In the lower solar atmosphere, magnetic reconnection is also very common in different forms of activity, such as microflares (e.g., Fang et al. 2006a; Brosius & Holman 2009; Gontikakis et al. 2013), the flux cancellation process in the photosphere (e.g., Zwaan 1987; van Ballegooijen & Martens 1989; Cameron et al. 2007), chromospheric jets (e.g., Chae et al. 2003; Shibata et al. 2007; Liu et al. 2009; Shen et al. 2011, 2012; Morton 2012; Bharti et al. 2013), and temperature enhancements in Ellerman bombs (e.g., Zachariadis et al. 1987; Denker 1997; Dara et al. 1997; Qiu et al. 2000; Georgoulis et al. 2002; Fang et al. 2006b). Magnetic reconnection probably plays an important role as the energy source for the heating of the chromosphere and corona (e.g., Parker 1972; Sturrock 1999; Klimchuk 2006). The source of the mass and energy in the solar wind may also be related to magnetic reconnection in the lower solar atmosphere.

Compared to the solar corona, the mass density and plasma β in the chromosphere and photosphere are higher, and the temperature is much lower, mostly below 10⁴ K, except in the upper chromosphere. Therefore, the plasma in this height range is only partially ionized with an ionization degree of about 10⁻⁴–10⁻¹ (e.g., Vernazza et al. 1981; Pneuman et al. 1986; Khomenko et al. 2008). The magnetic diffusivity is of order η ∼ (10³–10⁴) m² s⁻¹ in the range from the upper photosphere to the middle chromosphere, so that the Lundquist number (magnetic Reynolds number) for a length scale L ∼ 1 Mm and Alfvén velocity v_A ∼ (10–100) km s⁻¹ is of order S = Lv_A/η ∼ 10⁶–10⁸. The classical Sweet–Parker model for a current sheet of this length yields a reconnection rate γ_SP  ∼  S^−1/2 ∼ 10⁻⁴–10⁻³, far smaller than the reconnection rates inferred for events of chromospheric activity. For example, Nishizuka et al. (2011) estimated γ ∼ 0.02–0.1 for chromospheric anemone jets from their lifetime and the estimated Alfvén crossing time. The same estimate for type II spicules, which have lifetimes of 10–150 s, sizes of ∼200 km, and ambient Alfvén velocities based on the total density of ∼100 km s⁻¹ (de Pontieu et al. 2007; Moore et al. 2011), yields γ  ∼  0.01–0.2. Such and even larger disparity is typically found for weakly ionized astrophysical plasmas like, for example, the interstellar medium and protostellar and protoplanetary disks (Zweibel et al. 2011 and references therein).

Many numerical simulations applied to phenomena of solar activity employ anomalous resistivity enhancements to obtain fast reconnection. However, the exact forms of anomalous resistivity used in magnetohydrodynamic (MHD) models are not directly deducible from the kinetic theory and simulations of real physical systems, and therefore some simplifying assumptions are still required for the model of anomalous resistivity (see, e.g., Ni et al. 2012a and references therein).

Steady-state reconnection rates in weakly ionized plasma can be significantly enhanced above the Sweet–Parker value by radiative cooling and by ambipolar diffusion; however, both are efficient only if the guide field is weak. Radiative cooling reduces the pressure in the current sheet, thus allowing the sheet to thin and, correspondingly, the current density to rise. This implies higher Lorentz forces accelerating the reconnection outflow, resulting in an enhancement over the Sweet–Parker rate of order A^1/2, where A is the compression factor of the sheet due to the cooling (Dorman & Kulsrud 1995; Uzdensky & McKinney 2011). In the incompressible strong guide field case, the cooling leads to a weaker enhancement of the reconnection rate through the higher resistivity, but in a chromospheric environment this can hardly reach an order of magnitude. Ambipolar diffusion in the current sheet decouples the ions and neutrals, so that only the ion pressure is available to balance the Lorentz force, which results in a strong current sheet thinning for weak ionization, allowing rapid flux annihilation (Brandenburg & Zweibel 1994). The enhancement over the Sweet–Parker reconnection rate can be very large in the special case of antiparallel field (Heitsch & Zweibel 2003a); however, already a tiny guide field component suppresses this fast reconnection regime (Heitsch & Zweibel 2003b).

In partially ionized plasma, fast reconnection also results from the recombination effect under nonequilibrium ionization conditions. The numerical results from a two-fluid model in Leake et al. (2012, 2013) indicate that the reconnection rate can been increased to above 0.05 solely as a result of the decoupling of neutrals from ions. Therefore, when the neutral–ion collisional mean free path exceeds the current sheet width, the two-fluid model including the recombination effect should be applied in studying the reconnection process.

Reconnection also becomes fast when the current sheet thins down to the ion inertial scale d_i, where the Hall term is relevant (Mandt et al. 1994; Malyshkin & Zweibel 2011). However, a dynamic regime of reconnection that involves plasmoids (magnetic islands) realizes high reconnection rates already at current sheet widths typically several orders of magnitude larger than d_i; this is now commonly referred to as the plasmoid instability (Loureiro et al. 2007; Daughton et al. 2009; Bhattacharjee et al. 2009). The motion of the plasmoids in the current sheet is largely constrained by ideal MHD (Finn & Kaw 1977); hence, it is fast and largely independent of the resistivity, forming a highly efficient internal engine of fast reconnection. The instability cascades to smaller scales such that current sheet sections between newly formed islands break up as well; the resulting high current densities in the increasingly thin current sheet fragments also facilitate a high reconnection rate. Moreover, this mechanism operates for current sheets with and without a guide field. Its high relevance for reconnection in the solar corona, especially for flares, has recently been demonstrated in a number of simulation studies (e.g., Murphy et al. 2012; Mei et al. 2012; Guo et al. 2013). The plasmoid instability has recently been found to occur also in simulations of reconnection in weakly ionized plasma (Leake et al. 2012, 2013). Plasmoids indicating the action of the instability were detected in laboratory experiments (e.g., Dong et al. 2012), in the solar corona in the course of eruptions (e.g., Lin et al. 2008; Savage et al. 2010; Milligan et al. 2010; Takasao et al. 2012; Liu 2013), and in chromospheric jets (Shibata et al. 2007; Singh et al. 2012).

The plasmoid instability occurs only if the Lundquist number and the aspect ratio L/δ of the current sheet are sufficiently large, where δ is the current sheet width. Both are related through the Sweet–Parker scaling for the current sheet width, δ_SP ∼ S^−1/2L. For fully ionized plasma, the critical Lundquist number, S_cr, typically lies in the range of several 10³ to several 10⁴, decreasing for increasing plasma β in the inflow region (e.g., Bhattacharjee et al. 2009; Huang & Bhattacharjee 2010; Ni et al. 2010, 2012b, 2013). The minimum aspect ratio correspondingly is of order ∼10². Fast reconnection at a rate $\gamma \sim S_\mathrm{cr}^{-1/2}$ was found for S > S_cr. Leake et al. (2012, 2013) find the onset of the plasmoid instability in partially ionized plasma in basic agreement with this value.

Since the characteristic Lundquist number estimated for chromospheric reconnection events is much higher than the critical Lundquist number for onset of the plasmoid instability, and since this process yields high reconnection rates for weak and strong guide fields, we focus on the plasmoid instability in chromospheric reconnection in the present investigation, considering both the regimes of strong and of vanishing guide field. Parameters characteristic of the middle solar chromosphere and a low plasma β in the inflow region are considered. The energetically and dynamically important processes of radiative cooling and ambipolar diffusion are included, and their effect on the reconnection rate is compared with the effect of the cascading plasmoid instability. The formation of slow-mode shocks associated with the X-points between the plasmoids is demonstrated. Finally, we compare the obtained reconnection rates, upward reconnection outflow velocities, and temperature enhancements with the values typically seen in or inferred for events of chromospheric activity.

The model and numerical method are described in the following section. In Section 3 we present our numerical results and compare them with observations of chromospheric activity phenomena. A summary and discussion are given in Section 4.

We only consider the hydrogen gas, so that the plasma is composed of neutral hydrogen atoms, ions (protons), and electrons. In general, this system is described by the three-fluid MHD equations. However, a single-fluid approximation adding up all three equations can be used if the collisional coupling between ions and neutrals is strong, which is valid on length scales exceeding the neutral–ion collisional mean free path λ_ni. The relevant length scale here is the current sheet width δ. We find below that parameter values characteristic of the middle chromosphere yield critical current sheet widths for onset of the plasmoid instability several orders of magnitude above λ_ni, indicating that the onset and initial evolution of the instability in this part of the solar atmosphere can generally be described using the one-fluid equations. Since the instability tends to cascade to small scales comparable to the ion inertial length, which can be of similar magnitude to λ_ni or smaller (depending on the ion density), it is clear that a complete description of plasmoid-mediated reconnection in the chromosphere will require a multifluid treatment (such as recently developed by Leake et al. 2012, 2013). However, the present investigation demonstrates that the plasmoid instability occurs and yields high reconnection rates already within the validity range of the one-fluid approximation adopted here. Since the neutrals are coupled to the ions in this approximation, the kinematic behavior of the plasma bears a considerable degree of similarity to the fully ionized case, but radiative losses and ambipolar diffusion can nevertheless influence the dynamics strongly.

Additionally, we approximate the gas pressure tensor as a scalar, assume isothermal conditions (for simplicity and in the absence of detailed information to the opposite), T_i = T_n = T_e, and neglect heat conduction. The latter is based on the fact that heat conduction is energetically less important than radiative cooling for the high densities and low temperatures characteristic of the middle chromosphere (Lin et al. 1992). Additionally, heat conduction into and out of plasmoids is largely across the magnetic field, especially in a two-dimensional description; thus, it tends to be small (see Section 4 for a detailed discussion). This results in the set of basic equations (see, e.g., Khomenko & Collados 2012):

$$\begin{eqnarray} \partial _t \rho &=& -\nabla \cdot (\rho \boldsymbol {v}) \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} \partial _t \boldsymbol {B} &=& \nabla \times (\boldsymbol {v} \times \boldsymbol {B}-\eta \nabla \times \boldsymbol {B}+\boldsymbol {E}_\mathrm{AD}) \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \partial _t e &=& - \nabla \cdot \left[ \left(e+p+\frac{1}{2\mu _0 }\vert \boldsymbol {B} \vert ^2\right)\boldsymbol {v} \right] \nonumber \\ & &+\nabla \cdot \left[\frac{1}{\mu _0} (\boldsymbol {v} \cdot \boldsymbol {B})\boldsymbol {B}\right] \nonumber \\ & & + \nabla \cdot \left[ \frac{\eta }{\mu _0} \textbf {B} \times (\nabla \times \boldsymbol {B}) \right] \nonumber \\ & &-\nabla \cdot \left[\frac{1}{\mu _0}\boldsymbol {B} \times \boldsymbol {E}_\mathrm{AD}\right] \nonumber \\ & &+\rho \boldsymbol {g} \cdot \boldsymbol {v}+\mathcal {L}_\mathrm{rad}+\mathcal {H} \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} e &=& \frac{p}{\Gamma _0-1}+\frac{1}{2}\rho \vert \boldsymbol {v} \vert ^2+\frac{1}{2\mu _0}\vert \boldsymbol {B} \vert ^2 \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} p &=& \frac{(1+Y_\mathrm{i}) \rho }{m_\mathrm{i}} k_\mathrm{B}T. \end{eqnarray} \tag{ 6 }$$

As usual, ρ is the plasma mass density, $\boldsymbol {v}$ is the center-of-mass velocity, e is the total energy density, $\boldsymbol {B}$ is the magnetic field, η is the magnetic diffusivity, p is the plasma thermal pressure, and $\boldsymbol {g}=-273.9 \,{\rm m\,s}^{-2} \boldsymbol {e}_y$ is the gravitational acceleration of the Sun. The ambipolar electric field is given by $\boldsymbol {E}_\mathrm{AD} = \mu _0^{-1}\eta _\mathrm{AD}[(\nabla \times \boldsymbol {B}) \times \boldsymbol {B}]\times \boldsymbol {B}$, where η_AD is the ambipolar diffusion coefficient. $\mathcal {L}_\mathrm{rad}$ is the radiative cooling function, and $\mathcal {H}$ is the heating function. Y_i = n_e/n_H is the ionization degree of the plasma. In the above equations, the following definitions are used:

$$\begin{eqnarray} \rho &=& \rho _\mathrm{n}+\rho _\mathrm{i}+\rho _\mathrm{e} \simeq \rho _\mathrm{H}, \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} m_\mathrm{i} &\simeq & m_\mathrm{n}, \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} \boldsymbol {v} &=& \frac{\rho _\mathrm{n} \boldsymbol {v}_\mathrm{n} +\rho _\mathrm{i} \boldsymbol {v}_\mathrm{i} + \rho _\mathrm{e} \boldsymbol {v}_\mathrm{e}}{\rho }, \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} p &=& n_\mathrm{n}k_\mathrm{B}T_\mathrm{n}+n_\mathrm{i}k_\mathrm{B}T_\mathrm{i}+n_\mathrm{e}k_\mathrm{B}T_\mathrm{e} \nonumber \\ &=& (n_\mathrm{n}+n_\mathrm{i}+n_\mathrm{e})k_\mathrm{B}T \nonumber \\ &=& (1+Y_\mathrm{i})n_\mathrm{H}k_\mathrm{B}T \simeq \frac{(1+Y_\mathrm{i}) \rho }{m_\mathrm{i}} k_\mathrm{B}T. \end{eqnarray} \tag{ 10 }$$

The subscripts n, i, e, and H refer to neutral hydrogen atoms, ions, electrons, and the total number of hydrogen particles (neutral and ionized), respectively, such that the number densities are related by n_H = n_n + n_i = n_n + n_e.

In computing the radiative loss, we follow Gan & Fang (1990) and set

$$\begin{eqnarray} \mathcal {L}_\mathrm{rad} &=& -1.547\times 10^{-42}\, Y_\mathrm{i0} \left(\frac{\rho }{m_\mathrm{i}}\right)^2 \alpha T^{1.5}, \end{eqnarray} \tag{ 11 }$$

where the subscript 0 represents the initial value at time t = 0. The coefficient α in the above function is the same as the one in Equation (11) in Gan & Fang (1990). Their model is applicable up to T ∼ 10⁵ K and is well supported by recent investigations of the active solar chromosphere (e.g., Fang et al. 2003; Jiang et al. 2010; Xu et al. 2011). We assume the plasma to be in ionization equilibrium. The cooling processes considered in this model (Raymond et al. 1976) are permitted, forbidden, and semiforbidden line transitions, including contributions from dielectric recombination and bremsstrahlung, radiative recombination, and two-photon continua. The heating function is chosen to initially balance the radiative loss exactly,

$$\begin{eqnarray} \mathcal {H}&=& 1.547\times 10^{-42}\, Y_\mathrm{i0} \left(\frac{\rho }{m_\mathrm{i}}\right)^2 \alpha T_0^{1.5}. \end{eqnarray} \tag{ 12 }$$

The simulation domain extends from x = 0 to x = L₀ in the x-direction and from y = 0 to y = 2L₀ in the y-direction, with L₀ = 10⁶ m. Open boundary conditions are used in both x- and y-directions.

Three models of the initial current sheet are simulated in this work, differing in the guide field strength and in the inclusion of a density stratification and gravity. For each model we consider three cases, Cases A–C, D–F, and G–I in Models I, II, and II, respectively. In Cases A, D, and G radiative cooling, heating, and ambipolar diffusion are excluded. Cases B, E, and H include heating and radiative cooling but not ambipolar diffusion. Cases C, F, and I include ambipolar diffusion but not the heating and radiative cooling terms. The other initial and boundary conditions and the physical parameters are identical within the cases for each model.

The computations are performed using the MHD simulation code NIRVANA (ver. 3.6; Ziegler 2011). This code does not yet have the capability to describe ionization and recombination, so we work here with a non-time-dependent ionization ratio Y_i0 and defer the inclusion of these effects to a follow-up investigation.

Adaptive mesh refinement (AMR) is applied. We start the simulation from a base-level grid of 160 × 320. The highest refinement level is 10, which corresponds to a grid resolution Δx ≈ 6.1 m. This resolution is comparable to the neutral–ion collisional mean fee path λ_ni in the middle solar chromosphere (see Table 1). Convergence studies have been carried out by repeating some simulations with both a lower and a higher resolution, with the highest refinement level limited to 9 and 11, respectively. The numerical results in the lower-resolution case and in the higher-resolution case are very similar.

2.1. Models I and II

Since the effects of radiative cooling and ambipolar diffusion on the reconnection rate depend strongly on the guide field (Uzdensky & McKinney 2011; Heitsch & Zweibel 2003b), we consider the Harris sheet with strong guide field (Model I) and without a guide field (Model II). These are realized by employing two versions of the Harris current sheet equilibrium, both neglecting gravity and oriented vertically. A force-free Harris current sheet, which is appropriate in a low-β environment and has a strong guide field by its nature, is used as Model I,

$$\begin{eqnarray} B_{x0}&=&0 \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} B_{y0}&=&b_0\tanh [(x-L_0/2)/0.05L_0] \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} B_{z0}&=&b_0/\cosh [(x-L_0/2)/0.05L_0], \end{eqnarray} \tag{ 15 }$$

where b₀ = 0.01 T. The initial current sheet width thus is δ₀ = 0.1L₀. Owing to the force-freeness and neglect of gravity, the initial equilibrium thermal pressure is uniform. The initial plasma β is also uniform and set to β₀ = 0.1. This yields the initial plasma thermal pressure p₀ = 12.5π⁻¹ Pa.

To realize a model without guide field, we must deviate from force-freeness. Model II uses a standard Harris current sheet in equilibrium with a pressure gradient,

$$\begin{eqnarray} B_{x0}&=& 0 \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} B_{y0}&=& b_0\tanh [(x-L_0/2)/0.05L_0] \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} B_{z0}&=& 0\,. \end{eqnarray} \tag{ 18 }$$

From Equation (2), the initial plasma pressure in Model II is

$$\begin{eqnarray} p_0 &=& (1+\beta _0-(\tanh [(x-L_0/2)/0.05L_0])^2)\frac{b_0^2}{2\mu _0}, \end{eqnarray} \tag{ 19 }$$

and we set the initial asymptotic plasma β also to β₀(|x| → ∞) = 0.1.

In all cases we assume the initial equilibrium temperature as

$$\begin{eqnarray} T_0 &=& T_{0b}\left[1+(1+\tanh (12(y-L_0)/L_0))/2\right], \end{eqnarray} \tag{ 20 }$$

where T_0b = 3500 K. The functions p₀(x) and T₀(y) yield the total mass density ρ₀(x, y) and particle number density n_H0 + n_e0 for each case. Then we calculate the initial ionization degree Y_i0 = n_e0/n_H0 using Equation (12) in Gan & Fang (1990). Our settings yield rather similar values of the initial gas pressure, mass density, field strength, and Alfvén velocity v_A0 = b₀(μ₀ρ_H)^−1/2 in the inflow region (x = 0 and L₀, and y = L₀) for both models. However, in Model II the gas pressure and mass density in the current sheet region are higher by a factor of 11 compared to the inflow region. The resulting Alfvén velocity in the current sheet is lower than in Model I by a factor of ≈3. We should point out that the Alfvén velocity calculated in this work is based on the total plasma density, as is appropriate for strongly coupled partially ionized plasmas (at scales δ ≫ λ_ni); see, e.g., Zweibel (1989).

The basic plasma parameters in the current sheet region are compiled for the two models in Table 1 (along with the values of the resulting Lundquist number based on L₀, Alfvén velocity, and length scales); they correspond to the conditions in the middle chromosphere, i.e., at roughly 1 Mm above the photosphere, with the field strength being representative of active regions. The neutral–ion collisional mean free path is given by the ratio of the thermal velocity of the neutrals and the neutral–ion collision frequency, $\lambda _\mathrm{ni}=v_{T_\mathrm{n}}/\nu _\mathrm{ni}$, with $v_{T_\mathrm{n}}=\sqrt{2k_\mathrm{B}T_\mathrm{n}/m_\mathrm{n}}$, $\nu _\mathrm{kl}=n_\mathrm{l}\sigma _\mathrm{kl}\sqrt{8k_\mathrm{B}T_\mathrm{kl}/(\pi m_\mathrm{kl})}$, T_kl = (T_k + T_l)/2, m_kl = m_km_l/(m_k + m_l) (where the subscripts k and l denote the species), and the value of the collisional cross section σ_kl is taken from Khomenko & Collados (2012).

Initial perturbations of both magnetic field ($\boldsymbol {B}^{\prime }$, with $\nabla \cdot \textbf {B}^{\prime }=0$) and velocity ($\boldsymbol {v}^{\prime }$) are applied at t = 0 in all cases to trigger the reconnection process. The perturbations result in a thinning of the current sheet in two sections between a set of three primary islands, whose midpoints are located at y = 0, L₀, and 2L₀ (see Figure 1(a)). In this paper we focus only on the section in the domain L₀ < y < 2L₀, i.e., the bottom half of the box is used as an auxiliary part of the computation only. Its function is to yield a stationary primary plasmoid at the bottom of the height range of interest, which is not influenced by any effects of a numerical boundary; the primary plasmoid acts like a line-tied bottom of the current sheet of interest. The initial temperature T₀ varies from 5250 to 7000 K in the domain of interest, but it quickly increases to the upper value above the stationary plasmoid at y = L₀, such that essentially all of the upper current sheet initially is at T₀ ≈ 7000 K. The resulting ionization degree in the current sheet is 0.5% in Model I and 0.2% in Model II.

Zoom In Zoom Out Reset image size

We use a simple parameterization of the magnetic diffusion that relatively closely matches the diffusion computed from a model atmosphere in Khomenko & Collados (2012), η = 5 × 10⁴(3500/T)^1.5 m² s⁻¹. This yields the initial diffusion in the current sheet as 1.77 × 10⁴ m² s⁻¹. The Lundquist number based on this diffusivity, v_A0, and L₀, which corresponds to the global scale of the current sheet in the upper part of the box, is S₀ = 1.88 × 10⁶ in Model I and 5.66 × 10⁵ in Model II. η decreases as the temperature increases. From the numerical results presented in the following section, the value of η at the main X-point drops to ∼2000 m² s⁻¹ during the secondary instability process, and the Lundquist number then reaches ∼10⁷.

The ambipolar diffusion coefficient η_AD is defined as (Braginskii 1965)

$$\begin{eqnarray} \eta _\mathrm{AD}&=& \left(\frac{\rho _\mathrm{n}}{\rho }\right)^2(\rho _\mathrm{i}\nu _\mathrm{in}+\rho _\mathrm{e}\nu _\mathrm{en})^{-1}, \end{eqnarray} \tag{ 21 }$$

where ν_in and ν_en are, respectively, the ion–neutral and electron–neutral collision frequencies defined analogous to ν_ni above, with the cross sections σ_in = 5 × 10⁻¹⁹ m² and σ_en = 10⁻¹⁹ m² again taken from Khomenko & Collados (2012). This yields

$$\begin{eqnarray} \eta _\mathrm{AD}&=& 1.65 \times 10^{-11} \left(\frac{1}{Y_i}-1 \right) \frac{1}{\rho ^2\sqrt{T}}\, {\rm m^3\,s\, kg^{-1}}, \end{eqnarray} \tag{ 22 }$$

with the temperature inserted in kelvin. The ambipolar diffusion coefficient is a function of ionization degree, temperature, and plasma density. According to the model atmospheres by Khomenko et al. (2008) and Pneuman et al. (1986), the magnitude of η_AD varies from 10 to 10⁵ m³ s kg⁻¹ from the bottom to the top of the chromosphere. Here we set η_AD = 5 × 10⁴(T₀/T)^1/2(ρ₀/ρ)² m³ s kg⁻¹. This is representative of the ambipolar diffusion coefficient in the middle chromosphere.

For the radiative loss and heating function in Models I and II, we simplify the coefficient α compared to the expression in Gan & Fang (1990), which is slowly varying in the middle chromosphere (about 1000 km above the photosphere), the height range of interest in the present study. Since we do not investigate the height dependence of the reconnection rate, we here use a representative uniform value of α = 0.01, which is chosen such that the resulting magnitude of the initial radiative loss is of an order of 0.07–3.5 J m⁻³ s⁻¹, consistent with their results.

2.2. Model III

Since the active section of the current sheet in our computations, L₀ ≲ y < 2L₀, is of order 1 Mm, it encompasses a considerable range of density and, correspondingly, Alfvén velocity in the chromosphere. In Model III we therefore include gravity and the resulting stratification of thermal pressure and plasma density. It is not possible to construct an equilibrium consisting of a vertical Harris sheet without a guide field in the presence of gravity. Therefore, Model III combines the force-free Harris sheet of Model I with density stratification, gravity, and a balancing pressure gradient, given by

$$\begin{eqnarray} \partial _x p_0(y) &=& 0, \end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray} \partial _y p_0(y) &=& -273.9\rho _0(y), \end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray} p_0(y) &=& \frac{(1+Y_\mathrm{i0}) \rho _\mathrm{H0}(y)}{m_\mathrm{i}} k_\mathrm{B}T_0, \end{eqnarray} \tag{ 25 }$$

where, for simplicity, the initial temperature T₀ = 7000 K is chosen to be uniform, and ρ_H0 = n_H0m_i. According to Equation (12) in Gan & Fang (1990), the ionization degree is simplified as

$$\begin{eqnarray} &&Y_\mathrm{i0}=\frac{n_\mathrm{e0}}{n_\mathrm{H0}} \approx \sqrt{\frac{\phi _0}{n_\mathrm{H0}}}, \end{eqnarray} \tag{ 26 }$$

where ϕ₀ ≃ 1.21 × 10¹⁵ m⁻³ and ϕ₀ ≪ n_H0 in our model. From Equations (22)–(25), we obtain

$$\begin{eqnarray} &&n_\mathrm{H0} \simeq n_\mathrm{H00}\exp \left(-\frac{4.706y}{L_0}\right), \end{eqnarray} \tag{ 27 }$$

where we assume n_H00 = 10²³ m⁻³. The same perturbations as those in Models I and II are applied; we also only focus on the section in the domain L₀ < y < 2L₀. From Equations (24)–(26), we calculate the initial distributions of hydrogen number density n_H0, ion number density n_i0, ionization degree Y_i0, Alfvén velocity v_A0, ion inertial length d_i0, and neutral–ion collision mean free path λ_ni0. The values of these parameters in Model III at y = L₀, y = 1.5L₀ and y = 2L₀ are compiled in Table 1.

The initial Lundquist number S₀ = Lv_A0/η is of order 10⁶ (see Table 1). The initial ambipolar diffusion coefficient takes values of η_AD = 75.67 m³ s kg⁻¹ at y = L₀, 2.57 × 10³ m³ s kg⁻¹ at y = 1.5L₀, and 8.71 × 10⁴ m³ s kg⁻¹ at y = 2.0L₀.

Since the density stratification in the y-direction is included in Model III, we here use the full y-dependent expression for α from Gan & Fang (1990), so that the resulting height dependence of the radiative loss is similar to their results. Compared to their expression, we here increase the value of α by a factor of 10, such that the initial radiative loss is in the range ≈(0.2–1.3) J m⁻³ s⁻¹ for L₀ < y < 2L₀, considerably higher than the value of ≈0.07 J m⁻³ s⁻¹ in Model I. This choice is made to demonstrate that a main result we will obtain for Model I, namely, that radiative losses in the chromosphere do not significantly influence the evolution of current sheets with a strong guide field, remains valid for the much higher level of radiative losses. The levels chosen for the initial radiative losses in Models I and III both fall in the range deduced in Gan & Fang (1990).

From the above descriptions, we find that Model III is close to the conditions in the middle chromosphere, i.e., in the height range ≈(0.5–1.4) Mm above the photosphere. Figure 1(b) displays the initial magnetic field and plasma density.

3.1. Cascading Plasmoid Instability with Internal Slow-mode Shocks

The Lundquist number in our models exceeds the critical Lundquist number for onset of the plasmoid instability substantially. As the Harris sheet evolves to a thinner current sheet in response to the initial perturbation, it first passes through a Sweet–Parker-like phase (with the exception of Case F), as indicated by its smooth, very elongated structure (which contains an X-point at y  ≈  1.5L₀) and an average reconnection rate γ ∼ 10⁻³ in agreement with the Sweet–Parker scaling (see below for the method of quantitative analysis). When the current sheet sufficiently thins such that a critical aspect ratio is reached, the plasmoid instability sets in, forming multiple islands. Subsequently, the instability cascades to smaller scales. Throughout this process a high reconnection rate, rapidly fluctuating in the range ∼0.01–0.03, is realized (see Section 3.2). The plasmoid cascading process is found to be very similar in all nine cases run.

Figure 2 displays the plasmoid cascading process for Model III, Case I with ambipolar diffusion included. The left panel shows the primary island at y = L₀ produced by the initial perturbation and a number of secondary islands formed by the plasmoid instability (note that the primary island at y = 2L₀ was ejected by the slow reconnection outflows that evolved in the Sweet–Parker-like phase). As we zoom in (panels 2–3), several third-order plasmoids become visible in the current sheet, which has further narrowed between the secondary plasmoids. Zooming in further (panels 3–4), fourth-order plasmoids can be seen clearly. This plasmoid cascading process is very similar to the results of Bárta et al. (2011), who have simulated it for coronal parameters. It is found here for the first time in partially ionized plasma.

Zoom In Zoom Out Reset image size

The half-width of the thinnest fragment current sheet in Figure 2 is about 25 m, clearly larger than the neutral–ion collisional mean free path λ_ni (Table 1), thus consistent with the one-fluid approximation. The cascading plasmoid instability occurs in all nine runs considered in this paper, reaching comparable reconnection rates in the cases with and without a guide field, which differs strongly from the results for stationary reconnection in partially ionized plasma (Heitsch & Zweibel 2003b; Uzdensky & McKinney 2011). The gravity of the Sun also has only little effect on the time-dependent reconnection rate (see below).

Given the large body of experience from one-fluid, Hall-MHD and kinetic models of reconnection in the plasmoid-unstable range (e.g., Huang et al. 2011; Markidis et al. 2013; Daughton et al. 2009), one can reasonably expect that the plasmoid cascading process tends to continue to even smaller scales than resolved here, down to the ion inertial length d_i. In partially ionized plasma the process is modified by recombination within the plasmoids, reducing their current and flux (Leake et al. 2012), and it is not yet known how this influences the cascading process. A less dynamical evolution of the plasmoids is conceivable, but since the cascading occurs between previously formed plasmoids, it appears unlikely that recombination within plasmoids would suppress it completely. This is indeed indicated by the results in Leake et al. (2013), who find that the plasmoid instability can raise the reconnection rate also beyond the validity range of the one-fluid approximation. We thus expect that the multiple levels of the plasmoid instability process connect the global, large MHD-scale reconnection with the local, small kinetic-scale reconnection also in partially ionized plasma. It is well known that the plasmoid instability in fully ionized plasma yields a high reconnection rate throughout the cascading process, and our simulations demonstrate that the cascading plasmoid instability operates in weakly ionized plasma in the one-fluid limit in basically the same manner as in fully ionized plasma, reaching a high reconnection rate as well.

As is well known, reconnection at X-points always involves slow-mode shocks in the Petschek regime. In our simulations, we find that many small-scale slow-mode shocks form transiently at the edge of secondary plasmoids, attached to the neighboring X-point in a secondary fragment of the current sheet. Figure 3 shows a pair of such shocks. Sudden changes of the magnetic field direction and current density at the left and right edges of the plasmoid can clearly be seen; these are very similar to the slow-mode shock structure displayed, e.g., in Forbes (1988), Schumacher & Kliem (1997), and Mei et al. (2012). Figure 3(b) shows that the temperature increases strongly not only within the current sheet fragment but also downstream of the slow-mode shocks in the plasmoid. The angle between the shock front and the y-direction is ≈54. Figure 4 displays the magnetic field and the current density along a cut in the x-direction at y = 1.691L₀. One can see that the magnitude of the field component tangential to the shock, B_∥, decreases rapidly toward the downstream side and that the current density J_z has a peak where B_∥ changes rapidly, but the component normal to the shock, B_⊥, stays nearly uniform along the cut. This is exactly the behavior of a slow-mode shock. Similar structures exist in the outflow region of the current sheet. These slow-mode shocks could be a candidate mechanism to explain the chromospheric heating (Kumar et al. 2011).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2. Effects of Radiative Cooling and Ambipolar Diffusion

3.2.1. Strong guide field

We now compare the simulation results of the three cases that were run for Model I. The reconnection rate is computed as the rate of change of the magnetic flux accumulated between the O-point in the primary island at y = L₀ and the main reconnection X-point (see Ni et al. 2012b, 2013),

$$\begin{equation} \gamma (t)=\frac{\partial (\psi _X(t)-\psi _O(t))}{\partial t}\frac{1}{b_0 v_{A0}}. \end{equation} \tag{ 28 }$$

Here the flux function ψ is defined through the relations B_x = −∂ψ/∂y, B_y = ∂ψ/∂x, and the main reconnection X-point is determined as the X-point that has the highest ψ value of all X-points in the box. This analysis is performed on the relatively low grid refinement level 3, which ensures uniform coverage of the current sheet section from the primary O-point to the upper boundary throughout the runs. This choice still captures much of the temporal variability of the reconnection rate. Local quantities like the current density, temperature, and current sheet width at the main reconnection X-point are determined at the highest refinement level chosen by the code, to evaluate their full dynamic range; after the onset of the plasmoid instability this is level 10 in all nine cases. The current sheet width is determined as the FWHM of the current density profile in the x-direction.

The reconnection rate for Cases A–C is displayed in Figure 5(a). The three cases are very similar to each other and show a dynamical behavior that is basically analogous to reconnection in fully ionized plasma, as expected from our one-fluid approximation, especially for Case A in the absence of heating, radiative cooling, and ambipolar diffusion. Up to t ≈ 26 s the reconnection rate remains at a low level, reaching the Sweet–Parker value of γ  ∼  S^−1/2 ∼ 10⁻³ for t ≳ 20 s. At this time the decreasing current sheet width has also reached the Sweet–Parker value δx  ∼  S^−1/2L₀ ∼ 10³ m. This appears to be a Sweet–Parker reconnection phase, which is also supported by the elongated, smooth appearance of the current sheet (see Figure 1).

Zoom In Zoom Out Reset image size

Subsequently, a fast rise of the reconnection rate commences, which quickly saturates at a level γ ≈ 0.02. The first secondary islands due to the plasmoid instability have clearly developed by the end of the fast rise, which we thus interpret as the linear phase of the plasmoid instability. All quantities show the high variability characteristic of the plasmoid instability in the further evolution. The current sheet width, measured at the main X-point, gradually decreases in the Sweet–Parker phase and enters a rapid decrease down to δx ≈ 30 m simultaneously with the onset of the plasmoid instability (Figure 6(a)). This minimum current sheet width is comparable to the grid scale at refinement level 10 (δx ≈ 5Δx) but still much larger than the neutral–ion collisional mean free path and the ion inertial length (Table 1), so that our one-fluid approximation is valid throughout the computation. The current density at the main X-point shows a rather close inverse proportionality to the current sheet width. The average reconnection rate remains at a high level throughout the cascading phase of the instability. It shows a trend of gradual decrease to an average level of ≈0.015, which is most likely due to the shortening of the active section of the current sheet as the primary island at y = L₀ grows; the field lines in the inflow region must then bend increasingly before reconnecting at the X-point. The temperature at the main X-point increases rapidly from the onset of the plasmoid instability, saturating at ∼3 × 10⁴ K (Figure 7(a)). The corresponding decrease of the magnetic diffusivity raises the effective Lundquist number in the course of the plasmoid instability by nearly an order of magnitude. It should be noted that the peak temperature in the current sheet, located within the islands downstream of the slow-mode shocks, continues to rise, reaching ∼8 × 10⁴ K during the plasmoid instability.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The comparison of Cases A and B shows that the radiative cooling hardly has any effect on the evolution of the reconnecting current sheet at the chosen parameters of the middle chromosphere if the sheet includes a considerable guide field. This is fully in agreement with earlier investigations (e.g., Uzdensky & McKinney 2011), which had pointed out that a sufficiently strong guide field suppresses the current sheet thinning due to radiative cooling by making the sheet incompressible. The only remaining effect on the reconnection rate is given by the temperature dependence of the magnetic diffusivity. Uzdensky & McKinney (2011) considered the case that radiative cooling in a stationary Sweet–Parker sheet with negligible heating increases the reconnection rate in this way. However, heating is a key factor of the energy balance in the chromosphere, and so our model also includes a heating term that is adjusted to match the initial cooling rate and to have the same dependence on the density. Consequently, the radiative losses will not cool the plasma below the initial temperature T₀; they have a strong effect only when the temperature increases considerably above T₀. Lower temperatures can be reached by adiabatic cooling if the plasma expands (an effect seen below in Figure 9), but this is not the case here at the X-point, where the temperature has the same value as in Case A during the Sweet–Parker phase (Figure 7(a)). The subsequent fast reconnection dominated by the plasmoid instability is largely independent of the magnetic diffusivity (Huang & Bhattacharjee 2010). Temperatures of (4–5)T₀ are reached in this phase, as in Case A, which shows that the strong heating at the slow-mode shocks dominates the radiative losses in this phase.

Considering the effect of ambipolar diffusion in Case C, again only a weak influence on the reconnection rate is found. The only differences to Cases A and B are a more impulsive and slightly earlier onset of the plasmoid instability and a somewhat higher variability near the end of the Sweet–Parker phase. We interpret this behavior as follows. In the Sweet–Parker phase a strong thinning of the current sheet is inhibited by the guide field (Heitsch & Zweibel 2003b). The ambipolar electric field nevertheless increases strongly as the sheet thins in response to the initial perturbation. However, the area of enhanced E_AD does not reach the plane x = 0.5L₀, where the antiparallel field component reverses sign, but rather forms two layers on either side of the field reversal plane (Figures 8(a) and (c)). Thus, it does not amplify the annihilation of the field in this plane. The width of the layers, ∼500 m, agrees with the ion–neutral decoupling length scale L_AD = v_Ai/ν_in, where v_Ai is the Alfv$\acute{e}$n velocity based on the ion density (Zweibel et al. 2011). As a dominant X-point develops, the associated reconnection flows convect the two layers with enhanced ambipolar electric field to the field reversal plane (Figures 8(b) and (c)), so that it now contributes to the change of flux (Equation (2)). The comparison of Figures 8(a) and (b) shows that the rate of change of the magnetic field in the Sweet–Parker-like current sheet is dominated by the resistive contribution to the electric field, whereas it is dominated by the ambipolar electric field during the rapid rise of the reconnection rate at the first X-point of the commencing plasmoid instability. The layers of enhanced ambipolar field intermittently approach the field reversal plane already during the considerable thinning of the current sheet in the final part of the Sweet–Parker phase (Figure 6(a)), causing the intermittent moderate enhancements of the reconnection rate seen in Figure 5(a). Since the inclusion of the ambipolar diffusion terms causes a strong decrease of the adaptive time step, this run was terminated after the reconnection rate reached the saturation level.

Zoom In Zoom Out Reset image size

The three cases also show a very similar current sheet aspect ratio at the onset of the plasmoid instability. We refer to this quantity as the "onset aspect ratio" and determine it from the current sheet width δx and the similarly defined FWHM length of the current sheet in the y-direction, δy. We find δx ≈ 400 m and δy ≈ 3.5 × 10⁵ m, i.e., an onset aspect ratio δ_ons = δy/δx ≈ 875 in all three cases (see Figure 9 for Case B). This value lies considerably above the critical aspect ratio for onset of the plasmoid instability at the point of transition between the Sweet–Parker and plasmoid-mediated reconnection regimes, $S_\mathrm{cr}^{1/2}$, when the Lundquist number is varied to reach the critical value S_cr ∼ 10³–10⁴ (Huang & Bhattacharjee 2010; Ni et al. 2012b; Leake et al. 2012). Our simulations are performed with far higher Lundquist numbers in the range S ∼ 10⁶–10⁷ (varying in time as a result of the evolving δx and T), and correspondingly the current sheet aspect ratio rises considerably above the critical value $S_\mathrm{cr}^{1/2}$, roughly approaching S^1/2.

Zoom In Zoom Out Reset image size

In summary, in current sheets with a strong guide field, the plasmoid instability mediates a fast reconnection regime in partially ionized plasma, while radiative cooling and ambipolar diffusion have only very little effect. For the parameters of the middle chromosphere, the reconnection rate is enhanced by a factor of ∼20 above the Sweet–Parker rate, reaching the range of observationally inferred values.

3.2.2. Zero Guide Field

Next, we consider the Harris current sheet with vanishing guide field (Model II, Cases D–F). A major difference to all corresponding runs for Model I is the much slower evolution of the current sheet due to our choice to realize the necessary pressure enhancement in the sheet by a density enhancement and the resulting smaller Alfvén velocity in the sheet. The reconnection rate, which is normalized to the field strength and Alfvén velocity in the inflow region, thus takes a much longer time to rise. The density in the current sheet gradually decreases during our runs, and eventually the Alfvén velocity in the sheet becomes comparable to the external Alfvén velocity. The reconnection rate then reaches similar values to those in Model I.

In all three cases run for Model II, reconnection rates exceeding the Sweet–Parker value develop only from the onset of the plasmoid instability (Figure 9), which happens at t ≈ 92, 84, and 50 s for Cases D–F, respectively (Figures 5(b) and 9). A Sweet–Parker-like regime develops only transiently, for ∼10 s prior to the onset of the plasmoid instability, in Cases D and E as the current sheet width approaches the Sweet–Parker value (Figure 6(b)). As for Model I, the onset of the plasmoid instability is associated with a rapid decrease of the current sheet width. Owing to the high inertia of the plasmoids in the initially dense current sheet, the instability shows a far more gradual overall development compared to the strong guide field cases.

The reconnection rate in Case D without radiative cooling and ambipolar diffusion does not reach saturation but rather rises essentially monotonically (apart form the small timescale fluctuations characteristic of the instability) until the main X-point leaves the box through the upper boundary at t ≈ 230 s. A peak reconnection rate of γ ≈ 0.025 is then reached; it would be somewhat higher (lower) if the box were chosen somewhat larger (smaller). The current sheet width continues to decrease in the course of the cascading plasmoid instability until it reaches δx ≈ 30 m, as in Model I. The temperature at the X-point reached at the end of the run, T ∼ 3 × 10⁴ K, is also similar to the result for Model I.

The evolution in Case E with radiative cooling and heating is qualitatively similar to Case D, but it commences earlier and is more impulsive. Here the current sheet width drops to δx ≈ 30 m, and the reconnection rate rises to γ ≈ 0.025, in about half the time. This more dynamic behavior is related to a faster decrease of the density in the current sheet after t ≈ 100 s. It appears that the reconnection rate saturates at this level until it sharply drops as the main X-point leaves the box at t ≈ 170 s. The effect of the radiative cooling again remains minor prior to the onset of the plasmoid instability. The current sheet is then only slightly thinner than in Case D (Figures 6(b) and 8). The effect of radiative cooling becomes more obvious at the high temperatures in the course of the plasmoid instability, limiting the temperature to about half the peak value found in Case D (Figure 7(b)). Radiative cooling in Case D acts more efficiently than in Case B because here the densities are higher.

As expected in the absence of a guide field, the inclusion of ambipolar diffusion in Case F allows the current sheet to thin very efficiently (Brandenburg & Zweibel 1994). The evolution here begins with a rapid decrease of the current sheet width at the main X-point down to δx ≈ 30 m (within ≈10 s), which is near its smallest value observed in all runs at the AMR refinement levels 10 and 11. A thinning by a factor of 500 is expected from the ratio of ion and neutral pressure (the inverse ionization degree for our isothermal conditions) when only the ion pressure in the sheet balances the magnetic pressure of the external region. The actual thinning goes beyond that by a factor of ∼3, but this already includes the onset of the plasmoid instability at t ≈ 50 s and δx ≈ 80 m. The instability seamlessly grows out of the current sheet collapse mediated by the ambipolar diffusion, without leaving room for a Sweet–Parker-like phase, since the onset aspect ratio is obviously reached during the collapse. The reconnection rate grows in a manner similar to Cases D and E, at least to a level γ ∼ 0.02. It is possible that it would grow to higher values if the run could be continued. However, the main reconnection X-point has moved out of our simulation domain after t = 150 s. It is seen that the reconnection rate then suddenly decreases to a low value, similar to Cases E and D. Overall, the dynamical evolution of the plasmoid instability appears to be similar to Cases D and E, including the temperatures reached at the main X-point (Figure 7(b)), so that a saturation of the reconnection rate at similar levels is conceivable.

The current sheet widths and aspect ratios for onset of the plasmoid instability are, approximately, 1200 m and 500 in Case D, 950 m and 730 in Case E, and 80 m and 1250 in Case F (see Figure 9). The onset occurs at a higher aspect ratio when the sheet thins more dynamically.

To summarize the runs for Model II without guide field, we find a basic analogy to the strong guide field cases in that the plasmoid instability is the dominant mechanism to realize fast reconnection, with rates similar to the strong guide field cases. Although our simulations confirm the strong thinning of the current sheet allowed by ambipolar diffusion and anticipated to yield an accelerated Sweet–Parker-like reconnection regime (Heitsch & Zweibel 2003a), this regime does not occur because the aspect ratio for onset of the plasmoid instability is reached by the end of the ambipolar diffusion-driven current sheet thinning. An accelerated Sweet–Parker-like reconnection due to radiative cooling of the current sheet does not occur in our model because we include a background heating term representing the chromospheric heating.

3.3. Effects of the Density Stratification

We find that the inclusion of gravity and density stratification in Model III has surprisingly little effect on the main results found for Models I and II, although the density varies by two orders of magnitude in the relevant height range y > L₀. Specifically, (1) the plasmoid cascading process shown in Figure 2 is very similar to Models I and II; (2) the reconnection rate, current sheet width, and temperature at the main X-point for Cases G–I in Model III behave very similar to those for Cases A–C in Model I, respectively (Figures 5–7); (3) many fine slow-mode shock structures attached around the edges of the plasmoids appear, where the highest temperatures are produced; (4) the onset aspect ratios for Models I and III are the same, as shown in Figure 9 and Table 1; (5) radiative cooling and ambipolar diffusion have little effect on the plasmoid instability in Model III; (6) the maximum outflow velocity in Model III is the same as in Model I, about 40 km s⁻¹. The similarity of the reconnection rate, which is normalized by the same value of b₀v_A0 as in Model I, to other quantities at the main X-point can be understood from the rather similar values of the Alfvén velocity at that point in Models I and III, with the Alfvén velocity initially being higher in Model III by only a factor of ≈1.2; this implies a similar dynamic behavior in the vicinity of the main X-point.

There are some details in the reconnection process that are different in Models I and III. As shown in Figures 9 and 10, the current sheet gradually becomes inclined to the side, especially after secondary islands appear. This is due to gravity acting on the dense plasmoids, which form slightly asymmetrically about the current sheet as a consequence of round-off errors. Secondary islands in Model III start to appear later than in Model I for the same initial perturbation. In the period before secondary islands appear, the hottest plasma exists in the upflow region in Model III (Figure 11) but in the downflow region in Model I (not shown). Since low-density plasma can be heated more easily, the rapid decrease of the plasma density with increasing height in Model III is the main reason for this difference. For Case H, the hottest plasma usually appears around the edge of the plasmoids, but their center is still relatively much cooler (Figure 11), the reason being that the plasma density inside the plasmoids is higher and, therefore, radiative cooling is stronger. During the reconnection process, the main X-point in Model III gradually moves to a higher place compared to Model I, as seen in Figures 10 and 12. With radiative cooling, this results in a higher temperature for Model III (∼80, 000 K in Case H versus ∼70, 000 K in Case B).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.4. Comparison with Observations

Magnetic reconnection in the middle chromosphere, including radiative cooling, heating, and ambipolar diffusion effects, has been studied in the 2.5 dimensional MHD approximation. A Harris current sheet model with either a strong guide field (Models I and III) or a vanishing guide field (Model II) is considered. We have assumed that the plasma consisting of ions, electrons, and neutral hydrogen atoms is strongly coupled and in ionization equilibrium state, that the pressure is isotropic, and that the heat flux can be neglected. The ionization and radiative loss model by Gan & Fang (1990) is used, resulting in an initial ionization degree in the current sheet in the range ≈(0.2–0.5)%. A temperature-dependent magnetic diffusion coefficient adapted to the values in Khomenko & Collados (2012) is employed, leading to Lundquist numbers of S ∼ 10⁶–10⁷.

The main numerical results can be summarized as follows.

• 1.  
  The plasmoid instability mediates a fast reconnection regime under chromospheric conditions for a vanishing and a strong guide field. Reconnection rates of order 0.01–0.03, a factor of ∼20 above the Sweet–Parker rate and within the range of observationally inferred rates, are reached in either case, based on the classical Spitzer resistivity. The aspect ratio at the onset of the plasmoid instability at the supercritical values of the Lundquist number considered here is found to lie in the range 500–1250.
• 2.  
  AMR allowed us to resolve three levels of secondary islands and current sheet widths down to δx ≈ 30 m, where the one-fluid approximation is still valid.
• 3.  
  Ambipolar diffusion and radiative cooling have a significant influence on the reconnection process only in the model with a vanishing guide field. Ambipolar diffusion then leads to the expected strong current sheet thinning (Brandenburg & Zweibel 1994), but the current sheet thins very fast, reaching the onset aspect ratio, such that an accelerated Sweet–Parker reconnection regime does not develop, but rather the plasmoid instability sets in. The expected strong current sheet thinning in a Sweet–Parker regime due to radiative cooling (Uzdensky & McKinney 2011) does not occur in our simulations because we include a background heating term. This thinning can generally not be expected to be strong in the middle and lower chromosphere, where the temperature is hardly higher than twice the solar temperature minimum. Both ambipolar diffusion and radiative losses have an effect on the plasmoid instability for a vanishing guide field, implying, respectively, a much earlier onset and a faster development, but apparently no significant increase in the reconnection rate.
• 4.  
  Many slow-mode shocks are generated between secondary plasmoids and secondary fragments of the current sheet during the unstable reconnection process. Downstream from the shock front regions (in the plasmoids), the temperature is highest and reaches T  ∼  8 × 10⁴ K, i.e., lower transition region temperatures. This supports conjectures in the literature that fast magnetic reconnection is a candidate mechanism for the chromospheric heating.
• 5.  
  Upward outflow velocities in the course of the plasmoid instability reach ≈40 km s⁻¹ and ≈25 km s⁻¹ in the strong and zero guide field cases, respectively. Outflow velocities higher by a factor of ∼3 are expected in the zero guide field case if the initial current sheet was not chosen to be overdense. These values lie in the range observed for chromospheric jets. Dynamic plasmoids yield variations in the outflow, consistent with observations of plasmoids in chromospheric jets.

Recent work by Leake et al. (2012, 2013) has studied the effects of recombination on reconnection in the middle chromosphere in a two-fluid model (weak coupling between ions and neutrals). This also included the occurrence of the plasmoid instability. Reconnection rates of an order of 0.1 were found to result from recombination effects, while the plasmoid instability gave only a minor additional increase of ∼15%. Recombination within the plasmoids reduces their current, which may slow down the dynamics of the plasmoid instability. It thus appears that recombination must ultimately be included for a complete description of reconnection in partially ionized plasma. However, the small increase of the reconnection rate due to the plasmoid instability in Leake et al. (2012) results in part from their choice of a high magnetic diffusivity, η = 6 × 10⁴–2.4 × 10⁶ m² s⁻¹, which is about two orders of magnitude higher than a typical magnetic diffusivity in the middle chromosphere and did not change with temperature. Our diffusivity is η₀ = 1.8 × 10⁴ m² s⁻¹ initially and decreases to ∼2 × 10³ m² s⁻¹ as the current sheet is heated in the course of the plasmoid instability. This gives a difference in η of up to three orders of magnitude. In Leake et al. (2013) the diffusivity was computed for chromospheric temperatures and of order 3 × 10³ m² s⁻¹ before the onset of the plasmoid instability, much closer to our values. The reconnection rate then reached ≈0.05, only a factor of ∼2 above our value. Thus, the high reconnection rate found in Leake et al. (2012) results in part from the high diffusivity chosen in this paper. The reconnection rate of ≈0.05 in Leake et al. (2013) is due to recombination effects and higher than but comparable to the reconnection rate we find during the plasmoid instability. After the onset of the plasmoid instability in Leake et al. (2013), the reconnection rate shows a steep increase, which could only be followed up to ≈0.06, since the current sheet width then reached the resolution limit. The steep increase can be taken as an indication that the plasmoid instability is potentially very important in the weak coupling description of chromospheric reconnection as well. Additionally, we note that Leake et al. used the magnetic field and Alfvén velocity in the upstream region not very far from the current sheet to normalize the reconnection rate, while we used the asymptotic inflow values. Our reconnection rates increase by a factor of 2–3 if we use the magnetic field and Alfvén velocity at distances in the range of 5–20 current sheet half-widths ((2.5–10)δx). Since the initial plasma density at the initial main reconnection X-point (y ≈ 1.5L₀) in our models (∼10²⁰ m⁻³ in Models I and III) is more than 10 times higher than that in their models (6 × 10¹⁸ m⁻³), and the initial temperature (T₀ = 7000 K) is lower than in their models (T₀ = 8500 K), our models represent the conditions at somewhat smaller heights in the chromosphere.

Our simulations do not include thermal conduction, which can generally be expected to be energetically less important than radiative cooling at chromospheric temperatures of ∼10⁴ K, which dominate in the major part of our simulated volumes (Lin et al. 1992). Much higher temperatures and steep temperature gradients are built up in the plasmoids and at the slow-mode shocks and X-points. Conductive energy transport from these regions into the inflow regions and into the big plasmoid at the bottom end of the current sheet is largely directed across the magnetic field (conduction in the third direction, parallel to the field, is not relevant in our 2.5 dimensional simulations). Thermal conduction by the charged particles is very strongly suppressed in the direction perpendicular to the field; however, the contribution by the neutral component is not influenced. Deviations from local thermodynamic equilibrium (LTE) allow a neutral component to be present even at the high temperatures that we find in the plasmoids. To estimate the potential magnitude of thermal conduction for the conditions of our simulations, we follow Orrall & Zirker (1961). They find that the thermal conductivity across the field is provided by the neutral component in the whole temperature range relevant here and give the coefficient of thermal conductivity in a pure hydrogen gas as

$$\begin{equation} \kappa = \frac{75}{64}k_\mathrm{B} \left(\frac{k_\mathrm{B}T}{\pi m_\mathrm{i} M_\mathrm{H}} \right)^{1/2}\frac{n_\mathrm{n}}{\sigma _\mathrm{nn}(n_\mathrm{n}+n_\mathrm{i})}, \end{equation} \tag{ 29 }$$

where M_H = 1 is the molecular weight and σ_nn = 3.21 × 10⁻²⁰ m² is the neutral–neutral collisional cross section (note that Orrall & Zirker use $\sigma _\mathrm{nn}^2$ to designate the cross section). Since the ionization degree is small and does not change with time, n_n/(n_n + n_i) ≃ 1 in our models. Under these conditions, the coefficient of thermal conductivity is evaluated as

$$\begin{equation} \kappa \simeq 2.59\times 10^{-2} T^{1/2}, \end{equation} \tag{ 30 }$$

where T is measured in kelvin and the dimension of κ is J K⁻¹ m⁻¹ s⁻¹. The thermal conduction is given by

$$\begin{equation} F_C = \nabla \cdot (\kappa \nabla T). \end{equation} \tag{ 31 }$$

Our temperature gradient scale (slow-mode shock width) is typically ∼500 m. For the temperature as high as T = 80, 000 K, the resulting thermal conduction F_C  ∼  2.59 × 10⁻³T^3/2δ⁻² is about 2.3 J m⁻³ s⁻¹. This falls in the range of radiative cooling relevant in our models. Considering Model III, which encompasses densities n_H ∼ 10¹⁹–10²¹ m⁻³ and corresponding ionization degrees Y_i ∼ 10⁻³–10⁻², and using α ≃ 0.01 as a conservative estimate, plasma at 80, 000 K cools radiatively at rates of 0.35–350 J m⁻³ s⁻¹. This suggests that neutral thermal conduction could limit the rise of the temperature in Model I and in the upper part of the volume in Model III, where the densities are relatively low. However, if the evolution of the ionization degree were also included, neutral thermal conduction can be expected to be negligible. At T = 80, 000 K the ionization degree is very high, e.g., Orrall & Zirker (1961) give n_n/n_i ≈ 3 × 10⁻⁶ in their Table 5, so that the estimate of the neutral conductive heat flux drops by this factor and is thus lower than the radiative losses by at least four orders of magnitude throughout the range of parameters studied in this paper. Additionally, more recent considerations of the collisional cross sections (Khomenko & Collados 2012) indicate a value for σ_nn higher by an order of magnitude than given by Orrall & Zirker. The above estimates should be checked in future simulations of chromospheric reconnection when nonequilibrium ionization is also taken into account.

Effects of potential relevance but not yet included in our model are radiative transport and non-LTE effects, nonequilibrium ionization effects, and at very small scales the Hall effect. Additionally, the complex topology of the chromospheric magnetic field has not yet been addressed. It is currently impossible to include all of them in a single model. Perhaps the most relevant steps beyond our approximations are the inclusion of recombination effects in a two-fluid description, more complex magnetic topologies guided by the observations, and radiative transport. Regarding the latter extension, we note that several estimates of the reconnection rate in plasmas with radiative cooling in Uzdensky & McKinney (2011) were found to not depend sensitively on the specific properties of the cooling process. Also, the effect of radiative losses is only a weak to moderate one in the computations presented here. These facts may be taken as an indication that the inclusion of radiative transport in the optically thick case would not strongly influence most of our results.

We thank Kai Germaschewski for his help in analyzing and displaying multilevel AMR data and James Leake, Mark Linton, Udo Ziegler, and an anonymous referee for constructive, very helpful comments. This research is supported by NSFY (grant No. 11203069), the Yunnan Province (grant No. 2011FB113), the key Laboratory of Solar Activity grant KLSA201404, China Scholarship Council (CSC) No. 201404910269, Western Light of Chinese Academy of Sciences, Program 973 grants 2011CB811403 and 2013CBA01503, NSFC grant 11273055, NSFC grant 11333007, CAS grant KJCX2-EW-T07, and CAS grant XDB09000000. B.K. acknowledges the hospitality of the solar group at the Yunnan Observatories, where most of his work was carried out, and the associated support by the Chinese Academy of Sciences under grant no. 2012T1J0017. He also acknowledges support by the DFG. We have used the NIRVANA code v3.6 developed by Udo Ziegler at the Leibniz-Institut für Astrophysik Potsdam. The authors gratefully acknowledge the computing time granted by the Yunnan Observatories and provided on the facilities at the Yunnan Observatories Supercomputing Platform, as well as the help from all faculties of the Platform.