Tearing-mediated Reconnection in Magnetohydrodynamic Poorly Ionized Plasmas. I. Onset and Linear Evolution

In high-Lundquist-number plasmas, reconnection proceeds via the onset of tearing, followed by a nonlinear phase during which plasmoids continuously form, merge, and are ejected from the current sheet (CS). This process is understood in fully ionized, magnetohydrodynamic plasmas. However, many plasma environments, such as star-forming molecular clouds and the solar chromosphere, are poorly ionized. We use theory and computation to study tearing-mediated reconnection in such poorly ionized systems. In this paper, we focus on the onset and linear evolution of this process. In poorly ionized plasmas, magnetic nulls on scales below v A,n0/ν ni0, with v A,n0 the neutral Alfvén speed and ν ni0 the neutral–ion collision frequency, will self-sharpen via ambipolar diffusion. This sharpening occurs at an increasing rate, inhibiting the onset of reconnection. Once the CS becomes thin enough, however, ions decouple from neutrals and the thinning of the CS slows, allowing the onset of tearing in a time of order νni0−1 . We find that the wavelength and growth rate of the mode that first disrupts the forming sheet can be predicted from a poorly ionized tearing dispersion relation; as the plasma recombination rate increases and ionization fraction decreases, the growth rate becomes an increasing multiple of ν ni0 and the wavelength becomes a decreasing fraction of v A,n0/ν ni0.

Theoretical understanding of reconnection in high-Lundquist-number, fully ionized, magnetohydrodynamic (MHD) plasmas has recently approached maturity.Pioneering works have considered the linear growth of the tearing instability in static CSs (Furth et al. 1963;Coppi et al. 1976) and in dynamically forming CSs (Pucci & Velli 2013;Uzdensky & Loureiro 2016;Comisso et al. 2016;Huang et al. 2017;Tolman et al. 2018) 1 and the resulting nonlinear plasmoid chain (Shibata & Tanuma 2001;Samtaney et al. 2009;Uzdensky et al. 2010;Loureiro et al. 2012).Evidence supporting an emerging theoretical picture in which reconnection in a forming CS is initiated by strongly unstable tearing modes that lead to a nonlinear plasmoid chain has been found in astrophysical observations (Nishizuka et al. 2010;Takasao et al. 2011) as well as in targeted laboratory experiments (Hare et al. 2017).
However, many plasmas of current interest in the astrophysics and laboratory plasma physics communities are not fully ionized, and reconnection may proceed in a different manner in these plasmas.For example, star-forming molecular clouds have ionization fractions ranging from ∼10 −9 to ∼10 −4 (Caselli et al. 1998;Goicoechea et al. 2009), the diffuse interstellar medium has an ionization fraction that ranges from ∼10 −4 to ∼1 (Heiles & Haverkorn 2012;Fielding et al. 2023), protoplanetary disks have ionization fractions that vary widely about ∼10 −11 ( Öberg et al. 2011;Lesur 2021), the solar chromosphere plasma has an ionization fraction ∼10 −4 -1 (Leenaarts et al. 2007), and certain laboratory reconnection experiment plasmas have ionization fractions as low as almost 10 −2 (Lawrence et al. 2013).Seminal theoretical work on reconnection in poorly ionized environments is several decades old (Zweibel 1989;Brandenburg & Zweibel 1994, 1995).More recent theoretical studies (Heitsch & Zweibel 2003;Leake et al. 2013;Singh et al. 2019;Ni et al. 2020;Murtas et al. 2021) of poorly ionized reconnection have been less plentiful and have usually not made contact with the conceptual advances developed in recent fully-ionized studies (Pucci et al. (2020) does make this contact).
Thus, in this and an accompanying paper (Tolman et al. 2024), we present analytic theory, verified by highresolution simulations, of tearing-mediated reconnection in a magnetohydrodynamic, poorly ionized plasma; this theory is of the style used in recent fully-ionized studies.
In this paper we study the plasma motion that triggers reconnection, the onset of tearing, and its linear phase of evolution.Future work (Tolman et al. 2024) will investigate the nonlinear phase of plasmoid-mediated reconnection.We consider a scenario for triggering of reconnection that is expected to describe approximately the onset of reconnection in several poorly ionized plasmas.In particular, we assume that some bulk motion in the neutral gas creates a region where the magnetic field reverses over a sufficiently short lengthscale.The associated CS will rapidly self-sharpen further via ambipolar diffusion, a diffusive process existing in poorly ionized plasmas due to the frictional drag between the ionized and the neutral species.This self-sharpening ultimately triggers reconnection.
The physics of the self-sharpening process determines how reconnection onsets.Models assumed for fully ionized CS formation often have a constant formation rate (e.g., Uzdensky & Loureiro 2016).In contrast, ambipolar sharpening of the CS causes the formation rate to accelerate, such that the tendency to tear is never able to overwhelm ambipolar sharpening.As a result, tearing cannot onset during the main CS formation process.However, at sufficiently narrow CS widths, ions completely decouple from the neutrals and the thinning of the CS slows (Brandenburg & Zweibel 1995), allowing tearing onset.This onset occurs after a time of order the inverse neutral-ion collision frequency, ν −1 ni0 .Analysis of this onset stage shows that the wavelength and growth rate of the mode that first disrupts the forming sheet can be predicted from a poorly ionized tearing dispersion relation; as the plasma recombination rate increases and ionization fraction decreases, the growth rate becomes an increasing multiple of ν ni0 and wavelength becomes a decreasing fraction of v A,n0 /ν ni0 , where v A,n0 is the neutral Alfvén velocity.After the linear onset, the system ultimately settles into a nonlinear phase characterized by a stochastic plasmoid chain.This phase is analyzed in Tolman et al. (2024).
Our analysis begins in Section 2 with a specific statement of the plasma physics problem we consider and the analytic tools used to study this problem, some of which are discussed more fully in Appendix A. Section 3 considers how a CS would form in a poorly ionized plasma in the absence of tearing and what steady state it would approach.Certain characteristics of this steady state are derived in Appendix B. Then, Section 4 uses the results of these prior sections to understand how tearing onsets in a poorly ionized, dynamically thinning CS.Section 5 verifies this analysis using high-resolution numerical simulations.Section 6 summarizes the paper, discusses its implications for plasma physics and for astrophysics, and comments on avenues for future work.An additional appendix ( §C) provides useful numbers for placing our results in an astrophysical context.

FORMULATION OF THE PROBLEM
In this Section, we state the problem we solve in this paper, provide the equations we use to study this problem, describe a common approximation applied to these equations, and present some results from the stationary linear theory of the tearing instability that will be useful later in the paper.

Problem Statement
We study the tearing-mediated onset of magnetic reconnection in a plasma comprised of an ionized, quasineutral fluid of low density and a neutral fluid of much higher density.These fluids interact via frictional drag, ionization, and recombination.The bulk motions and magnetic field configurations that might exist in such a plasma, and the ways in which they might lead to reconnection, vary, but we show in Section 3 that the essential elements of the final stages of reconnection onset in any system must be similar.This implies that a general model for the reconnecting CS should have broad application to a wide variety of poorly ionized systems.
In particular, we consider a CS that forms in a poorly ionized plasma when motions in the bulk neutral fluid move the magnetic field in a way that creates a null about which the magnetic field reverses direction.To study this system, we orient the y axis along the initial magnetic field and place the origin at the null (we do not consider a guide field).That is, (1) This field reverses direction over a length scale a(t).We are interested in how Equation (1) evolves to become where δB(x, y, t) is a magnetic perturbation driven by the tearing instability.Throughout this process, the plasma far from the CS is taken to have constant magnetic field B 0 and constant ionized and neutral fluid densities, ρ i0 and ρ n0 respectively, with their ratio determined by ionization equilibrium.

Two-Fluid MHD Equations
We model the poorly ionized plasma that hosts the reversing magnetic field using a system of hydrodynamic and MHD equations that are coupled together by frictional drag from ion-neutral collisions, by ionization of the form that would be caused by cosmic rays, and by recombination.A discussion of this type of model can be found in Draine (1986).
The neutral fluid is governed by a momentum equation, ) and a continuity equation, where ρ n is the neutral fluid mass density, v n is the neutral fluid velocity, P n is the neutral fluid pressure, ρ i is the ionized fluid mass density, v i is the ionized fluid velocity, and ν ni is the collision frequency describing the rate of collisional momentum transfer from the neutral fluid to the ionized fluid.The latter results primarily from collisions between the ions and the neutrals, with where m i is the mass of the ion species, m n is mass of the neutral species, and ⟨σw⟩ in is the mean collisional rate (collisional momentum transfer from the neutral fluid to the electrons occurs at a rate ν ne that is ∼m e /m i slower than ν ni ).We take the neutrals to be ionized at a constant rate ξ (e.g., by cosmic rays or photoionization) and produced via recombination (dissociative and/or radiative) with constant coefficient α.As the neutral and ionized populations exchange mass through ionization and recombination, their momenta are exchanged as well, a process captured by the final two terms on the righthand side of Equation (3).The momentum and continuity equations for the ionized fluid reflect Newton's third law and total mass conservation, with the former making use of quasineutrality between the ions and electrons and including the Lorentz force associated with gradients in the magnetic field B. In units where the vacuum permeability µ 0 = 1, these equations are The magnetic field evolves according to Faraday's law, with an electric field E that accounts for flux freezing in the ionized fluid but for the diffusive effects of a constant Ohmic resistivity η: This version of Ohm's law neglects Hall drifts between the electron and ion fluids (cf., Malyshkin & Zweibel 2011).Together, Equations ( 8) and ( 9) provide the non-ideal MHD induction equation In the ideal-MHD limit, in which the fluid is perfectly conducting (η = 0) and the ions and neutrals are collisionally well coupled with negligible relative drift velocity (v i ≈ v n ), the right-hand side of Equation ( 10) may be dropped and the magnetic flux is effectively frozen into the bulk neutral fluid.We close Equations ( 3) and ( 6) by adopting isothermal equations of state for both the neutral and ionized fluids, with the isothermal sound speeds C n and C i taken to be comparable to one another.
From Equation (4), we find that the ionized and neutral fluid densities far from the CS (where ionization equilibrium holds) satisfy The Alfvén speed in this predominantly neutral fluid is given by likewise, the Alfvén speed defined using the ionized fluid density is denoted These Alfvén speeds may be compared to their corresponding isothermal sound speeds to form the plasma beta parameters and Note that β i0 ≪ β n0 .The neutral-ion collision frequency in the region far from the CS where the ionized fluid density is ρ i0 is ν ni0 .A particularly important length scale for assessing the importance of ambipolar diffusion is the distance traveled by an (undamped) Alfvén wave in the predominantly neutral fluid during one neutral-ion collision timescale, The importance of resistivity is quantified using the (time-dependent) Lundquist number defined with the CS width, with its value when a(t) = a 0 being denoted by S a0 (note that our definition of S a uses the Alfvén speed in the neutrals).All variables unadorned with a subscript "0", such as ν ni , v A,n , and v A,i , refer to those quantities evaluated at a particular location and time.
Typical values for some of the parameters presented in this Section can be found in Appendix C.

Strong-Coupling Approximation
At certain moments in our study of this CS, we adopt a simplification commonly used in many MHD treatments of poorly ionized plasmas.Namely, the ion-neutral drift velocity is related directly to the Lorentz force via This simplification, called a "strong-coupling approximation," assumes that a terminal velocity is quickly established as the Lorentz force pushes the ionized fluid through the collisional neutrals.From Equation (6), it is clear that this relation neglects not only the inertia of the ionized fluid, but also its pressure and the interspecies exchange of momentum due to ionization and recombination.
The ambipolar magnetic diffusivity implied by Equations (10) and ( 18) is nonlinear in the magnetic-field strength: where the subscript "⊥" denotes the component perpendicular to the local magnetic field.Indeed, for the CS profile described by Equation (1), the right-hand side of Equation ( 19) becomes In this case, the ion-neutral drift velocity is directly proportional to −dB 2 y /dx, so that the magnetic flux is pushed towards the magnetic null at a rate that increases with B 2 y .At the same time, the ambipolar diffusion coefficient B 2 y /ρ n ν ni vanishes at the center of the sheet, and so the field close to the null is unable to smooth itself out as it steepens.As a result, without finite ionized fluid pressure and Ohmic dissipation, the magnetic gradients sharpen into singularities (Brandenburg & Zweibel 1994).In particular, the steady-state solution of Equation ( 19) is with C 1 and C 2 being arbitrary constants.For the specific situation described by Equation (1), Equation ( 21) becomes The current J to which this magnetic field corresponds is given by J The singularity of this equation at x = 0 indicates the failure of the strong-coupling approximation; in Section 3.3 we show how this singularity is resolved by finite ionized fluid pressure and Ohmic resistivity.

Tearing Growth Rates in Poorly Ionized Plasma
To predict when a forming CS might become disrupted by the tearing instability, we first require a theory for this instability's linear growth rate and characteristic scale in a stationary poorly ionized plasma.That is, we require a theory for δB(x, y, t) in Equation ( 2) for the case of B y (x, t) = B y (x).This theory will serve in Section 4 as a tool to understand the stability of the time-dependent, forming sheet.
We develop this theory in Appendix A; key results are provided here.The perturbed quantities can all be Fourier analyzed, so that the magnetic-field fluctuation is decomposed as where γ is the tearing-mode growth rate and k is a wavenumber along the CS.The growth rate pertaining to a fully ionized CS is determined implicitly by the dispersion relation (Coppi et al. 1976) where a is the characteristic width of the CS.The parameter ∆ ′ is the tearing-mode stability index (Furth et al. 1963), which is a function of k and a that depends on the functional form of the CS's equilibrium magnetic field.In addition, is the characteristic timescale on which the CS diffuses resistively, and is the characteristic timescale for an (undamped) Alfvén wave propagating in the ionized fluid at speed v A,i to cross a fluctuation having wavenumber k.The eigenvalue Λ(γ is the growth rate of the tearing mode normalized by the rate of diffusion across the inner resistive boundary layer of width δ η .
To apply these expressions for a poorly ionized system, we introduce an effective Alfvén velocity, note that this quantity depends upon γ.It turns out that Equation ( 25) applies equally well to the poorly ionized case if we simply replace This fact was first pointed out by Zweibel (1989) and has also been used by Pucci et al. (2020).Together, Equations ( 25), ( 29), and (30) determine the linear growth rate of the tearing mode in a poorly ionized plasma.It will be helpful to have not only the full dispersion relation for the tearing growth rate (Equation ( 25) with Equation( 30)), but also various asymptotic expressions.For example, the asymptotic form of Equation ( 25) for values of ka such that ∆ ′ δ η ∼ Λ ≪ 1 provides the "FKR" growth rate, γ FKR (Furth et al. 1963); values of ka such that ∆ ′ δ η ≫ 1 (i.e., Λ → 1) provide the "Coppi" growth rate, γ Coppi (Coppi et al. 1966).These limits are discussed, for example, in Boldyrev & Loureiro (2018).To adapt these expressions for a poorly ionized plasma, we replace v A,i → v ′ A,i in Equation ( 25) and approximate the stability index using This index is the small-ka limit of the tearing stability index for both the classic Harris (1962) equilibrium and, as we will show in Section 3.3, the poorly ionized CS (Equation ( 56), with a = x 1 ).The resulting asymptotic expressions are used in §4 to describe cases in which the plasma densities and collision frequencies can be assumed to be almost constant throughout the CS, and the magnetic field can be assumed to have a characteristic strength given by its value far from the sheet.Accordingly, in the remainder of this subsection we adorn values of densities, Alfvén speeds, and collision frequencies with the subscript "0".When ionization equilibrium holds, the effective Alfvén speed (Equation ( 29)) may then be written as This expression reduces further in the strong-coupling approximation, for which ξ ≪ ν ni0 : Let us also assume that the plasma is sufficiently poorly ionized and sufficiently unstable that 1 ≪ γ/ν ni0 ≪ ρ n0 /ρ i0 .Then we can simplify Equation (33) all the way to in which case it becomes beneficial to introduce as the characteristic timescale for an (undamped) Alfvén wave propagating in the ion-neutral fluid at speed v A,n0 to cross a fluctuation having wavenumber k.The poorly ionized version of the FKR growth rate, good for S −1/6 a0 ≪ ka < 1, is then This matches equation ( 18) of Zweibel (1989) when ∆ ′ is given by Equation (31).In the Coppi limit, ka ≪ S −1/6 a0 and we find that Because γ Coppi ∝ k while γ FKR ∝ k −1/2 , we can obtain an approximate expression for the maximum tearing growth rate by setting Equations ( 36) and (37) equal to one another and finding the k at which they intersect.At this maximally growing wavenumber, the growth rate is This matches equation ( 14) of Pucci et al. (2020) for the fastest-growing tearing mode in the "intermediate regime" of 1 ≪ γ/ν ni0 ≪ ρ n0 /ρ i0 .The corresponding inner resistive layer thickness, defined in (28), satisfies . By way of comparison, the maximum tearing growth rate in the regime where the ions and neutrals are fully decoupled from one another (i.e., γ/ν ) is given by This approximate expression requires that the CS width satisfies a/a 0 ≪ (

MULTI-STAGE CS FORMATION
In this Section, we describe the evolution of the profile in Equation ( 1) in the absence of the tearing instability.This description will allow us to understand when in the formation process tearing onsets.The evolution of the profile comprises three stages: a stage in which neutral advection dominates CS evolution, a stage in which ambipolar diffusion dominates the evolution, and a final steady-state stage.These stages are diagrammed schematically in Figure 1; we describe each stage, and note relevant characteristics about it, in the Sections that follow.

Neutral-Advection-Dominated Formation
We begin by considering the first stage, in which the CS reverses over a large length scale a(t).By considering the relative sizes of the terms in Equation ( 19) pertaining to neutral advection and ambipolar diffusion, we see that, on length scales above a 0 , changes in the magnetic field are caused primarily by neutral advection.Thus, initially, the CS is formed by flows in the coupled neutral fluid, as depicted in the leftmost panel of Figure 1.How precisely this stage occurs and how long it lasts depend on the nature of those flows, which vary from context to context.However, we can describe the time-dependent rate of formation during this stage as This approximate expression will be helpful in Section 4 in determining if tearing can onset during this stage.

Ambipolar-Diffusion-Dominated Formation
Once motions in the neutral fluid thin the CS enough that it varies on a length scale a(t) ≲ a 0 , the ambipolar diffusion term in Equation ( 19) dominates, and the second stage of CS formation begins.As in Section 2.1, we assume that neutral advection constantly supplies magnetic flux to the CS, such that the magnetic field is pinned to a constant value at the edge of the region in which ambipolar diffusion is occurring, Likewise, we take the plasma at the edge of the ambipolar region to have constant ionized and neutral fluid densities, ρ i0 and ρ n0 respectively, with their ratio determined by ionization equilibrium, Equation ( 12).Because ambipolar diffusion is nonlinearly dependent upon the magnetic-field strength, it attempts to sharpen the CS towards the singular current configuration given by Equation ( 22) (Brandenburg & Zweibel 1994, 1995), as described in Section 2.3 and as depicted schematically in the center panel of Figure 1.
Initially, the rate of CS formation increases with time, reflecting that the ambipolar-diffusion rate, which can be estimated using Equation (20) as increases with decreasing CS width a.Eventually, after the magnetic field has diffused nearly all of the way to the center of the null, this sharpening process is curbed by contributions from resistivity and finite ionized fluid pressure (Brandenburg & Zweibel 1995).These cause the strong-coupling assumption to break down near the center of the sheet. 2The initial width of the magnetic null is a = a 0 , so that we can use Equation ( 42) to estimate the duration of the second stage of CS formation, t form , as i.e., the time for the CS to steepen nonlinearly by ambipolar diffusion is comparable to the initial neutral-ion collision time.

Steady-State CS
After the second stage of CS formation occurs, the CS enters its final stage, a steady-state CS.This stage can be thought of as an analog of the Sweet-Parker CS in a 2 In Section 3.3, we show that the pressure gradient of the ionized fluid becomes important when the CS width a approaches the value x 1 that characterizes the effective width of a poorly ionized Sweet-Parker CS, Equation ( 51).This width is much smaller than the width of the ambipolar diffusion region surrounding the magnetic null, i.e., x 1 /a 0 ≪ 1.
fully ionized plasma.This stage is depicted in the final panel of Figure 1.Here, we describe the structure of this sheet.
Close to the singularity at x = 0, Equation ( 21) is no longer valid.In this region, the current singularity suggested by the ambipolar diffusion term is resolved by resistivity and finite ionized fluid pressure (Brandenburg & Zweibel 1995), and we cannot apply the strongcoupling approximation, Equation (18).In steady state, we have from Equation (10) that near the magnetic null at x = 0, where the final approximate equality comes from the recognition that in steady state near a magnetic null, the electric field (determined from Equation ( 9)) should remain constant while the magnetic field and x-directed ionized fluid velocity go to zero, such that the resistive term dominates (Brandenburg & Zweibel 1995).Condition (44) implies that the magnetic-field strength should be linear near the center of the CS.Far from the center of the CS, the functional form enforced by ambipolar diffusion, described by Equation ( 21), continues to hold.Altogether then, the magnetic field within the ambipolar-diffusion-dominated region can be described by where x 1 is the point where the strong-coupling assumption breaks and B 1 is the magnetic-field strength at that point.In the region where the magnetic field is linear, the magnetic pressure gradient is supported by a sharp gradient in the pressure of the ionized fluid.This magnetic configuration is shown schematically in the rightmost panel of Figure 1.
Precise values for the parameters in Equation ( 45) must be found numerically.However, we can motivate scaling laws for how they depend on the plasma properties by adopting physical arguments supported by a more detailed study of a one-dimensional version of our setup (Heitsch & Zweibel 2003).These arguments also provide scaling laws for the ionized fluid inflow velocity to the sheet, v i,AD , which transports magnetic flux inwards and supports the magnetic profile in the outer region; the characteristic velocity in the linear region −x 1 < x < x 1 , v i,η (which is smaller than v i,AD ); and the ionized fluid over-density at the center of the sheet, ρ i,CS .
First, the continuity equation for the ionized fluid, Equation (4), implies that plasma moving inside the inner linear-magnetic-field region recombines, yielding the following balance: Second, the induction equation ( 10) states that the magnetic field in the inner region is diffused away by resistivity, giving Third, approximate pressure balance throughout the CS provides where ρ n,CS is the neutral fluid density at the center of the CS, ρ n1 is the neutral fluid density at x 1 , and we have neglected the contribution of the ionized fluid pressure at x 1 and far from the CS.To simplify this expression further, we approximate that the neutral fluid pressure is approximately constant throughout the CS, so that it is primarily the gradient in the ionized fluid pressure (governed by recombination at the center of the CS) that balances the magnetic pressure gradient.As a result, and The range of validity of this assumption is checked a posteriori.From these balances, Equations ( 46)-( 50), we find the following scaling laws for the CS width4 and the reconnecting ionized fluid velocity From mass flux continuity we can estimate or, using Equations ( 49), (50), and (52), Under these scalings, a larger recombination coefficient implies a larger ionized fluid inflow velocity.
To determine the range of validity of these scalings, we can use Equation ( 54) to check the relation (50).In particular, in the case when B 1 differs noticeably from B 0 , we can relate the inflow velocity to B 1 by using Equation ( 45) to evaluate lim x→x + 1 ∇ B 2 /2 and then applying the strong-coupling relation ( 18) with ).Using Equation ( 54), we find that B 1 /B 0 ≪ 1 when C i β i0 / √ αηρ n0 ≪ 1, and our assumption that B 1 ∼ B 0 then fails.This case corresponds to the "Z ≪ 1" case in Heitsch & Zweibel (2003), in which scalings different from those above must be used.Thus, our results are good only for i.e., the recombination rate must be sufficiently small to allow the ionized fluid pressure in the inner region x 1 to grow large enough to support the CS's magnetic pressure.This criterion is satisfied in all of our numerical simulations ( §5), where α must be small so that x 1 is a resolvable length, but it breaks down in systems with low ionization fractions (such as star-forming molecular clouds; see Appendix C).In this case (i.e., Z ≪ 1), our expression for x 1 (Equation ( 51)) simply acquires a multiplicative factor of Z −1/4 (Heitsch & Zweibel 2003).
In the next section ( §4.3), these expressions are used to identify the tearing mode that disrupts the CS.To do this, it is also necessary to have a quantitative expression for the tearing-mode stability index ∆ ′ in Equation ( 25) for the profile (45).This is derived in Appendix B, with the result being that For small kx 1 , this quantity approaches ∆ ′ ≈ (kx 2 1 ) −1 , which agrees with the trend found numerically by Brandenburg & Zweibel (1995).Note that this longwavelength limit matches the tearing-mode stability index corresponding to a Harris (1962) equilibrium having width x 1 .

TEARING ONSET IN THE FORMING CS
The goal of the paper is to determine when and how in the formation process outlined in Section 3 tearing onsets.We conduct this task in this Section, starting by determining at what time tearing onsets and concluding by determining the wavelength of the mode that disrupts the sheet.

Time of Tearing Onset
Tearing onset in time-dependent systems has been considered in several recent works (Pucci & Velli 2013;Uzdensky & Loureiro 2016;Comisso et al. 2016;Huang et al. 2017;Tolman et al. 2018).A rough consensus of these works is as follows.
First, a necessary requirement for tearing growth is that the tearing stability index ∆ ′ (k) be positive for a wavenumber k that fits within the length of the CS.The value of ∆ ′ increases as k and the width decrease.The minimum possible value of k is k = 2π/L CS with L CS the CS length.We assume that our CS is sufficiently long that a k small enough for ∆ ′ to be positive exists.
Second, the growth rate of the unstable tearing mode, given implicitly by Equation ( 25) and approximately by Equation ( 38), must also be large enough to exceed significantly the rate of CS formation, allowing for an asymptotic separation between the evolution of the tearing perturbation and the evolution of the background CS (Tolman et al. 2018).Let us consider when and if these conditions are met during the formation process outlined in Section 3.
In this consideration, we require an expression for the tearing stability index ∆ ′ .The form of this index depends on the functional form of the magnetic field, and the evolving CS experiences a variety of magnetic configurations.For our analysis, we choose to use the stability which was used to derive the approximate expression (38) for the maximum tearing growth rate.This index is the small-ka limit of the tearing stability index for both the classic Harris (1962) equilibrium and the poorly ionized CS (Equation ( 56), with a = x 1 ).5 Thus, Equation ( 57) is expected to be sufficient to describe overall trends occurring in a generic narrowing magnetic configuration.
During the neutral-advection-dominated formation stage discussed in Section 3.1, we have from Equations ( 38) and ( 40) that where in the final inequality we have recognized that the minimum value of a(t) during this stage is a 0 .The value of S −2/3 a0 is small in high-Lundquist number plasmas, showing that the tearing growth rate is always smaller than the CS formation rate.As a result, tearing cannot onset during this stage.
Next we consider if an asymptotic separation of timescales between the formation rate and the tearing growth rate can occur during the ambipolar-diffusion stage discussed in Section 3.2.We compare the fastest tearing growth rate (Equation ( 38)) to the rate of CS formation via ambipolar diffusion, Equation (42): which is again a small number that remains nearly constant throughout CS formation.Asymptotic separation between the tearing growth rate and CS formation rate cannot be achieved during this second stage.6However, in the third stage identified in Section 3, when the strong-coupling approximation breaks down and the CS settles into a steady state with a steep inner profile, tearing can easily onset.At this time, given by Equation ( 43), the CS formation rate slows to zero, and the unstable tearing modes can grow.Thus, we predict that the tearing onset time is given by We conclude that tearing onsets at the end of the second stage, while the breaking of strong coupling arrests CS formation.

Dimensionless Variables
We have identified in Equation ( 16) the characteristic length scale of the problem, a 0 , and in Equation ( 60) the characteristic timescale of the problem, ν −1 ni0 .Henceforth, we refer to quantities normalized to these values and to the background plasma parameters, with normalized quantities being adorned with a tilde.That is, we replace and introduce the dimensionless variables Note that η is equivalent to the inverse Lundquist number defined using a 0 , S −1 a0 .The dimensionless neutral-ion collision frequency, is a function of the normalized ionized fluid density, recombination coefficient, and ionization coefficient.The cross section and mass that enter in Equation ( 5) affect the lengthscales and timescales of the problem but no other element of the physics.When normalized, the ionized momentum equation ( 6) becomes with similar forms holding for the other equations we solve.

Tearing Mode That Disrupts the CS
Let us now consider how tearing onsets when the strong-coupling approximation breaks.The strong coupling approximation does not break instantaneously; after the ambipolar diffusion stage, CS formation slows down continuously until the CS reaches its final state, in which the magnetic configuration does not change.Tearing onsets at some point during this slowing-down process, before CS formation has ceased completely.However, expressions for the fully formed CS parameters, given in Section 3.3, are sufficient to describe the CS at the time when tearing onsets, and the onset time is given approximately by the formation time, Equation (43).
When tearing onsets, a wide range of modes could be destabilized.These modes have growth rates determined by Equation ( 25) modified according to Equation (30), with the tearing stability index given by Equation (56). 7(Full dispersion relations, rather than asymptotic expressions, allow a deeper understanding of this stage.)A plot of such growth rates for three 7 Technically, tearing instability occurs in a background that includes the flows estimated by Equations ( 52) and ( 54).These flows do not influence tearing if the tearing growth rate is much larger than the flow shear rate (Tolman et al. 2018), i.e., where γt is the maximum growth rate and a flow is the length scale of the flow speed v flow .The flow shear is highest in the inner region of the CS, and so we may estimate the denominator of Equation ( 64) using Equations (51) (i.e., a flow ∼ x 1 ) and ( 52) (i.e., v flow ∼ v i,η ).Then we find that, for the tearing growth rate given by Equation ( 38), condition (64) holds if This condition holds in all of the high-Lundquist-number plasmas that we study.25), with ∆ ′ (k) given by Equation ( 56).In all cases, ξ = 0.1, η = 10 −4 , and Ci = Cn = 1.The proportionality coefficient for x1 in Equation ( 51) and the values for ionized density and magnetic field are set using our simulation results (see Equations ( 66)-( 68)).
different recombination parameters α is shown in Figure 2. The mode with the largest growth rate grows fastest, and quickly overwhelms any other modes that might grow.Thus, for a given set of plasma parameters, we can predict the mode that dominates linear onset by finding the value of k that maximizes Equation (25) modified according to Equation (30), with ã set to Equation ( 51) and the ionized fluid Alfvén speed defined using Equations ( 49) and ( 50).
The resulting values of k for a set of ionization fractions are shown in Figure 3.This figure illustrates a central finding of this paper: as the plasma recombination rate increases and ionization fraction decreases, the wavelength of the mode that dominates the linear stage becomes a decreasing fraction of a 0 .

ILLUSTRATION AND VERIFICATION BY DIRECT NUMERICAL SIMULATIONS
In this Section, we test the conclusions developed analytically in previous Sections of the paper numerically.Namely, we use nonlinear two-fluid numerical simulations to demonstrate that ambipolar diffusion can cause a CS to evolve from tearing stable to tearing unstable, that this process occurs in a time t ∼ 1, and that the wavenumber that dominates the linear onset can be found as described in Section 4.3.We also verify multiple aspects of the tearing physics developed in Appendix A, including the tearing-mode eigenfunction and its growth rate.

Numerical Approach
We perform a suite of simulations demonstrating CS formation, tearing onset, linear tearing growth, and the production of a chain of nonlinear plasmoids (studied further in Tolman et al. (2024)).Our numerical simulations employ the GPU-enabled astrophysical magnetohydrodynamic code suite AthenaK, which allows the simultaneous, coupled simulation of magnetohydrodynamic and hydrodynamic species.This suite also implements novel implicit-explicit methods (Pareschi & Russo 2005) that enable the efficient treatment of the different timescales of these species, which has historically been an impediment to the study of reconnection in poorly ionized plasmas (Zweibel 2015).
All of our simulations focus on the portion of the CS in which ambipolar diffusion is important, with the active mesh covering the range x ∈ [−1, 1].The length of the active mesh is ỹ ∈ [−2, 2]; the longer ỹ domain helps the plasmoids escape in the nonlinear portion of the simulation.Following Equation (1), the magnetic field in the active mesh lies primarily in the ŷ direction.However, we introduce a slight pinch in this field towards x = 0 to aid the escape of any plasmoids produced from the simulation domain.Specifically, the initial magnetic configuration is chosen to reflect the magnetic-field structure between the boundary of two adjacent magnetic islands (depicted in Fig. 1 of Huang & Bhattacharjee (2010)), i.e., we adopt a vector potential given by Ã(x, ỹ) = − 5 π tanh(x) sin(πx/5) cos(π ỹ/9) ẑ. (65) The coefficients 5 and 9 are selected such that the magnetic-field strength at (x, ỹ) = (±1, 0) is close to ±1 and the pinch is sufficient to help plasmoids escape but is not too pronounced.
The neutral fluid pressure is taken to balance the magnetic pressure approximately, following the method of Huang & Bhattacharjee (2010); the ionized fluid pressure is set such that the ionized fluid density is in ionization balance with the neutral fluid density.As in other Parameter Value ξ 0.10 α 10., 16, 22, 28, 34, 40., 60.
simulations in the literature demonstrating sharpening of a CS caused by ambipolar diffusion (e.g., Brandenburg & Zweibel 1995), this is not a pressure-balanced equilibrium, but it very quickly settles into one: the ionized fluid starts to move, and then drags on the neutrals, such that the neutral fluid pressure begins to support the magnetic pressure via the drag term.We verify that our simulations approximately satisfy the strong coupling condition (18) soon after the CS starts to evolve.
To introduce perturbations that can grow into tearing eigenmodes and to introduce asymmetry in the ŷ direction (necessary to prevent the artificial lingering of plasmoids near the center of the sheet), the ionized fluid is seeded with Gaussian-random momentum perturbations in the x and ŷ directions having amplitude 10 −7 .
Near the center of the CS, these simulations all show tearing behavior, which has very fine-scale structure that is only resolvable with high spatial resolution.Far from the CS, where no instability exists, these simulations exhibit only ambipolar diffusion, which results in relatively large-scale structure and therefore is resolvable with less resolution.Accordingly, we run our simulations using static mesh refinement (SMR): far from the center of the CS, the resolution of the simulation is 256 × 512, while the inner part of the CS is resolved using several SMR levels.The number of SMR levels for each simulation is given in Table 1.The simulations refined to a maximum SMR level of 6 are refined over a range of x values given by x ∈ [−0.03, 0.03], at a resolution equivalent to a resolution of 16384 × 32768 over the entire mesh.The simulations refined to a maximum SMR level of 7 are refined over a range of x values given by x ∈ [−0.015, 0.015], at a resolution equivalent to a resolution of 32768 × 65536 over the entire mesh.
Boundary conditions in ŷ are free outflow.Boundary conditions in x have the magnetic field fixed to the value determined by Equation ( 65), the density of both fluids fixed to their initial values at the edge of the sheet, and the momentum of both fluids set such that its values in the ghost zones are equal to its values in the last cell of the active mesh.
The range of parameters that can be effectively simulated is limited.At high ionization fractions (equivalent for fixed ionization coefficient ξ to low values of recombination coefficient α), the strong-coupling assumption does not hold.At low ionization fractions (high values of α), the scale of the reconnecting CS is very small, making resolved simulations prohibitively expensive.Thus, we complete a scan over a limited set of ionization and recombination rates, namely, ξ = 0.1 and α = 10, 16, 22, 28, 34, 40, and 60; the resulting ionization fractions at the edge of the CS are given in Table 1.

Evolution of Simulations
At the start of all our simulations, the magnetic-field profile begins to sharpen via ambipolar diffusion at an increasing rate.This evolution is highlighted in Figure 4 for the simulation with α = 10.As this process occurs, the ionized fluid pressure gradient at the center of the CS increases in magnitude, as shown in Figure 5.
Eventually, the increasing ionized fluid pressure gradient near the center of the CS breaks the strong-coupling assumption and a CS of the kind described in Section 3.3 forms (Figure 6): an outer region where the strongcoupling approximation holds meets an inner region where it is broken and the profile steepens to become linear.In the α = 10 simulation, the boundary between these two regions is measured to occur at x1 = 0.0035, where the reconnecting magnetic field B1 = 0.84.
Figure 7 demonstrates how the magnetic profile changes in time during the CS formation process by displaying the rate of change of the reconnecting field, ∂ t log Bỹ , evaluated at x = x1 .Initially, this quantity increases as the width of the CS decreases nonlinearly due to ambipolar diffusion.However, when the ionized fluid pressure gradient starts to become important (in the α simulation, around t = 0.25), the rate of change of Bỹ stops increasing and subsequently decreases towards zero.
When the rate of change is slow enough, the CS goes unstable to the tearing instability, with tearing modes clearly visible in the profile of Bx (Figure 8).For each of our simulations, we can estimate the tearing onset time by noting when the tearing eigenfunction is visible in Bx evaluated at ỹ = 0, and the CS formation time by the time at which the value of ∂ t log Bỹ decreases below unity.These values are shown in Figure 9, demonstrating that the predicted condition for onset, Equation (60), approximately holds.
The onset time and formation time both decrease slightly with α.The decrease in the formation time is likely due to the physics of the transient phase between the second and third stages of CS formation, when the strong-coupling assumption starts to break and the final stationary CS equilibrium is approached.The decrease in the onset time is likely due to both the decrease in the formation time and the increase in the linear growth rate that occurs with α.A higher linear growth rate means that tearing growth can overwhelm the CS formation earlier during the transient phase between the second and third stages.After tearing onsets, the simulations evolve to produce a nonlinear stochastic plasmoid chain, whose structure and evolution will be discussed in Tolman et al. (2024).
We can compare the wavenumber, growth rate, and eigenfunction of the linear modes observed in our simulations to the predictions obtained via the methods we  have laid out in our theory.This process starts by using our simulations to verify the predicted scalings for x 1 , ρ i,CS , and B 1 given in Section 3.3 and to obtain the corresponding numerical proportionality coefficients.For example, we fit the scaling for x 1 , Equation (51), to the values of the CS width at the time when a tearing eigenfunction is first visible in the profile of Bx at ŷ = 0 (the value of x 1 is found by identifying the point where the nearly flat magnetic-field profile associated with the strong-coupling approximation intersects the inner re- gion where the profile is linear).Figure 10 in agreement with our predictions.Using these expressions, we can make a comparison of the most unstable mode predicted via the spectra described in Section 4.3 to those observed in the simulation.This comparison is shown in Figure 11.The comparison seems quite favorable given that the higherα simulation may not be fully resolved, that the higher-α simulation approaches parameters where Equation (55) fails and other scalings for CS parameters should be used, and that our methods for have been approximate.
We can also compare growth rates measured in the simulations (found by fitting an exponential function in time to the linear growth of the peak value of the eigenfunction of tearing, seen in Bx evaluated at ỹ = 0), with the analytic predictions described in Section 4.3.Figure 12 demonstrates favorable results in this regard.
The results for mode wavenumber demonstrate less agreement than those for growth rate.This may reflect that the spectra, shown in Figure 2, are somewhat flat near the maximum growth rate, such that any mode near the mode with the maximum growth rate might onset.
Finally, the eigenfunctions of the modes observed in the simulations also appear similar to those predicted by the theory (derived in Appendix B).An example eigenfunction, taken from the α = 10 simulation at t = 1.6, is provided in Figure 13 alongside the analytically predicted eigenfunction.

CONCLUSION
The recent plasma physics literature provides an extensive discussion on the tearing-induced onset of magnetic reconnection in fully ionized plasmas.In this paper, we have presented the analogous discussion for poorly ionized plasmas, developed using analytical arguments and verified by nonlinear two-fluid MHD simulations.Our analysis yields the following important conclusions.Once the width of the forming CS shrinks below the scale a 0 (Equation ( 16)), ambipolar diffusion causes its profile to steepen nonlinearly at an increasing rate.The consequent current singularity breaks the strong-coupling approximation within a distance x 1 from the null (Equation ( 66)) and the CS formation rate stalls at a time t ∼ ν −1 ni0 .Only then is the tearing instability able to onset, grow exponentially, and eventually disrupt the CS.As the plasma recombination rate increases and ionization fraction decreases, the growth rate becomes an increasing multiple of ν ni0 and the wavelength becomes a decreasing fraction of v A,n0 /ν ni0 .
Our study has neglected several effects likely to be important in poorly ionized plasmas.Specifically, CSs in many poorly ionized plasmas are likely to have a guide field (that is, they will have a null in one component of the magnetic field only).The presence of such a guide field will inhibit CS formation through ambipolar diffusion (Zweibel & Brandenburg 1997), perhaps preventing the triggering of the tearing instability and certainly modifying the nature of its onset.In addition, CS formation occurring in a system where significant turbulent motion exists might occur in a different manner than the one we have described.Tearing onset in these situations is a fruitful avenue for future work.
In addition, poorly ionized plasmas can be subject to the Hall effect, which becomes important on scales larger than in a fully ionized plasma because the ions are weighed down relative to the electrons not only by their own inertia but also by the inertia of the neutral species to which they are collisionally coupled (e.g., Pandey & Wardle 2008).This effect may be important in the systems that we study and deserves further investigation.Finally, we have taken an isothermal equation of state, neglecting any heating that may occur during the reconnection process (both ambipolar and Ohmic).We expect that this final assumption will have limited effect on the linear stage, though it may become important in the nonlinear stage, e.g., if the heating becomes large enough to thermally ionize the plasma and thereby change the degree of ionization appreciably.This will be checked in subsequent work (Tolman et al. 2024).
Nonetheless, we can identify areas where our work is likely to advance our understanding of astrophysical plasmas.For example, one frontier area in the research of fully ionized plasmas concerns the impact of tearing modes on MHD turbulence.Arguments similar in spirit to those employed here have been used to predict a range of wavenumbers where the turbulent cascade of energy from large to small scales is mediated by reconnection (Loureiro & Boldyrev 2017).We anticipate that similar physics could be at play in poorly ionized turbulent plasmas, the most immediate applications being to protoplanetary disks and the colder phases of the interstellar medium.Concerning the latter, the structures produced by tearing onset could affect the transport of cosmic rays and the scintillation of pulsar signals (Stinebring et al. 2022;Fielding et al. 2023;Hopkins et al. 2022;Kempski & Quataert 2022;Reardon & Coles 2023).Finally, our work complements multi-fluid numerical studies of magnetic reconnection in the low solar atmosphere (e.g., Leake et al. 2013;Ni et al. 2018), in which the ionization fraction varies from ∼10 −4 in the photosphere to ∼1 at the top of the chromosphere and explosive events such as jets and flares occur that have been attributed to magnetic reconnection.With our analysis of the onset and linear evolution of tearing-mediated reconnection in poorly ionized CSs complete, understanding its nonlinear evolution and the impact on turbulent cascades and chromospheric dynamics is the clear next step.

C. TYPICAL VALUES FOR KEY PARAMETERS
In this Appendix we provide typical values for various key parameters that appear in our calculations.These values are based on fiducial conditions in star-forming molecular clouds, an example poorly ionized system in which magnetic reconnection might be important (though they may be readily scaled for application to other systems).Molecular clouds contain gas that is cold (T ≈ 10 K) and composed mostly of neutral molecular hydrogen H 2 with 20% He by number, trace amounts of electrons and molecular and atomic ions (e.g., HCO + , Na + , Mg + ), and ∼1% by mass of dust grains (both charged and neutral).The neutral fluid mass density is ρ n = 1.2 × 10 −20 n n 3 × 10 3 cm −3 g cm −3 , (C17) where m p is the proton mass and m n = 2.33m p is the mean mass per neutral particle.A typical Alfvén velocity in the neutrals is then (C20) Here we have taken the dominant ion mass to be that of HCO + (29m p ); other relevant ion species include Na + (23m p ) and Mg + (24m p ), for which radiative recombination rates should be used.Using these fiducial values, ionization equilibrium (Equation ( 12)) then provides The mean collisional rate ⟨σw⟩ in is equal to 1.69 × 10 −9 cm 3 s −1 for HCO + -H 2 collisions, and almost identical to this value for Na + -H 2 and Mg + -H 2 collisions; the presence of He effectively reduces this value by a factor of 0.81 by lengthening the slowing-down time relative to the value it would have if only HCO + -H 2 collisions were considered (Mouschovias 1996).Substituting these values into Equation ( 5), adopting ionization equilibrium, and setting m i = 29m p yields ν ni = 2.9 × 10 −13 n n 3 × 10 3 cm −3 (C24)

Figure 1 .
Figure 1.Schematic diagrams of the three stages of CS formation, from left to right: (1) an initial neutral-advection driven stage, (2) an ambipolar diffusion stage, and (3) a final steady-state CS.

Figure 3 .
Figure 3. Wavenumber k of the tearing mode predicted to dominate linear onset in a poorly ionized CS for different α and, therefore, different equilibrium ionization fractions ρi0.Other parameters are as in Figure 2.

Figure 7 .
Figure 7. Rate of change of the magnetic field at the location x = x1 vs. time, illustrating how the strong-coupling approximation breaks down in the α = 10 simulation.
Figure 8. Pseudo-color plots of Bx from the simulations with (top) α = 10 and (bottom) α = 60, with the former taken at the same time shown in Figure 6.

Figure 9 .
Figure 9. CS-formation and tearing-onset times measured in simulations with different recombination coefficients α.

Figure 10 .
Figure 10.Values of the CS width when strong coupling breaks, x1, measured in the simulations (red) and their leastsquares fit using the predicted scaling law, Equation (66).

Figure 11 .Figure 12 .Figure 13 .
Figure11.Predicted wavenumbers of the onset tearing mode (blue) compared to the dominant wavenumbers measured in the simulations (red), for each α listed in Table1.