Coalescence Instability in Chromospheric Partially Ionised Plasmas
PPhys. Plasmas , 032901 (2021); https://doi.org/10.1063/5.0032236 , 032901© 2021 Author(s). Coalescence instability in chromosphericpartially ionized plasmas
Cite as: Phys. Plasmas , 032901 (2021); https://doi.org/10.1063/5.0032236Submitted: 07 October 2020 . Accepted: 30 January 2021 . Published Online: 02 March 2021 Giulia Murtas, Andrew Hillier, and Ben Snow
COLLECTIONS
This paper was selected as an Editor’s Pick oalescence instability in chromospheric partiallyionized plasmas
Cite as: Phys. Plasmas , 032901 (2021); doi: 10.1063/5.0032236Submitted: 7 October 2020 . Accepted: 30 January 2021 . Published Online: 2 March 2021Giulia Murtas, a) Andrew Hillier, and Ben Snow
AFFILIATIONS
College of Engineering, Mathematics and Physical Sciences, Harrison Building, Streatham Campus, University of Exeter,North Park Road, Exeter EX4 4QF, United Kingdom a) Author to whom correspondence should be addressed: [email protected]
ABSTRACT
Fast magnetic reconnection plays a fundamental role in driving explosive dynamics and heating in the solar chromosphere. Thereconnection time scale of traditional models is shortened at the onset of the coalescence instability, which forms a turbulent reconnectingcurrent sheet through plasmoid interaction. In this work, we aim to investigate the role of partial ionization in the development of fastreconnection through the study of the coalescence instability of plasmoids. Unlike the processes occurring in fully ionized coronal plasmas,relatively little is known about how fast reconnection develops in partially ionized plasmas (PIPs) of the chromosphere. We present 2.5Dnumerical simulations of coalescing plasmoids in a single fluid magnetohydrodynamic (MHD) model and a two-fluid model of a partiallyionized plasma (PIP). We find that in the PIP model, which has the same total density as the MHD model but an initial plasma densitytwo orders of magnitude smaller, plasmoid coalescence is faster than the MHD case, following the faster thinning of the current sheet andsecondary plasmoid dynamics. Secondary plasmoids form in the PIP model where the effective Lundquist number S ¼ : (cid:2) , but areabsent from the MHD case where S ¼ : (cid:2) : these are responsible for a more violent reconnection. Secondary plasmoids also form inlinearly stable conditions as a consequence of the nonlinear dynamics of the neutrals in the inflow. In the light of these results, we can affirmthat two-fluid effects play a major role in the processes occurring in the solar chromosphere. Published under license by AIP Publishing. https://doi.org/10.1063/5.0032236
I. INTRODUCTION
Magnetic reconnection is a process responsible for explosiveevents in astrophysical, space, and laboratory plasmas, allowing plas-mas to break free from the frozen-in constraint imposed by ideal mag-netohydrodynamics (MHD). The magnetic field lines reconnectthrough a narrow diffusion region, enabling the conversion of storedmagnetic energy into kinetic and thermal energy, and particleacceleration.
A classical description for reconnection is provided by theSweet–Parker model, which describes steady-state reconnection in along diffusion region of length L and width d . The reconnection timescale goes as s SP (cid:3) L = d (cid:3) S = , where S ¼ Lv A = g is the Lundquistnumber, g is the diffusivity, and v A is the Alfv (cid:2) en speed.The presence of localized, transient outflows in the solar chromo-sphere (such as the chromospheric jets ) is believed to be the result ofthe onset of magnetic reconnection. Further evidence for magneticreconnection is provided by the identification of bubbles of plasma inthe outflow of chromospheric jets and UV bursts, which are gener-ally interpreted to be plasmoids. Plasmoids are commonly present in reconnecting systems, and they are believed to play a major role in speeding up reconnec-tion by having an influence on the variation of the current sheetsize. Under the formation of plasmoids, the current sheet breaksinto fragments, and the resulting high current densities in each ofthese sections facilitate a high reconnection rate. Plasmoid for-mation due to the instability of Sweet–Parker current sheets hasbeen extensively examined through numerical studies andfound in direct imaging observations of solar flares. Severalworks proved that in fully ionized plasmas it is possible to triggerplasmoid formation in thin current sheets ( d = L (cid:4)
1) for a criticalLundquist number (cid:3)ð (cid:5) Þ . Numericalstudies found an upper limit for the critical Lundquist number (cid:3) (cid:2) . Above this critical value of S , the current sheet becomesrapidly unstable to the resistive tearing instability, forming a chainof plasmoids along the current sheet. In astrophysical plasmas with very large S , the onset of the tearingmode takes place in current sheets that collapse to a thickness on theorder of S (cid:5) = , larger than the Sweet–Parker thickness of S (cid:5) = . Phys. Plasmas , 032901 (2021); doi: 10.1063/5.0032236 , 032901-1Published under license by AIP Publishing Physics of Plasmas
ARTICLE scitation.org/journal/php he trigger of fast reconnection before a Sweet–Parker-type configura-tion can form was examined in detail by several studies.
Many works proved that the critical Lundquist number andaspect ratio for the onset of the tearing instability are affected bychanges in the initial setup, which might result in a discrepancy oforders of magnitude in the critical Lundquist number.
A majorrole is played by the initial current sheet configuration as well asby the amplitudes of viscosity and perturbation noises, and theplasma b . The variation of the critical S as a function of theinitial noise, investigated in some works, covers several orders ofmagnitude (from S (cid:3) to S (cid:3) ). The role of the configurationof the simulation domain in affecting the threshold of the criticalLundquist number has been pointed out in studies of magnetic recon-nection in fully ionized plasmas. In these works, it was shown thata 2.5D simulation of magnetic reconnection with a force-free currentsheet and uniform plasma pressure as initial conditions lead to a muchlower critical S than those obtained in 2D cases with an initial Harriscurrent sheet and nonuniform plasma pressure. Their results are con-sistent with the findings of a recent study, where magnetic reconnec-tion is examined on different spatial scales in weakly ionized plasmasby using a reactive 2.5D multi-fluid plasma–neutral model. The nonlinear tearing mode shows the development of an impor-tant secondary instability called the coalescence instability. This insta-bility is driven by the coalescence of neighboring plasmoids sharing anX-point and results from the attractive forces between parallel cur-rents. The coalescence instability is characterized by two differentphases: ideal and resistive. The ideal phase has a growth rate that isalmost independent of g . The resistive phase is driven by the currentsheet reconnection. In a single-fluid MHD approach, the coalescenceinstability produces a reconnecting current sheet by driving plasmoidinteraction and is a key process that might explain fast reconnectionwithout the need of anomalous resistivity terms to be added into thesystem. In the solar chromosphere plasmas, as well as many other plas-mas found in the universe, are partially ionized, and their ionizationdegree falls in the range of 10 (cid:5) (cid:5) (cid:5) . Multi-fluid effectslinked to the different behavior of the particle species must be takeninto account for a correct physical description of this atmosphericlayer. The low chromospheric densities do not allow a complete colli-sional coupling between ions and neutral species: the low ion fractionallows the fewer charged particles to be coupled to the neutrals, butthe neutrals may not be completely coupled to the ions. A partial cou-pling between the two species results in the presence of relativemotions. While ideal MHD equations are applicable to fully ionizedplasmas, two-fluid effects should be described by taking into accountthe relative motions between plasma and neutral components in thesolar chromosphere. The role of partial ionization in the onset of mag-netic reconnection and development of the resistive tearing instabilitywas investigated in many studies.
In a system where thetwo fluids are coupled through elastic collisions and subject to ioniza-tion and recombination, the reconnection rate in the coalescence pro-cess depends on the ion fraction.
The role of fast magnetic reconnection in triggering events in thesolar chromosphere is undoubtedly fundamental. An additional com-plexity comes from the partially ionized nature of the solar chromo-sphere. In this paper, we discuss the role of partial ionization inmagnetic reconnection through the study of plasmoid coalescence, with the aim of understanding to what extent the two fluid effectsinfluence such process. In order to do this, we first compare two refer-ence MHD and partially ionized plasma (PIP) simulations (Sec. III)and then investigate in more detail how two-fluid properties affect thecoalescence instability through a parameter survey (Sec. IV). In Sec. V,our results are connected back to the physical scales of reconnection inthe solar chromosphere.
II. METHODS
We perform simulations of the coalescence instability in fully(MHD) and partially ionized plasmas (PIP), using the (PIP) code. The code makes use of a four-step Runge-Kutta and a fourth-ordercentral difference scheme. The fully ionized plasma consists of asingle-fluid model of a hydrogen plasma. The partially ionized plasmaenvironment is simulated through a two-fluid model consisting of twoseparate sets of equations, describing a neutral fluid and a charge-neutral ion–electron plasma which are collisionally coupled. The equa-tions are derived from those found in previous models.
All sets of equations are non-dimensionalized. The choice of thisparticular normalization is performed in order to have a dependencyon characteristic length scales, which are comparable to the size of theplasmoids involved in the merging. This allows the model to beapplied to plasmoids at different scales in the solar chromosphere,from a few meters to a few hundred kilometers, as all the quantitiesscale together with the plasmoid size. Such normalization also affectsthe collisional coupling between the two fluids. As a characteristicdimensional time scale s can be derived from the physical plasmoidsize and the characteristic speed in the solar chromosphere (in ourcase this is the sound speed), the non-dimensional collisional fre-quency is easily re-scaled back to physical quantities by dividing it by s . Further details on the non-dimensionalization are provided at theend of this section and in Sec. II A.The neutral fluid is described by non-dimensional inviscidhydrodynamics equations @ q n @ t þ r (cid:6) ð q n v n Þ ¼ ; (1) @@ t ð q n v n Þ þ r (cid:6) ð q n v n v n þ p n I Þ ¼ (cid:5) a c ð T n ; T p ; v D Þ q n q p ð v n (cid:5) v p Þ ; (2) @ e n @ t þ r (cid:6) v n ð e n þ p n Þ½ (cid:7) ¼ (cid:5) a c ð T n ; T p ; v D Þ q n q p (cid:2) ð v n (cid:5) v p Þ þ p n q n (cid:5) p p q p (cid:2) (cid:3)" ; (3) e n ¼ p n c (cid:5) þ q n v n ; (4) T n ¼ c p n q n ; (5)while inviscid resistive magnetohydrodynamics relations govern theplasma, here stated in non-dimensional form @ q p @ t þ r (cid:6) ð q p v p Þ ¼ ; (6) Physics of Plasmas
ARTICLE scitation.org/journal/php
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.0032236 , 032901-2Published under license by AIP Publishing @ t ð q p v p Þ þ r (cid:6) q p v p v p þ p p I (cid:5) BB þ B I (cid:2) (cid:3) ¼ a c ð T n ; T p ; v D Þ q n q p ð v n (cid:5) v p Þ ; (7) @@ t e p þ B (cid:2) (cid:3) þ r (cid:6) v p ð e p þ p p Þ þ (cid:5)ð v p (cid:2) B Þ (cid:2) B þ g ðr (cid:2) B Þ (cid:2) B (cid:4) (cid:5) ¼ a c ð T n ; T p ; v D Þ q n q p ð v n (cid:5) v p Þ þ p n q n (cid:5) p p q p (cid:2) (cid:3)" ; (8) @ B @ t (cid:5) r (cid:2) ð v p (cid:2) B (cid:5) g r (cid:2) B Þ ¼ ; (9) e p ¼ p p c (cid:5) þ q p v p ; (10) r (cid:6) B ¼ ; (11) T p ¼ c p p q p : (12)In the equations above, the subscripts p and n refer, respectively, to theion–electron plasma and the neutral fluid, v , p , q , T , and e are thevelocity, gas pressure, density, temperature, and internal energy ofeach species, c ¼ = B is the magneticfield. Both fluids follow the ideal gas law. The factor 2 in Eq. (12) is totake into account the electron pressure. The parameter a c , given by Eq.(13), is associated with the two fluids collisional coupling. This particu-lar formulation of a c is new to models of magnetic reconnection inpartially ionized plasmas, and it accounts for the increased amount ofcollisions at supersonic drift velocities. The non-dimensional expres-sion for a c including charge exchange is as follows: a c ¼ a c ð Þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi T n þ T p r ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ p c ð T n þ T p Þ v D s ; (13)where a c ð Þ is the initial coupling and v D ¼ j v n - v p j is the magnitudeof the drift velocity between the neutral components and the hydrogenplasma. When the drift velocity becomes bigger than the thermalvelocity, the particles are subject to a higher number of collisions asthey are drifting past each other. The collisional coupling between ionsand electrons is represented by setting a small finite diffusivity g that isassumed to be spatially uniform and not varying with time. In thiswork, we are not including the Hall effect.The two systems of equations are non-dimensionalized by areference density q and the total sound speed c s ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c ð p n þ p p Þ = ð q n þ q p Þ q , initially set equal to 1. For the MHD sim-ulation, where the plasma is fully ionized, the initial density and pres-sure are constant and equal to q p ¼ n p q ¼ ; (14) p p ¼ p ¼ c (cid:5) : (15)For the PIP simulations, the bulk density and pressure are equalto the MHD values q n þ q p ¼ n n q þ n p q ¼ ; (16) p n þ p p ¼ n n ð n n þ n p Þ p þ n p ð n n þ n p Þ p ¼ c (cid:5) ; (17) and they are uniform in all the domain. Initially the two fluids are inthermal equilibrium. A. Initial conditions
The bulk initial conditions of the PIP case are equal to the initialconditions of the MHD case. The initial setup is provided by a force-free modified Fadeev equilibrium.
The magnetic scalar potential ofthe classical 2D Fadeev equilibrium in the xy (cid:5) plane is given by w ð x ; y Þ ¼ B k ln cosh ð ky Þ þ (cid:2) cos ð kx Þ½ (cid:7) ; (18)where B is the field intensity for the limit j y j ! 1 . In our simula-tions, B is equal to ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c (cid:5) b (cid:5) p , where plasma b ¼ p = B ; k ¼ p ,and we set (cid:2) ¼ :
5, which corresponds to a moderately peaked currentlocalization at the plasmoid center, shown in Fig. 1. As (cid:2) ! (cid:2) ! (cid:2) ¼ b may becomevery small. Although the photospheric magnetic field emerging fromthe convection zone is not force-free, its structure is rearranged by thetime it reaches the corona as the non-force-free components decaydue to the action of chromospheric neutrals. It is hence of interest tostudy the coalescence instability in a regime that is initially force-free.The magnetic field B x and B y components from the classic Fadeevequilibrium are not sufficient to satisfy the condition J (cid:2) B ¼ B z , making it force-free. The magneticfield components are shown in the following equations: B x ¼ (cid:5) ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c (cid:5) b (cid:5) p (cid:2) sin ð kx Þ cosh ð ky Þ þ (cid:2) cos ð kx Þ ; (19) B y ¼ (cid:5) ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c (cid:5) b (cid:5) p sinh ð ky Þ cosh ð ky Þ þ (cid:2) cos ð kx Þ ; (20) B z ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c (cid:5) b (cid:5) p ffiffiffiffiffiffiffiffiffiffiffiffi (cid:5) (cid:2) p cosh ð ky Þ þ (cid:2) cos ð kx Þ : (21) FIG. 1.
Initial distribution for the current density J z ( t ¼
0) for both the MHD and thePIP simulations. Two initial plasmoids (blue spots), which are concentrations of pos-itive current, are positioned with their center on the x -axis. The magnetic field linesare displayed in black. The left boundary of the domain for the reference MHD sim-ulation and the cases discussed in the parameter survey is shown in green. Physics of Plasmas
ARTICLE scitation.org/journal/php
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.0032236 , 032901-3Published under license by AIP Publishing etting (cid:2) ¼ B y component leads to acurrent sheet with the characteristic tanh profile of the well knownHarris sheet. Our initial conditions for the current density are dis-played in Fig. 1 and are the same for both MHD and PIP cases.The Fadeev equilibrium is unstable to the coalescence instability. We hence choose a small perturbation in the velocity of both plasmaand neutral components to break the initial equilibrium by pushingneighboring plasmoids toward each other. The velocity perturbation isgiven by v x ; p ¼ v x ; n ¼ (cid:5) :
05 sin kx (cid:2) (cid:3) e (cid:5) y þ v noise ; (22)where v noise is a white noise component simulating small environmen-tal perturbations. The sine term dependent on x in the main perturba-tion results in a push on the pair of plasmoids, so they move closer toeach other. As the domain has periodic boundary conditions on thesides, there is an effective chain of plasmoids moving along the x -axis.The perturbation causes the coalescence to take place for each pair ofplasmoids separately, while moving the other plasmoids away. Theterm dependent on y localizes the perturbation to a small regionaround the plasmoids center. The white noise perturbation, which isset equal for plasma and neutrals, has a magnitude of 0.0005, twoorders of magnitude smaller than the main perturbation in Eq. (22).Choosing such value prevents the noise from dominating the motionof the two plasmoids during coalescence, but allows the developmentof dynamics at a smaller scale by breaking the symmetry of the system.The same random noise seed was used in all simulations.The reference MHD simulation in Sec. III is resolved by2062 (cid:2) D x ¼ : (cid:2) (cid:5) and D y ¼ : (cid:2) (cid:5) . In the PIP case, the partial ionizationeffects lead to the thinning of the current sheet and the developmentof sharp small-scale magnetic structures as a result of the neutralsdecoupling from ions and leaving the current sheet. Our sim-ulations show the formation of a thinner current sheet in between theplasmoids merging due to this two-fluid effect. The reference PIP casein Sec. III was hence run at the higher resolution of D x ¼ : (cid:2) (cid:5) and D y ¼ : (cid:2) (cid:5) to ensure the current sheet is resolved by ourgrid. The grid in the PIP case is composed of 6478 points in the x direction and 4862 points in the y direction. In Sec. IV, we present aparameter study of the coalescence process. The resolution used foreach simulation and the total number of grid points are detailed inSec. IV.The initial separation between the plasmoids (calculated from O (cid:5) point to O (cid:5) point, identified as blue spots in Fig. 1) is equal to 4 L ,where L is resolved by 515 grid points in the MHD case and by 809grid points in the PIP case. The plasmoid width, calculated as the dis-tance between top and bottom edges of the separatrix and which initialvalue is 1 : L (resolved, respectively, by 638 grid points in the MHDcase and by 1037 grid points in the PIP case), is determined by theconditions of the Fadeev equilibrium for the magnetic field.The non-dimensional diffusion length scale is calculated as L diff ¼ ffiffiffiffiffiffiffi gs p for both the reference cases and the simulations of theparameter survey. Taking s ¼
1, the diffusion length scale for theMHD and PIP reference cases is 4 : (cid:2) (cid:5) for a diffusivity g ¼ (cid:2) (cid:5) , while for the cases in the parameter survey, which arecharacterized by g ¼ : (cid:2) (cid:5) ; L diff ¼ (cid:2) (cid:5) . The approximatenumber of grid cells per diffusion length scale at the lower resolution of the MHD case is (23, 17) in ( x , y ). These, respectively, increase to 41grid cells along x and 31 grid cells along y for L diff in the parametersurvey. Therefore, the diffusion scale is resolved in all simulations. Thecollisional ion–neutral time scale s col ; pn ¼ ð a c q n Þ (cid:5) for the referencePIP case is 10 (cid:5) , while the collisional neutral–ion time scale, definedas s col ; np ¼ ð a c q p Þ (cid:5) , is 1. These values lead to the non-dimensionalcoupling length scales L pn ¼ (cid:5) and L np ¼
1, which are both wellresolved by our grid.
B. Boundary conditions
While coalescing, the plasmoids move toward each other alongthe x -axis. Because of the symmetry of the problem, in the referenceMHD simulation (Sec. III) and in the set of simulations performed forthe parameter survey (Sec. IV), we cut the computational domain at x ¼ v x and B y change signacross each boundary and v y and B x remain the same. The left bound-ary is shown in Fig. 1 as a dashed green line. The computationaldomain size is chosen equal to x ¼ ½ ; (cid:7) and y ¼ ½(cid:5) ; (cid:7) .The dynamics of the plasmoids merging in the reference PIP case(Sec. III) is evaluated in a full computational domain, with x ¼ ½(cid:5) ; (cid:7) and y ¼ ½(cid:5) ; (cid:7) . This arrangement was made to be able to betterexamine the dynamics in the region of the current sheet at higher reso-lution. In this case, the top and bottom boundaries are kept symmetric,while the side boundaries are chosen to be periodic. III. RESULTS
First, we explore the coalescence instability in both a single-fluidfully ionized plasma and a two-fluid partially ionized plasma by com-paring two simulations, an MHD case and a PIP case. The single-fluidcase acts as a reference case for the more complex two-fluidsimulation.The initial parameters for both simulations are constant diffusivity g ¼ : b ¼ :
1. In the PIP simulation, we set the col-lisional coefficient a c ¼
100 and the ion fraction n p ¼ q p = ð q p þ q n Þ¼ :
01, while the effective ion fraction in the MHD case is n p ¼
1. Theeffects of the parameters variation on the coalescence in PIP simula-tions are investigated later on in Sec. IV.Figure 2 displays a sequence of the evolution of the current den-sity J z , which is directed out of the plane. For a better comparison, theframes show times where similar physics takes place in both the MHDand PIP cases. As g
0, the magnetic flux reconnects and leaves thecurrent sheet that is formed in between the coalescing plasmoids: thisis the region of strong negative current between the two plasmoids inpanels (a), (b), (c) and (d) of Fig. 2. In a first phase, the current sheetlength rapidly increases when the plasmoids approach, and then itprogressively reduces with the size of the coalescing plasmoids [panels(c) and (d) of Fig. 2]. The reconnection results in the formation of asingle large plasmoid, as shown in panels (e) and (f) of Fig. 2.The left–right symmetry in the PIP case is broken during the coa-lescence, as evident, in particular, from panels (d) and (f) of Fig. 2. Theasymmetry arises from the initial noise perturbation in Eq. (22), allow-ing the onset of small-scale dynamics in the central current sheet. Thesymmetry in the MHD case is also reinforced from the presence of acentral boundary at x ¼
0, as introduced in Sec. II B.Figure 3 displays in blue the separation of the two plasmoids inthe MHD case, calculated as the distance between O (cid:5) points. Thesquares along the curve identify the times of panels (a), (c), and (e) in Physics of Plasmas
ARTICLE scitation.org/journal/php
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.0032236 , 032901-4Published under license by AIP Publishing ig. 2. This distance fluctuates in time, with peaks that appear regularlyduring the reconnection phase. The merging plasmoids acceleratetoward each other, move slightly apart as they bounce on the currentsheet, and accelerate back again toward the center. Such movement can be associated with the high gas pressure generated inside the cur-rent sheet.As observed from the black curve in Fig. 3, the distance betweenthe two plasmoids reduces rapidly in the PIP case, following the fasterreconnection process and without displaying the same oscillations thatare remarkable in the MHD case. The reason could be that the recon-nection is happening fast enough that the plasmoids do not need torebound off each other. The squares here identify the times of panels(b), (d), and (f) in Fig. 2.The shortening of the coalescence time scale in the PIP case canbe associated with the decoupling of ions and neutrals in the reconnec-tion region. This results in a faster thinning of the current sheet as thelower ion density allows a stronger compression than in the MHDcase. The decoupling is discussed in more detail in Sec. III A. Thesharpening of the magnetic field profile has already been observed inmany studies as a result of the ambipolar diffusion. The reducedsingle-fluid model for two-fluid effects provided by the ambipolar dif-fusion provides a good approximation for strongly coupled systems,and it is consistent with the effects that we record in our simulations.However, it might fail in describing the complexity of weakly andintermediate coupled systems, and could not be used to explain thedifferent time scale of the first phase of coalescence when the plas-moids move toward each other. The increased complexity of our sys-tem is, therefore, better investigated through a full two-fluid modelsuch as the one used in this work. FIG. 2.
Comparison of J z between the MHD case (left column) and the reference PIP case (right column). The frames identify different steps of the coalescence instability.Panels (a) and (b) show the initiation of the reconnection process. In panels (c) and (d) the evolution of coalescence is displayed at later stages. The final stage of coalescenceis shown in panels (e) and (f), with the formation of the resulting plasmoid. The same magnetic field lines are displayed in black for all frames. Times are given in the samenon-dimensional unit. FIG. 3.
Time variation of the distance between the merging plasmoids, calculated asthe distance between the O (cid:5) points for the MHD case (solid blue) and the PIP case(solid black). The blue squares refer to the times shown in panels (a), (c), and (e) ofFig. 2. The black squares identify the times shown in panels (b), (d), and (f) of Fig. 2. Physics of Plasmas
ARTICLE scitation.org/journal/php
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.0032236 , 032901-5Published under license by AIP Publishing t a first qualitative view, several differences are present betweenthe fully ionized and the partially ionized cases. First, the plasmoidmerging occurs faster in the PIP simulation, ending at t (cid:3)
4, whilecoalescence takes a longer time in the MHD case, ending at t (cid:3) : A. Ideal phase of coalescence
The mutual attraction of parallel currents pulls the two initial plas-moids together, and a current sheet forms as a result of the anti-parallelmagnetic field being pushed together. In the MHD simulation, the idealphase of coalescence is characterized by a plasma inflow forming alongthe x -axis that contributes to the formation of the central current sheet.The charged species in the PIP simulation are pulled to the center of thedomain by the Lorentz force, and the neutrals are dragged by collisionsresulting in a small drift velocity that can be seen in Fig. 4 for t ¼ v n (cid:5) v p ) differs from zero in the inflow, which indi-cates that the two species are weakly coupled in this first phase of coalescence, and increases steadily in time with the acceleration of theplasmoids motion toward each other. After t ¼ g J (shown in Fig. 5) has an effect on the neutral temperature, which FIG. 4.
Plots of the velocity (top left), pressure (top right), density (bottom left) and temperature (bottom right) of ions (blue) and neutrals (red) at t ¼ t ¼ FIG. 5.
Frictional heating at t ¼ t ¼ x -axis in the PIP simulation, compared to the Ohmic heating (reddashed line for t ¼ t ¼ Physics of Plasmas
ARTICLE scitation.org/journal/php
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.0032236 , 032901-6Published under license by AIP Publishing ncreases due to the thermal coupling between the species. As the neu-trals are not completely coupled to the ions, they are expelled from thecurrent sheet, as shown by the drop in the neutral density (bottom leftpanel of Fig. 4).Both fluids are also heated up in the inflow region through fric-tional heating, the non-dimensional definition of which is ð = Þ a c ð T n ; T p ; v D Þ q n q p ð v n (cid:5) v p Þ : (23)The frictional heating at y ¼ t ¼ t ¼ B. Current sheet and reconnection rate
Once the current sheet is generated, the MHD and PIP casesshow very different reconnection processes. Laminar reconnectiontakes place in the MHD case, independent of the initial white noiseperturbation without the onset of further instabilities. The long thincurrent sheet can be compared to the steady-state Sweet–Parkermodel. Both the length D MHD and the width d MHD are estimated bytaking the full width at 1/8 of the maximum current density J z , respec-tively, along the y -axis and the x -axis. We choose this ratio as it accu-rately represents the termination of the reconnection region. Thecurrent sheet width and length are d MHD (cid:3) :
02 and D MHD (cid:3) : D MHD as characteristic length of the system and the maximumvalue of the Alfv (cid:2) en speed that occurs at the boundary of the currentsheet, v A (cid:3) :
26, it is possible to calculate the Lundquist number. Wefind that S ¼ D MHD v A = g (cid:3) : (cid:6) .In the case of Sweet–Parker-like reconnection, the current sheetaspect ratio scales as S (cid:5) = . From the value of S we obtain that theexpected aspect ratio for the MHD case is d = D (cid:3) : (cid:2) (cid:5) , anestimate comparable to the measure obtained by d MHD = D MHD (cid:3) : (cid:2) (cid:5) . Having laminar reconnection in a longthin current sheet which does not develop instabilities nor break intosmaller parts, we may suggest that the MHD case is subject to a recon-nection process that is Sweet–Parker-like.Unlike the MHD case, reconnection in the PIP case is nonlami-nar. The presence of plasmoids breaks the current sheet into multiplethinner current sheets, which reconnect faster than the original struc-ture. The dimensions of the current sheet at the time immediatelybefore the generation of the first plasmoid ( t ¼ d PIP (cid:3) : D PIP (cid:3) :
46 in length, estimated as J z ; max = x -axis and the y -axis, respectively. The onset of the plasmoid insta-bility in a fully ionized plasma might take place below a criticalaspect ratio of 1 = ¼ (cid:2) (cid:5) : the same result was obtainedfor a multi-fluid plasma. Our current sheet aspect ratio is (cid:3) : (cid:2) (cid:5) , about four times larger than the predicted aspectratio implying further physics may be involved in the onset of plas-moid formation. This is investigated in Sec. IV E.In a partially ionized plasma, three different Alfv (cid:2) en speeds can beidentified: a total Alfv (cid:2) en speed v A ; t related to the total density q n þ q p , an ion Alfv (cid:2) en speed v A ; p , based on the density of the chargedparticles, and an effective Alfv (cid:2) en speed v A ; e , based on the combineddensity of charged particles and neutrals that are coupled through col-lisions. The expression for v A ; e might be non-trivial; however, it is possible to provide a close estimate for it from the plasma outflowvelocity ( v out (cid:3) v A ; e ). We chose to use the plasma velocity as the ion-ized fluid is the one directly accelerated by the reconnected magneticfield lines. The Alfv (cid:2) en speed is inversely proportional to the density: as v A increases at the decrease in density, the ion Alfv (cid:2) en speed is biggerthan the total Alfv (cid:2) en speed. The Lundquist number calculated byusing v A ; t (cid:3) :
03 is S t ¼ v A ; t D PIP = g (cid:3) : (cid:2) . If we consider theion Alfv (cid:2) en speed only, which is v A ; p (cid:3) :
29, the Lundquist numberbecomes S p (cid:3) : (cid:2) , which is consistent with the thresholdvalue S ¼ (cid:2) for the onset of the tearing instability and plasmoidformation. In the presence of collisional coupling, reconnection scales withthe Lundquist number associated with the effective Alfv (cid:2) en speed v A ; e .Estimating v A ; e (cid:3) :
09 from the maximum outflow velocity, we find S e (cid:3) : (cid:2) , which is below the threshold value for the onset ofthe tearing instability. The discrepancy suggests that effects due to par-tial ionization might affect the dynamics of reconnection, allowing theformation of secondary plasmoids in the presence of a lowerLundquist number. The answer to this discrepancy between the mod-els can be sought in the modifications due to the ion–neutralinteraction.The reconnection rate M is defined by M ¼ g J max v A B up ; (24)where J max is the absolute maximum value of the current density insidethe current sheet, v A is the initial maximum bulk Alfv (cid:2) en speed ( v A ; t inthe PIP case), and B up is the initial maximum value of B y in theinflow. The MHD mean reconnection rate is M ¼ : : t ¼ ½ : ; (cid:7) for the PIP reference case. Both are FIG. 6.
Temporal evolution of the drift velocity magnitude, compared to J z at thecenter of the current sheet (in red) for the PIP case. The solid black lines indicatethe maximum values of drift velocity and the solid blue lines indicate the medianvalues. The vertical dashed lines indicate the times displayed in panels (b), (d), and(f) of Fig. 2. Physics of Plasmas
ARTICLE scitation.org/journal/php
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.0032236 , 032901-7Published under license by AIP Publishing ompared to the evolution of current density at the center of the cur-rent sheet ( x ¼ ; y ¼ (cid:3) :
1, increasing smoothlyin the last phases of the coalescence after t ¼ J z decreases. The peak in the median drift velocity is reached at thecomplete merging, where it reaches a value j v D j (cid:3) x directionduring its collapse, and it takes place during the formation and expul-sion of secondary plasmoids (for more details on secondary plasmoidssee Sec. III D). This can be seen from the peak in the maximum driftvelocity that is reached in the interval between the two central verticaldashed lines. These vertical lines represent the times identified in pan-els (b) and (d) of Fig. 2. The drop in the maximum j v D j occurringapproximately between t ¼ t ¼ J z . Such smooth variation of J z is linkedto the formation and growth of the first secondary plasmoid ant thecenter of the current sheet. More details about the investigation of thecurrent density are presented in Sec. IV B. C. Shocks
During reconnection and in the final phase of the merging in theMHD case [panels (c) and (e) of Fig. 2], there is evidence of shocks vis-ible as thin elongated lines corresponding to both positive and negativepeaks of the current density magnitude. The structures that can be dis-tinguished in the current density are enhanced in the divergence of v p ,shown in Fig. 7, where they are identified as thin red lines.The minimum in the divergence of the plasma velocity field iden-tifies a region in which the flow is highly compressed, i.e., a shock.Across the shock, the magnetic field components B y and B z drop, whileplasma density and pressure rise steeply. The behavior of magneticfield and pressure identifies this as a slow-shock. The presence of slow-mode shocks is expected as they are a common feature in recon-necting systems, being part of a huge variety of fine structures that canbe identified when plasmoid dynamics takes place. Comparing panels (c) and (d) of Fig. 2, the PIP case shows farfewer clear shock structures. Here we analyze the mechanisms sup-pressing slow-mode shocks in the PIP simulations. We examine thedivergence of the plasma and neutral velocity fields (Fig. 8) at t ¼ r (cid:6) v are a signature of shocks.Comparing it to the divergence of the velocity in the MHD case (Fig.7), there are no structures that can be associated with slow-modeshocks. We present more information later on in Sec. IV B.However, a wide range of structures appear in both r (cid:6) v n and r (cid:6) v p , and they are particularly enhanced in the neutrals. In theplasma velocity divergence, the most prominent structure is associatedwith the neutral jet discussed in Sec. III E, but other structures cut the x -axis symmetrically at both sides of the reconnection region. Thesestructures, which form in the neutrals and later couple to the plasma,are hydrodynamic shocks generated by the motion of the neutral spe-cies in the inflow and not slow-mode shocks as found in the MHDcase. FIG. 7.
Divergence of the plasma velocity at t ¼ FIG. 8.
Divergence of the plasma velocity (top) and the neutral velocity (bottom) at t ¼ Physics of Plasmas
ARTICLE scitation.org/journal/php
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.0032236 , 032901-8Published under license by AIP Publishing uring the merging, the neutrals are expelled and travel awayfrom the reconnection region with a flow that is denser in the directionperpendicular to the current sheet. This can be seen from panel (a) ofFig. 14. In their motion, the neutrals interact with the dense plasmaflow [panel (b) of Fig. 14], and they are halted by the collisions. Thecompression of the neutral flow leads to the formation of multipleshock fronts that are perpendicular to the x -axis. The neutral shockscoupling to the plasma manifest as the lines in r (cid:6) v p . The lines, whosefront moves away from the current sheet center, are present in theneutral v x component, pressure and density but do not display a coun-terpart in the plasma variables.Other hydrodynamic shocks are visible along the y -axis. Theseshocks are formed by the material accelerated inside a neutral jet,which will be examined in Sec. III E. D. Secondary plasmoids
During coalescence in the PIP case the central current sheet issubject to the tearing instability, and secondary plasmoids are pro-duced as evident in panel (b) of Fig. 2. Figure 9 shows three secondaryplasmoids forming, moving along the current sheet and beingexpelled. The motion along the current sheet is triggered by the whitenoise perturbation that breaks the symmetry of the system.In the framework of plasmoid dynamics, it is interesting to inves-tigate whether the secondary plasmoids have any characteristic incommon with the initial plasmoids. We look at the force balancebetween the total pressure gradient and the Lorentz force( J (cid:2) B (cid:5) r p ), shown in Fig. 10 along the y -axis at t ¼ t ¼ t ¼ y ¼½(cid:5) : ; : (cid:7) and moves to the right at t ¼ y ¼ ½(cid:5) : ; (cid:5) : (cid:7) at the later time. The force components canceleach other at the plasmoids location, while the current sheet around isstill out of balance. Inside the plasmoids, the major contributions tothe total force are provided by the magnetic pressure B = y component of the magnetic tension ð B (cid:6) rÞ (cid:6) B (red), whilethe gradient of the neutral pressure (green) is negligible across thewhole region. This suggests that the secondary plasmoids are in analmost force-free condition other than a small region at their centerwhere the plasma pressure ( (cid:5)r p p is shown in blue in Fig. 10)becomes significant.The thermal coupling between ions and neutrals, whose termsare displayed on the right hand side of Eqs. (3) and (8), however, con-tributes to change the plasma pressure gradient with time. The effectof the thermal coupling is shown in Fig. 11, where the plasma pressure,the plasma temperature, and the two terms associated with the cou-pling with the neutrals (frictional heating and thermal damping) aredisplayed at the center of the secondary plasmoid 1 in the time interval t ¼ ½ : ; : (cid:7) . The thermal damping term, whose non-dimensionaldefinition is ð = c Þ a c ð T n ; T p ; v D Þ q n q p ð T n (cid:5) T p Þ , drives the thermalequilibrium. When negative, the thermal damping indicates thatenergy is transferred from the hotter plasma to the neutrals. Theplasma pressure decreases, together with the plasma temperature,under the effect of the thermal damping that reaches more negativevalues in time, a trend that is associated with energy passing from theplasma to the neutrals. The trend shown by the thermal dampingreflects the neutral temperature, which tends to the plasma tempera-ture until t ¼ FIG. 9.
Formation and expulsion of plasmoids from the central current sheet at t ¼ FIG. 10.
Force balance J (cid:2) B (cid:5) r p (black solid line) calculated along the currentsheet in the y -axis at t ¼ t ¼ (cid:5)r p p (blue), (cid:5)r p n (green), the magnetic pressure (magenta), and the magnetictension (red). Physics of Plasmas
ARTICLE scitation.org/journal/php
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.0032236 , 032901-9Published under license by AIP Publishing s also seen to increase remarkably with respect to time, as a result ofthe combined effect of the neutral pressure increasing from the inflowand the plasmoid motion. The coupling with the neutrals acts on r p p by sharpening its peak at the plasmoid center, but the gradient contin-ues to be a relevant contribution to the total force. For this reason, thesecondary plasmoids do not become completely force-free before theyare expelled from the current sheet.The detail of one of the secondary plasmoids reconnecting at oneend of the current sheet is shown in Fig. 12, where J z , the drift velocitymagnitude, and the plasma and neutral velocity components are dis-played. The plasma flow is faster than the neutrals, and this leads to anon-negligible drift velocity between the two fluids. Looking at boththe current density and the drift velocity magnitude, a thin elongatedvertical structure is observed from the current sheet to the center ofthe secondary plasmoid. This structure, visible in red in the bottomright panel of Fig. 12 (neutral v y ), is a jet in the neutral flow, whichextends in the direction opposite to the plasmoid motion, going backto the current sheet. This jet, accelerated by the neutral pressure gradi-ent inside the plasmoid, is present only in correspondence of second-ary plasmoids: there is no similar structure forming in the biggercoalescing plasmoids. The absence of this feature might depend on thefact that the bigger plasmoids are initially in a force-free condition andthe pressure gradients are too small to expel the neutrals through a jet. E. Extended neutral reconnection jet
A prominent feature developing in the PIP simulation is the for-mation of a jet-like structure that extends asymmetrically along the y -axis during coalescence. This large structure must be distinguishedfrom the small-scale neutral jets discussed in Sec. III D. In standard FIG. 11.
Evolution of pressure, temperature, frictional heating and thermal damping at the center of secondary plasmoid 1 in the time interval t ¼ ½ : ; : (cid:7) . In the top panelsplasma properties are displayed in blue while neutral properties are shown in red. FIG. 12.
Detail of the secondary plasmoid leaving the current sheet at t ¼ J z (top left) and drift velocity magnitude (top right). Theplasma velocity components are displayed in the central panels ( v x ; p on the left, v y ; p on the right), while the neutral velocity components are shown in the bottompanels ( v x ; n on the left, v y ; n on the right). Physics of Plasmas
ARTICLE scitation.org/journal/php
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.003223628
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.003223628 , 032901-10Published under license by AIP Publishing econnection models, magnetic energy can be released to form aplasma jet, a feature that is also found in many observations.From Fig. 13, however, we can clearly see that the neutral jet issignificantly longer than the plasma jet. The ion velocity increases tosupersonic values along the reconnection region, but the enhancementis localized near the center of the domain. The velocity of the extendedneutral jet is supersonic, and the neutral Mach number reaches valuesof (cid:3) :
6. The larger drift velocity in Fig. 13 shows that the species aresignificantly decoupled in the jet.A more detailed picture of the interaction of the two fluids alongthe jet structure can be provided by looking at the physical propertiesin the smaller region where the jet develops. Such region is identifiedin the top panel of Fig. 13. Plots of neutral and plasma densities andtemperatures, frictional heating, and temperature difference betweenthe two species are shown in Fig. 14.The decoupling of neutrals and plasma along the jet is favored bythe very low density of both species, as shown in panels (a) and (b) ofFig. 14. The two species reach similar peaks in temperature, as shownin panels (d) and (e) of Fig. 14, but the heating distribution is differentfor each fluid. The species are thermally decoupled, as shown by thedifference between the neutral and ion temperatures in panel (f) ofFig. 14. During its evolution, the jet appears to be very turbulent.There is presence of many coherent vortices mostly concentrated at the jet truncation that are particularly evident in panels (c) and (f) ofFig. 14. Along the jet the ion temperature is the highest and reaches itsmaximum in correspondence of the current sheet, while the neutraltemperature is higher than the ion temperature at the center of thevortices. Neutrals and ions are, however, heated up in the current sheetand along the jet by the thermal energy. The thermal energy is releasedthrough the frictional heating, defined in Eq. (23), which is associatedwith collisions between the two fluids, and it is shown in panel (c) ofFig. 14.The velocity difference at the interfaces between the jet and theenvironment leads to the onset of shear flow instabilities. The sinusoi-dal shape of the jet is characteristic of the Kelvin–Helmholtz instability(KHI), a classical shear flow instability that tears apart vorticity sheetsat the surface of separation of the two fluids.
In order to confirmwhether the system is KH unstable, we compare the neutral jet to thesimple Bickley jet, a steady two-dimensional laminar jet which isunstable to the sinusoidal-mode of the KHI. Under the action of theKHI, the Bickley jet develops a sinusoidal structure at a preferredwavelength of (cid:3) : k can be calculated at t ¼ y ¼ k (cid:3) : d (cid:3) : k = d of6.6, which is similar to the value predicted for the Bickley jet undergo-ing KHI. We can conclude that the jet is subject to the KHI.For our jet, the KHI is seen to evolve to a turbulent state. At thetermination point, in this location shocks are generated and can beseen as weak structures in the neutral and drift velocities in Fig. 13(top and bottom panels). These structures are the hydrodynamicshocks discussed in Sec. III C. IV. PARAMETER SURVEY
We investigate the changes in the coalescence process due to thediffusivity ( g ), the collisional coupling ( a c ), the ion fraction ( n p ) andthe plasma b . In this section, we present a survey over these four keyparameters of our physical system. The simulations are identified bynumbers, and the respective physical parameters and the spatial reso-lution are listed in Table I.The PIP simulations in this section have the same resolution asthe MHD cases or an even lower resolution( D x ¼ : (cid:2) (cid:5) ; D y ¼ : (cid:2) (cid:5) ). Due to the parameter variationchanging the size of the central current sheet, it was possible to use alower resolution without losing the possibility to resolve the currentsheet. A. Variation of resistivity
We begin the investigation of partial ionization effects on thecoalescence instability by considering the role of varying the resistivityin PIP simulations. The seven cases that are examined in this section( g ¼ J z at the center ofthe current sheet ( x ¼ y ¼
0) is displayed in Fig. 15, and it is seen todecrease as result of the increasing diffusion. The chosen locationallows us to identify the beginning and the end of reconnection withthe formation of the final plasmoid.
FIG. 13.
Top: drift velocity magnitude j v D j between neutral and charged fluids.Center: ion Mach number v p = c s ; p . The magnetic field lines are shown in black.Bottom: neutral Mach number v n = c s ; n . All the plots display the quantities distribu-tion at t ¼ Physics of Plasmas
ARTICLE scitation.org/journal/php
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.0032236 , 032901-11Published under license by AIP Publishing or g ¼ :
005 (in black in Fig. 15) reconnection happens inthe central current sheet with the formation of some secondaryplasmoids, although fewer with respect to the less diffusive cases( g ¼ : g ¼ : j J z j as soon as it moves away from thecenter of the current sheet. In the simulation with g ¼ :
015 thecurrent is highly diffused. During the plasmoid merger there is noobserved formation of secondary plasmoids. This is the first case inwhich the diffusion is high enough to prevent the formation of sec-ondary plasmoids during coalescence in a PIP case. A similar resultis obtained for the case with g ¼ :
05. The cases having the highestresistivity ( g ¼ : ; :
5) do not develop the coalescence instabil-ity, as the initial plasmoids are quickly diffused. As g is increased,the Lundquist number S decreases: when S becomes sufficiently low, the tearing instability stops taking place in the central currentsheet, thus explaining the lack of secondary plasmoids.The variation of resistivity also plays a role in changing the timescale of coalescence. For higher g , coalescence starts at a later time ineach simulation. This is shown in Fig. 15 by the position in time of thefirst drops in J z . Such effect depends on how efficiently the currentsheet is generated in the phase where the two plasmoids are approach-ing, as for smaller g it is easier to buildup current. B. Variation of collisional frequency
Here we investigate the effects of the collisional coupling betweenions and neutrals. We compare the simulations listed in Table I withthe numbers 3, 9, 10, 11, 12, 13, and 14. An MHD case with total den-sity equal to q p ¼ a c ! 1 .A second MHD case (simulation 9) having initial density equal to FIG. 14.
Profiles of neutral (panel a) andion density (panel b), frictional heating(panel c), neutral (panel d), and plasmatemperature (panel e) and differencebetween neutral and ion temperatures(panel f ) are displayed at t ¼ b , the separatrices are shown inblack. The selected region is shown inFig. 13. Physics of Plasmas
ARTICLE scitation.org/journal/php
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.003223628
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.003223628 , 032901-12Published under license by AIP Publishing p ¼ :
01 is considered as the limit for a c ¼ a c ¼ 1 theneutral and plasma species are entirely coupled hence the systembehaves like a single fluid MHD model, with the density and pressure being the bulk (ion þ neutral) values. For a c ¼
0, the species arecompletely decoupled, i.e., the plasma evolves independently from theneutrals and can be considered to be a single fluid MHD system withthe density/pressure based on the plasma values only. An equivalentMHD simulation can be performed by changing the initial plasmabeta. While the magnetic field strength is unchanged, the variation ofplasma beta modifies the pressure. To maintain the same initial tem-perature of the calculation, we use a lower density (the same as theplasma density of our two-fluid calculations). Therefore, the differencein the plasma density between the two limit MHD cases result in aneffective difference in the plasma b of the two simulations.In Fig. 16, we display the variation of current density J z at x ; y ¼ a c the time scale for the plasmoid coalescence becomesconsiderably shorter. The time scale does not vary linearly, but itshows the presence of two accumulation points, corresponding to thetwo limits identified by the MHD simulations. At lower a c ( a c ¼ ; a c ! a c ( a c ¼ ; a c ! 1 .The more ions and neutrals are coupled, the later reconnectionstarts. The slowing down of the first phase prior to the onset of recon-nection can be associated with the damping effect that neutrals haveon the ions in the inflow, which increases at the increase in a c as theions (which fraction is much smaller than the neutral fraction) interact TABLE I.
List of the simulation parameters.
Nr. Type H a c n p b No. x grid points No. y grid points D x D y a a : (cid:2) (cid:5) : (cid:2) (cid:5) : (cid:2) (cid:5) : (cid:2) (cid:5) : (cid:2) (cid:5) : (cid:2) (cid:5) : (cid:2) (cid:5) : (cid:2) (cid:5) : (cid:2) (cid:5) : (cid:2) (cid:5) : (cid:2) (cid:5) : (cid:2) (cid:5) : (cid:2) (cid:5) : (cid:2) (cid:5) : (cid:2) (cid:5) : (cid:2) (cid:5) a a : (cid:2) (cid:5) : (cid:2) (cid:5)
10 PIP 0.0015 1 0.01 0.1 2062 3086 1 : (cid:2) (cid:5) : (cid:2) (cid:5)
11 PIP 0.0015 10 0.01 0.1 2062 3086 1 : (cid:2) (cid:5) : (cid:2) (cid:5)
12 PIP 0.0015 1000 0.01 0.1 2062 3086 1 : (cid:2) (cid:5) : (cid:2) (cid:5)
13 PIP 0.0015 3000 0.01 0.1 2062 3086 1 : (cid:2) (cid:5) : (cid:2) (cid:5)
14 MHD 0.0015 a a : (cid:2) (cid:5) : (cid:2) (cid:5)
15 PIP 0.0015 100 0.5 0.1 2062 3086 1 : (cid:2) (cid:5) : (cid:2) (cid:5)
16 PIP 0.0015 100 0.1 0.1 2062 3086 1 : (cid:2) (cid:5) : (cid:2) (cid:5)
17 PIP 0.0015 100 0.001 0.1 2062 3086 1 : (cid:2) (cid:5) : (cid:2) (cid:5)
18 PIP 0.0015 100 0.01 1 1038 1550 3 : (cid:2) (cid:5) : (cid:2) (cid:5)
19 PIP 0.0015 100 0.01 0.01 1038 1550 3 : (cid:2) (cid:5) : (cid:2) (cid:5) a These data are the effective values of the two-fluid parameters a c and n p for the single-fluid cases, which are chosen as limits for the PIP simulations. FIG. 15.
Evolution over time of the current density J z at the central point of the cur-rent sheet formed during coalescence ð x ¼ ; y ¼ Þ at the variation of the initialresistivity. The cases displayed are for g ¼ : g ¼ : g ¼ :
005 (black), g ¼ :
015 (blue), g ¼ :
05 (brown), g ¼ :
15 (magenta), and g ¼ : Physics of Plasmas
ARTICLE scitation.org/journal/php
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.003223628
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.003223628 , 032901-13Published under license by AIP Publishing ith a higher number of neutrals. This result is consistent with previ-ous simulations. We can understand this result by comparing theion Alfv (cid:2) en time s A ; p with both the ion and the neutral collisionaltimes. This leads to the identification of two values for a c that definewhen the species couple with each other through collisions. For a c (cid:3)
20, the plasma–neutral collisional time s col ; pn ¼ ð a c q n Þ (cid:5) becomes smaller than s A ; p and the ions couple with the neutrals. Thecoupling of the neutrals on the ions, that happens when the neutral–-plasma collisional time ð a c q p Þ (cid:5) becomes smaller than s A ; p , takesplace for a c (cid:3) a c increases, reconnection takes a longer time before the plas-moids are completely merged. This results in a decrease in the averagereconnection rate as a c increases, as shown in the top panel of Fig. 17,where the error bar is given by the standard deviation. The reducedreconnection rate can be explained by the variation of the effectiveAlfv (cid:2) en speed, as introduced for the two-fluid case in Sec. III B. If thecollisional frequency is increased, the effective density of the magneticfluid increases as more neutrals interact with the ions and conse-quently v A ; e decreases. The variation of the Lundquist number isshown in the central panel of Fig. 17, as calculated from v A ; e .After reconnection begins, the PIP cases (Fig. 16) show a strongfluctuation of the current toward less negative values, behavior, that is,not present for a full coupling (MHD case for a c ! 1 ). This fluctua-tion is a consequence of plasmoid formation in the current sheet. Thegrowth and expulsion of plasmoids contributes initially to slow downthe reconnection by saturating the negative current in a blob betweenthe merging plasmoids and then to leave the current sheet unstable,leading to the formation of further plasmoids. For the simulationswith secondary plasmoids, S is calculated at the time before the forma-tion of the first plasmoid in the central current sheet, while in the caseof simulations without secondary plasmoid S is evaluated at the firstminimum of J z at the center of the current sheet. The value of theLundquist number decreases at the increase in the collisional coupling,following the variation of both the effective Alfv (cid:2) en speed and the char-acteristic length of the current sheet. S is more sensitive to the variationof v A ; e , as such factor varies over an interval covering orders of magnitude, while D , being about the same size of the plasmoids, dis-plays a lesser variation in length across the cases. As shown in Fig. 17, S appears to be above the threshold of 10 for all the cases showing for-mation of secondary plasmoids, with the exception of the case with a c ¼ a c ¼ Therefore, wewant to investigate whether slow-mode shocks are produced at ahigher or lower collisional coupling that approach the MHD cases,and whether they form but are dissipated at a later time of coalescence.Slow-mode shocks are indeed generated in the two-fluid simulations,and their propagation is clearly visible particularly at low a c , as shownin the first two panels (simulations 10 and 11) of Fig. 18. At theincrease in collisional coupling ( a c (cid:8) FIG. 17.
Top: average reconnection rate as a function of a c . The error bar is associ-ated with the standard deviation on each measure. Center: Lundquist number S with respect to a c calculated using the effective Alfv (cid:2) en speed v A ; e . The horizontaldashed line indicates the threshold limit of S ¼ for the onset of the tearinginstability. Bottom: current sheet aspect ratio d = D as a function of a c . In both cen-tral and bottom panels, diamonds are associated with simulations with secondaryplasmoids, crosses are associated with cases without secondary plasmoids. FIG. 16.
Evolution over time of the current density J z at the central point of the cur-rent sheet formed during coalescence ð x ¼ ; y ¼ Þ at the variation of the initialcollisional coupling. The cases displayed are the MHD limit for a ¼ a ¼ a ¼
10 (blue), a ¼ a ¼ a ¼ a ! 1 (black, solid line). Physics of Plasmas
ARTICLE scitation.org/journal/php
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.003223628
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.003223628 , 032901-14Published under license by AIP Publishing requencies ( a c ¼ ; a c (cid:9)
1, increases with the coupling between the two speciesaround the inflow region as the two species that move in oppositedirections are subject to an increasing interaction. For a c ¼ a c ¼ C. Variation of ion fraction
In this section, we compare five cases with a different ion fraction n p in the range 10 (cid:5) (cid:5) (cid:2) (cid:5) , corresponding to the numbers 3, 14,15, 16, and 17 listed in Table I.The variation in time of J z at the center of the current sheet is dis-played in Fig. 19, with the speed in the process of coalescence drasti-cally increased at the decrease in n p . Such behavior, which shows avariation in the time scale of both ideal phase (when the plasmoidsattract each other and the current sheet is formed) and reconnectionphase, might be explained by similar arguments to those presented inSec. IV B.At the variation of the ion fraction, the initial ion Alfv (cid:2) en timescale s A ; p and the ion collisional time scale increase with n p , while theneutral collisional time scales decreases. We compare s A ; p with the ionand the neutral collisional times to find the values of n p for which thetwo species couple with each other. In our range, the ion collisionaltime ð a c q n Þ (cid:5) is always smaller than s A ; p with the only exception ofthe MHD case (simulation 14) which is taken into account as the limitvalue for n p ¼
1. Therefore, the ions are always coupled to the neutralsin all the PIP simulations in this section. The value of n p below whichthe ion dynamics becomes fast enough to decouple from the neutralsis (cid:3) (cid:2) (cid:5) . Conversely, the neutrals coupling on the ions takesplace for n p ¼ : n p ¼ : ; : n p the effectiveAlfv (cid:2) en speed v A ; e decreases, following the increase in theplasma density. This affects the Lundquist number (central panelof Fig. 20), which in turn extends the time scale of the reconnec-tion phase. Sub-critical plasmoid formation, discussed in Sec.IV E, is observed for the case at the lowest ion fraction( n p ¼ : S below the threshold limit of 10 at the onsetof the tearing instability. D. Variation of plasma beta
In this subsection, three cases at different plasma b (simulations3, 18 and 19 in Table I) are investigated. The speed of coalescence isgreatly affected by the change of b , and increases at the decrease in FIG. 18.
Divergence of the plasma velocity for the PIP cases at a ¼ ; ; FIG. 19.
Evolution over time of the current density J z at the central point of the cur-rent sheet formed during coalescence ð x ¼ ; y ¼ Þ at the variation of the initialion fraction. The cases displayed are for the MHD case ( n p ¼
1, black) n p ¼ : n p ¼ : n p ¼ :
01 (green), and n p ¼ :
001 (red).
Physics of Plasmas
ARTICLE scitation.org/journal/php
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.003223628
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.003223628 , 032901-15Published under license by AIP Publishing his parameter as shown in Fig. 21. This is explained by the fact that atsmaller b the magnetic field get stronger with respect to the plasmapressure: the Alfv (cid:2) en speed increases (as its numerator increases), andso does the Lundquist number.The consequences of plasma b variations are not directly associ-ated with a two-fluid effect, as the different dynamics dependsuniquely on the variation of B and can be purely reduced to an MHDeffect. However, it provides important context for how changing dif-ferent parameters in the two-fluid case, which effectively alter v A , canincrease the merger rate.The presence of sub-critical plasmoid formation is observed forthe case with the largest plasma b ( b ¼ . The physics behind this case is discussed in Sec. IV E. E. Sub-critical plasmoid formation
Several PIP simulations presented in Secs. IV A–IV D (identifiedin Table I by number 2, 4, 12, 13, 17, and 18) show the signs of second-ary plasmoids formation in the central current sheet below theLundquist number threshold S ¼ . The values of the critical S found in these cases are consistent with the results obtained by recentstudies of magnetic reconnection in a multi-fluid partially ionizedplasma, already discussed in Sec. I. One might expect that the lowercritical Lundquist number in our simulations could depend on the ini-tial setup, and, in particular, the low ion fraction. In terms of the initialsetup, a comparison with previous works is difficult: many studies are performed with a static current sheet setup, while our simulationspresent a driven reconnection, which leads to a different current sheetstructure and dynamics. A more detailed evaluation of the role playedby the initial perturbation and the random white noise in Eq. (22) onthe onset of the tearing instability must be investigated in future devel-opments of this research. Several studies already found that the role ofthe amplitude of perturbation noises is major in determiningthe critical Lundquist number and current sheet aspect ratio in a rangeof initial configurations. However, tests performed on the variation ofthe white noise perturbations on our simulations proved that the sec-ondary plasmoids are always generated at the same time in every sim-ulation that has the same initial set of parameters.In this section, we often refer to the effective Lundquist numberas the critical Lundquist number for simulations displaying secondaryplasmoid formation. The critical Lundquist number for the onset ofthe tearing instability was not isolated in previous studies, but itwas considered to fall in an interval of values whose limits are set bythe absence (lower boundary) and presence (upper boundary) of plas-moid formation. Such intervals had been determined by varying twomain parameters of the current sheet: characteristic length and resis-tivity. In our study, this interval is determined by changing the charac-teristic length scale of the current sheet. We do not change L directly,but it varies in time as the system evolves. The time interval betweentwo outputs is small enough that the current sheet length does not dis-play a large variation: this reduces to have Lundquist numbers that arevery close to each other before and after the first secondary plasmoid FIG. 20.
Top: average reconnection rate as a function of n p . The error bar is associ-ated with the standard deviation on each measure. Center: Lundquist number S with respect to n p calculated using v A ; e . The horizontal dashed line indicates thethreshold limit of S ¼ for the onset of the tearing instability. Diamonds areassociated with simulations with secondary plasmoids, and crosses are associatedwith cases without secondary plasmoids. Bottom left: current sheet aspect ratio d = D as a function of n p . FIG. 21.
Evolution over time of the current density J z at the central point of the cur-rent sheet formed during coalescence ( x ¼ ; y ¼
0) at the variation of the initialplasma b . The cases displayed are for b ¼ :
01 (red), b ¼ : b ¼ J z was normalized with respectto the case with plasma b ¼ : ffiffiffiffiffiffiffiffiffiffiffi b = : p . Physics of Plasmas
ARTICLE scitation.org/journal/php
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.003223628
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.003223628 , 032901-16Published under license by AIP Publishing ppear, which we identify by the formation of an O -point in the cur-rent sheet magnetic field. For such reason, we choose to identify thecritical Lundquist number as a single value rather than specifying aninterval.We might expect the sub-critical plasmoid formation to betriggered in all the PIP cases because of the small ion fraction.However, the critical Lundquist number obtained for the MHDcase at lower ion density (matching that of the plasma density ofthe PIP simulations), shown in Fig. 17, is well above the thresholdof 10 (limit case for a c ¼
0, listed with number 9 in Table I). Thisproves that the change in the plasma density itself, along with thevariation of the plasma b of the plasma that this change implies,does not affect the system so that it develops sub-critical plas-moids. We also see that for the MHD calculation, that is, equiva-lent to a c ¼ 1 the current sheet is stable for S ¼ : (cid:2) , againimplying a critical S > . At the same time, PIP cases character-ized by an S and an aspect ratio d = D incredibly similar to thisMHD calculation (see cases 12 and 13, where a c is, respectively,10 and 3 (cid:2) , in Fig. 17), show sub-critical plasmoid formation.The inclusion of the coupling between the ion and neutral fluids(effectively looking at systems between those two MHD simula-tions) allows plasmoids to form below S ¼ , even withoutchanging the plasma b , the initial conditions or the perturbationmagnitude. Therefore, such sub-critical plasmoid formation doesnot depend on the characteristics of the plasma itself, but it mightbe ascribed to the combined result of two-fluid effects.From previous studies of the onset of the tearing instability inpartially ionized plasmas, the critical aspect ratio for the initiation ofa growth rate independent of the Lundquist number and neutral toionic density was derived for three regimes (coupled by collisions,intermediate and uncoupled). In the intermediate regime, where ionsand neutral are partly coupled, the critical aspect ratio for a genericequilibrium configuration is d D (cid:3) S (cid:5) = p ð (cid:3) pn s A ; p Þ (cid:5) = ; (25)where the collisional frequency (cid:3) pn is equal to a c q n and s A ; p ¼ D = v A ; p is the ion Alfv (cid:2) en time scale. In general, this critical aspect ratio is big-ger than both the threshold S (cid:5) = obtained from classical argumentsthat involve the Lundquist number and the value of 1/200 found inseveral works. Comparing our current sheet aspect ratios at the time of the onsetof the tearing instability in the PIP simulations with the critical aspectratio, we find that in most of the cases the aspect ratio is still smallerthan the critical value obtained by Eq. (25). Therefore, the sub-criticalplasmoid formation cannot be explained uniquely by applying suchcondition for the onset of the tearing instability.Focusing on a physical explanation for the sub-critical plasmoidformation, we take a close look at what happens in the current sheetbefore the onset of tearing instability. We examine in detail the PIPcase corresponding to simulation 13 in Table I and evaluate how thetwo-fluid parameters vary up to the time immediately before the for-mation of a central plasmoid, which signature appears in the currentdensity and magnetic field at t ¼ a c ð Þ isset equal to 3000, and the Lundquist number before the onset of thetearing instability is S ¼ (cid:2) . The plasma velocity along x is shown in row (a) of Fig. 22. Asenhanced by the contour lines, before the formation of the secondaryplasmoid the plasma in the inflow slows down around the center ofthe current sheet where a stagnation point in the flow exists, forming a FIG. 22.
Evolution of v p ; x (row a), x component of the drift force (row b) and p n (row c) in a PIP simulation with an initial strong coupling between the two species( a c ¼ Physics of Plasmas
ARTICLE scitation.org/journal/php
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.003223628
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.003223628 , 032901-17Published under license by AIP Publishing inching flow that promotes reconnection in two points above andbelow the current sheet central point.In order to understand how the plasma motion is slowed down,we evaluate the force balance in the inflow region. All the force contri-butions are in equilibrium with each other, with the exception of thedrift force, stated on the right hand side of Eqs. (2) and (7), x compo-nent of which for the plasma is shown in row (b) of Fig. 22. The driftforce magnitude appears to be quite large; however, this does notmean that the drift velocity between the two species is also big in theinflow region. The drift force has a strong dependency on the colli-sional term a c , which is itself dependent on both the drift velocity andthe temperatures of the two fluids. Because of the strong dependenceon the a c term, when the species are weakly coupled (and consequentlythe drift velocities are large), the drift force is small.The drift force acts by slowing down the flow as it approaches thecurrent sheet: as its sign depends on the term ð v n (cid:5) v p Þ , its directionsuggests that the plasma is slowed down by the interaction with theneutrals, whose velocity in the inflow is slower than the plasma. Theneutrals are naturally slower as they move toward the center of thecurrent sheet as they are being pulled in by the plasma motions.However, as shown in panel (c), the neutral pressure increases at thecenter slowing the inflow motions and increasing the drag. The neutralpressure increases under the effect of the Ohmic heating g J and ofthe adiabatic compression of the neutrals in the inflow. As already dis-cussed for the reference PIP case in Sec. III A, the Ohmic heating playsa major role in heating the plasma inside the current sheet (see Fig. 5),which results in the increase in the plasma pressure and temperature.As the two species are thermally coupled, the Ohmic heating acts byindirectly heating the neutrals, driving their pressure to increase. Thenon-adiabatic contribution from the Ohmic heating is responsible forthe increase in the plasma temperature inside the current sheet, whichaffect the neutral temperature. The heated neutrals expand in theinflow region, leading to an adiabatic compression directed out of thecurrent sheet, which contributes to increase the neutral pressure.The sub-critical plasmoid formation can therefore be triggered bythe two-fluid interaction between plasma and neutrals. The PIP casesat lower collisional coupling and higher ion fraction, however, developthe formation of critical plasmoids following the linear condition forthe tearing instability before this non-linear mechanism becomesimportant. Therefore, the two-fluid interaction plays a lesser role inthe formation of plasmoids when the effect of neutrals is weak, i.e., atsmaller a c and higher n p . V. DISCUSSION
Among the processes promoting the development of fast mag-netic reconnection, the coalescence instability can play an importantrole by driving the interaction of plasmoids (and their subsequentreconnection) on dynamic time scales. We have investigated how plas-moid coalescence behaves in a partially ionized plasma, a situationreflected in a range of solar atmospheric layers and, in particular, thesolar chromosphere. Through the comparison of a fully ionized case(MHD) and a partially ionized case (PIP), we find that partial ioniza-tion noticeably shortens the coalescence time scale and creates newdynamics, producing neutral jets and secondary plasmoids, suppress-ing slow-mode shocks while promoting hydrodynamic shocks, andleading to sub-critical plasmoid formation. Observations of chromospheric anemone jets performed with theCa II H filtergram of the Solar Optical Telescope onboard Hinode showed the presence of recurrent plasmoids expulsion with a size ofabout a few hundred km from the jets. Assuming an approximatediameter MAX (cid:3)
500 km our characteristic length is (cid:3)
100 km, as inour system the initial plasmoid length is equal to 4 L . Knowing that thesound speed is about 10 km s (cid:5) , we identify a time scale of (cid:3)
10 s.Taking the appropriate ion–neutral collisional cross section for elasticscattering and charge exchange and the characteristic temperature,density and ion fraction of the chromosphere the ion–neutralcollisional frequency is between 10 and 10 s (cid:5) . In non-dimensional form, to compare with our simulations, an equivalent a c to the one observed for the Hinode plasmoids is in the range of10 (cid:5) .A prediction on the behavior at these higher a c can be made bystudying the simulations presented in Sec. IV B, whose trend is shownin Figs. 16 and 17. When a c increases, the evolution tends to the one ofthe upper limiting MHD case ( a c ! 1 ): such case can be consideredas an accumulation point. In the evolution of the current density, allthe simulations at the chromospheric a c would reach the beginning ofreconnection in a time interval between the current minimum ofthe PIP case with a c ¼ (cid:2) and the one of the MHD case, so itis sensible to say that the coalescence evolution and time scale canbe predicted. These values suggest that the coalescence between thebiggest observed plasmoids would take place in a regime, that is,almost MHD. At the specific values for ion fraction, plasma b , andresistivity used in this study, faster coalescence would become rele-vant for plasmoids with a diameter of (cid:3) (cid:5) , we obtain a time scale s ¼ : (cid:3) s (cid:5) , consistent with the collisional fre-quency observed in the chromosphere and with our case at a non-dimensional a c ¼ (cid:2) .This length scale is, however, dependent on the parameters of themedium studied. In many regions of the solar chromosphere the resis-tivity is often smaller, which results in the formation of more plasmoiddynamics, and both n p (equal to 0.01 in the a c section of the parametersurvey) and the plasma b (set equal to 0.1) can be lower in the chro-mosphere, leading to an enhancement of two-fluid effects. At lower n p and b , for the range of a c considered in this study, PIP effects thatinclude a faster coalescence and the sub-critical plasmoid formationwould become important for larger plasmoids and are potentiallyobservable on scales that are currently resolved. The 3 km plasmoidlength scale must therefore be considered as a lower limit for PIPeffects to become observable.Partial ionization plays a role in changing the way many effects,such as the Hall diffusion, develop and act on magnetic reconnection.The Hall effect results from the different velocities of electrons andions, is dependent on the ion fraction and is seen to play a majorrole in both highly ionized and weakly ionized plasmas. In a partiallyionized environment, the physical scales over which the Hall diffusionbecomes important drastically change from the ones found for fullyionized plasmas. Several works show that in magnetic fluxtubes, of which plasmoids are to be considered 2D sections, the Halldiffusion is less important than the contribution from ion–neutral
Physics of Plasmas
ARTICLE scitation.org/journal/php
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.003223628
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.003223628 , 032901-18Published under license by AIP Publishing ffects at the chromospheric heights. As such we viewed that includingthe Hall effect is beyond the scope of this paper.A very important effect in evidence in our results is the sub-critical plasmoid formation discussed in Sec. IV E. Many studies onthe onset of the tearing instability in partially ionized plasmas oftenfocused on linear stability criteria, while in our study plas-moid formation is also promoted by a nonlinear effect linked to thetwo-fluid collisional and thermal coupling. The role played by the cou-pling between ions and neutrals in determining the dynamics of plas-moid formation has already been acknowledged in a previous study looking at the onset of the tearing mode at several levels down to thekinetic scale. In that study the authors only use a linear criterion forthe onset of the tearing instability, while in this work we find that non-linear effects are able to play a major role in plasmoid formation.Including the nonlinear physics may result in new physical parameterregimes that are able to tear down to kinetic scales.The physics that leads to sub-critical plasmoid formation in oursimulations is expected to be largely affected by the non-equilibriumionization–recombination processes. These are currently not includedin our model as the interaction between the two species is provideduniquely by elastic collisions and charge-exchange, but future develop-ments of this research strictly require their inclusion. As pointed outfrom previous studies of magnetic reconnection in a multi-fluid par-tially ionized plasma at low plasma b , the non-equilibriumionization–recombination effect leads to a strong ionization of thematerial in the reconnection region. In a recent paper magneticreconnection has been examined through 2.5D simulations in weaklyionized plasmas with an initial n p ( n p ¼ (cid:5) ) lower than our refer-ence cases. Their results suggest that, for a plasma b larger than 1 andweak reconnection magnetic fields, the non-equilibrium ionization–re-combination effect is responsible for a strong ionization of the neutralfluid in the reconnection region and a faster reconnection rate occur-ring before the onset of plasmoid instabilities. These are consistentwith a previous paper, where an increase by an order of magnitudewas recorded for the ionization degree within the current sheet duringthe reconnection process. The strong ionization is responsible for abigger interaction between the neutral fluid and the plasma, which willbe better coupled both in the inflow and outflow region. These verysame effects are to be expected in the case discussed in Sec. IV E: thedrastic increase in temperature due to the Ohmic heating in the recon-nection region would promote the ionization of the neutral fluid thatin the absence of such process is forced to expand outward, halting theplasma inflow. In the case of a small plasma b smaller than 1 like inour study, however, plasmoid instability is still the main process pro-moting fast magnetic reconnection. The ionization/recombination effects are largely affected by theaction of the ionization potential. When collisional ionization takesplace, the work done against the ionization potential to ionize theatom removes energy from the electrons and acts as a cooling term inthe plasma. As the recombination process is associated with photonsbeing released, this overall effect can be modeled as a radiative loss.Researches investigating the role of radiative cooling in magneticreconnection proved that such process, linked to the collisionalionization, thins the reconnection layer by decreasing the plasma pres-sure and density inside the current sheet. Therefore, the inclusion ofradiative losses speeds up reconnection to rates that are bigger thanthe ones found in models without radiation and might lead to time scales and outflows that are consistent with those found in spiculesand chromospheric jets. In a recent study, it was found that astrong recombination process in the reconnection region, combinedwith Alfv (cid:2) enic outflows, can lead to a fast reconnection rate indepen-dent of Lundquist number. While the decoupling of neutrals andplasma has been recorded in the inflow region, these findings showthat the two fluids are well coupled in the outflows, which is oppositeto what we find for our intermediate coupled case (see Sec. III E). We,therefore, expect that including the radiative losses by adding an ioni-zation potential would lead to a better coupling of the two speciesaround the reconnection region as well as lower temperatures and afaster reconnection rate.Radiation is not only important in terms of the radiative lossesfrom the ionization/recombination processes in the chromosphere.Beyond the action of the ionization potential, effects of ionization andexcitation could be triggered by an external radiation field. This is afundamental factor in the chromosphere, where the ionization degreeis largely determined by the incident external radiation. In a recent study it was demonstrated that the plasmoid cascad-ing process, for which the current sheet breaks into smaller sectionsfollowing the formation of multiple plasmoids, is terminated in theMHD scale. The progressive reduction of secondary plasmoids in thePIP cases as an MHD-like regime is approached is an aspect that hasbeen already marginally observed in our simulations, specifically in theparameter survey linked to the variation of the collisional coupling(Sec. IV B). Multiple plasmoids are seen to form in the simulationshaving an initial lower a c , while at higher collisional coupling only twoand one secondary plasmoids are produced, respectively, for a c ¼ a c are generated by the pinching action ofthe neutrals in the inflow, rather than the onset of instabilities in thecurrent sheet. Therefore, we expect the plasmoid cascading process tobe interrupted in these cases approaching the MHD regime. A moredetailed study needs to be performed on these secondary plasmoidsnumber and characteristics, however, we can already confirm a goodagreement with the trend showed by previous studies. ACKNOWLEDGMENTS
The authors are grateful to Dr. N. Nakamura for theinspiration of his original work on this problem that led to thisstudy. A.H. and B.S. are supported by STFC Research Grant No.ST/R000891/1. A.H. is also supported by his STFC ErnestRutherford Fellowship Grant No. ST/L00397X/1.
DATA AVAILABILITY
The data that support the findings of this study are availablefrom the corresponding author upon reasonable request. The (PIP)code is available at the following url: https://github.com/AstroSnow/PIP.
REFERENCES E. Priest and T. Forbes,
Magnetic Reconnection: MHD Theory and Applications (Cambridge University Press, New York, 2000). D. Biskamp, in
Magnetic Reconnection in Plasmas , Cambridge Monographs onPlasma Physics (Cambridge University Press, Cambridge, UK, 2000), Vol. 3, p.387.
Physics of Plasmas
ARTICLE scitation.org/journal/php
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.003223628
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.003223628 , 032901-19Published under license by AIP Publishing E. N. Parker, “Sweet’s mechanism for merging magnetic fields in conductingfluids,” J. Geophys. Res. , 509–520, https://doi.org/10.1029/JZ062i004p00509(1957). P. A. Sweet, “The neutral point theory of solar flares,” in
ElectromagneticPhenomena in Cosmical Physics , Proceedings from International AstronomicalUnion Symposium No. 6 , edited by B. Lehnert (Cambridge University Press,1958), p. 123. K. Shibata, T. Nakamura, T. Matsumoto, K. Otsuji, T. J. Okamoto, N.Nishizuka, T. Kawate, H. Watanabe, S. Nagata, S. UeNo, R. Kitai, S. Nozawa, S.Tsuneta, Y. Suematsu, K. Ichimoto, T. Shimizu, Y. Katsukawa, T. D. Tarbell, T.E. Berger, B. W. Lites, R. A. Shine, and A. M. Title, “Chromospheric anemonejets as evidence of ubiquitous reconnection,” Science , 1591 (2007). N. Nishizuka, T. Nakamura, T. Kawate, K. A. P. Singh, and K. Shibata,“Statistical study of chromospheric anemone jets observed with hinode/SOT,”Astrophys. J. , 43 (2011). K. A. P. Singh, K. Shibata, N. Nishizuka, and H. Isobe, “Chromospheric anem-one jets and magnetic reconnection in partially ionized solar atmosphere,”Phys. Plasmas , 111210 (2011). K. A. P. Singh, H. Isobe, N. Nishizuka, K. Nishida, and K. Shibata, “Multipleplasma ejections and intermittent nature of magnetic reconnection in solarchromospheric anemone jets,” Astrophys. J. , 33 (2012). L. J. Guo, B. De Pontieu, Y. M. Huang, H. Peter, and A. Bhattacharjee,“Observations and modeling of the onset of fast reconnection in the solar tran-sition region,” arXiv:2009.11475 (2020). H. P. Furth, J. Killeen, and M. N. Rosenbluth, “Finite-resistivity instabilities ofa sheet pinch,” Phys. Fluids , 459–484 (1963). S. Tanuma, T. Yokoyama, T. Kudoh, and K. Shibata, “Two-dimensional magne-tohydrodynamic numerical simulations of magnetic reconnection triggered by asupernova shock in the interstellar medium: Generation of x-ray gas in the gal-axy,” Astrophys. J. , 312–332 (2001). K. Shibata and S. Tanuma, “Plasmoid-induced-reconnection and fractalreconnection,” Earth Planets Space , 473–482 (2001). R. Samtaney, N. F. Loureiro, D. A. Uzdensky, A. A. Schekochihin, and S. C.Cowley, “Formation of plasmoid chains in magnetic reconnection,” Phys. Rev.Lett. , 105004 (2009). N. F. Loureiro, A. A. Schekochihin, and S. C. Cowley, “Instability ofcurrent sheets and formation of plasmoid chains,” Phys. Plasmas , 100703(2007). N. F. Loureiro, R. Samtaney, A. A. Schekochihin, and D. A. Uzdensky,“Magnetic reconnection and stochastic plasmoid chains in high-Lundquist-number plasmas,” Phys. Plasmas , 042303 (2012). N. F. Loureiro and D. A. Uzdensky, “Magnetic reconnection: From the Sweet-Parker model to stochastic plasmoid chains,” Plasma Phys. Controlled Fusion , 014021 (2016). E. G. Zweibel, “Magnetic reconnection in partially ionized gases,” Astrophys. J. , 550 (1989). L. Ni, B. Kliem, J. Lin, and N. Wu, “Fast magnetic reconnection in the solarchromosphere mediated by the plasmoid instability,” Astrophys. J. , 79(2015). W. Park, D. A. Monticello, and R. B. White, “Reconnection rates of magneticfields including the effects of viscosity,” Phys. Fluids , 137–149 (1984). R. S. Steinolfson and G. van Hoven, “Nonlinear evolution of the resistive tear-ing mode,” Phys. Fluids , 1207–1214 (1984). D. Biskamp, “Magnetic reconnection via current sheets,” Phys. Fluids ,1520–1531 (1986). L. C. Lee and Z. F. Fu, “Multiple X line reconnection 1. A criterion for the tran-sition from a single X line to a multiple X line reconnection,” J. Geophys. Res. , 6807–6815, https://doi.org/10.1029/JA091iA06p06807 (1986). S. P. Jin and W. H. Ip, “Two-dimensional compressible magnetohydrodynamicsimulation of the driven reconnection process,” Phys. Fluids B , 1927–1936(1991). M. Ugai, “Computer studies on powerful magnetic energy conversion by thespontaneous fast reconnection mechanism,” Phys. Plasmas , 388–397 (1995). N. F. Loureiro, S. C. Cowley, W. D. Dorland, M. G. Haines, and A. A.Schekochihin, “X-point collapse and saturation in the nonlinear tearing modereconnection,” Phys. Rev. Lett. , 235003 (2005). S. Takasao, A. Asai, H. Isobe, and K. Shibata, “Simultaneous observation ofreconnection inflow and outflow associated with the 2010 August 18 solarflare,” ApJ Lett. , L6 (2012). A. Bhattacharjee, Y.-M. Huang, H. Yang, and B. Rogers, “Fast reconnection inhigh-Lundquist-number plasmas due to the plasmoid instability,” Phys.Plasmas , 112102 (2009). P. A. Cassak, M. A. Shay, and J. F. Drake, “Scaling of Sweet-Parker reconnec-tion with secondary islands,” Phys. Plasmas , 120702 (2009). Y.-M. Huang and A. Bhattacharjee, “Scaling laws of resistive magnetohydrody-namic reconnection in the high-Lundquist-number, plasmoid-unstableregime,” Phys. Plasmas , 062104 (2010). L. Ni, K. Germaschewski, Y.-M. Huang, B. P. Sullivan, H. Yang, and A.Bhattacharjee, “Linear plasmoid instability of thin current sheets with shearflow,” Phys. Plasmas , 052109 (2010). L. Ni, U. Ziegler, Y.-M. Huang, J. Lin, and Z. Mei, “Effects of plasma b on theplasmoid instability,” Phys. Plasmas , 072902 (2012). L. Ni, J. Lin, and N. A. Murphy, “Effects of the non-uniform initial environ-ment and the guide field on the plasmoid instability,” Phys. Plasmas ,061206 (2013). F. Pucci and M. Velli, “Reconnection of quasi-singular current sheets: The‘ideal’ tearing mode,” ApJ Lett. , L19 (2014). A. Tenerani, A. F. Rappazzo, M. Velli, and F. Pucci, “The tearing mode insta-bility of thin current sheets: The transition to fast reconnection in the presenceof viscosity,” Astrophys. J. , 145 (2015). A. Tenerani, M. Velli, A. F. Rappazzo, and F. Pucci, “Magnetic reconnection:Recursive current sheet collapse triggered by ‘Ideal’ tearing,” ApJ Lett. , L32(2015). F. Pucci, M. Velli, A. Tenerani, and D. Del Sarto, “Onset of fast ‘ideal’ tearingin thin current sheets: Dependence on the equilibrium current profile,” Phys.Plasmas , 032113 (2018). Y. M. Huang, L. Comisso, and A. Bhattacharjee, “Plasmoid instability in evolv-ing current sheets and onset of fast reconnection,” Astrophys. J. , 75 (2017). L. Ni and V. S. Lukin, “Onset of secondary instabilities and plasma heatingduring magnetic reconnection in strongly magnetized regions of the low solaratmosphere,” Astrophys. J. , 144 (2018). L. Comisso and D. Grasso, “Visco-resistive plasmoid instability,” Phys.Plasmas , 032111 (2016). L. Comisso, M. Lingam, Y. M. Huang, and A. Bhattacharjee, “General theoryof the plasmoid instability,” Phys. Plasmas , 100702 (2016). L. Comisso, M. Lingam, Y. M. Huang, and A. Bhattacharjee, “Plasmoid instabil-ity in forming current sheets,” Astrophys. J. , 142 (2017). J. E. Leake, V. S. Lukin, M. G. Linton, and E. T. Meier, “Multi-fluid simulationsof chromospheric magnetic reconnection in a weakly ionized reacting plasma,”Astrophys. J. , 109 (2012). J. E. Leake, V. S. Lukin, and M. G. Linton, “Magnetic reconnection in a weaklyionized plasma,” Phys. Plasmas , 061202 (2013). T. Tajima and J. I. Sakai, “Explosive coalescence of magnetic islands,” IEEETrans. Plasma Sci. , 929–933 (1986). J. E. Vernazza, E. H. Avrett, and R. Loeser, “Structure of the solar chromo-sphere. III. Models of the EUV brightness components of the quiet sun,” ApJSuppl. , 635–725 (1981). G. W. Pneuman, S. K. Solanki, and J. O. Stenflo, “Structure and merging ofsolar magnetic fluxtubes,” Astron. Astrophys. , 231–242 (1986), https://ui.adsabs.harvard.edu/abs/1986A%26A...154..231P/abstract. E. Khomenko, M. Collados, and T. Felipe, “Nonlinear numerical simulations ofmagneto-acoustic wave propagation in small-scale flux tubes,” Sol. Phys. ,589–611 (2008). K. A. P. Singh, A. Hillier, H. Isobe, and K. Shibata, “Nonlinear instability andintermittent nature of magnetic reconnection in solar chromosphere,” PASJ , 96 (2015). E. G. Zweibel, E. Lawrence, J. Yoo, H. Ji, M. Yamada, and L. M. Malyshkin,“Magnetic reconnection in partially ionized plasmas,” Phys. Plasmas ,111211 (2011). P. D. Smith and J. I. Sakai, “Chromospheric magnetic reconnection: Two-fluidsimulations of coalescing current loops,” Astron. Astrophys. , 569–575(2008).
Physics of Plasmas
ARTICLE scitation.org/journal/php
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.003223628
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.003223628 , 032901-20Published under license by AIP Publishing J. I. Sakai and P. D. Smith, “Two-fluid simulations of coalescing penumbra fila-ments driven by neutral-hydrogen flows,” ApJ Lett. , L45–L48 (2009). A. Hillier, S. Takasao, and N. Nakamura, “The formation and evolution ofreconnection-driven, slow-mode shocks in a partially ionised plasma,” Astron.Astrophys. , A112 (2016). S. I. Braginskii, “Transport processes in a plasma,” Rev. Plasma Phys. , 205(1965), https://ui.adsabs.harvard.edu/abs/1965RvPP....1..205B/abstract. E. T. Meier and U. Shumlak, “A general nonlinear fluid model for reactingplasma-neutral mixtures,” Phys. Plasmas , 072508 (2012). B. T. Draine, “Multicomponent, reacting MHD flows,” MNRAS , 133–148(1986). G. P. Zank, L. Adhikari, L. L. Zhao, P. Mostafavi, E. J. Zirnstein, and D. J.McComas, “The pickup ion-mediated solar wind,” Astrophys. J. , 23(2018). A. Hillier, “Ion-neutral decoupling in the nonlinear Kelvin-Helmholtz instabil-ity: Case of field-aligned flow,” Phys. Plasmas , 082902 (2019). V. M. Fadeev, I. F. Kvabtskhava, and N. N. Komarov, “Self-focusing of localplasma currents,” Nucl. Fusion , 202–209 (1965). T. D. Arber, G. J. J. Botha, and C. S. Brady, “Effect of solar chromospheric neu-trals on equilibrium field structures,” Astrophys. J. , 1183–1188 (2009). E. G. Harris, “On a plasma sheath separating regions of oppositely directedmagnetic field,” Il Nuovo Cimento , 115–121 (1962). A. Brandenburg and E. G. Zweibel, “The formation of sharp structures byambipolar diffusion,” Astrophys. J. Lett. , L91–L94 (1994). A. Brandenburg and E. G. Zweibel, “Effects of pressure and resistivity on the ambi-polar diffusion singularity: Too little, too late,” Astrophys. J. , 734 (1995). E. N. Parker, “The solar-flare phenomenon and the theory of reconnection andannihiliation of magnetic fields,” ApJ Suppl. , 177 (1963). E. T. Vishniac and A. Lazarian, “Reconnection in the interstellar medium,”Astrophys. J. , 193–203 (1999). F. Heitsch and E. G. Zweibel, “Fast reconnection in a two-stage process,”Astrophys. J. , 229–244 (2003). A. Alvarez Laguna, A. Lani, N. N. Mansour, H. Deconinck, and S. Poedts,“Effect of radiation on chromospheric magnetic reconnection: reactive and col-lisional multi-fluid simulations,” Astrophys. J. , 117 (2017). T. Shibayama, K. Kusano, T. Miyoshi, T. Nakabou, and G. Vekstein, “Fastmagnetic reconnection supported by sporadic small-scale Petschek-typeshocks,” Phys. Plasmas , 100706 (2015). S. Zenitani and T. Miyoshi, “Magnetohydrodynamic structure of a plasmoid infast reconnection in low-beta plasmas,” Phys. Plasmas , 022105 (2011). P. G. Drazin and W. H. Reid, “Hydrodynamic stability,” NASA STI/ReconTech. Rep. A , 17950 (1981). P. Drazin, “Dynamical meteorology—Kelvin–Helmholtz instability,” in
Encyclopedia of Atmospheric Sciences , edited by G. R. North, J. Pyle, and F.Zhang, 2nd ed. (Academic Press, Oxford, 2015), pp. 343–346. A. Hillier and V. Polito, “Observations of the Kelvin-Helmholtz instabilitydriven by dynamic motions in a solar prominence,” ApJ Lett. , L10 (2018). B. Snow and A. Hillier, “Intermediate shock sub-structures within a slow-modeshock occurring in partially ionised plasma,” Astron. Astrophys. , A46(2019). F. Pucci, K. A. P. Singh, A. Tenerani, and M. Velli, “Tearing modes in partiallyionized plasmas,” arXiv:2006.03957 (2020). J. Vranjes and P. S. Krstic, “Collisions, magnetization, and transport coeffi-cients in the lower solar atmosphere,” Astron. Astrophys. , A22 (2013). J. E. Leake, T. D. Arber, and M. L. Khodachenko, “Collisional dissipation ofAlfv (cid:2) en waves in a partially ionised solar chromosphere,” Astron. Astrophys. , 1091–1098 (2005). E. Khomenko and M. Collados, “Simulations of chromospheric heating byambipolar diffusion,” in Proceedings of the Second ATST-EAST Meeting:Magnetic Fields from the Photosphere to the Corona, Astronomical Society ofthe Pacific Conference Series, Vol. 463, edited by T. R. Rimmele, A. Tritschler,F. W € oger, M. Collados Vera, H. Socas-Navarro, R. Schlichenmaier, M.Carlsson, T. Berger, A. Cadavid, P. R. Gilbert, P. R. Goode, and M. Kn € olker,2012, p. 281. E. Khomenko and M. Collados, “Heating of the magnetized solar chromo-sphere by partial ionization effects,” Astrophys. J. , 87 (2012). B. P. Pandey and M. Wardle, “Hall magnetohydrodynamics of partially ionizedplasmas,” MNRAS , 2269–2278 (2008). N. A. Murphy and V. S. Lukin, “Asymmetric magnetic reconnection in weaklyionized chromospheric plasmas,” Astrophys. J. , 134 (2015). L. Ni, V. S. Lukin, N. A. Murphy, and J. Lin, “Magnetic reconnection instrongly magnetized regions of the low solar chromosphere,” Astrophys. J. , 95 (2018). L. Ni, V. S. Lukin, N. A. Murphy, and J. Lin, “Magnetic reconnection in the lowsolar chromosphere with a more realistic radiative cooling model,” Phys.Plasmas , 042903 (2018). B. Snow and A. Hillier, “Collisional ionisation, recombination and ionisationpotential in two-fluid slow-mode shocks: Analytical and numerical results,”arXiv:2010.06303 (2020). V. L. Dorman and R. M. Kulsrud, “One-dimensional merging of magneticfields with cooling,” Astrophys. J. , 777 (1995). D. A. Uzdensky and J. C. McKinney, “Magnetic reconnection with radiativecooling. I. Optically thin regime,” Phys. Plasmas , 042105 (2011). J. Leenaarts, M. Carlsson, V. Hansteen, and R. J. Rutten, “Non-equilibriumhydrogen ionization in 2D simulations of the solar atmosphere,” Astron.Astrophys. , 625–632 (2007). T. P. Golding, M. Carlsson, and J. Leenaarts, “Detailed and simplified nonequi-librium helium ionization in the solar atmosphere,” Astrophys. J. , 30 (2014). J. Mart (cid:2) ınez-Sykora, B. De Pontieu, M. Carlsson, V. H. Hansteen, D. N (cid:2) obrega-Siverio, and B. V. Gudiksen, “Two-dimensional radiative magnetohydrody-namic simulations of partial ionization in the chromosphere. II. Dynamics andenergetics of the low solar atmosphere,” Astrophys. J. , 36 (2017). R. J. Rutten, “Non-equilibrium spectrum formation affecting solar irradiance,”Sol. Phys. , 165 (2019).
Physics of Plasmas
ARTICLE scitation.org/journal/php
Phys. Plasmas , 032901 (2021); doi: 10.1063/5.003223628