Parametric Study of Resistive Plasmoid Instability

By using 2.5-dimensional resistive MHD simulations, dynamics of the plasmoid instability in a Harris current sheet has been studied with taking into account two main controlling parameters: the plasma-β in the range (0 < β < 1) and the amplitude ratio of magnetic guide field to the reconnection plane field in three different cases with zero, uniform, and non-uniform guide field. Varying the plasma-β changes the plasma compressibility which affects significantly on the linear and nonlinear growth rates of the plasmoid instability. For each of three cases, some associated scaling relations between the instability growth rate, the plasma-β and the magnitude of guide field are obtained.


INTRODUCTION
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 4 ) (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 12 -10 14 ) and hence, the growth of plasmoid instability in elongated current sheets of solar flares are expected 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 , 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 midcorona (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 4 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  Frontiers in Astronomy and Space Sciences | www.frontiersin.org October 2021 | Volume 8 | Article 768965 3 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 Frontiers in Astronomy and Space Sciences | www.frontiersin.org October 2021 | Volume 8 | Article 768965 6 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 . 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 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 andGrasso, 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.

EQUATIONS AND NUMERICAL MODEL
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.5dimensional cartesian coordinate system. The basic equations are 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 /2 is total pressure and ϵ p/(Γ − 1) + ρv 2 /2 + B 2 /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 (ρ 0 ), the magnetic field to the initial magnetic field (B 0 ), the flow velocity to the Alfvén velocity (V A ), the electric current density to L 0 B 0 /μ 0 and resistivity to L 0 v A . Note that, L 0 is the length scale of the system and all spatial variables are normalized to L 0 . 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/zz 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 ]): where B 0 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: Assuming the isothermal equation of state (P 2ρT) and solving the above Equation 7, the initial plasma mass density is satisfied as follows: The initial plasma pressure (Eq. 7) and plasma density (Eq. 8) depend on the plasma-β( 2ρT/B 2 0 ). 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 η η 0 exp [ − (x 2 + y 2 )] where η 0 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.

SIMULATION RESULTS
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 FIGURE 13 | The maximum reconnection rate for three cases.

Instability With Zero Guide Field (B g = 0)
In the first step, we consider the case without the presence of any guide field and discuss the growth of plasmoid instability with different values of plasma-β. The subsequent variation of main plasma parameters in the inflow region, such as (p, ρ, v A , S) with plasma-β is shown in Table 1. For all cases in Table 1, the Lundquist number is larger than the critical value (S > S c ), so the plasmoid instability is expected to be triggered in the current sheet. We diagnose the magnetic reconnection rate during plasmoid instability with taking the average of the electric field component perpendicular to the reconnection plane, E z , on the y 0 line: Where ψ is the reconnected flux function defined by B ∇ψ(x, y) ×ẑ. The reconnection rate is normalized to B 0 v A . The time evolution of the reconnection rate for different values of β( 0.1, 0.2, 0.4, 0.6, 0.8, 1.0) is shown in Figure 1, in which different linear and non-linear phases of instability can be distinguished.
As can be seen, in both linear and non-linear phases, the reconnection rate decreases as plasma-β increases. This means that it will take a much longer time for the system to develop plasmoid instability and the subsequent transition to the nonlinear stage occurs at longer time scales. Here, larger β corresponds to higher density plasmas (see Table 1) with more flow incompressibility which results in a slow reconnection. For smaller β cases (β 0.1, 0.2) where the system transits faster into the non-linear phase, we find that the fast reconnection evolution becomes more drastic. To better understand the growth of plasmoid instability in the current layer, we consider a case that includes linear, transition, and non-linear phases, i.e. β 0.2 (case b in Figure 1). Figure 2 shows the contour plots of the magnetic field lines at different phases according to plot b in Figure 1. From Figure 2 we can see that how the primary current sheet becomes unstable to plasmoid instability which is eventually filled with multiple X-points and O-points. Figures 2A,B are plotted in the linear phase. The initial X-point of Figure 2A is transformed into an elongated Sweet-Parker current sheet, while Figure 2C (at time t 45) represents the transition state to the non-linear stage. New X-points are formed around the original O-point and in continue a small magnetic island, which its size is small compared to the width of current sheet (a B 0.7) grows ( Figure 2D). Figure 2E represents the field lines in the non-linear regime. Following the transition into the non-linear phase at t > 60 (for β 0.2) the secondary magnetic islands are generated. As a result, they merge and form a larger magnetic island. As the merging of magnetic islands continues, a large size plasmoid, named "monster plasmoid", can be produced whose size even exceeds the scale length of magnetic shear. Note that, inside plasmoid, high density and so high-pressure plasma is enclosed.
We now discuss the effect of different values of β on the dynamics of plasmoid instability by using contour plots of magnetic field lines. These contour plots are coordinate with the reconnection rate plots in Figure 1. Figure 3 shows the magnetic field lines for different values of β at the same time, t 95. Figures 3A,B correspond to β 0.1 and 0.2 show that the secondary magnetic islands form inside the current sheet and so the plasmoid instability developed at there. Since the β is low in these cases, the small plasma compressibility leads to increased reconnection rate, and thus the current sheet becomes unstable. The oscillations of reconnection rate decrease as β increases in the nonlinear phase. Each peak in the reconnection rate plot is related to the generation of new X-points where the dissipation of magnetic energy is significant. As the plasma-β increases the plasmoid instability growth rate decreases and the plasmoid generation turns out to be slower. This point can be seen from the size of the monster plasmoid in different panels of Figure 3. According to Figure 3, we see that in cases with smaller β, the magnetic field lines are piled up around the outflow regions. Therefore, the plasma outflow jets are ejected more rapidly due to the high pressure behind the plasma outflow jet, inside the islands. Figure 4 shows the variations of the logarithmic average of reconnection rate (Ln dER dt ) with β in the linear phase of the plasmoid instability in the time range (∼ 40 < t < ∼ 60) where the growth rate increases almost linearly with time for all β cases. In fact, measurements have been done at t ∼ 50. The blue squares correspond to the numerical values at the beginning of the transition phase (where the primary plasmoid is born in the origin of the current sheet and its size is smaller than the width of the initial current layer). The red dashed line also represents the quadratic fitted curve (scaling) of the numerical values obtained in this simulation, which is a scaling of the growth rate of the plasmoid instability in the linear phase. At the beginning of the transition phase where the primary plasmoid appears, the scaling is found to be: where the coefficients are a 0.51, b − 2.74 and c − 12.14. As a result, it can be seen that in plasmas with lower compressibility (a case with the higher β), the plasmoid instability growth rate also decreases following a reduction in the reconnection rate. Also, at the longer times (when the primary magnetic island grows), the growth rate of instability decreases faster.

Instability With Uniform Guide Field
In this section, we consider the presence of the uniform guide field and discuss the effect of this magnetic field on the dynamics of plasmoid instability. Different values of uniform guide field assumed in our study: B z B g 0.1, 0.3, 0.5B 0 . Due to the dependence of the initial plasma pressure to the guide field in our model (Eq. 7), plasma pressure decreases in the current sheet and upstream region with increasing the magnitude of the guide field. For a better understanding of the uniform guide field effect on the physical parameters, we consider β 0.4, 1.0 in Table 2. In all of cases shown in Table 2, according to Eq. 7, and Eq. 8, increasing the magnitude of the guide field leads to increase in the Alfvén Frontiers in Astronomy and Space Sciences | www.frontiersin.org October 2021 | Volume 8 | Article 768965 velocity and the Lundquist number. Therefore, it is expected that as the Lundquist number increases, the magnetic reconnection rate will increase. Figure 5 shows the magnetic reconnection rate for β 0.4 with different values of the guide field. In a weak guide field (B g 0.1B 0 ), the system behaves similar to the case of zero guide field (B g 0), but in the medium one the rate of magnetic reconnection increases, and for B g 0.5B 0 the nonlinear stage appears faster, which means that the plasmoid instability developed and the secondary magnetic islands are formed in the current layer. For all cases (β 0.1, 1.0 and B g 0.1, 0.3, 0.5B 0 ), the simulations are carried out. The results show that increasing the strength of the uniform guide field not only leads to the rise in the maximum reconnection rate, but also the linear phase of the system becomes shorter. Also, for all cases β > 0.4 and B g 0.5B 0 , although, the Lundquist number is greater than the critical value (S > S c ), the nonlinear phase was not observed, and secondary plasmoids do not form in the current layer, because in high β systems the plasma pressure is high at the center of the current sheet and the plasma behaves like an incompressible plasma which suppresses up to some extent the plasmoid instability.
We now present a scaling for the growth rate of instability. Figure 6 shows the effect of the uniform guide field on the growth rate of instability for β 0.4, 1.0 and B g 0, 0.1, 0.3, 0.5B 0 . Therefore, it can be seen that the presence of a guide field leads to an increase in the growth rate of instability for both β values. The blue stars and the blue diamonds represent the numerical values obtained from the simulations, and the red dashed lines are the quadratic fitted curve as a scaling of the instability growth rate as: with a − 3, b − 0.55, c − 7.9 for β 0.4 and a − 2.5, b − 0.64, c − 8.9 for β 1.0. The variation of magnetic energy with uniform and zero guide fields are shown in Figure 7 which is normalized to the initial magnetic energy. As we know, the magnetic reconnection process is a mechanism to convert the stored magnetic energy to plasma kinetic and thermal energies. With an increased reconnection rate by applying the uniform guide field, the conversion of the magnetic energy becomes faster. This point is clear from Figure 7. Figure 8 shows the time variation of the maximum outflow velocity at the X-points (V x,out ) in the presence of different values of the uniform guide field. In all three cases of the uniform guide field, the plasma outflow velocity from X-points increases until the system reaches the saturation state. After that, due to the increase of pressure in the outflow regions, the velocity slows down and becomes uniform. In general, as shown in Figure 8, the maximum outflow velocity increases with increasing guide field and around X-points is of the order of the Alfvén velocity. Now, we are looking to find a scaling relation for variation of the maximum outflow velocity in different values of the guide field. Figure 9 shows the variation of the V x,out with uniform B g in our simulations. Plots are considered for two cases. The first case (red squares) is for the maximum outflow velocity between t 1 65 and t 2 75, while the second case (blue stars) is the maximum outflow velocity between t 2 75 and t 3 85. In both cases, the maximum outflow velocity increases as the uniform guide field increases. The dashed lines in Figure 9 represent the V x,out scaling with the uniform B g , which has been obtained from the quadratic fitted curve of the numerical values: with a 8.9, b − 0.13, c 2.15 for t 1 < t < t 2 curve and a 8.7, b − 0.13, c 2.3 for t 2 < t < t 3 curve.

Instability With the Non-uniform Guide Field
Now, we discuss the effect of non-uniform guide field. Note that the general structure of three components of the initial equilibrium of the magnetic field is in the form B 2 x + B 2 y + B 2 z B 2 0 , where B y 0 and B 0 1.0. Therefore, in our study we consider the non-uniform guide field as: As we have described in Section 2, Eq. 7, we obtain the initial plasma pressure as: Hence, for the case of non-uniform guide field the initial plasma pressure and so the plasma density are uniform everywhere.
The initial profiles of the plasma pressure (p), magnetic pressures (B 2 x /2, B 2 z /2), reconnection plane magnetic field (B x ), guide field (B z ) and total pressure (P t P + B 2 x /2 + B 2 z /2) are shown in Figure 10 for β 0.1 and 1. The total pressure remains constant inside and outside of the current layer. In our study, the uniform temperature was assumed. Accordingly, the variation of β results from the change in density, so in the presence of the non-uniform guide field the plasma pressure and the plasma density only dependent on the β. The magnetic pressure due to the non-uniform guide field reaches its maximum value in the center of the diffusion region.
The reconnection rate plots are shows in the presence of the non-uniform guide field for different β in Figure 11. As seen, similar to the previous cases, the reconnection rate decreases by increasing the plasm-β. The presence of a non-uniform guide field makes much shorter the linear phase of the instability compared to previous cases (with zero and uniform guide field) for all β cases. Therefore, the instability rapidly enters the transition phase, and the initial plasmoid is generated and grows in the center of the current sheet. Unlike the previous cases, when a non-uniform guide field is added to the reconnection plane, the system cannot enter the nonlinear phase due to the increased magnetic pressure caused by the guide field in the center of the current layer, except β 0.1 case, so secondary magnetic islands do not form in the current layer. Over time, the primary plasmoid grows, and when the width of the primary plasmoid is much larger than the width of the current layer, the system becomes saturated, and then the reconnection rates decrease. Figure 12 shows the linear growth rate of the plasmoid instability with β with a non-uniform guide field. As expected, the growth rate of the instability decreases as β increases. By comparing this case with the case A (without the guide field), we find that due to the presence of the non-uniform guide field, the decreasing slope of the linear instability growth rate is sharper. From the numerical values obtained in this simulation (blue stars) and quadratic fitted curve (red dashed line) in Figure 12, we can calculate the scaling of the linear growth rate as: where the coefficients are a 1.5, b − 1.4 and c − 10.14. The coefficient a in the quadratic equations expresses the slope of the parabolic graph, and its sign is related to its direction. According to the scaling equations for the instability growth rate in cases A and C, this coefficient increases with the presence of the guide field. Also, in case B, the presence of the medium guide field has led to an increase in this coefficient.

SUMMARY AND CONCLUSION
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 0 ) and case C: instability with the non-uniform guide field 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 4 ).
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 0 ) (case B), we saw that the weak guide field (0.1B 0 ) does not much affect on the reconnection process, but for the larger guide fields (0.3B 0 and 0.5B 0 ) the reconnection rate increases and the system becomes unstable faster. In particular, for the β 0.4 case, the medium guide filed (0.5B 0 ) 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 nonuniform 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).

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
HL: Doing the numerical simulations, visualizing the data and preliminary writing. MH: Analyzing the data and writing and revising some part of manuscript.