MULTI-FLUID SIMULATIONS OF CHROMOSPHERIC MAGNETIC RECONNECTION IN A WEAKLY IONIZED REACTING PLASMA

A partially ionized reacting multi-fluid hydrogen plasma model is used to simulate reconnection in the solar chromosphere. For a detailed description and derivation of the model the reader is directed to Meier (2011) and Meier & Shumlak (2012). The model is implemented in the implicit, adaptive high order finite (spectral) element code framework, HiFi (Lukin 2008).

The full model consists of three fluids, ion (i), electron (e), and neutral (n). The fluids can undergo recombination, ionization, and charge exchange (CX) interactions, with Γ^r_α denoting the reaction rate for interaction r affecting fluid α.

The recombination reaction rate for ions (the rate of change of ion number density due to recombination), Γ^rec_i, is defined as

$$\begin{equation} \Gamma _{i}^{{\rm rec}} \equiv -n_{i}\nu ^{{\rm rec}}, \end{equation} \tag{ 4 }$$

where the recombination frequency

$$\begin{equation} \nu ^{{\rm rec}} = n_{e}\frac{1}{\sqrt{T_{e}^{*}}} 2.6\times 10^{-19}\ \hbox{m}^{3}\ \hbox{s}^{-1} \end{equation} \tag{ 5 }$$

is obtained from Smirnov (2003) and T*_e is the electron temperature T_e specified in eV.

The ionization reaction rate for neutrals, Γ^ion_n, is defined as

$$\begin{equation} \Gamma _{n}^{{\rm ion}} \equiv -n_{n}\nu ^{{\rm ion}}, \end{equation} \tag{ 6 }$$

where the ionization frequency

$$\begin{equation} \nu ^{{\rm ion}} = n_{e}A \frac{1}{X+\phi _{{\rm ion}}/T_{e}^{*}}\left(\frac{\phi _{{\rm ion}}}{T_{e}^{*}}\right)^{K}e^{-\phi _{{\rm ion}}/T_{e}^{*}}\ \hbox{m}^{3}\ \hbox{s}^{-1} \end{equation} \tag{ 7 }$$

is given by the practical fit from Voronov (1997), using the values A = 2.91 × 10⁻¹⁴, K = 0.39, X = 0.232, and the hydrogen ionization potential ϕ_ion = 13.6 eV. Note that Γ^ion_i = −Γ_n^ion, and Γ^rec_n = −Γ_i^rec.

The CX reaction rate, Γ^cx is defined as

$$\begin{equation} \Gamma ^{{\rm cx}} \equiv \sigma _{{\rm cx}}(V_{{\rm cx}})n_{i}n_{n}V_{{\rm cx}}, \end{equation} \tag{ 8 }$$

where

$$\begin{equation} V_{{\rm cx}}\equiv \sqrt{\frac{4}{\pi }v_{Ti}^{2}+\frac{4}{\pi }v_{Tn}^{2}+v_{{\rm in}}^{2}} \end{equation} \tag{ 9 }$$

is the representative speed of the interaction and $v_{{\rm in}}^{2} \equiv {|\mathbf {v}_{i}-\mathbf {v}_{n}|}^2$ with $\mathbf {v}_{\alpha }$ denoting the velocity of species α. The thermal speed of species α is given by $v_{T\alpha } = \sqrt{{2k_{B}T_{\alpha }}/{m_{\alpha }}}$, where T_α is the temperature, m_α is the corresponding particle's mass, and k_B is Boltzmann's constant. Functional forms for the CX cross-section σ_cx(V_cx) can be found in Meier (2011).

The multi-fluid model can be expressed by taking the neutral continuity equation, momentum equation, and energy (or pressure) equation to obtain a set of equations for the neutral fluid, and combining the electron and ion versions of these equations to obtain a set of equations for the "ionized" fluid. These equations are supplemented with Faraday's law, and the required transport closure equations. Assuming charge neutrality in a hydrogen plasma, the electron and ion number densities are set equal to each other (n_i = n_e). Also the ionized and neutral atom masses are set equal to the proton mass (m_i = m_n = m_p). The resulting system of partial differential equations (PDEs) is given below.

Continuity. Due to charge neutrality, only the ion and neutral continuity equations are required:

$$\begin{eqnarray} \frac{\partial n_{i}}{\partial t} + \nabla \cdot (n_{i} \mathbf {v}_i) & = & \Gamma _i^{{\rm ion}} + \Gamma _i^{{\rm rec}}, \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} \frac{\partial n_n}{\partial t} + \nabla \cdot (n_n \mathbf {v}_n) & = & \Gamma _n^{{\rm rec}} + \Gamma _n^{{\rm ion}}. \end{eqnarray} \tag{ 11 }$$

Momentum. The electron and ion momentum equations are summed and terms of order (m_e/m_p)^1/2 and higher are neglected to give

$$\begin{eqnarray} && \frac{\partial }{\partial t} (m_i n_{i} \mathbf {v}_i) + \nabla \cdot (m_i n_{i} \mathbf {v}_i \mathbf {v}_i + \mathbb {P}_i+ \mathbb {P}_e) = \mathbf {j} \times \mathbf {B} + \mathbf {R}_i^{{\rm in}} \nonumber\\ &&\quad + \Gamma _i^{{\rm ion}} m_i \mathbf {v}_n - \Gamma _n^{{\rm rec}} m_i \mathbf {v}_i + \Gamma ^{{\rm cx}} m_i (\mathbf {v}_n - \mathbf {v}_i) \nonumber\\ &&\quad + \mathbf {R}_{{\rm in}}^{{\rm cx}} - \mathbf {R}_{{\rm ni}}^{{\rm cx}}. \end{eqnarray} \tag{ 12 }$$

The neutral momentum equation is

$$\begin{eqnarray} && \frac{\partial }{\partial t} (m_i n_n \mathbf {v}_n) + \nabla \cdot (m_i n_n \mathbf {v}_n \mathbf {v}_n + \mathbb {P}_n) = - \mathbf {R}_i^{{\rm in}} + \Gamma _n^{{\rm rec}} m_i \mathbf {v}_i \nonumber\\ &&\quad -\Gamma _i^{{\rm ion}} m_i \mathbf {v}_n + \Gamma ^{{\rm cx}} m_i (\mathbf {v}_i - \mathbf {v}_n) - \mathbf {R}_{{\rm in}}^{{\rm cx}} + \mathbf {R}_{{\rm ni}}^{{\rm cx}}. \end{eqnarray} \tag{ 13 }$$

Here the current density is

$$\begin{equation} \mathbf {j} = en_{i}(\mathbf {v}_{i}-\mathbf {v}_{e})=\frac{\nabla \times \mathbf {B}}{\mu _{0}}. \end{equation} \tag{ 14 }$$

The momentum transfer R^αβ_α is the transfer of momentum to species α due to identity-preserving collisions with species β:

$$\begin{equation} \mathbf {R}_\alpha ^{\alpha \beta } = m_{\alpha \beta }n_{\alpha }\nu _{\alpha \beta }(\mathbf {v}_{\beta }-\mathbf {v}_{\alpha }), \end{equation} \tag{ 15 }$$

where m_αβ = m_αm_β/(m_α + m_β). The collision frequency ν_αβ is given by

$$\begin{equation} \nu _{\alpha \beta } = n_{\beta }\Sigma _{\alpha \beta }\sqrt{\frac{8k_{B}T_{\alpha \beta }}{\pi m_{\alpha \beta }}}, \end{equation} \tag{ 16 }$$

with T_αβ = (T_α + T_β/2). The cross-section Σ_in = Σ_ni is 1.41 × 10⁻¹⁹ m², and the cross-section Σ_en = Σ_ne is 1 × 10⁻¹⁹ m², assuming solid sphere elastic collisions (Draine et al. 1983). The cross-section for ion–electron collisions, Σ_ei = Σ_ie, is assumed to be πr²_d where r_d is the distance of closest approach (e²/(4π₀k_BT_ei) with ₀ the permittivity of free space and e the elementary charge).

The momentum transfer from species β to species α due to CX is R^cx_αβ. As found by Pauls et al. (1995) and detailed in Meier (2011) and Meier & Shumlak (2012), appropriate approximations for these terms are

$$\begin{eqnarray} \mathbf {R}_{{\rm in}}^{{\rm cx}} &\, {\approx} \, & -m_{i}\sigma _{{\rm cx}}(V_{{\rm cx}})n_{i} n_{n}\mathbf {v}_{{\rm in}}v_{Tn}^{2} \left[\! 4\left(\!\frac{4}{\pi }v_{Ti}^{2}\,{+}\,v_{{\rm in}}^2 \!\right)+\frac{9\pi}{4 }v_{Tn}^{2}\!\right]^{-1/2} \nonumber\\ \end{eqnarray} \tag{ 17 }$$

and

$$\begin{eqnarray} &&\mathbf {R}_{{\rm ni}}^{{\rm cx}} \approx m_{i}\sigma _{{\rm cx}}(V_{{\rm cx}})n_{i} n_{n}\mathbf {v}_{{\rm in}}v_{Ti}^{2}\left[4\left(\frac{4}{\pi }v_{Tn}^{2}+v_{{\rm in}}^2 \right)+\frac{9\pi}{4 }v_{Ti}^{2}\right]^{-1/2}. \nonumber\\ \end{eqnarray} \tag{ 18 }$$

The pressure tensor is $\mathbb {P}_\alpha = P_\alpha \mathbb {I} + \pi _\alpha$ where P_α is the scalar pressure and π_α is the viscous stress tensor, given by $\pi _{\alpha } = -\xi _{\alpha }[\nabla \mathbf {v}_{\alpha } + (\nabla \mathbf {v}_{\alpha })^{\top }]$ where ξ_α is the isotropic dynamic viscosity coefficient for the fluid α.

Internal energy. Again, combining the electron and ion energy equations together and neglecting terms of the order (m_e/m_p)^1/2 and higher gives:

$$\begin{eqnarray} && \frac{\partial }{\partial t}\left(\varepsilon _i + \frac{P_e}{\gamma -1}\right) \nonumber\\ &&\quad + \nabla \cdot \left(\varepsilon _i \mathbf {v}_i + \frac{P_e\mathbf {v}_e}{\gamma -1} + \mathbf {v}_i\cdot \mathbb {P}_i + \mathbf {v}_e\cdot \mathbb {P}_e + \mathbf {h}_i + \mathbf {h}_e\right) = \mathbf {j}\cdot \mathbf {E} \nonumber \\ &&\quad + \mathbf {v}_i\cdot \mathbf {R}_i^{{\rm in}} + Q_i^{{\rm in}} - \Gamma _n^{{\rm rec}}\frac{1}{2} m_i v_i^2 - Q_{n}^{{\rm rec}} + \Gamma _i^{{\rm ion}} \left(\frac{1}{2} m_i {\bf v}_n^2 - \phi _{{\rm eff}} \!\right) \nonumber\\ &&\quad + Q_i^{{\rm ion}} + \Gamma ^{{\rm cx}}\frac{1}{2}m_i\left({\bf v}_{n}^{2} - {\bf v}_{i}^{2}\right) + \mathbf {v}_n \cdot \mathbf {R}_{{\rm in}}^{{\rm cx}} - \mathbf {v}_i \cdot \mathbf {R}_{{\rm ni}}^{{\rm cx}} \nonumber \\ &&\quad + Q_{{\rm in}}^{{\rm cx}} - Q_{{\rm ni}}^{{\rm cx}}. \end{eqnarray} \tag{ 19 }$$

The neutral energy equation is

$$\begin{eqnarray} && \frac{\partial \varepsilon _n}{\partial t} + \nabla \cdot (\varepsilon _n \mathbf {v}_n + \mathbf {v}_n \cdot \mathbb {P}_n + \mathbf {h}_n)\nonumber \\ &&\quad = -\mathbf {v}_n \cdot \mathbf {R}_i^{{\rm in}} + Q_n^{{\rm ni}} - \Gamma _i^{{\rm ion}}\frac{1}{2} m_i\mathbf {v}_n^2 - Q_i^{{\rm ion}} + \Gamma _n^{{\rm rec}}\frac{1}{2} m_i\mathbf {v}_i^2 \nonumber \\ &&\qquad + Q_n^{{\rm rec}} + \Gamma ^{{\rm cx}} \frac{1}{2} m_i \big(\mathbf {v}_i^2 - \mathbf {v}_n^2 \big) + \mathbf {v}_i \cdot \mathbf {R}_{{\rm ni}}^{{\rm cx}} - \mathbf {v}_n\cdot \mathbf {R}_{{\rm in}}^{{\rm cx}} \nonumber \\ &&\qquad + Q_{{\rm ni}}^{{\rm cx}} - Q_{{\rm in}}^{{\rm cx}}. \end{eqnarray} \tag{ 20 }$$

Here ε_α ≡ m_αn_αv²_α/2 + P_α/(γ − 1) is the internal energy density of fluid α. The term Γ^ion_iϕ_eff represents optically thin radiative losses, and ϕ_eff = 33 eV (Meier 2011). The ratio of specific heats is denoted by γ. Q^αβ_α is the heating of species α due to interaction with species β, which is a combination of frictional heating and a thermal transfer between the two populations: Q^αβ_α = R^αβ_α · (v_β − v_α) +  3m_αβn_αν_αβ(T_β − T_α). The heat fluxes h_e, h_i are calculated using the Braginskii closure for a magnetized plasma, and can be written as

$$\begin{equation} \mathbf {h}_{\alpha } = [ \kappa _{\Vert ,\alpha }\hat{\mathbf {b}}\hat{\mathbf {b}} + \kappa _{\bot ,\alpha }(\mathbb {I}-\hat{\mathbf {b}}\hat{\mathbf {b}}) ]\cdot \nabla k_{B} T_\alpha , \end{equation} \tag{ 21 }$$

where κ_{||, α}(T_α) and κ_⊥,α(n_α, T_α, |B|) account for the effects of thermal diffusion parallel to and perpendicular to the magnetic field direction ($\hat{\mathbf {b}}$), respectively, and whose functional forms can be found in Braginskii (1965).

The neutral thermal conduction is isotropic h_n = −κ_n∇k_BT_n. Q^r_α denotes thermal energy gain of species α due to a reaction r, with Q^ion_i = Γ_i^ion(3/2)k_BT_n and Q^rec_n = Γ_n^rec(3/2)k_BT_i. Q^cx_αβ denotes heat flow from species β to species α due to CX (Meier 2011; Meier & Shumlak 2012). The temperature of the neutrals is given by T_n = P_n/(n_nk_B); and the ion and electron temperatures are assumed to be equal such that T_i = P_i/(n_ik_B) = P_e/(n_ek_B) = T_e.

Ohm's law. The generalized Ohm's law is given by the electron momentum equation:

$$\begin{eqnarray} && \frac{\partial }{\partial t} (m_e n_e \mathbf {v}_e) + \nabla \cdot (m_e n_e \mathbf {v}_e \mathbf {v}_e) - \Gamma _n^{{\rm ion}} m_e \mathbf {v}_n + \Gamma _i^{{\rm rec}} m_e \mathbf {v}_e \nonumber \\ &&\quad = - e n_e (\mathbf {E} + \mathbf {v}_e \times \mathbf {B}) - \nabla \cdot \mathbb {P}_e + \mathbf {R}_{e}^{{\rm ei}} + \mathbf {R}_{e}^{{\rm en}}. \end{eqnarray} \tag{ 22 }$$

In the HiFi implementation all terms on the left-hand side of Equation (22) are neglected, as in the chromosphere of the Sun, the electron inertial scale c/ω_pe is likely much smaller than magnetic diffusion length scales. However, the electron viscous stress tensor is preserved in order to represent the effects of microturbulence and three-dimensional instabilities (Che et al. 2011), and also to damp the dispersive Whistler and kinetic Alfvén waves at the shortest resolvable wavelengths. Note that the electron–neutral collision term R^en_e cannot be neglected. Using the identity v_i = v_e +  j/en_i, and the definition w ≡ v_i − v_n, Equation (22) can then be written as

$$\begin{equation} \mathbf {E} + (\mathbf {v}_i \times \mathbf {B}) = \eta \mathbf {j} + \frac{\mathbf {j}\times \mathbf {B}}{en_{i}} -\frac{1}{en_{i}}\nabla \cdot \mathbb {P}_{e} - \frac{m_{e}\nu _{{\rm en}}}{e}\mathbf {w}, \end{equation} \tag{ 23 }$$

where

$$\begin{equation} \eta =\frac{m_{e}n_{e}(\nu _{{\rm ei}}+\nu _{{\rm en}})}{(e n_{e})^2} \end{equation} \tag{ 24 }$$

is the electron, or Spitzer, resistivity, which includes electron–ion and electron–neutral collisions. The system is closed by the use of Faraday's law (∂B/∂t) = −∇ × E.

It is worth noting how the model presented here differs from the partially ionized single-fluid model used in Leake & Arber (2006) and Arber et al. (2007). The single-fluid approach implicitly assumes that the ions and neutrals are in ionization balance, which is determined by the average density and temperature, and the center of mass velocity used in Ohm's law is an average over ions and neutrals. The Pedersen resistivity, present in this partially ionized single-fluid approach, is a consequence of this center of mass velocity. The multi-fluid model presented here follows the ions and neutrals separately, thus self-consistently including the interactions between ions and neutrals which are represented by the Pedersen resistivity in the single-fluid approach. A key advantage of the multi-fluid model used in this paper is that ions and electrons are allowed to be out of local thermodynamic equilibrium (LTE). This will be shown to be vital for the onset of fast reconnection in chromospheric plasmas.

The equations are non-dimensionalized by dividing each variable (C) by its normalizing value (C₀). The set of equations requires a choice of three normalizing values. Normalizing values for the length (L₀ = 1 × 10⁵ m), number density (n₀ = 3.3 × 10¹⁶ m⁻³), and magnetic field (B₀ = 1 × 10⁻³ T) are chosen. From these values the normalizing values for the velocity ($v_{0}=B_{0}/\sqrt{\mu _{0}m_{p}n_{0}}=1.20\times 10^{5}\ \hbox{m}\,\hbox{s}^{-1}$), time (t₀ = L₀/v₀ = 0.83 s), temperature (T₀ = m_pB²₀/k_Bμ₀m_pn₀ = 1.75 × 10⁶ K), pressure (P₀ = B²₀/μ₀ = 0.80 Pa), and resistivity (η₀ = μ₀L₀v₀ = 1.5 × 10⁴ Ωm) can be derived.

The simulation domain extends from −36 L₀ to 36 L₀ in the x-direction and −6 L₀ to 6 L₀ in the y-direction, with a periodic boundary condition in the x-direction and perfectly conducting boundary conditions in the y-direction. The size of the domain has been chosen so that the boundaries and the particular boundary conditions do not affect any properties of the reconnection region centered and localized near x = 0, y = 0.

The initial neutral fluid number density is 200n₀ (6.6 × 10¹⁸ m⁻³), and the ion fluid number density is n₀ (3.3 × 10¹⁶ m⁻³). This gives a total (ion+neutral) mass density of 1.11 × 10⁻⁸ kg m⁻³, and an initial ionization level (ψ_i ≡ n_i/(n_i + n_n)) of 0.5%. The initial electron, ion, and neutral temperatures are set to 0.005T₀ (8750 K). These initial conditions are consistent with lower to middle chromospheric conditions, based on 1D semi-empirical models of the quiet Sun (Vernazza et al. 1981).

In this plasma parameter regime, ion–neutral identity-preserving collisions and CX interactions are equally important and have a very similar effect of collisionally coupling neutral and ion fluids with a neutral–ion collision mean free path of L_ni = v_{T, n}/ν_{n, i} = 140 m. However, the detailed CX physics is substantially more complicated and so for simplicity of interpretation CX interactions are neglected in this initial study and will be considered in future work. CX terms (terms with the superscript cx) in Equations (12)–(20) are thus dropped.

The ion inertial scale for these plasma parameters is c/ω_pi ≈ 1 m, which is much smaller than any scale of interest. Consequently, we neglect the Hall ($\mathbf {j}\times \mathbf {B}$) and the electron pressure tensor ($\nabla \cdot \mathbb {P}_{e})$ terms in Equation (23), the electron viscous stress tensor (∇ · π_e) in Equation (12) and Equation (19), and set electron velocity $\mathbf {v}_e$ equal to ion velocity $\mathbf {v}_i$ in Equation (19). Similarly, Malyshkin & Zweibel (2011) showed that in the geometry of a reconnection current sheet electron–neutral collisions are only important in calculating resistivity, and so the term $({m_{e}\nu _{{\rm en}}}/{e})\mathbf {w}$ in Equation (23) is also dropped. With these simplifications, we solve the following set of governing PDEs:

Continuity.

$$\begin{eqnarray} \frac{\partial n_i}{\partial t} + \nabla \cdot (n_i \mathbf {v}_i) & = & \Gamma _i^{{\rm ion}} + \Gamma _i^{{\rm rec}}, \end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray} \frac{\partial n_n}{\partial t} + \nabla \cdot (n_n \mathbf {v}_n) & = & \Gamma _n^{{\rm rec}} + \Gamma _n^{{\rm ion}}. \end{eqnarray} \tag{ 26 }$$

Momentum.

$$\begin{eqnarray} && \frac{\partial }{\partial t} (m_i n_i \mathbf {v}_i) + \nabla \cdot (m_i n_i \mathbf {v}_i \mathbf {v}_i + \mathbb {P}_i+ P_e) = \mathbf {j} \times \mathbf {B} + \mathbf {R}_i^{{\rm in}} \nonumber\\ &&\quad + \Gamma _i^{{\rm ion}} m_i \mathbf {v}_n - \Gamma _n^{{\rm rec}} m_i \mathbf {v}_i, \end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray} && \frac{\partial }{\partial t} (m_i n_n \mathbf {v}_n) + \nabla \cdot (m_i n_n \mathbf {v}_n \mathbf {v}_n + \mathbb {P}_n) \nonumber\\ &&\quad = {-} \mathbf {R}_i^{{\rm in}} + \Gamma _n^{{\rm rec}} m_i \mathbf {v}_i -\Gamma _i^{{\rm ion}} m_i \mathbf {v}_n. \end{eqnarray} \tag{ 28 }$$

Internal energy.

$$\begin{eqnarray} \frac{\partial }{\partial t}\left(\varepsilon _i + \frac{P_e}{\gamma -1}\right) & + & \nabla \cdot \left(\varepsilon _i \mathbf {v}_i + \frac{\gamma P_e\mathbf {v}_i}{\gamma -1} + \mathbf {v}_i\cdot \mathbb {P}_i + \mathbf {h}_i + \mathbf {h}_e\right) \nonumber\\ \end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray} & = & \mathbf {j}\cdot \mathbf {E} + \mathbf {v}_i\cdot \mathbf {R}_i^{{\rm in}} + Q_i^{{\rm in}} - \Gamma _n^{{\rm rec}}\frac{1}{2} m_i v_i^2 - Q_{n}^{{\rm rec}} \end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray} & +& \Gamma _i^{{\rm ion}} \left(\frac{1}{2} m_i v_n^2 - \phi _{{\rm eff}} \right) + Q_i^{{\rm ion}}, \end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray} && \frac{\partial \varepsilon _n}{\partial t} + \nabla \cdot (\varepsilon _n \mathbf {v}_n \mathbf {v}_n \cdot \mathbb {P}_n + \mathbf {h}_n) = -\mathbf {v}_n \cdot \mathbf {R}_i^{{\rm in}} + Q_n^{{\rm ni}} \nonumber \\ &&\quad - \Gamma _i^{{\rm ion}}\frac{1}{2} m_i\mathbf {v}_n^2 - Q_i^{{\rm ion}} + \Gamma _n^{{\rm rec}}\frac{1}{2} m_i\mathbf {v}_i^2 + Q_n^{{\rm rec}}. \end{eqnarray} \tag{ 32 }$$

Ohm's law.

$$\begin{equation} \mathbf {E} + (\mathbf {v}_i \times \mathbf {B}) = \eta \mathbf {j}. \end{equation} \tag{ 33 }$$

Note that to investigate the scaling of the reconnection rate with the Lundquist number (S), the Spitzer resistivity η is made a parameter of the simulations. The values used for η are [0.5, 1, 2, 4, 8, 20] × 10⁻⁵ η₀, and so η lies in the range [0.08,3] Ωm.

The viscosity coefficients for neutrals and ions are set to ξ_i = ξ_n = 10⁻³ ξ₀ and the neutral thermal conduction coefficient is κ_n = 4 × 10⁻³ κ₀. The normalizing constants are ξ₀ = m_pn₀L₀v₀ and κ₀ = m_pn₀L₀v³₀/T₀. Note that for these plasma parameters the isotropic neutral heat conduction is much faster than any of the anisotropic heat conduction tensor components for either ions or electrons, with collisional ion–neutral heat exchange Qⁱⁿ_i ≫ ∇ · (h_i + h_e). Thus, thermal diffusion for all species is dominated by neutral heat conduction and is primarily isotropic.

A Harris current sheet is used for the initial magnetic configuration, and is given here in terms of the in-plane magnetic flux A_z:

$$\begin{equation} A_{z} = - B_{0}\lambda _{\psi }\ln {\cosh {(y/\lambda _{\psi })}}, \end{equation} \tag{ 34 }$$

where $\mathbf {B} = \nabla \times A_{z}\hat{\mathbf {e}}_{z}$ and λ_ψ = 0.5 L₀ is the initial width of the current sheet.

To provide the outward force to maintain this current sheet, both the ionized pressure P_p = P_i + P_e and the neutral pressure P_n are increased in the sheet:

$$\begin{eqnarray} P_{p}(y) & = & P_{p} + \frac{1}{2}\frac{F}{\cosh ^{2}{(y/\lambda _{\psi })}}, \end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray} P_{n}(y) & = & P_{n} + \frac{1}{2}\frac{1-F}{\cosh ^2{(y/\lambda _{\psi })}}, \end{eqnarray} \tag{ 36 }$$

where F = 0.01 P₀ is chosen to maintain an approximate ionization balance of 0.5% inside the current sheet. These perturbations ensure that the total (ion and neutral) pressure perturbation balances the Lorentz force of the magnetic field in the current sheet:

$$\begin{equation} 0 = -\nabla (P_{n}+P_{p}) + \mathbf {j}\times \mathbf {B}. \end{equation} \tag{ 37 }$$

At the same time, a relative velocity between ions and neutrals is required, so that the frictional force Rⁱⁿ_i can couple the ionized and neutral fluids and keep the fluids individually in approximate force balance:

$$\begin{eqnarray} 0 & = & {-}\nabla P_{p} + \mathbf {j}\times \mathbf {B} + \mathbf {R}_{i}^{{\rm in}}, \ \hbox{and} \end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray} 0 & = & {-}\nabla P_{n} - \mathbf {R}_i^{{\rm in}}. \end{eqnarray} \tag{ 39 }$$

To give this balance, we choose an initial ion velocity

$$\begin{equation} {v_{i}}_{y}(y) = \frac{(F-1)}{n_i n_n\nu _{{\rm in}}}\frac{\tanh {(y/\lambda _{\psi })}}{\lambda _{\psi }\cosh ^{2}{(y/\lambda _{\psi })}}\frac{n_{0}^2v_{0}^2}{P_{0}}. \end{equation} \tag{ 40 }$$

To this initial steady state a small, local perturbation of the flux, A¹_z, is added to initiate the primary tearing instability and start the magnetic reconnection:

$$\begin{equation} A_{z}^{1}= - \epsilon e^{-{\left(\frac{x}{4\lambda _{\psi }}\right)}^{2}}e^{-{\left(\frac{y}{\lambda _{\psi }}\right)}^{2}}, \end{equation} \tag{ 41 }$$

where = 0.01 B₀L₀. The initial state is shown in Figure 1. Only a subset of the domain is shown, and the y-axis is stretched by a factor of 40, so that thin structures can be clearly seen in the plots. Using the symmetry of the initial conditions, only the top right quadrant of the domain is simulated with the use of appropriate boundary conditions.

Zoom In Zoom Out Reset image size

As mentioned in the previous section, η is made a parameter of the simulations, and is not dependent on the local plasma parameters, so that a general scaling law of reconnection rate with resistivity can be derived. Six simulations are performed which have values of the resistivity of [0.5, 1, 2, 4, 8, 20] × 10⁻⁵ η₀. First, the generic processes that are evident in the simulations are highlighted by focusing on one particular value of η = 0.5 × 10⁻⁵ η₀. Later, a relationship between reconnection rate and η is determined for the range of η in these simulations.

Figure 2 shows the early evolution of the reconnection region resulting from the initial perturbation. The Harris current sheet undergoes the tearing instability, magnetic field reconnects at the X-point and flux is ejected. By t = 537.5 t₀ a Sweet–Parker-like reconnection region has been formed. At this stage, ions are pulled in by the magnetic field and drag the neutrals with them.

Zoom In Zoom Out Reset image size

Figure 3 shows the ion and neutral flow fields, on a smaller subdomain of the simulation, centered on the reconnection region. Panel (a) shows the vertical velocities at time t = 537.5 t₀, with ions on the top left quadrant and neutrals on the top right quadrant. Panel (b) shows the magnitude of the horizontal velocity, again with ions on the top left quadrant and neutrals on the top right quadrant. The bottom two quadrants of panels (a) and (b) show the flow vectors for ions, on the left, and for neutrals, on the right, both scaled to the same magnitude. The vertical velocities in panel (a) show that the ions and neutrals are decoupled, with the ions flowing faster into the reconnection region than the neutrals. The difference between neutral and ion velocity is approximately 50% of the ion velocity. The horizontal velocities in panel (b) show that the ion and neutral outflows (the strong horizontal velocity at x = ±2.5 L₀) are coupled. The difference between the ion and neutral outflow is negligible compared to the actual flow. The Alfvén speed can be estimated using the upstream magnetic field of B_up = 0.6 B₀, taken at the point in the current sheet where the current amplitude reaches half its maximum value, and the total number density at the center of the sheet of 280 n₀ to give v_A = 0.035 v₀. The coupled outflow of ions and neutrals, which has a maximum of 0.015 v₀, is thus approximately half the Alfvén speed (based on total number density, as the outflow is coupled) at this time. Later in time the outflow increases to the Alfvén speed.

Zoom In Zoom Out Reset image size

This feature of the η = 5 × 10⁻⁶ η₀ simulation is common to all six simulations. The ion and neutral inflows are decoupled, as ions are pulled in by reconnecting magnetic field and the neutrals are dragged in via collisions. The timescale of inflow is fast enough that the collisions cannot keep the neutrals completely coupled to the ions and an excess of ions builds up in the reconnection region, creating an ionization imbalance. This is the situation considered for astrophysical plasmas by both Vishniac & Lazarian (1999), who treated weakly ionized reconnection with an analytic approach, and Heitsch & Zweibel (2003a), who calculated analytic and numerical solutions of 1D steady state models of weakly ionized reconnection.

As described in the previous section, the plasma in the reconnection region is out of ionization balance as the neutrals are largely left behind by the ions. Figure 4 shows the four components contributing to (∂n_i/∂t) in the ion continuity equation. The top left quadrant shows the loss due to recombination, the bottom left quadrant shows the loss due to outflow (horizontal gradient in horizontal momentum of ions), the top right quadrant shows the gain due to inflow (vertical gradient in vertical momentum), and the bottom right shows the gain due to ionization. The color scheme is based on a log-scale, and shows that within the reconnection region gains due to ionization are negligible relative to the other three terms. Looking at the largest values for the remaining three terms, the losses due to recombination and outflow are comparable, and add up to equal the gain due to inflow. This shows that the reconnection region is close to a steady state, with inflow of ions balanced by comparable contributions from recombination and outflow. Recall that the 1D models of Vishniac & Lazarian (1999) and Heitsch & Zweibel (2003a) assumed that recombination was fast enough to dominate over the outflow, and thus the horizontal direction (along the current sheet) was ignorable. Instead, Figure 4 shows that for the self-consistently created reconnection region in this parameter range the recombination and outflow are comparable and so 2D effects cannot be neglected.

Zoom In Zoom Out Reset image size

The nature of the steady state balance in these reconnection simulations has important consequences for the scaling of the reconnection rate. In the standard Sweet–Parker reconnection scenario, inflow of ions into the reconnection region is assumed to be balanced by outflow of ions. This, along with other assumptions, gives the standard scaling of the normalized inflow rate (reconnection rate) of $M\propto \sqrt{\eta } \propto 1/\sqrt{S}$. For the weakly ionized plasma in these simulations, the situation is very different due to the ionization imbalance.

The standard steady state argument of the Sweet–Parker model can be modified to include the reacting multi-fluid equations, in order to derive a scaling relationship for the reconnection rate. Figure 5 shows a schematic of the reconnection region in these simulations. The reconnection region, inside which the frozen-in constraint is broken by resistivity η, has width δ and length L. There is an external ion density of n_{i, ext} and a current sheet ion density of n_{i, CS}. The inflow of ions into the reconnection region is v_in, the outflow of ions is v_out, and the upstream magnetic field is B_up. Assuming that the system is in steady state, i.e., that (∂n_i/∂t) = 0, and integrating around the reconnection region gives

$$\begin{equation} n_{i,{\rm ext}}v_{{\rm in}}L = n_{i,{\rm CS}}(\delta v_{{\rm out}} + \delta L \nu ^{{\rm rec}} - \delta L \nu ^{{\rm ion}}), \end{equation} \tag{ 42 }$$

where ν^rec and ν^ion are the recombination and ionization frequencies, defined in Equations (5) and (6). Defining ν_out ≡ v_out/L, and equating Ohm's law evaluated inside (E = ηj) and outside (E = v_inB_up) of the reconnection region,

$$\begin{equation} v_{{\rm in}}B_{{\rm up}} = \eta j \approx \frac{\eta B_{{\rm up}}}{\delta \mu _{0}}, \end{equation} \tag{ 43 }$$

to eliminate δ, this steady state equation can be rewritten as

$$\begin{equation} v_{{\rm in}} \approx \sqrt{ \frac{\eta }{\mu _{0}}\frac{n_{i,{\rm CS}}}{n_{i,{\rm ext}}} (\nu _{{\rm out}} + \nu ^{{\rm rec}} - \nu ^{{\rm ion}})}. \end{equation} \tag{ 44 }$$

Zoom In Zoom Out Reset image size

For a plasma in ionization balance, ν^rec = ν^ion, and n_{i, ext} ≈ n_{i, CS}. Using these relationships and the total momentum equation to derive v_out = v_A (where $v_{A} = ({B_{{\rm up}}}/{\sqrt{\mu _{0}\rho _{{\rm total}}}})$), recovers the standard Sweet–Parker scaling law

$$\begin{equation} M\equiv \frac{v_{{\rm in}}}{v_{A}} \approx \sqrt{\frac{\eta }{v_{A}L\mu _{0}}} = \sqrt{\frac{1}{S}}, \end{equation} \tag{ 45 }$$

where B_up is the magnetic field upstream of the current sheet and ρ_total is the total (ion + neutral) density in the reconnection region.

In contrast, in these reacting two-fluid simulations, the system nonlinearly and self-consistently forms a current layer where plasma is out of ionization balance, where the recombination is comparable to the outflow, and where ionization is negligible compared to both recombination and outflow. Hence the inflow rate can be approximated by

$$\begin{equation} v_{{\rm in}} \approx \sqrt{ \frac{\eta }{\mu _{0}}\frac{n_{i,{\rm CS}}}{n_{i,{\rm ext}}}(\nu _{{\rm out}} + \nu ^{{\rm rec}}) }. \end{equation} \tag{ 46 }$$

Note that this equation does not by itself indicate how the reconnection rate will depend on S, as it is not clear how ν^rec, n_{i, CS}, and n_{i, ext} depend on S from the equations. We will therefore use our series of simulations performed over a range of η to determine the reconnection rate dependence on S.

For all six simulations the effective Lundquist number is defined by

$$\begin{equation} S_{{\rm sim}} \equiv \frac{v_{A,0}L_{{\rm sim}}\mu _{0}}{\eta }, \end{equation} \tag{ 47 }$$

where L_sim = 2.75 L₀. Note that L_sim is much smaller than the horizontal extent of the domain. We define v_{A, 0} as the Alfvén velocity given a field strength of B₀ and number density 201 n₀, i.e., based on the initial background magnetic field and number density. Note that S_sim varies over simulations due to η only. The reconnection rate is defined by

$$\begin{equation} M_{{\rm sim}} \equiv \frac{\eta j_{{\rm max}}}{v_{{\rm out}}B_{{\rm up}}}. \end{equation} \tag{ 48 }$$

Here, j_max is the maximum value of the current density, located at (x, y) = (x_j, 0), within the reconnection region, and B_up is evaluated at (x_j, δ_sim) where δ_sim is the y location on the line x = x_j at which the current density reaches j_max/2, i.e., half-width at half-maximum of the current sheet.

For each simulation, the value of M_sim is taken at a time when the length of the current sheet, defined by the distance from the original X-point to the location of maximum outflow v_out equals L_sim. This happens at a different time for each simulation, but in all cases prior to the onset of any secondary instabilities of the current sheet. This instantaneous value of M_sim is plotted for each simulation against the effective Lundquist number S_sim in Figure 6. The dashed line shows the single-fluid Sweet–Parker predicted scaling $M \propto 1/\sqrt{S}$. Over the range of S which these six simulations cover, the reconnection rate is only weakly dependent on the Lundquist number, with a slow decrease in M_sim with increasing S_sim. The reconnection rate is much faster than the Sweet–Parker reconnection rate, as a direct result of the decoupling of ions from neutrals and the enhanced recombination in the reconnection region. As discussed in Section 3.3, plasmoid formation increases the reconnection rate for S_sim > S_c. We also note that each of these simulations have been performed with the same viscosity coefficient, and it is known that the reconnection rate dependency on resistivity weakens when viscous effects in the reconnection layer become important (Park et al. 1984).

Zoom In Zoom Out Reset image size

Figure 7 shows the scaling of the reconnection region aspect ratio σ_sim ≡ (δ_sim/L_sim) with S_sim for the six simulations. Also shown is the Sweet–Parker scaling of $1/\sqrt{S}$. The simulations do not exhibit the $\sigma \propto 1/\sqrt{S}$ scaling but show an approximate 1/S scaling (the power-law fit through the simulation points has an exponent of −1.1 ± 0.17). This 1/S scaling is consistent with the result of Figure 6, that M_sim is approximately independent of S_sim. This is shown as follows: assume that v_out does not substantially vary over different simulations, and note that S_sim only changes over simulations due to η. Then using the definition of M_sim and j_max ≈ B_up/δ_sim gives M_sim∝ηj_max/B_up ≈ η/δ_sim∝1/S_simσ_sim ≈ const.

Zoom In Zoom Out Reset image size

As discussed above, at a critical aspect ratio of σ_c = 1/200, a resistive current sheet can become unstable to a secondary tearing instability known as the plasmoid instability. Of the six two-fluid simulations in this paper, three cases show evidence of secondary tearing. Figure 8 shows the onset of the plasmoid instability for the simulation with η = 0.5 × 10⁻⁵ η₀, with the two contours of magnetic flux at A_z = −0.0681 B₀L₀ (white line) and A_z = −0.0687 B₀L₀ (black line) following the evolution of two particular fieldlines in time. The initial laminar current sheet breaks up into a number of plasmoids, with thinner current sheets between them. The current density and recombination rate within this fragmented reconnection region are shown in Figure 8 on a pseudo-color scale. During this phase of the reconnection the formation of a plasmoid chain redistributes the ionized plasma within the current sheet, which is otherwise approximately uniformly distributed throughout the reconnection region, into regions of higher and lower electron number density. Since Γ^rec_i∝n_e², the recombination rate is increased within the plasmoids with respect to the sub-layers between them. This leads to the plasmoid magnetic flux collapsing on itself (i.e., disappearing as it would in vacuum) at a rate comparable to or faster than the plasmoids are exhausted out of the reconnection region. This collapse can affect the expected distribution of plasmoid sizes and the ionization fraction within the plasmoids relative to the background medium in the reconnection exhaust.

Zoom In Zoom Out Reset image size

Huang & Bhattacharjee (2010) found that in a fully ionized plasma, the onset of the plasmoid instability occurred at a critical value of Lundquist number of 4 × 10⁴ which, with the Sweet–Parker scaling of current sheet width with Lundquist number, corresponds to a critical current sheet aspect ratio of 1/200, as shown by the dashed line in Figure 7. The three simulations that undergo secondary tearing are shown as diamonds in Figure 7. The intersection of a power-law fit to the σ_sim(S_sim) data intersects the σ_sim = 1/200 line at a value of S_sim = 10⁴. The simulation with S_sim = 10⁴ exhibits the plasmoid instability, as do the two simulations with higher S_sim (shown as diamonds in the figure). These two facts support the postulation of Loureiro et al. (2007) that the criterion for the plasmoid instability is that the aspect ratio decreases below a critical value. For our simulations, as for Huang & Bhattacharjee (2010), this critical value is 1/200.

To demonstrate the change in reconnection rate due to the plasmoid instability, M_sim is evaluated again later in each of the three plasmoid-unstable simulations. It is difficult to measure M_sim during the plasmoid instability, as the reconnection rate varies with time, depending on the plasmoid evolution. To give an idea of how the plasmoid instability is affecting the reconnection rate we plot in red in Figure 6 the range of M_sim observed after plasmoid formation in each simulation until the current sheet width becomes spatially unresolved.

It is apparent, particularly in higher S_sim simulations, that the reconnection rate increases as the plasmoids develop. Thus, the reconnection rate in this regime is determined by both the fast recombination of ions, and the effect of the secondary tearing instability.