Magnetic reconnection is a fundamental phenomenon in highly conductive magnetized plasmas such as space and astrophysical plasmas in which the magnetic energy is abruptly released and converted to the heating and kinetic energies as well as the acceleration of particles. Magnetic reconnection takes place in narrow regions where the frozen-in-flow constraint breaks down. As a result the field lines cut and reconnect to each other continuously and the large scale topology of magnetic field lines changes significantly (Birn and Priest, 2007; Yamada et al., 2010; Priest and Forbes, 2000; Gonzalez and Parker 2016).

The initial steady-state classical models of the magnetic reconnection process were discussed in MHD framework in two-dimensional by Sweet-Parker and Petschek models. In high Lundquist numbers, S = Lv _A /η (L is the system size, v A ( = B 0 / μ 0 ρ 0 ) is the Alfvén velocity and η is the magnetic diffusion) when the Lundquist number exceeds a critical value (S _c ≃ 10⁴) (Biskamp, 1986) the elongated thin current sheet is fragmented and new X-points (reconnection site) are generated. Therefore, multiple X-points and secondary magnetic islands fill the current layer. These magnetic islands merge with each other and form a larger islands. This type of MHD instability is known as “Plasmoid instability” (Bhattacharjee et al., 2009). The formation of plasmoids leads to a reconnection rate much faster than that predicted by previous models (Huang and Bhattacharjee, 2010; Uzdensky et al., 2010; Comisso et al., 2015; Comisso and Grasso, 2016). For weakly collisional systems such as solar corona, the Lundquist number is very large (∽ 10¹²–10¹⁴) and hence, the growth of plasmoid instability in elongated current sheets of solar flares are expected (Comisso et al., 2016; Comisso et al., 2017).

Numerical simulation of the magnetic reconnection event plays a crucial role in the conception of plasmoid instability dynamics and what happens inside the current sheets. Hence, in the last decade, the study of numerical simulation of the plasmoid instability was interested (Huang and Bhattacharjee, 2013; Loureiro and Uzdensky, 2015). Some of these studies are MHD simulation in two (Yu et al., 2010; Murphy 2010; Bárta et al., 2008) and three dimensions (Ugai 2008; MacTaggart and Fletcher, 2019). These simulations have uncovered some new aspects of the non-linear evolution of magnetic reconnection and plasmoid instability. In addition, the plasmoid instability has also been observed and discussed in many kinetic simulations (Stanier et al., 2019; Daughton and Karimabadi, 2007) and also there is tentative observational evidence (Bemporad 2008) that plasmoid might play a key role in the dynamics of magnetic reconnection. Theoretically, plasmoid formation has been proposed as a mechanism of fast reconnection (Lapenta 2008; Daughton et al., 2006) and non-thermal particle acceleration in reconnection related events (Drake et al., 2006).

In MHD numerical simulations, some physical parameters play a major role in the dynamics of plasmoid instability and magnetic reconnection. Some publications have considered the effects of some main parameters on the plasmoid instability such as plasma viscosity (Comisso, et al., 2015; Comisso and Grasso, 2016), asymmetric magnetic field (Murphy et al., 2013), shear flow (Hosseinpour et al., 2018) and non-uniform plasma mass density in two sides of the current layer (Doss et al., 2015; Birn et al., 2008). Besides, magnetic and plasma pressure play important roles in MHD studies. The ratio of plasma pressure to the magnetic pressure is described by the plasma-β(≡ p _plasma /p _mag ). This parameter in the upstream region is a key parameter in the reconnection dynamics and might vary from β > 1 in the photosphere at the base of the field line to β ≪ 1 in the mid-corona (Gary 2001).

In many of the researches on the plasmoid instability, the plasma-β effect has been ignored (Hosseinpour et al., 2018; Dong et al., 2018; Huang et al., 2017) and it is usually considered a constant and large. Hence, the β effect is not clear yet on the plasmoid dynamics. So, it is important to understand the physics of the plasmoid instability with various β values. The β effect not only shows itself in the linear phase of the instability but also in the nonlinear phase. The Lundquist number depends on the β and in the linear stage is proportional to β as S ∝ β ^−1/2.

Ni et al. (2012) and Baty (2014) with considering Harris current sheet configuration investigated the effect of plasma-β on the critical Lundquist number, S _c , for the onset of the plasmoid instability. In an initially uniform temperature configuration, Ni et al. (2012) concluded that the critical Lundquist number is strongly dependent on the plasma-β and claimed that the S _c varies between 10⁴ for β = 0.2 and 2000 for β = 50. Also, Baty (2014) with changing plasma resistivity (η), investigated the effect of varying plasma-β (0.2–15) and the equilibrium structure (uniform temperature or uniform density) on the plasmoid instability in a Sweet-Parker like layer. They found that the critical Lundquist number depends on the β, and indicated that for higher β, the critical Lundquist number is smaller. Also, the simulation results of Ni et al. (2012) showed that the reconnection rate for larger β is higher than smaller β, and the current sheet becomes turbulent earlier in larger β. To reach this conclusion they used β = 0.2, 1, 5, 50. Recently Zenitani and Miyoshi (2020) investigated the properties of plasmoid instability for β = 0.2, 1, 5. They showed that the reconnection rate increases as β decreases.

Many physical systems observed in laboratory experiments, space and astrophysical plasmas can be modeled as a simple current sheet with a magnetic field reversal (Harris sheet). A current sheet in the solar corona is one of the few concrete examples of solar applications where force-free state is possible, which is not necessary the case in other applications. However, often additional component of the magnetic field need to be considered. In the Earth’s magnetopause, there is an out-of-plane magnetic field (referred to as guide field) of strength comparable to the lobe field on either side of the Harris sheet. The magnetic reconnection in the presence of a guide field is an important issue in magnetospheric physics. In astrophysical systems, such as jets from acceleration disk, the magnetic field is believed to be primarily aligned with the current and can be represented by Harris sheet with a very strong guide field much larger than the reconnection plane field. Additionally, in magnetic confinement fusion devices, reconnection develops in the presence of strong toroidal fields and the configuration can be represented by a Harris sheet in the poloidal plane with a strong guide field in the toroidal direction.

The presence of guide field has been shown to be important for particle acceleration issue (Zenitani and Hoshino, 2008; Hamilton et al., 2003; Li and Lin, 2012) and also is widely considered in kinetic simulations (Inoue et al., 2015; Fu et al., 2007; Huang et al., 2011), which have shown that the strength of the guide field controls the growth of secondary magnetic islands and can evidently alter not only the trajectory of the particles but also the structure of the electric and velocity fields in the vicinity of the reconnection region. These simulations showed that the reconnection rate decreases with the guide field. The dependence of Hall mediated magnetic reconnection dynamics on the guide field is also investigated (Yang et al., 2008; Huba 2005), and found that the reconnection rate and plasma energization are reduced by increasing the guide field strength. In the MHD scale, the guide field effect on the plasmoid instability is very different from the results of the kinetic simulations (Ni et al., 2013). However, in 3D MHD simulations it has been shown that the guide field makes the reconnected field lines, as well as the jet, inclined from the current sheet normal direction, and also the guide field acts as an obstacle to the jet formation (Huang and Bhattacharjee, 2016; Hashimoto et al., 2005).

In spite of many related simulation works, there are still some unambiguities regarding the detailed effects of the plasma-β and the guide field on the plasmoid instability. Therefore, in this work, we investigate the effect of plasma-β in the range (β ≤ 1) with and without a guide field on the dynamics of plasmoid instability. Three different cases are considered: A: with a zero guide field, B: with uniform guide field and C: with non-uniform guide field. It should be mentioned that the effect of plasma viscosity is ignored in our study by ignoring the respective term in MHD equations. The effect of plasma viscosity on the dynamics of plasmoid instability has already been discussed in detail by Comisso, et al., 2015, Comisso and Grasso, 2016. They have shown that plasma viscosity has the effect of decreasing the linear growth rate and the wavenumber of the instability. However, despite its damping effect, for very high Lundquist numbers the plasmoid instability turns out to be very rapid in the linear regime. The initial temperature is assumed to be uniform. Therefore, the initial plasma density and pressure are strongly dependent on the plasma-β. The paper is organized as follows. In Section 2, the basic MHD equations, numerical setup, and initial condition are described. In Section 3, numerical results are presented. A summary and a short discussion of the results are then given in Section 4.

The compressible single-fluid resistive MHD equations are solved by using the OpenMHD code being developed by Zenitani (2016) to investigate the dynamics of the plasmoid instability in 2.5-dimensional cartesian coordinate system. The basic equations are ∂ t ρ = − ∇ . ( ρ . V ) , (1) ∂ t ( ρ V ) = − ∇ . [ ρ V V + p t o t I − B B ] , (2) ∂ t ϵ = − ∇ . [ ( ϵ + p t o t ) V − ( V . B ) B + η j × B ] , (3) ∂ t B = − ∇ . ( V B − B V ) − ∇ × ( η j ) , (4) E + V × B = η j , (5) which are written in conservation-law form. Here, ρ, V, B, I, j, η are the plasma mass density, the flow velocity, the magnetic field, the unitary tensor, the electric current density and the magnetic diffusivity, respectively. The p _tot = p + B ²/2 is total pressure and ϵ = p/(Γ − 1) + ρv ²/2 + B ²/2 is total energy density. For the convenience of numerical calculations, all variables are normalized. For this purpose, the plasma mass density is normalized to the initial plasma mass density (ρ ₀), the magnetic field to the initial magnetic field (B ₀), the flow velocity to the Alfvén velocity (V _A ), the electric current density to L ₀ B ₀/μ ₀ and resistivity to L ₀ v _A . Note that, L ₀ is the length scale of the system and all spatial variables are normalized to L ₀. We set the adiabatic index Γ = 5/3. All variables are function of space (x, y) and time (t) and the variation of variables in the z-direction is ignored, ∂/∂z = 0. Equations 1–5, are solved by a Godunov-type code. The code employs the HLLD scheme (Miyoshi and Kusano, 2005) to calculate numerical fluxes. The second-order Runge-Kutta method is used for time marching. Also, the hyperbolic divergence cleaning method (∇.B = 0) is employed for the solenoidal condition (Dedner et al., 2002). Moreover, the time step is based on convective and diffusive CFL condition. Our simulations are carried out in the x − y plane. The size of simulation box is set to be x = [ − 90, 90] and y = [ − 9, 9]. The number of grid points are N _x = 3,000 and N _y = 300 so that the grid sizes are Δx = 2L _x /N _x = 0.06, Δy = 2L _y /N _y = 0.06. We set open boundary condition at x = ±L _x direction, so that the reconnected field lines can leave boundaries freely. On the other hand, conducting boundary condition is considered for the bottom y = − L _y and top y = + L _y boundaries.

An initial Harris current sheet is used in the form (B = [B _x , 0, B _z ]): B x ( y ) = B 0 t a n h y a B , B y = 0 , B z ( y ) = B g , (6) where B ₀ = 1.0 is initial asymptotic magnetic field strength and a _B = 0.7 is the half width of the current sheet in the initial geometry. B _g is the guide field that is perpendicular to the reconnection plane. Initial plasma velocity is zero.

The initial static pressure, P, is obtained by solving the equilibrium equation and is given by: P ( y ) = 1 2 [ B 0 2 ( 1 + β ) − ( B x 2 + B z 2 ) ] . (7)

Assuming the isothermal equation of state (P = 2ρT) and solving the above Equation 7, the initial plasma mass density is satisfied as follows: ρ ( y ) = P ( y ) 1 + β . (8)

The initial plasma pressure (Eq. 7) and plasma density (Eq. 8) depend on the plasma- β ( = 2 ρ T / B 0 2 ) . In these simulations, the initially isothermal condition is used so that the variation of plasma-β is associated with the variation of plasma mass density in the upstream region. Also, in the MHD scale, the guide field appears as a term in the pressure equation. Hence, the Alfvén velocity is determined by both the plasma-β and the strength of the guide field. In all our simulations, the time is normalized with the Alfvén transit time scale τ _A = a _B /v _A . On the other hand, to achieve fast reconnection the initial resistive disturbance is set during 0 < t < 5 as η = η ₀  exp [ − (x ² + y ²)] where η ₀ = 0.025. As a result, an X-point is formed at the origin (x = y = 0). However, for t > 5, a uniform resistivity η = 2 × 10^–3 is assumed.

In this section, we discuss three different cases: A: instability with a zero guide field (B _g = 0) where the effect of plasma-β in the range (0.1 ≤ β ≤ 1) is addressed. B: instability with a uniform guide field where B _g = 0.1, 0.3, 0.5B ₀ and β = 0.4, 1.0 is assumed. C: instability with the non-uniform guide field with β = 0.1, … , 1.0.

In this study, the plasma-β effect in the range of β(= 0.1, 0.2, 0.4, 0.6, 0.8, 1.0) is investigated on the dynamics of plasmoid instability as a fundamental parameter in the magnetic reconnection with and without guide field. We used 2.5-dimensional MHD simulations and considered a standard Harris current sheet profile to establish an initial equilibrium. We added a magnetic field component perpendicular to the reconnection plane (so-called “guide field”) to the standard Harris current sheet. Three different cases was considered: case A: instability with zero guide field, case B: instability with the uniform guide field (B _g = 0.1, 0.3, 0.5B ₀) and case C: instability with the non-uniform guide field ( B g = B 0 2 − B x 2 ) . In the early times, a non-uniform localized resistivity was applied to trigger a Petschek type reconnection at the origin. Latter, a uniform resistivity was set which corresponds to a Lundquist number sufficient for the plasmoid instability. The Petschek structure is shortly converted to an elongated thin Sweet-Parker layer, which is subsequently fragmented to multiple X- and O-points (magnetic islands). Then, secondary magnetic islands merge and produce larger plasmoids. In these simulations, the temperature is assumed to be uniform and the plasma mass density profile was obtained from the initial equilibrium condition. The plasma mass density profile has a gradient from the upstream region to the center of the current sheet, thus the plasma mass density is non-uniform and β dependence could be largely attributed to the density variation. Therefore, as β increases the Alfvén speed and the Lundquist number decrease, but in all of our simulations, the Lundquist number is greater than the critical Lundquist number (S _c ⋍10⁴).

The simulation results can be listed as follows: 1. From case A, we found that the plasma-β leads to change in reconnection rate and growth rate of the plasmoid instability by affecting the plasma compressibility. Plasma compressibility decreases with increasing β, so the reconnection rate slows down, and the linear phase of the system becomes longer, which means a longer time is needed to achieve instability. For β > 0.2 cases, in which the Lundquist number is greater than the critical value (S > S _c ), the nonlinear stage is not observed, and secondary magnetic islands do not form. Also, the growth rate of instability (Eq. 10) decreases with increasing β. 2. By adding the uniform guide field (B _g = 0.1, 0.3, 0.5B ₀) (case B), we saw that the weak guide field (0.1B ₀) does not much affect on the reconnection process, but for the larger guide fields (0.3B ₀ and 0.5B ₀) the reconnection rate increases and the system becomes unstable faster. In particular, for the β = 0.4 case, the medium guide filed (0.5B ₀) was investigated in detail. Comparing this case with the case of zero guide field showed that in the presence of the guide field the reconnection rate increases, and the system reaches the nonlinear phase, so the secondary magnetic islands are generated. For β > 0.4 cases, although the Lundquist number is larger than the critical value, but due to the significant flow incompressibility, the nonlinear stage is not observed at time scales on the order of previous cases. Moreover, according to Figure 8, by comparing the maximum outflow velocity (V _x,out ), we found that V _x,out increases with increasing the uniform guide field, and we obtained the respective scalings. In the last case (case C), the non-uniform guide field was added. This guide field leads to a constant pressure and density everywhere. The magnetic pressure due to the guide field in the current sheet is an obstacle to the reconnection process. However, the linear stage for all β cases is very shorter than the cases of A and B, so the primary plasmoid in the center of the current sheet grows rapidly. In this state, the system cannot enter the non-linear phase and the secondary magnetic islands do not form. The reconnection rate decreases with increasing β in all cases. The presence of the uniform guide field has more impact on the reconnection rate than a non-uniform guide field. In cases B and C for β > 0.6, the reconnection rate shows almost equality, and for β < 0.4 cases the guide field has a more effect on the reconnection process (see Figure 13).