{'fraction_non_alphanumeric': 0.05086932301287126, 'fraction_numerical': 0.027706922332658986, 'mean_word_length': 3.93379622932639, 'pattern_counts': {'":': 0, '<': 194, '<?xml version=': 0, '>': 52, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A New Approach of Linear Theory of Tearing Instability in Uniform Resistivity

The linear perturbation equation of the tearing instability derived in LSC theory (Loureiro, Schekochihin, and Cowley, PoP2007) is numerically examined as an initial value problem, where the inner and outer regions are seamlessly solved under uniform resistivity. Hence, all regions are solved as the resistive MHD (magnetohydrodynamics). To comprehensively study physically acceptable perturbation solutions, the behaviors of the local maximum points required for physically acceptable solutions and zero-crossing points, at which \phi=0 and \psi=0, are examined. Eventually, the uniform resistivity assumed in the outer region is shown to play an important role in improving some conclusions derived from the theory. In conclusion, the upper limit \lambda_{up} of the growth rate obtained in the improved (modified) LSC theory is shown to be regulated by the Alfven speed measured in the outer region. It is also shown to be partially consistent with the growth rate in the linear developing stage of the impulsive tearing instability observed in the compressible MHD simulation of the plasmoid instability (PI) based on uniform resistivity.


Introduction
The magnetic reconnection process is an energy conversion mechanism from magnetic energy to plasma kinetic energy. To explain the explosive energy conversion observed in the solar flares and geomagnetic substorms, the magnetic reconnection process must be fast. The fast magnetic reconnection process has been studied during the past exceeding 60 years.
Some theories of the magnetic reconnection process were proposed in the 1960s (Vasyliunas 1975), some of which have survived until today. In particular, two famous theoretical MHD (magnetohydrodynamic) models, i.e., the Sweet-Parker (SP) model and the Petschek (PK) model, are well known steady-state models (Parker 1957;Petschek 1964). In general, the PK model is much more efficient than the SP model (Vasyliunas 1975); hence, the PK model is believed to be a candidate for explaining the fast magnetic reconnection process required for solar flares and geomagnetic substorms.
FKR (Furth, Killeen, and Rosenbluth (Furth, et.al. 1963)) theory was also proposed in the 1960s to explore the tearing instability caused by the magnetic reconnection process. Since FKR theory is a linear theory, the magnetic reconnection process based on FKR theory may be developed to some nonlinear instability. If the nonlinear instability finally settles to a steady state, either the SP model or PK model will be applicable. Meanwhile, if the nonlinear instability maintains any non-steady state, the resulting magnetic reconnection process may be characterized by neither the SP model nor the PK model. Recently, plasmoid instability (PI) has been actively studied as such a nonsteady state model (Loureiro, et.al. 2007;Baalrud, et.al. 2012;Samtaney, et.al. 2009;Loureiro, et.al. 2009;Bhattacharjee,et.al. 2009;Cassak & Drake 2009;Ni,et.al. 2012;Baty 2014;Shibayama,et.al. 2015;Tenerani,et.al. 2015a,b;Landi,et.al. 2015;Loureiro & Uzdensky 2016;Huang,et.al. 2017).
In the 1970s, two theoretical models, i.e., the spontaneous and externally driven models, were proposed to numerically reproduce the PK model by means of MHD simulations. In the spontaneous model (Ugai & Tsuda 1977;Ugai 1984;Ugai,et.al. 2005;Shimizu & Ugai 2003;Shimizu,et.al. 2009a,b;Shimizu & Kondoh 2013;Shimizu,et.al. 2016) proposed by Ugai, the locally enhanced nonuniform resistivity at X-point results in the PK model. Such nonuniform resistivity is called anomalous resistivity. Ugai proposed current-driven anomalous resistivity, which is self-consistently enhanced by a positive feedback mechanism between the plasma inflow and outflow to and from the magnetic diffusion region, leading to the appearance of the PK model. The spontaneous model does not require any additional mechanism to establish the PK model, with the exception of the current-driven anomalous resistivity. Moreover, the externally driven model (Hayashi & Sato 1978) proposed by Sato also requires current-driven anomalous resistivity, also resulting in the PK model. In contrast to the spontaneous model, some externally driven mechanism in the upstream is required to maintain the PK model. If the externally driven mechanism is removed, the fast magnetic reconnection process terminates. Hence, the assistance of the externally driven mechanism is essential to maintain the PK model in the externally driven model, and whether the spontaneous model or externally driven model causes the fast magnetic reconnection process has become a controversial topic.
The magnetic reconnection process requires some form of resistivity to reconnect magnetic field lines at X-point. As shown by many numerical MHD studies, nonuniform resistivity, such as anomalous resistivity, results in the PK model. Moreover, when the resistivity is assumed to be uniform in time and space rather than nonuniform, whether the PK model can be generated or not is a delicate problem (Shibayama,et.al. 2015;Biskamp 1986;Kulsrud 2001;Malyshkin,et.al. 2005;Baty,et.al. 2006;Kulsrud 2011). If the uniform resistivity can reproduce the PK model, nonuniform resistivity is not necessarily required for the fast magnetic reconnection process. If so, one may say that the fast magnetic reconnection needs some form of resistivity but is independent of the resistivity. Then, the uniform versus nonuniform resistivity for the PK model is also currently a controversial topic.
Separate from steady-state models, such as the PK and SP models, non-steady state models may also be explored. The plasmoid instability (PI) is a non-steady state and multiple-tearing instability model developed from FKR and LSC theories. In fact, LSC theory (Loureiro, et.al. 2007) predicts that when the uniform resistivity is extremely weak, the linear growth rate can be extremely high. When the Lundquist number exceeds a critical value S c , such a high-speed non-steady and multiple-tearing instability is caused, which is called PI. In fact, some recent highly sophisticated numerical MHD studies reported that, once PI occurs, the reconnection rate exceeds the value predicted by the steady-state SP model (Loureiro, et.al. 2007;Baalrud, et.al. 2012;Loureiro, et.al. 2009;Bhattacharjee,et.al. 2009;Ni,et.al. 2012;Shibayama,et.al. 2015;Loureiro & Uzdensky 2016), suggesting that PI is another potential cause of the fast magnetic reconnection process of solar flares and substorms.
In contrast, recent numerical MHD studies (Ng & Ragunathan 2010;Shimizu,et.al. 2017) suggest that when numerical dissipations are sufficiently removed, the tearing instability tends to be less active; hence, it suggests that PI under uniform resistivity cannot cause the fast magnetic reconnection process, at least in the spontaneous model.
In other words, it suggests that S c does not exist, except in the externally driven model. From another perspective, PI may be associated with the non-steady state PK model rather than the non-steady state SP model. If so, the controversial topic of the uniform versus nonuniform resistivity for the PK model (Shibayama,et.al. 2015) may be extended to the non-steady state PK model.
The main theme of this paper is to propose the modified LSC theory as a new linear theory of tearing instability, improving the original LSC theory (Loureiro, et.al. 2007). In Section 2, LSC theory is numerically examined, where the inner and outer regions are seamlessly solved under uniform resistivity as an initial value problem (Shimizu 2018a,b). Then, the upper limit λ up of the linear growth rate is obtained for physically acceptable perturbation solutions. In addition, a candidate of the physically acceptable perturbation solutions, which is proposed in this paper as zero-crossing solution, is studied. In Section 3, the MHD simulation of PI under uniform resistivity is shown, where how the modified LSC theory can be applied to the tearing instability in PI is demonstrated. The modified LSC theory is then shown to be partially consistent with the MHD simulation of PI. In Section 4, the following topics are discussed. First, it is shown that the existence of S c is not directly supported by the modified LSC theories. Second, the original and modified LSC theories are compared for PI application. Third, the viscosity effect, which is not considered in LSC theory but is employed in the MHD simulation for numerical stabilization, is briefly discussed. Fourth, FKR theory viewed from the perspective of the modified LSC theory is discussed. Section 5 provides a summary of this study.

Basic equation
In this section, LSC theory (Loureiro, et.al. 2007) is numerically examined as an initial value problem (Shimizu 2018a,b), where the inner and outer regions are seamlessly solved on the assumption of uniform resistivity. LSC theory is based on incompressible onecomponent MHD equations. The largest difference between the FKR and LSC theories is that the zero-order flow field is, respectively, null and nonzero. By assuming a nonzero flow field, the thickness of the current sheet in LSC theory can be finite, thus maintaining the rigorous equilibrium. LSC theory can examine how the tearing instability grows in the nonzero flow field generated by the steady-state SP current sheet. In other words, FKR theory is not related to the SP model, because of the null flow field. At this point, LSC theory is better and more realistic than FKR theory. The basic equations derived in LSC theory are as follows.
φ − κ 2 2 φ = −κf (ξ)ψ + κf 2 (ξ)φ/λ + f (ξ)ψ/λ, (2.1) ψ − κ 2 2 ψ = κλψ − κf (ξ)φ. (2.2) Eqs.(2.1) and (2.2) basically follow LSC's notations, where φ and ψ are the scalar potential functions of the 2D incompressible velocity and 2D magnetic fields, respectively. The prime of φ and ψ indicates the derivative with respect to ξ, which is normalized by the characteristic width δ cs of the current sheet, as ξ = y/δ cs . In contrast to the LSC's notations, x and y are exchanged with each other, to compare LSC theory with the MHD simulation examined in Section 3. δ cs is close to the half thickness of the current sheet but is not exactly the distance between the center (i.e. origin ξ = 0) and edge of the current sheet. The edge is slightly separated from ξ = 1, and rather, is located at ξ = 1.307, i.e., y = 1.307δ cs .
Following the LSC's notation, λ is the growth rate normalized by l cs /V A , where l cs is the wave length of the plasmoid chain and V A is Alfven speed measured in the upstream magnetic field region. κ = πL cs /l cs is the wave number along the current sheet, where L cs is the total length of the steady state SP sheet along the sheet. In Section 3, L cs is calculated from the spatial gradient of the outflow velocity u x0 measured at X-point, as du x0 /dx = 2V A /L cs , where u x0 is defined, as below. Following LSC's procedures, the zero-order nonzero flow field and antiparallel magnetic field in the inner region of the current sheet, i.e., y < 1.307δ cs , are defined as described below.
Here, Γ 0 = 2V A /L cs is set. In addition, = 2δ cs /L cs = 2/ √ S corresponds to the uniform resistivity, where S is the Lundquist number and the steady state SP model is assumed. Then, f (ξ) is defined as shown below.
where ξ 0 = 1.307. Next, the field in the outside of the current sheet, i.e., y > 1.307δ cs = y 0 , is defined as described below.
It is important that Eqs.(2.1) and (2.2) can be solved as an initial value problem when φ(0), ψ(0), φ (0) and ψ (0) are given as the initial values for a set of κ, , and λ. Here, V A = 1 is assumed, and φ(0) = 0, ψ(0) = 1, and ψ (0) = 0 are set without the lack of generality of solutions, which are based on the symmetricity at ξ = 0. Then, φ (0) is a control parameter to uniquely determine a solution, in addition to κ, , and λ. Then, ψ (0) = 0 means when ∆ index employed in FKR theory is rigorously zero. The tearing instability in FKR theory is caused in ∆ > 0. At this point, ∆ index defined in LSC theory is different from that of FKR theory but a discontinuity of ∆ > 0 defined in FKR theory is assumed at the origin. Meanwhile, that in this paper is caused in ψ (0) > 0, as explained in the third paragraph of Section 2.2. It means that any discontinuity at ξ = 0 is not assumed in this paper.
To numerically solve φ and ψ, the forward Euler method is employed. The numerical resolution is set in 0.001 > ∆ξ > 0.00025 to suppress numerical errors. For a much higher numerical resolution, i.e., 0.00025 > ∆ξ, the convergence test of the numerical results was done in some typical cases. The following results are extremely sensitive for changing some control parameters but basically not changed, evenwhen the Runge-Kutta method is employed instead of the forward Euler method.

Physically acceptable φ and ψ
Depending on the control parameters of φ (0), κ, , and λ, a lot of various φ and ψ solutions can be obtained by numerically solving the initial value problem. We must select "physically acceptable" sets of φ and ψ from those various solutions. Such selection means to specify the physically acceptable upstream condition of φ and ψ. Simply, φ and ψ, which converge to zero at ξ = +∞, will be the most preferable candidate of physically acceptable solutions. Let us call it the zero-converging solution. In addition, φ and ψ, which reach zero in 0 < ξ +∞, will be another preferable candidate. Let us focus on when φ = 0 and ψ = 0 are simultaneously attained at a finite ξ value. Let us call it the zero-crossing solution as we define ξ = ξ c at the zero-crossing point of φ = ψ = 0.
In general, there will be many other physically acceptable solutions. In other words, the zero-converging and zero-crossing solutions are not the whole physically acceptable solutions. Rather, if φ and ψ do not diverge to infinity in the finite ξ range, they also may be a candidate for physically acceptable solutions in the range. It should be noted that they include extreme cases, in which φ and ψ diverge to infinity at ξ = +∞. In such extreme cases, the linear growth rate may exceed a unity, leading to super-Alfvenic plasmoid (tearing) instability (Huang,et.al. 2013). It will be fairly difficult to numerically explore the total behaviors of such physically acceptable solutions. Hence, in this paper, let us focus only on the zero-converging solutions and zero-crossing solutions introduced above, where the linear growth rate dose not exceed a unity, resulting in sub-Alfvenic tearing instability, as shown from Figs.6 to 11(b).
In the zero-converging and zero-crossing solutions, it is important that φ and ψ must have local maximum points in 0 < ξ < +∞, because, φ must start with φ(0) = 0 and φ (0) > 0 and must reach zero in 0 < ξ +∞. Note that φ (0) > 0 means the plasma inflow toward X-point and the outflow from O-point. Additionally, ψ must start with ψ(0) = 1 and ψ (0) = 0 and must reach zero in 0 < ξ +∞. At this point, ψ (0) > 0 is additionally required to cause the tearing instability because ψ (0) > 0 means when the current sheet becomes thin in the vicinity of X-point, and hence, corresponds to when ∆ -index is positive. However, ∆ = 0 is resolved through this paper, because of ψ (0) = 0. In other words, the unstable mode of tearing instability requires that ψ always has a downward convex feature at origin, and hence, a region of ψ > 0, at least, in the vicinity of origin. This is why ψ must have a local maximum point in 0 < ξ < +∞.
As a result, we must find φ and ψ, which have local maximum points in 0 < ξ < +∞. However, in numerical studies, such as the initial value problem started from ξ = 0, it is impossible to precisely explore the behaviors of φ and ψ at ξ = +∞. Instead, as the first step of the exploration, let us observe the behaviors of φ and ψ observed in the finite region, i.e., 0 < ξ < ξ b , where ξ b is simply the end of the numerical calculation. In other words, ξ b is an upper limiter of ξ to explore the local maximum points of φ and ψ in the range of 0 < ξ < ξ b . Increasing ξ b , the behaviors at ξ = +∞ can be deduced, followed by the upper limit of the growth rate λ. Let us call the upper limit λ up which is not necessarily equal to λ but will be useful, as shown below. Next, as the second step of the exploration, let us focus on the zero-crossing solution. Increasing ξ c , i.e., the crossing point, let us explore how the growth rate λ of the zero-crossing solutions depends on ξ c . Finally, the behaviors of zero-crossing solutions of ξ c = +∞ can be deduced, which will coincide with the zero-converging solution.
In addition, as will be discussed in Section 2.3.6, the zero-converging solutions and zero-crossing solutions observed in this paper may be classified into the "inner-triggered" and "outer-triggered" tearing instability. Let us define the former as when both local maximum points of φ and ψ are located in the inner region of the current sheet. Then, let us define the latter as when either of the local maximum points is at least located in the outer region of the current sheet. The latter is different from the externally-driven reconnection process proposed by Sato, et al. (Hayashi & Sato 1978). At this point, since φ and ψ are perturbation terms and hence have much weaker intensity than the zeroorder equilibrium defined as Eqs.(2.3)-(2.11), the latter does not mean that the tearing instability is caused by some apparently-strong externally-driven force in the upstream region. Note that tearing instability is essentially driven by the Alfven waves propagating in the inner or outer region of current sheet. Thus, the inner-triggered tearing instability is triggered by the weak Alfven waves in the inner region and the outer-triggered tearing instability is all other cases. At this point, we may assume that the main part of the trigger, i.e., the maximum amplitude of the weak Alfven wave, is located at the local maximum point of φ and ψ.

Basic features
Figure 1(a) shows the behaviors of φ and ψ numerically obtained from Eqs.(2.1) and (2.2) for φ (0) = 2.5 and κ = = 0.2 as λ varies between 0.1 and 2.0. This case is typical for the higher φ (0) range and is limited in the lower λ range. In the upper panel of Figure 1(a), as λ increases from label a to c, the local maximum point of φ shifts away from ξ = 0. Then, since φ of labels d-f exceeds the vertical axis scale range, the existence of a local maximum point is unclear. At this point, Figure 3(a) confirms the lack of local maximum point of φ in 0 < ξ < 5. In the middle panel of Figure 1(a), the behavior of ψ tends to be inverted with respect to that of φ. In fact, ψ of labels a-b does not have a local maximum point except at ξ = 0, while ψ of labels c-f does have a local maximum point slightly separated from the origin. Of note, ψ of labels c-d appears to be a flat-top function but has a local maximum point slightly separated from the origin, as shown in 5ψ (ξ = 0) = 0 in the middle panel of Figure 5. It means ψ (0) > 0 for every ψ shown in this figure. Eventually, for 0.2 < λ < 0.5, i.e., between labels b and d, φ and ψ may simultaneously have local maximum points in 0 < ξ < 5. Such φ and ψ may be physically acceptable. Figure 1(b) shows the behavior of φ and ψ for φ (0) = 2.5 and κ = = 0.2 while λ varies between 2 and 100. This is typical for the higher λ range than Figure 1(a), where label f in Figure 1(a) is the same as label a in figure 1(b), excluding the vertical axis scale. As λ increases, the local maximum point of φ goes to ξ = 0, while that of ψ shifts away from ξ = 0. As a result, Figure 1(b) does not have φ and ψ that simultaneously have local maximum points in 0 < ξ < 5. At this point, labels a and b are unclear for the existence of the local maximum point of φ, which is confirmed in Figure 3(a). Also, that of ψ is unclear for labels c-f, which is confirmed in Figure 3(a). Eventually, no physically acceptable solution exists in Figure 1 Figire 2 shows the behaviors of φ and ψ for φ (0) = 1.0 and κ = = 0.2 while λ is varied between 0.1 and 100. This is typical for the lower φ (0) range than Figures 1(a) and (b). In Figure 2, it seems that φ and ψ do not simultaneously have local maximum   Fig.1(a). Label a is the same as label f in Fig.1(a) but the vertical axis scale is different. points in 0 < ξ < 5. At this point, the existence of the local maximum point is confirmed in Figure 3   and yellow regions indicate where it does not exist. Comparing with the upper and lower panels, the blue region in the lower panel represents where the local maximum points of φ and ψ simultaneously exist in 0 < ξ < ξ b , and the yellow region represents where both do not exist there. Hence, the blue region may include physically acceptable φ and ψ in term of 0 < ξ < ξ b .
The pink regions in the upper and lower panels are completely separated by the blue or yellow region in the lower panel. In fact, the blue and yellow regions are contacted around λ = 0.9 and φ (0) = 1.52. The contact point indicates the location of the upper limit of λ for the physically acceptable φ and ψ. Let us define the upper limit as λ up . The location of the contact point is difficult to specify exactly because the numerically obtained φ and ψ are extremely sensitive to λ and φ (0). However, by carefully changing λ and φ (0), we can specify λ up , as shown below. Figure 3(b) shows the existence of the local maximum points of φ and ψ for 0.90 < λ < 0.92 and 1.52 < φ (0) < 1.54 for κ = = 0.2 and 0 < ξ < ξ b = 5. Accordingly, Figure 3 shows only a limited region of Figure 3(a). The contact point between the blue and yellow regions is clearly observed between 0.9116 < λ < 0.9120 and 1.5278 < φ (0) = 1.5281. Hence, λ up = 0.9116 is roughly the upper limit of the growth rate for κ = = 0.2 and 0 < ξ < ξ b = 5. Figure 4 shows the existence of the local maximum points of φ and ψ for the same parameter range as in Figure 3(a), with the exception that 0 < ξ < ξ b = 1.307, which corresponds to the case of the inner-triggered tearing instability. In other words, this figure shows when the local maximum points of φ and ψ are found inside of the current sheet. In the same manner as discussed in Figures 3(a) and (b), 0.7001 < λ up < 0.7101 is obtained around φ (0) = 1.84. In comparison to Figure 3(a), λ up tends to decrease as ξ b decreases from 5.0 to 1.307.

Zero-crossing solutions
We would ideally like to find φ and ψ that converge to zero at ξ = +∞, but as long as we rely on only a numerical analysis, such values are impossible to determine. Instead, Figure 3(a). The existence map of the local maximum points of φ (upper panel) and ψ (lower panel) for κ = = 0.2 and ξ b = 5. The horizontal axis is shown for 0.0001 < λ < 7, and the vertical axis is for 0.0001 < φ (0) < 5. Labels A, B, and C, respectively, indicate the parameter range examined in Figs.1(a), (b), and 2. The slight color gradation shows how far the local maximum point is separated from the origin.
by modifying the zero-order equilibrium, such a zero-converging solution can be found, as shown in Appendix A. However, such zero-converging solutions are not necessarily required for physically acceptable φ and ψ. Rather, if both φ and ψ reach zero in a limited finite ξ range, they may also represent a physically acceptable solution in the finite ξ range. In this section, let us define such solutions as "zero-crossing solution" and try to find it through the parameter survey executed in Figures 1-3. In fact, such zero-crossing solutions are shown as label c in Figures 1(a) and 5. Figure 5 shows a zero-crossing solution in which φ = 0 and ψ = 0 are simultaneously established at ξ = (ξ c =)7.36525, which is located outside of the current sheet, i.e., ξ > 1.307. This solution is obtained for κ = = 0.2, λ = 0.540043, and φ (0) = 1.9. To precisely determine this zero-crossing solution, a primitive pinching method is employed, as shown in Appendix B. ψ is also plotted in the middle panel to clearly illustrate the local maximum point of ψ, which is separated from the origin. The local maximum points of φ and ψ are, respectively, ξ = 1.071 and 1.035, which are located inside of the current sheet. Hence, as will be discussed in Section 2,3,6, this is the case for the inner-triggered tearing instability. It may be noted that label c of Figures 1(a) and 5 are obtained for the same κ and but, respectively, for λ = 0.349 and 0.540. However, those crossing points ξ c , i.e. the locations of upstream boundary, are different. It suggests that the growth rate significantly depends on the upstream boundary condition.
2.3.4. κ and dependences of λ up Figure 6 shows the ξ b dependence of the upper limit λ up for varying κ and , where 0 < ξ < ξ b is the parameter survey range of when the local maximum points of φ and ψ Figure 3(b). The existence map of the local maximum points of φ (upper panel) and ψ (lower panel) for κ = = 0.2 and ξ b = 5. The horizontal axis is shown for 0.9 < λ < 0.92 and the vertical axis is for 1.52 < φ (0) < 1.54. The right edge of the blue region is detected at λ = 0.9116 and the left edge of the yellow region is detected at λ = 0.9120. Hence, the contact point between the blue and yellow regions is located in 0.9116 < λ < 0.9120 and 1.5278 < φ (0) < 1.5281. are simultaneously found. Hence, ξ b may be considered to be another control parameter to explore the local maximum points, and then, measure λ up . As ξ b increases, λ up monotonically increases. Hence, this figure shows that either of the local maximum points of φ and ψ, at least, shifts to a larger ξ value for a larger growth rate λ. Additionally, Figure 6 suggests that λ up does not exceed unity, i.e., the Alfvenic rate V A /l cs ; in other words, the tearing instability always proceeds at sub-Alfvenic speed. In addition, even for = 0, i.e., the zero resistivity limit, φ and ψ can be obtained, and λ up appears to converge to unity as ξ b increases, which may be strange but not surprising, as will be discussed below. Figure 7 shows how λ up appears to converge to a value less than unity as is close to zero, where ξ b = 1.307 is set. As shown in figure 6, since λ up increases as ξ b increases, λ up for ξ b = 1.307 is considered to be the highest growth rate of the inner-triggered tearing instability. In addition, Figure 7 shows that, when κ is close to zero, λ up is not sensitive to while λ up for larger κ rapidly decreases as increases from zero.
Next, let us discuss why λ up can be obtained even for = 0, i.e., the zero resistivity limit. As simply predicted, the reconnection process, i.e., tearing instability, must stop in = 0, i.e. ideal MHD. However, the limit of = 0 examined in this paper does not directly mean = 0. In fact, the spatial scale ξ is normalized by the thickness δ cs of the current sheet, where δ cs is defined for the steady-state SP sheet. Let us consider when L cs is constant. As = 2δ cs /L cs approaches zero, δ cs is unlimitedly close to zero. At this time, since the current density unlimitedly increases, the instability will be promoted, i.e., not stop. In other words, evenwhen the current sheet unlimitedly becomes thinner in the limit of = 0, tearing instability can rapidly occur in the thin current sheet. In that case, λ up normalized by l cs /V A can become a value less than unity, which means    the tearing instability can proceed at sub-Alfvenic speed in the zero resistivity limit, i.e., the infinite Lundquist number S = +∞. In addition, Figure 7 suggests that the tearing instability cannot rapidly proceed beyond Alfven speed when local maximum points of φ and ψ exist, i.e. for zero-converging and zero-crossing solutions. It may be noted that λ up and λ studied in this paper is different from the growth rate γ/Γ 0 studied in the original LSC theory (e.g. FIG.4 (Loureiro, et.al. 2007)), as will be discussed in Section 2.3.7.
It may be noted that = 2δ cs /L cs = 0 can be established by either of δ cs = 0 or L cs = +∞. However, to establish the = 0 limit, δ cs = 0 is not equivalent to L cs = +∞. Because, the δ cs = 0 limit directly means that current sheet is unlimitedly thin but the L cs = +∞ limit changes κ = πL cs /l cs . Then, to keep a constant κ value, l cs = +∞ is required, leading to the change of the λ scale, i.e., l cs /V A .
2.3.5. κ, and ξ c dependences of λ for zero-crossing solutions Figure 8 shows the κ dependence of λ up for ξ b = 5.0, 1.307, and 0.5 for = 0. Regardless of ξ b , as κ increases from 0.002 to 20, three lines of λ up plotted for ξ b = 0.5, 1.307 and 5 tend to monotonically decrease. This result of λ up is applied to MHD simulation in the Figure 8. The κ dependence of λup for = 0 and ξ b = 0.5, 1.307, and 5. In addition, the one-dot, two-dot, and three-dot chain lines show the κ dependence of λ for the zero-crossing solution of = 0, in which the crossing point for each line is respectively fixed at ξc = 3.6, 7.2 and 10.8.
next section. Moreover, three lines of dot chain show the growth rate λ of zero-crossing solutions for = 0, in which the crossing point is respectively fixed at ξ c = 3.6, 7.2 and 10.8. In contrast to λ up , these three lines of the dot chain have a maximum point around κ = 1. More exactly, as ξ c increases from 3.6 to 10.8, the maximum point gradually shifts to a smaller κ. Thus, the growth rate λ of zero-crossing solutions depends on the crossing point ξ c , i.e., the upstream boundary condition, as mentioned in end of Section 2.3.3. The case of the ξ c = +∞ limit will correspond to the zero-converging solutions which may coincide with λ up in the ξ b = +∞ limit. As shown in Figure 6, λ up in the ξ b = +∞ limit is predicted to be unity.
Note that Figures 3 and 4 obtained for = 0.2 is directly inapplicable for Figure 8 for = 0.0. In addition, note that these three lines of dot chain are the growth rate λ, itself, but not the upper limit λ up . However, as ξ c increases, λ shown in Figure 8 tends to increase. This tendency is consistent with the λ up increase by increasing ξ b , which is shown in Figure 6. Figure 9 shows how λ of zero-crossing solutions depend on ξ c for κ = 0.001, 0.01, 0.05, 0.1, 0.2, 1, 3, and 5 of = 0. In every κ line, as ξ c increases, λ monotonically increases. Inversely observing, as ξ c decreases toward 1.307, λ seems to converge to zero. Due to the lack of numerical precision, the detail of λ behaviors around ξ c = 1.307 cannot be observed. Hence, it is unclear whether λ = 0 is attained exactly at ξ c = 1.307. As shown in Figure 8, Figure 9 also shows how the maximum point of λ shifts in κ space, depending on ξ c . In fact, as ξ c increases, the maximum point tends to shift from higher κ line to lower κ line. E.g., the highest λ value at ξ c = 5 is on the κ = 1 line, while the highest λ value at ξ c = 10 is on the κ = 0.2 line. Figure 10(a) shows how λ of zero-crossing solutions depends on ξ c for κ = 0.1 of = 0, 0.2, 1, 2, and 5. Evenwhen is non-zero, the monotonical increase in λ for the larger ξ c observed in Figure 9 seems to be retained. As increases from zero, λ tends to monotonically decrease. However, since λ lines in = 0 and 0.2 almost coincides, λ is not sensitive to for 0 < < 0.2. For this reason, in Section 3, the growth rate measured at = 0 is applied to the MHD simulation of the tearing instability which is not the case of = 0 but close to = 0. Figure 10(b) is similar to Figure 10(a) but shows the case of higher κ values, i.e., κ = 1  for = 0, 0.2 and 0.5 and the cases of κ = 3 for = 0 and 0.2. Evenwhen is non-zero, the monotonical increase in λ for larger ξ c observed in Figures 9 and 10(a) seems to be retained. However, in contrast to Figure 10(a), Figure 10(b) shows that λ for κ > 1 rapidly decreases for increasing .

Inner and outer-triggered tearing instability
As discussed in the end of Section 2.2, these numerical results shown in this section may be classified into the inner-triggered and outer-triggered tearing instability. The innertriggered tearing instability has the local maximum points of φ and ψ located in the inner region of the current sheet. Since such a tearing instability is initiated inside the current sheet, this will be called spontaneous tearing instability. Meanwhile, the outer-triggered tearing instability is all other cases, i.e., when either of the local maximum points of φ and ψ is located in the outer region. Figure 11(a) shows how the growth rate λ of zero-crossing solutions observed in Figures  8 and 9 changes for the movements of local maximum point of φ in ξ space, where the horizontal axis is the location ξ of the local maximum point of φ. Due to the lack of numerical precision, Figure 11(a) is unclear for the vicinity of λ = 0 and larger λ. The unclear range (the end of each line) largely depends on κ. As the location of the local maximum point of φ is separated from zero, λ increases. This feature is consistent with Figure 6 in which, as ξ b increases, λ up monotonically increases. In addition, as κ is close to zero, those curves tend to get lower. However, for all κ cases, it seems that the locations do not exceed ξ = 1.307, i.e., the boundary point of the inner and outer regions of the current sheet. Figure 11(b) is similar to Figure 11(a) but shows the movements of the local maximum point of ψ, where the horizontal axis is the location ξ of the local maximum point of ψ. As the location of the local maximum point of ψ is separated from zero, λ increases. This feature is consistent with Figures 6 and 11(a). In contrast to Figure 11(a), the local maximum point of ψ can shift beyond ξ = 1.307. Hence, when the location of the local maximum point of ψ is located in ξ < 1.307 is the case of the inner-triggered tearing instability. Figure 11(b) shows that the growth rate λ in the outer-triggered tearing instability is higher than that of the inner-triggered tearing instability.
2.3.7. Comparison of the growth rate with original LSC theory.
As mentioned in Section 2.1, the growth rate in this paper is normalized by l cs /V A . When we compare the growth rate λ and upper limit λ up studied in this paper with the growth rate in the original LSC theory (Loureiro, et.al. 2007), note that their growth rate is measured to be κλ = γ/Γ 0 , which is not λ, itself. In addition, their growth rate is obtained by separating inner and outer regions with the ∆ -index. For example, FIG.4 in Loureiro,et al. shows that κλ increases as κ 2/3 for a lower κ regime and decreases as κ −2/5 for higher κ regime. The finding means that λ decreases as κ −1/3 for a lower κ regime and decreases as κ −7/5 for a higher κ regime. Hence, λ monotonically decreases as κ increases. This feature is qualitatively consistent with Figure 8 for λ of κ > 1 and  λ up of the entire κ range. However, those behaviors of λ are difficult to quantitatively compare, because, λ and λ up in Figure 8 is very difficult to explore in κ > 5 due to the lack of numerical precision. Moreover, FIG.4 in Loureiro,et al. (Loureiro, et.al. 2007) appears to be unclear in κ << 1, where λ may be higher because it varies as κ −1/3 . To explore extremely high κ range of Figure 8, the λ = 0 limit solutions of Eqs.(2.1) and (2.2) will be helpful, which are extensively discussed in Appendix C. Figure 8 is similar to what was reported in some papers (Landi,et.al. 2015;Pucci & Velli 2014;Tenerani,et.al. 2015;Zanna,et.al. 2016;Papini,et.al. 2018), where Eqs.(2.1) and (2.2) were solved in the upstream boundary conditions which are different from this paper. At this point, this paper shows that the growth rate significantly depends on the upstream boundary condition, i.e., the locations of the zero-crossing point ξ c . In addition, their equilibrium seems to be different from that of this paper. At this point, most of them took the Γ 0 = 0 limit in the equilibrium of the LSC theory but this paper does not. In other words, they did not directly apply Eqs.(2.3)-(2.11) introduced by Loureiro. In fact, instead of Eqs.(2.4) and (2.7), the Harris type of current sheet as the magnetic field equilibrium was employed, and, instead of Eqs.(2.5), (2.6) and (2.10), the null-flow field was assumed. In the Γ 0 = 0 limit, since L cs = 2V A /Γ 0 becomes infinity, their results cannot be applied to the MHD simulation in the next section.
More recently, Shi,et.al. also numerically studied perturbation equations similar to Eqs.(2.1) and (2.2) (Shi,et.al. (2018)), where the inner and outer regions were seamlessly solved under uniform resistivity. Hence, their approach may be similar to what is shown in this section. However, there are some differences between their results and what is shown in this section. For example, they reported that tearing instability is stabilized by the background flow, i.e., Eq.(2.6), when Lundquist number S is less than 70. Meanwhile, as shown in Appendix C, this paper suggests that the growth rate λ can keep to be positive even when resistivity is extremely large, as long as the criterion of κ < 1.0915 is satified. Since low S corresponds to large , the criterion suggests that tearing instability can occur at small κ even in extremely low S.
It is unclear why the difference is caused but let us put an intuitive explanation for the tearing instability in small κ and low S, as below. First, note that tearing instability can occur when the current sheet at X-point becomes thin. Such a thinning steadily occurs if the outflow from the X-point is stronger than the inflow to the X-point. Because, at the time, the X-point becomes close to vacuum. The equilibrium in linear theory takes the balance between the outflow and inflow even in extremely low S. Hence, if the perturbed outflow is stronger than the perturbed inflow, i.e., φ (0) is large, the tearing instability will start. Next, let us consider about meanings of small κ. As defined in Fig.17, when dV x/dx measured at the X-point is large, L cs = 2V A /(dV x/dx) is small, leading to small κ = πL cs /l cs . At the time, resistivity will prevent the thinning at the diffusion speed. However, if we focus only on the sufficiently vicinity of the X-point, the thinning will be able to overcome the diffusion even in extremely low S. Hence, we cannot say that sufficiently low S steadily stabilizes the tearing instability.
However, it may be noted that the WKB approximation employed in LSC theory (Loureiro, et.al. 2007) fails for the small κ. Hence, the intuitive explanation mentioned above must be carefully re-examined in the future, improving the WKB approximation.

Spontaneous plasmoid instability
In this section, we apply the modified LSC theory, which resulted in Figure 8, to the MHD simulation of the spontaneous PI. At the end of this section, the modified LSC theory is shown to be partially consistent with the MHD simulation results.
The word "spontaneous" means that the subsequent tearing instability is spontaneously driven by the preceding tearing instability after the first tearing instability is externally initiated by an initial disturbance. Let us call such multiple tearing instabilities spontaneous PI. The spontaneous PI is enhanced and developed by a kind of positive feedback mechanism, including the PI itself. In addition to the spontaneous PI model, an externally driven PI model may be possible, where the tearing instability is maintained by an externally driven mechanism which may be a weak noise (Ng & Ragunathan 2010). Hence, if the driven mechanism is removed, the PI stops. However, the existence of such an externally driven mechanism in solar flares and substorms is unclear. In this paper, we focus on the spontaneous PI model.
In our previous paper (Shimizu,et.al. 2017), we found that when the current sheet becomes extremely thin, the numerical error may fatally affect the simulation results of the spontaneous PI model. Similar result has also been reported by Ng. et al. (Ng & Ragunathan 2010). In fact, the active PI observed at a lower numerical resolution tends to be less active at a higher numerical resolution. At a higher numerical resolution, the magnetic reconnection rate may momentarily exceed the value predicted by the steady state Sweet-Parker (SP) theory but may not constantly reach the level independent of the resistivity, which suggests that the critical Lundquist number S c does not exist (Shimizu,et.al. 2017).
In this paper, to avoid the numerical error due to lower numerical resolution, the case of relatively large resistivity, i.e. a relatively low Lundquist number S, is examined. Then, on the basis of modified LSC theory, the S = +∞ limit is discussed in Sections 4.1.

Simulation procedures
The one-component compressible 2D MHD eqs. are numerically solved. The reconnection process studied herein is essentially the same as that of our previous studies (Ugai 1984;Shimizu & Ugai 2003;Shimizu,et.al. 2017). In this paper, 2-step Lax-Wendroff scheme (Ugai & Tsuda 1977) is employed, but the results shown in this section are independent of the numerical scheme, because the numerical resolution is kept to be sufficiently high. The first tearing instability is initiated by a small resistive disturbance induced around the origin of the 1D current sheet. After the resistive disturbance is removed at t = 4, uniform resistivity is maintained for t > 4. Then, the thinning of the current sheet spontaneously starts around the most intensive region, i.e. the origin. Accordingly, we can examine whether the 1D current sheet is spontaneously destabilized by the tearing instability. In fact, after the first tearing instability is terminated by nonlinear saturation, subsequent tearing instabilities intermittently start in the current sheet. At t > 4, since no externally driven mechanism is applied to compress the current sheet to maintain the reconnection process, this situation represents the spontaneous PI model under uniform resistivity.
The simulation box size is limited to 0 < x < L x and 0 < y < L y , where (L x , L y ) = (400, 200). All boundary conditions of the simulation box are set to symmetric boundary conditions. The initial current sheet has a 1D structure, i.e., the magnetic field B = [B x (y), 0, 0] is assumed to be B x (y) = B x0 tanh y for 0 < y < L y . In contrast to previous studies (Shimizu,et.al. 2017), an inversed current sheet does not exist in the upstream region, but instead, L y = 200 is much larger than the initial thickness of the current sheet. Accordingly, the upstream boundary at y = L y will have a minimal numerical effect on the tearing instability in the current sheet near y = 0.
The plasma static pressure P initially satisfies the pressure-balance condition, i.e., 2P + B 2 x = 1 + β 0 , where β 0 is the ratio of the plasma pressure to the magnetic pressure in the magnetic field region of B x = 1.0. In this paper, β 0 = 0.15. The initial fluid velocity is u = (0, 0, 0) throughout the simulation box. Hence, the initial state is not equilibrium for uniform resistivity, leading to the magnetic annihilation. However, the magnetic reconnection, i.e., tearing instability, overcomes the magnetic annihilation and grows. The initial plasma density ρ = (P/2)/(1 + β 0 ), where the initial plasma temperature is uniform. The Alfven speed V A = B x / √ ρ measured in the upstream region of B x = 1.0 is approximately 5.5 with ρ = 0.033. The local speed of sound C s = γP/ρ in the center of the current sheet (y = 0) is approximately 2.0 with P = 0.575, γ = 5/3, and ρ = 0.25. The initial resistive disturbance for 0 < t < 4 is as follows.
η a = η a0 e −(x 2 +y 2 ) (3.1) where η a0 = 0.005. This resistive disturbance initiates the first tearing instability around the origin, i.e., (x, y) = (0, 0). The corresponding Lundquist number S for this resistive disturbance in 0 < t < 4 is estimated to be S = L x V A /η a0 = 440000. After t = 4, η a , which is defined by Eq.(3.1), is removed and, instead, uniform resistivity η b = 0.016 in time and space is assumed, which drives the first tearing instability, leading to PI. The corresponding Lundquist number S in t > 4 is estimated to be S = L x V A /η b = 137500. Since the first tearing instability initially occurs in a much smaller region than L x , a realistic S will be much smaller than 137500 and, hence, may be smaller than the critical Lundquist number S c = 10 5 ∼ 10 6 predicted previously (Loureiro, et.al. 2007). Because of the relatively low S, we can sufficiently suppress the numerical errors that may fatally affect the tearing instability. Hence, the numerical errors are not discussed in this paper, in contrast to our previous paper (Shimizu,et.al. 2017). The S = +∞ limit is then discussed in Section 4.1, using Figure 8. The simulation box (L x , L y ) is divided by numerical grids of ∆x(= ∆y) = 0.02 that are constant in time and space, where the time step ∆t = 0.004 is set to maintain the CFL conditions. Accordingly, the grid size is (N x , N y ) = (20000, 10000).
Unfortunately, the 2-step Lax-Wendroff scheme requires artificial viscosity to maintain numerical stability. In this paper, artificial viscosity that is uniform in time and space is applied to the mass, momentum, and energy conservation equations in the MHD eqs. but not for the magnetic flux conservation equation, i.e., Faraday's law, because uniform resistivity η b has been applied. The intensity of the artificial viscosity corresponds to 0.5 when translated as the magnetic Prandtl number. By contrast, LSC theory does not include any viscosity, i.e., the magnetic Prandtl number is zero. This discrepancy in viscosity between the theories and MHD simulations is briefly discussed in Section 4.3.

Overview of the numerical simulations
In this section, the first and second tearing instabilities are examined. After them, the third and subsequent tearing instabilities are impulsively repeated. After the third tearing instability ends, the numerical result starts to be gradually affected by numerical error, because of thinning of the current sheet. Hence, the numerical results after the third tearing instability is unreliable, which are not shown in this paper. Figure 12(a) shows the magnetic field lines and current density J z at t = 200. The initial resistive disturbance η a has been removed until this time, and the uniform resistivity η b is assumed. Red and yellow represent negative J z , and blue, which does not appear in this figure and is observed in Figures 12(b) and (c), represents positive. The J z intensity is indicated by the darkness of each color. More exactly, as shown in the color scale bar, the color intensities between red, white, and blue are respectively normalized by the maximum, zero, and minimum of J z in this figure. This figure shows the beginning of the first tearing instability directly resulting from the initial resistive disturbance, i.e., Eq.(3.1). The red colored current sheet is strongly localized around the origin. The red current sheet extending toward the surroundings of the plasmoid located at 20 < x < 50 is similar to the slow shock layer often observed in PK model. Hence, the first tearing instability may be expected to develop to the PK model but is immediately changed to SP-like sheet, as will be shown later. In addition, the color of the initial 1D current sheet observed in 100 < x < 400 is yellow and hence the current intensity is weakened because of magnetic annihilation, by which the current sheet is gradually diffused by uniform resistivity η b . Note that the initial 1D current sheet is not exactly the equilibrium under uniform resistivity. Eventually, the first tearing instability grows, overcoming the magnetic annihilation. The details of the first tearing instability are briefly presented in Appendix D. Figure 12(b) shows the magnetic field lines and the current density J z at t = 308. The Figure 12(a). Magnetic field lines and current density Jz color contour map at t = 200. This figure is shown for 0 < x < Lx = 400 and 0 < y < 8(<< Ly = 200), which is a limited part of the simulation box. Accordingly, note that the aspect ratio of this figure is extremely distorted to show the vicinity of the neutral sheet. The first tearing instability has been slowly developing around the origin. As in the PK-like current sheet, the red current sheet slightly splits from the neutral sheet, i.e., y = 0. Figure 12(b). Magnetic field lines and current density Jz color contour map at t = 308. The first tearing instability has been saturated, and an SP-like sheet is formed around 0 < x < 80. At this time, the second tearing instability has started near x = 75 but is still invisible.
SP-like sheet, which is shown in red, is formed around 0 < x < 75. Additionally, a largescale plasmoid, generated by the first tearing instability, grows around 75 < x < 200. In addition, the second tearing instability is starting, and a new X-point appears around x = 75. However, the new X-point is still invisible at this time. Figure 12(c) shows the magnetic field lines and the current density J z at t = 360. At this time, the second tearing instability develops around x = 140, and the red current sheet is strongly localized around this point. As the second tearing instability proceeds, a new plasmoid appears around 80 < x < 130. Since the second tearing instability is caused by the first tearing instability, it is a spontaneous tearing instability. In comparison to Figure 12(a), the second tearing instability develops more rapidly than the first tearing instability because the first tearing instability develops slowly during 4 < t < 300, whereas the second tearing instability develops rapidly during 308 < t < 380. After Figure 12(c), the third and fourth tearing instability start, leading to fully developed PI. Figure 13(a) shows the magnetic field lines and B y contour map, which is the recon-  nected field intensity at t = 308. The red region represents B y > 0, and the blue region represents B y < 0. The B y intensity is indicated by the darkness of each color. In fact, the left-side region of the large-scale plasmoid formed around 75 < x < 150 is deep red, and the right side is deep blue. These deep colors are associated with steady growth of the large-scale plasmoid. At this moment, the second tearing instability is starting, but the associated X-point and plasmoid are still invisible. However, a small light-blue region appears in the SP-like sheet, as indicated by a black arrow. The appearance of the light-blue region indicates the appearance of the second plasmoid and, hence, the beginning of the second tearing instability. Figure 13(b) shows the magnetic field lines and B y contour map at t = 360. In addition to the first (large-scale) plasmoid observed at 150 < x < 300, the second plasmoid generated by the second tearing instability grows in 80 < x < 120. Accordingly, the small light-blue region observed in Figure 13(a) develops into a vertically-wider deepblue region in 120 < x < 140 in Figure 13(b). Note that the aspect ratio of this figure is extremely distorted to show the vicinity of the neutral sheet. Figure 14(a) shows the magnetic field lines and V y contour map at t = 308. The blue  region, i.e., V y < 0, represents the inflow toward the current sheet, and the red region, i.e., V y > 0, represents the outflow from the current sheet. A negative high-intensity region of V y indicated by deep blue (i.e., labels A and B) is observed in the surroundings of the plasmoid at 75 < x < 150. Meanwhile, the color around the X-point at the origin, i.e., x = 0, is light blue or no color. This result is inconsistent with LSC theory because the highest negative V y (i.e., deep-blue) region assumed in the theory must be located around the X-point rather than around the plasmoid. This inconsistency will be associated with the compressibility and nonlinearity of this MHD simulation, while the theory is based on an incompressible and linearized MHD. Figure 14(b) shows the magnetic field lines and V y contour map at t = 360. As in Figure 14(a), a negative high-intensity region of V y (i.e., deep-blue region) is observed in the surroundings of the large scale plasmoid rather than the two X-points, i.e., at x = 0 and x = 140. This result suggests that the strong inflow to the SP-like sheet is driven by the movement of the plasmoids. In addition, it suggests that the plasma inflow around the X-point weakens as the plasmoid moves away from the X-point. To study how the inflow speed V y weakens, let us observe the time dependence of V y .   Figure 15 shows the magnetic field lines and ∆V y contour map at t = 308, where ∆V y = V y (t, x, y) − V y (t − 4, x, y) is defined as a velocity increment. The red region (∆V y > 0) indicates either deceleration of the inflow speed (V y < 0) or acceleration of the outflow speed (V y > 0). Inversely, the blue region indicates either acceleration of the inflow speed or deceleration of the outflow speed. Note that the locations of label A in Figure 14(a) and 15 are exactly the same. Since label A represents V y < 0 in Figure  14(a) and ∆V y > 0 in Figure 15, the plasma inflow toward the SP-like sheet is decelerated around label A. Furthermore, the locations of label B in Figures 14(a) and 15 are also exactly the same. Since label B is V y < 0 in Figure 14(a) and ∆V y < 0 in Figure 15, the plasma inflow toward the SP-like sheet is accelerated around label B. Hence, as the SP-like sheet elongates and the plasmoid moves away from the X-point, the plasma inflow around the plasmoid accelerates but the plasma inflow around the X-point decelerates. This deceleration suggests that the linear growth of the tearing instability is terminated as the SP-like sheet elongates, as will be discussed in Section 3.3.7.

x-directional profiles of J z and V x
In this section, we focus mainly on the x-directional behaviors of the second tearing instability observed during 308 < t < 380. Figure 16 shows the −J z profile in the xdirection at y = 0 at t = 308, 320, 340, and 360. As shown in Figures 12(a) to (c), Figure 16. The x-directional profile of the current density −Jz along the neutral sheet (y = 0) at t = 308, 320, 340, and 360. The first X-point is located at the −Jz peak of x = 0. The second X-point corresponds to another −Jz peak moving from x = 75 to 150.
keeping the first X-point at the origin, the second X-point has appeared at t = 360. The −J z peak observed at the origin in Figure 16 corresponds to the first X-point, and another −J z peak corresponding to the second X-point appears around x = 80 at t = 308. Figure  16 shows that the second peak gradually grows, moving in the +x direction. In FKR and LSC theories, ψ (0) > 0, which corresponds to ∆ > 0, predicts that this −J z peak appears at X-point. At this point, this MHD simulation is consistent with those theories. Figure 17 shows the reconnection outflow speed V x profile in the x-direction at y = 0 at t = 308, 320, 340, and 360. The V x peak, moving from x = 90 at t = 308 to x = 160 at t = 360, corresponds to where the reconnection jet collides with the plasmoid. Accordingly, the SP-like sheet is formed between the origin and the V x peak. This figure shows that V x does not yet reach the Alfven speed at these times, which is initially V A = 5.5 and almost unchanged over time. Hence, the first and second tearing instabilities do not yet reach the steady-state SP model. In addition, at t = 308 and 320, as x increases, V x monotonically increases in the SP-like sheet, while at t = 340 and 360, the second domed V x peak newly appears in the SP-like sheet, which are, respectively, located at x = 75 and 90. The appearance of the second domed V x peak is associated with the interaction between the first and second tearing instabilities. In addition, since V x at the second X-point is not zero and takes a positive value, the second X-point itself is moving in the +x direction. Because of the movement of X-point, the tearing instability in this MHD simulation cannot be directly applied to LSC theory, as will be discussed later in Section 3.3.4.
Additionally, Figure 17 explains how to measure L cs in the modified LSC theory. Once ∂V x /∂x is measured at X-point, L cs is calculated from L cs = 2V A /(∂V x /∂x), where V A is measured in the upstream magnetic field. In the same manner, L cs can be measured in each tearing instability, such as first, second, third and so on. Figure 17. The x-directional profile of the outflow speed Vx along the neutral sheet (y = 0) at t = 308, 320, 340, and 360. As shown, Lcs is measured from the Alfven speed VA and the spatial gradient ∂V x/∂x at the X-point, which is located around x = 75 at t = 308. Figure 18(a) shows the −J z profile in the y-direction at t = 308, 320, 340, and 360 at the second X-point. Note that since the symmetry boundary condition is set at y = 0, −J z profile for y < 0 is mirrored by this −J z profile for y > 0. In 308 < t < 340, the height and width of the −J z peak located at y = 0 are almost unchanged, but the −J z peak rapidly becomes higher during 340 < t < 360. Thus, the rapid increase of the −J z peak is delayed by the growth of the second tearing instability started at t = 308. In addition, Figure 18(a) shows that the −J z peak has a double current sheet structure during 308 < t < 340. In other words, the thin current sheet of 0 < y < 0.5 is embedded in the thick current sheet of 0 < y < 2. These two sheets are separated by the dent indicated by label A in Figure 18(a). This double current sheet structure is similar to what was observed in a previous study (Papini,et.al. 2019). This double current sheet structure cannot be directly applied to LSC theory, as will be discussed later in Section 3.3.6. Figure 18(b) shows the −J z profile in the y-direction at t = 308, 320, 340, and 360 in the second plasmoid, where these profiles are plotted at the local maximum point of B y at each time, e.g., which is located in the blue contour region around x = 75 of Figure  13(a) and x = 120 of Figure 13(b). The −J z peak at t = 360, which is indicated by label A, is separated from y = 0. This separation is caused by the plasmoid formation, where label A corresponds to the outer edge of the plasmoid. Hence, as the tearing instability is developed, the current sheet around the plasmoid shown in Figure 18(b) gradually thickens due to the growth of the plasmoid, while the thickness around the X-point shown in Figure 18(a) gradually thins. As a result, the current sheet thickness gradually becomes nonuniform along the sheet. This nonuniformity means that the linear theory, such as the modified LSC theory, is inapplicable at this moment. It has entered into a nonlinear phase until t = 360. Figure 18(a). The y-directional profile of the current density −Jz at the second X-point shown in Fig.12(c). The X-points at t = 308, 320, 340, and 360, respectively, are located at x = 78.9, 97.5, 127.8, and 155.1. The dent in the −Jz profile observed at t = 320 and 340 is indicated by label A, which is the boundary between the thin (inner) and thick (outer) current sheets. Figure 18(b). The y-directional profile of the current density −Jz in the second plasmoid. These profiles are plotted at the local maximum points of By at t = 308, 320, 340, and 360, which are, respectively, located at x = 67.5, 84.6, 114.0, and 134.1. Figure 19(a) shows the −B y profile in the y-direction at t = 308, 320, 340, and 360 in the second plasmoid. This B y profile corresponds to the ψ in LSC theory. In fact, as shown in Figure 19(a), these profiles have a local maximum point indicated by thick blue arrows. These profiles are plotted at the local maximum point of B y in the second Figure 19(a). The y-directional profiles of −By plotted at the same location as that in Fig.18(b). The local maximum points of By in the second plasmoid are indicated by thick blue arrows in each profile. The local maximum point of By at t = 308 is barely observed as the small light-blue region in Fig.13(a). plasmoid at each time. Hence, the x locations plotted at each time are exactly the same as those in Figure 18(b). First, the local maximum point of B y appears around y = 0.6 at t = 308, which corresponds to the small light-blue region in Figure 13(a). Then, the local maximum point gradually shifts to a larger y value from around y = 0.6. Finally, it grows and reaches y = 1.56 at t = 360, which is located around the outer edge of the growing plasmoid, i.e., the outer edge of the current sheet of t = 360 of Figure 18(b). Since LSC theory assumes that the local maximum point does not move, this MHD simulation is inconsistent with the theory at this point. Figure 19(b) shows how the local maximum point of −B y observed in Figure 19(a) moves in the y-direction. The purple solid line shows the y-directional movement with respect to time. The green and blue solid lines, respectively, show the movements of the locations of the one-half and one-quarter value of the maximum value of the −J z peaks measured in Figure 18(b). The dashed line shown only for 304 < t < 308 is the prediction line extrapolated from the purple solid line by which the generation point of the local maximum point of B y can be deduced. Since the beginning of the second tearing instability is during 304 < t < 308, the local maximum point of B y is deduced to be generated at y = 0.42 at t = 304. Hence, at the beginning of the tearing instability, the local maximum point of B y appears to be separated from the origin, i.e. y = 0.

The local maximum point of ∆V y around the X-point
Figure 20(a) shows the −V y profile in the y-direction at t = 308, 320, 340, and 360 at the second X-point. As the tearing instability grows, the plasma inflow toward the X-point is accelerated in 308 < t < 340. In fact, since the inflow speed V y takes a negative value, i.e., V y < 0, the −V y profile shown in Figure 20(a) tends to rise over time until t = 340. Then, the overall −V y profile for y > 1.5 starts to fall between t = 340 and 360. This fall indicates that the second tearing instability starts to terminate. Nevertheless, the second plasmoid still continues to grow beyond t = 340, as shown in Figure 19(a). Figure 19(b). The y-directional movements of the local maximum point of By. The movements of the half (green solid line) and quarter (blue solid line) thickness locations of the current sheet are indicated at the same x locations as measured in Fig.18(b). The local maximum point appears somewhere in 304 < t < 308 because it is not yet observed at t = 304 but is observed at t = 308, as shown in Figs.13(a) and 19(a). Hence, the generation point of the local maximum point predicted from the movement in t > 308 is around y = 0.42 at t = 304.
Thus, the reconnection process is maintained in the SP-like sheet beyond t = 340. In addition, the −V y profiles in 308 < t < 340 are slightly distorted around 0.4 < y < 1.0. This distortion is associated with the double current sheet structure observed in Figure  18(a). Then, the distortion disappears until t = 360; this disappearance appears to be associated with the disappearance of the double current sheet structure in Figure 18(a).
To adapt the V y profile observed in Figure 20(a) to φ in the modified LSC theory, at least, three problems must be considered. First, the largest problem is that V y does not have only the perturbed component φ but also the zero-order component defined by Eqs.(2.5) and (2.10). The perturbed component φ is difficult to rigorously extract from V y because at the beginning of the second tearing instability, i.e., t = 304 ∼ 308, the current sheet has already been in an unsteady state that is not the equilibrium. However, the growth of φ can be approximately measured by the increment of V y , i.e., ∆V y = V y (t, x, y)−V y (t−∆T, x, y), where ∆T = 4 is the data sampling time in the MHD simulation. Then, the growth rate Λ normalized by the MHD simulation time scale is calculated from ∆V y (t)/∆V y (t − ∆T ) = e Λ∆T , and then, translated to the growth rate λ ∆Vy normalized in the time scale of the modified LSC theory, as will be shown in Table  1.
The second and third problems, respectively, originate in the compressibility and moving of the second X-point in the MHD simulation. Because of the compressibility, the intensive region of −∆V y (∝ ∂V y /∂t) in Figure 15 is located around the plasmoid rather than around the X-point. This result is inconsistent with the assumption in incompressible LSC theory that the local maximum point of φ is located around the X-point. Adapting to the theory, we measure the local maximum value of ∆V y (t) around the X-point by using the −V y profiles shown in Figure 20(a). At this point, −∆V y around the X-point will not be seriously affected by the compressibility because the measurement point is close to the X-point. In addition, the third problem is that the X-point gradually moves from x = 78.9 to x = 155.1 during 308 < t < 360. At this point, ∆V y must be measured by an observer attached to the moving X-point. As a result, the ∆V y value is Figure 20(a). The y-directional profile of −Vy at the second X-point, where the x location is the same as that in Fig.18(a). The dent in the profiles observed around 0.4 < x < 1 and 308 < t < 340 is associated with that of the double current sheet structure shown in Fig.18(a). different from ∆V y in Figure 15, which was measured by an observer standing on the ground. Figure 20(b) shows the y-directional profiles of −∆V y , which is redefined as ∆V y = V y (t, x 1 (t), y)−V y (t−∆T, x 1 (t−∆T ), y), where x 1 is the x location of the moving X-point and ∆T = 4. Hence, this −∆V y is measured by an observer attached to the moving Xpoint. In Figure 20(b), the local maximum point of the redefined −∆V y is barely observed at t = 308. Then, the point for 312 < t < 360 moves in 0.3 < y < 0.6 and is always in the inner region of the current sheet shown in Figure 18(a). Then, during 340 < t < 360, the peak height rapidly decreases, indicating that the growth of the second tearing instability is terminated around t = 340.
Figure 20(c) shows how the local maximum point of −∆V y redefined in Figure 20(b) moves in the y-direction with respect to the thickness of the current sheet observed in Figure 18(a). The dashed line at 304 < t < 312 predicts the generation location of the local maximum point, which is extrapolated from the purple solid line observed at 312 < t. The local maximum point of −∆V y is predicted to be generated around x = 0.3 at t = 304. Hence, the local maximum point stagnates near the outer edge of the inner current sheet of the double sheet structure observed in Figure 18(a). At the beginning of the tearing instability, the local maximum point of −∆V y appears to be separated from the origin, i.e. y = 0.
3.3.5. The growth rates, λ By and λ ∆Vy , measured in the MHD simulation Figure 21 shows the time variations of B ymax and ∆V ymax , which are defined as the local maximum values of −B y and −∆V y , respectively, observed in Figures 19(a) and 20(b). Both B ymax and ∆V ymax monotonically increase at the beginning of the tearing instability. Then, B ymax continues to increase until t = 360, but ∆V ymax is saturated around t = 340 ∼ 350. The former growth is still maintained by the reconnection process based on the SP model, but the latter saturation indicates the termination of the second Figure 20(b). The y-directional profile of −∂Vy/∂t at the second X-point. This time differential is approximately measured as ∆Vy/∆T = (Vy(t, x1(t))−Vy(t−4, x1(t−4)))/4 with time interval ∆T = 4, where x1(t) is the location of the X-point in Fig.18(a). Hence, this figure corresponds to −∂Vy/∂t measured by an observer located at the moving X-point. In comparison to Fig.18(a), the local maximum points are always located inside the current sheet, i.e., y < 1.  Fig.20(b). The movements of the half and quarter thickness locations of the current sheet, which are measured in Fig.18(a), are also indicated. The local maximum point is not yet observed at t = 304. The generation point of the local maximum point predicted from the movement observed at t > 312 is around y = 0.3 at t = 304.  Fig.19(a), and the latter is measured in Fig.20(b). The former (i.e., Bymax) grows smoothly until the end of the second tearing instability, but the latter (i.e., ∆Vymax) is saturated at t = 340, which is around the end of the linear phase of the instability. tearing instability. In other words, around t = 340 ∼ 350, the second tearing instability switches from the linear phase to the nonlinear phase. Table 1 shows the numerical data plotted in Figure 21 and also the growth rates λ By and λ ∆Vy calculated from the numerical data. For example, λ By = 0.21 shown at t = 312 in Table 1 is calculated from B ymax (t = 312)/B ymax (t = 308) = −0.0014/ − 0.0007 = e 2πλ By (312−308)V A /lcs , where V A = 5.5 is the Alfven speed measured in the upstream region and l cs is the average value between t = 308 and 312, i.e., (l cs (t = 308) + l cs (t = 312))/2 = 41.5. As shown in Figure 13(b), l cs is measured as 4 times of the x-directional distance between the second X-point and the local maximum point of B y in the plasmoid left side of the second X-point, e.g., which is also indicated by the black arrow in Figure 13(a). At this point, since the plasmoid chain shown in Figure 13(b) is not exactly the sinusoidal shape assumed by linear theory, this l cs measurement is just an approximation, where B ymax is assumed to be located between X-point and O-point and the distance between the X-point and O-point is assumed to be 1/2 of the wave length of the plasmoid chain. Table 1 shows that l cs almost monotonically increases from 41 to 162 during 312 t 380, and λ By varies between 0.11 and 0.23 during 312 t 360. In 364 t 380, −J z at the X-point and λ By start to slowly decrease due to overelongation of the SP-like sheet, which is characterized by simultaneous increases in l cs and L cs . The λ By values can be compared with the upper limit λ up obtained in the modified LSC theory, as discussed below.
Similarly, λ ∆Vy shown in Table 1 is approximately calculated in the same manner as λ By . For example, λ ∆Vy = 0.27 shown at t = 316 in Table 1 is calculated from ∆V ymax (t = 316)/∆V ymax (t = 312) = −0.0287/ − 0.0124 = e 2πλ ∆Vy (316−312)V A /lcs . As shown in Table 1, λ ∆Vy takes its highest value at t = 312, i.e., the beginning of the tearing instability, and then monotonically decreases over time. Notably, λ ∆Vy takes negative values after t = 352; thus, the exponential growth of the tearing instability is  Table 1. The growth rates, λB y and λ∆V y , respectively, are measured from the time variations in the local maximum values, Bymax and ∆Vymax, shown in Fig.21. Jz is measured at X-point. How to measure lcs and Lcs is explained in Figs.13(b) and 17, respectively. Then, κ = πLcs/lcs is obtained. The upper limit λup is obtained from Fig.8, for either of the inner and outer ξ b , where = 0 is assumed. The inner ξ b (=inξ b ) and outer ξ b (=outξ b ) are measured for the inner and outer current sheet, respectively shown as labels B and C of Fig.18(a). terminated at this time. Since λ By > 0 is maintained even at t > 352, the reconnection process still continues to follow the SP model, maintaining the growth of the plasmoid. The λ ∆Vy values can be compared with λ up in the modified LSC theory, as shown next.
3.3.6. The upper limit λ up of growth rate predicted by the modified LSC theory This section presents how to compare the growth rates obtained from the MHD simulation and modified LSC theory, where it is shown that both growth rates, to some extent, appear to be consistent. First, note that λ up obtained from Figure 8 is just the upper limit of the growth rate. In the theory, the exact growth rate is unknown but generally depends on the upstream condition, such as zero-converging, zero-crossing and so on. The upstream condition in the theory is not easy to be matched to that of MHD simulation. Rather, let us focus on the fact that the exact growth rate is always smaller than λ up . As mentioned in Sections 2.3.4 and 5, λ up is determined from κ, , and ξ b , which can be measured in the MHD simulation result. In fact, κ, which is defined as the ratio of L cs and l cs , is shown in Table 1. The measurement of l cs is shown in Figure  13(b). The assumption of = 0 and measurements of L cs and ξ b are presented below.
As shown in Figure 17, L cs is obtained from ∂V x /∂x at the X-point and Alfven speed V A . This definition of L cs is essentially the same as that in the original LSC theory but the measurement method is evidently different. Because, LSC theory assumes that the SP sheet is in the steady state, but the MHD simulation in this section has not yet reached a steady state. In fact, the V x profile shown in Figure 17 does not reach V A = 5.5 at the maximum peak. Even in such a unsteady state, L cs can be measured. According to our measurement, Table 1 shows that L cs decreases from 150 to 50 during 308 t 356, which is associated with the fact that the magnetic diffusion region is localized by the tearing instability until t = 356. After t = 356, L cs starts to increase, which means that the magnetic diffusion region starts to elongate along the current sheet, resulting in the SP-like current sheet.
The current sheet thickness δ cs must be measured to determine = 2δ cs /L cs . As shown in Figure 18(a), the thickness is measured at approximately y = 2, that is, between y = 0 and the outer edge of the current sheet, and it is much smaller than the L cs shown in Table 1, resulting in < 0.1. According to Figure 7, λ up is not sensitive to when < 0.1 and κ < 2; hence, below we assume = 0. Figure 7 shows that λ up rapidly decreases in κ > 2. However, since λ up tends to decrease for larger , λ up obtained for = 0 will be available even for κ > 2. Then, λ up can simply be obtained from Figure 8 when κ and ξ b are known. By contrast, λ up is sensitive to ξ b .
The measurement of ξ b , where ξ b is equal to either the local maximum point's location of −∆V y or B y , is associated with some controversial problems to be studied in the future. In this paper, we employ −∆V y to measure ξ b . At this point, we can employ B y . However, as shown in Figures 19(b) and 20(c), the ξ b values obtained from B y and −∆V y are not drastically different; thus, we employ −∆V y to measure ξ b in Table 1.
Another problem to be studied in the future is that the current sheet consists of double sheets, as shown in Figure 18(a). Therefore, two choices, i.e., the inner or outer current sheet, are available to measure the current sheet thickness. Both cases are examined below.
First, we employ the thickness of the outer current sheet defined by label C in Figure  18(a), which is measured to be approximately y = 2. In addition, the local maximum point of −∆V y is located near y = 0.34 at t = 312, as shown in Figure 20(c). The resulting ξ b is roughly estimated to be 0.2 = 1.307 * 0.34/2. Note that the factor "1.307" originates from the fact that the outer edge of the current sheet in Eq.(2.4) is located at ξ = 1.307 on the ξ scale. The ξ b value listed as the outer ξ b in Table 1 varies between 0.2 and 0.4. Because ξ b < 1.307 at all times, this result represents the case of the inner-triggered tearing instability. Considering the κ listed in Table 1 and = 0, λ up < 0.1 is roughly estimated from Figure 8. Because λ up is the upper limit, λ up < 0.1 is too small to be consistent with the λ By and λ ∆Vy measured in the MHD simulation.
Second, we employ the thickness of the inner current sheet defined by label B in Figure  18(a), which is measured to be approximately y = 0.5. In the same manner as for the outer ξ b , the resulting ξ b varies between 0.7 and 1.4, as listed in the inner ξ b of Table 1. Since ξ b = 1.4 at t = 336 is slightly larger than 1.307, this scenario may be classified as the outer-triggered case. Considering that λ up is the upper limit, λ up for the inner ξ b is consistent with the λ By and λ ∆Vy measured in the MHD simulation. Finally, none of the growth rates obtained in the MHD simulation or modified LSC theory exceed unity, i.e., they are sub-Alfvenic. This result is reasonable for tearing instabilities driven by Alfven waves.
3.3.7. Relation of t s (linear phase) and t v (linear + nonlinear phases). Figures 12(a)-(c), the tearing instability is impulsively repeated, where each tearing instability grows and then slows during elongation of the SP-like sheet. Then, the next tearing instability starts in the over-elongated SP-like sheet. Evidently, there is a time interval t v for repeating each tearing instability and a duration time t s (< t v ) for which the modified LSC theory is applicable, i.e., a linear phase. For the second tearing instability, t v is roughly measured to be 380 − 308 = 72 because the second tearing instability starts at t = 308 and the third tearing instability starts at t = 380 which is not shown in this paper to reduce the number of Figures. Furthermore, t s is less than 348 − 308 = 40 because λ ∆Vy changes from positive to negative during 348 < t < 352, as shown in Table 1. Thus, the linear growth of the tearing instability has been terminated until t = 348. The factors that dominate t v are unclear, but t s will be regulated by two time scales, i.e., l cs /V A and l cs /C s , where C s is the local speed of sound in the current sheet. In other words, as Alfven and sound waves spread over the wavelength of the plasmoid chain (Shibata & Tanuma 2001), the tearing instability will shift from the linear phase to the nonlinear phase, where the linear phase is explained by the modified LSC theory and the nonlinear phase may be explained by the SP model. In this MHD simulation, the PK model is not observed in the nonlinear phase. According to Table 1, t s = l cs /V A varies between 41/5.5 = 7.5 at t = 312 and 84/5.5 = 15 at t = 360. In fact, as shown in the time variation of λ ∆Vy in Table 1, the second tearing instability, which started around t = 308, rapidly slows during 308 + 7.5 < t < 308 + 15. Then, as the SP-like sheet elongates along the current sheet, λ up for the inner ξ b and λ ∆Vy gradually separate. In contrast to λ ∆Vy , λ By does not drastically decrease because the reconnection process itself is maintained by SP-like sheet formation. This discussion is further continued in Section 4.1. Finally, beyond t = 364, λ By and −J z start to decrease, leading to the start of the third tearing instability. Since the MHD simulation studied in this paper is for the low beta plasma, the discussion based on t s = l cs /C s is not substantially different from that of t s = l cs /V A mentioned above.

As shown in
3.3.8. The first and third tearing instabilities.
Finally, let us briefly examine the first and third tearing instabilities observed in the MHD simulation. As shown in Figures 12(a)-(c), the growth of the first tearing instability at the origin in 4 < t < 280 is much slower than that of the second tearing instability in 308 < t < 348. According to the modified LSC theory proposed in this paper, this slow growth can be explained by the extremely large L cs measured at the origin. In fact, the ∂V x /∂x at the origin is smaller than that of the second tearing instability. The resulting κ is then larger, leading to the slower growth rate predicted from Figure 8. The details of the first tearing instability are summarized in Appendix D.
Conversely, the growth of the third tearing instability at 380 < t < 420 is faster than that of the second tearing instability, which is not shown in this paper to suppress the number of figures. However, as reported in many MHD simulations of PI (Samtaney, et.al. 2009;Bhattacharjee,et.al. 2009;Cassak & Drake 2009;Landi,et.al. 2015;Shimizu,et.al. 2017;Papini,et.al. 2019), as the tearing instability is repeated, the wave length l cs of the plasmoid chain gradually becomes shorter and the growth becomes faster. According to the modified LSC theory, this rapid growth can be explained by the small L cs . In addition, the relatively small l cs enhances the growth because the realistic growth rate Λ defined as B ymax (t + ∆T )/B ymax (t) = e 2πλ∆T V A /lcs = e Λ∆T tends to be larger for smaller l cs evenwhen λ, V A and ∆T are constant (Tajima & Shibata 2002). Then, as PI nonlinearly proceeds, l cs tends to be much smaller in the subsequently repeated tearing instabilities. Hence, the subsequent tearing instabilities tend to grow at a higher Λ evenwhen λ is constant. This is, as PI proceeds, why each single event of the tearing instability tends to be gradually faster in the realistic time scale of the MHD simulation.

Application of the modified LSC theory for PI
This section shows how the modified LSC theory proposed in Section 2 can be applied to PI. At this point, the modified LSC theory cannot be directly applied to PI. Furthermore, the original LSC theory (Loureiro, et.al. 2007) is also inapplicable; thus, the existence of the critical Lundquist number S c is not directly supported by those two theories. If PI actually can occur, it may have to be studied as a nonlinear process, which unfortunately is not considered in this paper. Figure 22 schematically shows a basic image of PI based on LSC theory. This image shows how a plasmoid chain appears in the SP sheet of length L cs /2, which is not necessarily equal to L cs /2 defined in this paper, i.e., Figure 17. At this point, note that Γ in Figure 22 is not necessarily equal to Γ 0 defined in Eqs.(2.5) and (2.6) because L cs defined in the macroscopic SP sheet may be different from L cs defined in each single event of tearing instability. As shown in Figure 22(a), the tiny plasmoid (magnetic island) generated around the origin gradually grows and propagates downstream, i.e., in the +x direction. This gradual growth is accompanied with the plasma outflow, which is shown as the solid oblique straight line in Figure 22 In general, as long as the amplitude of the perturbation solution is kept to be sufficiently weaker than that of the zero-order equilibrium, each trajectory x(t) of the X-points in the plasmoid chain will be traced as dx(t)/dt = V x0 = Γ x(t). This process results in x(t) = l cs e Γ t , where l cs is defined as the wavelength of the plasmoid chain generated around the origin. For example, Figure 22(c) shows four trajectories, i.e., for l cs = x 0 /4, x 0 /2, x 0 , and 2x 0 , where x 0 is defined as x 0 = L cs /(2e Γ T cs ) = L cs /(2e) = L cs /5.436, where T cs is defined as T cs = L cs /(2V A ). In other words, T cs is defined as the travel time from the origin to x = L cs /2 when a tiny plasmoid is assumed to move constantly at V A , as shown in the dashed oblique straight line in Figure 22(c).
According to Loureiro's paper (Loureiro, et.al. 2007), L cs = L cs , i.e. Γ = Γ 0 , is assumed. Then, Eqs.(2.1) and (2.2) are established when Γ = 2V A /L cs is sufficiently smaller than the growth time and l cs is almost constant in time. These conditions correspond to the case of l cs << L cs , e.g., l cs = x 0 /4 and x 0 /2 shown in Figure 22(c), i.e., l cs < x 0 = L cs /5.436. As shown in the x 0 /4 trajectory of Figure 22(c), since the movement of the X-point is very slow, the plasmoid remains in the SP sheet for a long time until it reaches x = L cs /2. In fact, as shown in Section 3, the X-point of every tearing instability remains almost still until each tearing instability is terminated by nonlinear saturation, i.e., the over-elongation of the SP sheet. Thus, t s = l cs /V A << T cs . Note that, as l cs becomes shorter, t s also tends to become shorter.
Two remarkable points are noted for the reality of the basic image of PI shown in Figure 22. First, as shown in Section 3, L cs and l cs strongly depend on each single event of tearing instability. For example, when the first tearing instability is terminated, the perturbed amplitude V x1 , which is shown as the dashed line in Figure 22(b) and corresponds to φ, becomes almost as large as the zero-order amplitude of V x0 , which is shown as the solid oblique straight line in Figure 22(b). In fact, as will be shown in Appendix D, the V x profiles in Figure 17 are largely different from those in Figure 27. Hence, L cs for the first tearing instability may be equal to L cs defined in Figure 22(a) but is evidently not applicable for the second tearing instability: L cs for the second and subsequent tearing instabilities must be individually measured. Furthermore, l cs must also be individually measured for each tearing instability. Hence, after the first tearing instability is terminated, the total length L cs of the macroscopic SP sheet does not affect how each subsequent tearing instability occurs. Hence, L cs is more important than L cs .
The second remarkable point is that, as mentioned in Section 2, λ up obtained for < 0.1 and κ < 5 is almost independent of , i.e., the resistivity. In addition, even for = 0, λ up takes a finite value less than unity. Note that the ξ scale is normalized by the current sheet thickness; thus, the tearing instability in the real time and space of the MHD simulations follows the similarity law. In other words, the MHD simulation based on finite resistivity is directly applicable to the much smaller resistivity case, where l cs and L cs are much smaller for the smaller resistivity, and the time intervals of t s and t v will also become much smaller. Moreover, note that the Alfven speed V A in the upstream magnetic field region and the local speed of sound C s in the current sheet do not change, even in the limit of = 0. Hence, the case of much smaller resistivity, i.e., much higher Lundquist number, will not essentially change how PI occurs, at least on the basis of the modified LSC theory. Thus, the existence of S c is not directly supported by the modified LSC theory, i.e., linear theory.

Comparison of the original and modified LSC theories for PI
The growth rate κλ(= γ/Γ 0 ) of the original LSC theory (Loureiro, et.al. 2007) is normalized by L cs defined in Figure 22, as shown in FIG.4 in the paper, while λ and λ up of the modified LSC theory, shown in Figure 8, are normalized by l cs , which will be much smaller than L cs . Therefore, the growth rate derived in the original LSC theory may be assumed to be constant during 0 < t < T cs , i.e., until the plasmoid generated around the origin reaches the exit of the macroscopic SP sheet at x = L cs /2. However, as shown in Table 1, the growth rate measured in the MHD simulation immediately slows at t > t s , where t s is much shorter than T cs . Thus, the growth rate κλ derived in the original LSC theory (Loureiro, et.al. 2007) appears to be overestimated for the PI application, at least, in terms of the MHD simulation shown in this paper.

Viscosity effect
Here, we discuss the viscosity effect, which is the largest difference between the theory and MHD simulation shown in this paper. Unfortunately, viscosity cannot be easily embedded in the modified LSC theory because Eqs.(2.1) and (2.2) are fairly complicated by viscosity. Moreover, the nonviscous MHD simulation results in a fatal numerical explosion due to the appearance of an extremely thin current sheet. Hence, rigorously, the modified LSC theory shown in this paper will be insufficient to be applied to the MHD simulations of PI.
It may be additionally noted that the nonviscous tearing instability is similar to the nonviscous turbulence in normal fluid dynamics, in which the simulation numerically fails because of the appearance of unlimitedly short wave length. Note that, since the thickness of the current sheet in LSC theory is normalized by δ cs , it can have unlimitedly thin current sheet in real space of (x, y) even in = 0. In other words, the real thickness of the steady state SP sheet can be unlimitedly thin even in = 0, i.e., the double limits of δ cs = 0 and L cs = 0. In addition, Figures 7, 10(a) and (b) suggest that such a unlimitedly thin current sheet always has a unstable κ range, in which λ > 0 is established (also, see Appendix C for the stability criterion), with exception of that LSC theory fails at κ = 0. Then, the finite viscosity can suppress the appearance of short wave length and will contribute to stabilize the tearing instability.
Hence, the introduction of viscosity is very important for the realistic linear theory of the tearing instability, and also MHD simulations of PI. There are some studies for the visco-resistive tearing instability (Tenerani,et.al. 2015;Betar,et.al. 2020), where the growth rate is observed to tend to decrease for a larger viscosity. Viscosity will be considered also in our next work.

FKR theory viewed from the perspective of modified LSC theory
The upper limit λ up of the growth rate in FKR theory can be explored in consideration of the limit of L cs = +∞ in the modified LSC theory, because, L cs = +∞ results in ∂V x /∂x = 0 at the X-point, which corresponds to the null-flow field of the equilibrium. Then, L cs = +∞ leads to κ = πL cs /l cs = +∞. According to Figure 8, λ up tends to decrease as κ increases; hence, λ up predicted in FKR theory will be close to zero. Thus, the tearing instability in FKR theory may be too slow in real cases. In other words, it may be important that the growth of the tearing instability is essentially dominated by the nonzero gradient ∂V x /∂x at the X-point, i.e. the nonzero gradient of outflow speed along the neutral sheet. Further discussion of FKR theory is provided at the end of Appendices C and E.

Summary
In this paper, LSC theory (Loureiro, et.al. 2007) under rigorously uniform resistivity has been examined, where inner and outer regions were seamlessly solved as an initial value problem. In other words, we did not use the ∆ index which is widely used to connect the inner and outer regions (Appendix E). Then, it was shown that the upper limit λ up of the growth rate can be estimated given the wave number κ, resistivity , and ξ b , which is related to the local maximum point location of all physically acceptable solutions, including the zero-converging and zero-crossing solutions. In addition, it was shown that the growth rate λ of the zero crossing solutions can be obtained from κ, , and ξ c , which represents the upstream condition, i.e. the zero-crossing point. These results for λ up and λ are summarized in Figure 8.
The MHD simulation presented in this paper is not new, i.e. it has been widely studied in many previous MHD studies. Rather, as shown in Table 1, here we showed that the modified LSC theory is applicable only to the beginning of each tearing instability, which exists in the linear phase. Then, when t s = l cs /V A is exceeded, the linear growth of the tearing instability terminates and enters the non-linear phase. In the MHD simulation presented in this paper, the nonlinear phase always resulted in the SP model, where the magnetic diffusion region was elongated along the current sheet. In other words, the nonlinear phase did not appear to result in the PK model, which is consistent with previous studies (Biskamp 1986). Then, after an interval of time t v from the beginning, the next tearing instability started in the over-elongated SP-like sheet. How t v is determined is unclear, but the MHD simulation in this paper suggests t v >> t s . If t v >> t s is always satisfied, the modified LSC theory is directly inapplicable for PI because PI should be fully nonlinear.
Most importantly, the modified LSC theory can predict tearing instabilities in the zero resistivity limit, i.e., S = +∞, as shown in Figure 8. In the limit, the current sheet becomes unlimitedly thin, but the tearing instability can still occur at a sub-Alfvenic growth rate, i.e., λ up < 1. Apparently, the zero resistivity "limit" case is different from the "perfect" zero resistivity case, in which the reconnection process is completely stopped. Additionally, it should be noted that the tearing instability essentially follows the similarity law in realistic time and space. Hence, the tearing instability in the zero resistivity limit will be similar to that in the finite resistivity case shown in Section 3. In other words, evenwhen S increases to +∞, no drastic changes appear in the tearing instabilities of PI, at least, in term of linear theory.

Appendix A. Zero-converging solutions in modified LSC theory
Moreover, let us extend the modified LSC theory. For φ and ψ to be physically acceptable, it is preferable that they converge to zero at ξ = +∞. However, in Section 2, we avoided finding such φ and ψ. In this section, let us attempt to find such zeroconverging solutions by modifying f (ξ). In LSC theory, f (ξ) is assumed to be constant for ξ > ξ 0 , but in this section, f (ξ) is modified to be constant for ξ > ξ 1 , where ξ 1 = ξ 0 (= 1.307) is not necessarily assumed. In this case, we can easily obtain the following analytical solutions for Eqs.(2.1) and (2.2) in ξ > ξ 1 .
Condition 1: Condition 2: Condition 3: where G 1 and G 2 are as follows.
G 1 = κ 2 2 + (κ/λ)(λ 2 + f 2 (ξ)), (A 6) We recognize that there are two types of solutions as for G 1 and G 2 . Hence, the general solutions will have the form of linear combinations of those types. However, let us only examine either type to confirm the existence of the zero-converging solutions. For C = E = 0, either of G 1 or G 2 for Eqs.(A3), (A4), and (A5) must be satisfied by adjusting λ, ξ 1 , and φ (0) for a given set of κ and . Figure 23 shows an example for G 2 . Hence, φ and ψ shown in Figure 23 simultaneously converges to zero at ξ = +∞, where f (ξ) is set to be constant in ξ > ξ 1 = 2.55625, as shown at the bottom panel of Figure  23. Inevitably, f (ξ) is discontinuous at ξ = 2.55625. Since Conditions 1, 2, and 3 can be exactly satisfied by adjusting λ, φ (0) and ξ 1 , the zero-converging perturbation solution actually exists.
Note that Figures 1-11 shown in Section 2 are inapplicable for this zero-converging solution because f (ξ) has been modified in ξ > ξ 0 . However, the growth rate λ = 0.868631 in Figure 23 is less than 1.0, i.e., it is sub-Alfvenic. In addition, the local maximum points of φ and ψ are, respectively, located at ξ = 1.6005 and 1.67275, i.e., between ξ 0 and ξ 1 , which is located in the reversed current sheet. This condition is the case for the outertriggered tearing instability.

Appendix B. Pinching method to find zero-crossing solution
This is a very primitive and steady technique, as shown in Figure 24. The format of Figure 24 is similar to the upper and middle panels of Figure 1(a) but this figure schematically shows how to numerically find the zero-crossing solution. In each panel of Figure 24, three lines labeled as "1", "2", and "3" correspond to φ and ψ solved for three different values of a control parameter, e.g., λ 1 , λ 2 , and λ 3 . In this case, the zero-crossing solution is established at λ = λ 2 , i.e., line "2", but it is impossible to exactly find line "2" by solving the numerical initial value problem because of inevitable numerical errors of the numerical study.
Rather, if lines "1" and "3" are found, and the transition from line "1" to "3" is assumed to be continuous, line "2" between lines "1" and "3" should actually exist, Because, as λ changes from λ 1 to λ 3 , the zero-crossing point of φ shifts from ξ φ1 to ξ φ3 and the zero-crossing point of ψ shifts from ξ ψ1 to ξ ψ3 . It is important that those zero-crossing points of φ and ψ shift in the opposite direction.
Appendix C. λ = 0 solutions of Eqs.(2.1) and (2.2) If there is a critical condition between stable and unstable modes for tearing instability, it will be worth finding φ and ψ solutions in λ = 0 limit. Taking the limit of λ = 0, Eqs.(2.1) and (2.2) are modified, as follows. ψ = (κ 2 2 + f (ξ)/f (ξ))ψ, (C 1) φ = −(f (ξ)/κf 2 (ξ))ψ (C 2) Eq.(C1) can be numerically resolved as the initial value problem with ψ(0) = 1 and ψ (0) = 0. Then, φ is directly obtained from Eq.(C2). Figure 25 shows the result for wave number κ = 5.4575 and resistivity = 0.2, which is a zero-crossing solution, because ψ = 0 at ξ = 5 and φ = 0 is constantly kept in ξ > 1.307. However, Eq.(C2) does not satisfy φ(0) = 0, because φ(0) diverges at the origin, as shown in the upper panel of Figure 25. In addition, φ generally has a discontinuity at ξ = 1.307, because f (ξ) has a discontinuity, which disappears for κ = 3.215 because ξ = 1.307 becomes the zero-crossing point of φ and ψ. It is important that, as long as κ = 1.0915 is kept, the zero-crossing solution is kept for any κ and , because Eq.(C1) is not changed. Hence, in = 0 limit, the zero-crossing solution shifts to κ = +∞. To compare with Figure Figures 26(a) and (b) are respectively obtained for φ (0) = 10 and 50. It may be expected that the zero-crossing solutions and zeroconverging solution are found in the extremely narrow range of λ, if they exist. At this point, it may be remembered that zero-converging solution was predicted between labels b and d in Figure 1(a), i.e., 0.2 < λ < 0.5. Unfortunately, due to the numerical error, we cannot deeply explore the extremely narrow range. The zero-crossing solutions deduced in Figures 26(a) and (b) appear to be close to that in Figure 25, suggesting that φ (0) = +∞ limit of Eqs.(2.1) and (2.2) corresponds to Figure 25. Those three figures are equally obtained for κ = 1.0915.
For κ < 1.0915, as shown in Section 2, we can easily find zero-crossing solutions while it is unclear whether zero-crossing solutions for κ > 1.0915 exist. Hence, it is unclear whether the tearing instability in κ > 1.0915 is stable or unstable. If it is stable, the critical condition of κ (= 2πδ cs /l cs ) = 1.0915 seems to be similar to that of FKR theory, Figure 25. φ (upper panel) and ψ (middle panel) of λ = 0 limit solution for κ = 5.4575 and = 0.2, where ψ has a zero-crossing point at ξ = 5 and φ = 0 is kept in ξ > 1.3069. Because of φ(0) = +∞ in the upper panel, this is not the physically acceptable solution that must be started from φ(0) = 0. In addition, φ has a discontinuity at ξ = 1.307 but is almost invisible in this figure. although the growth rate in FKR theory may be zero. Because, as discussed in Section 4.4, FKR theory corresponds to the L cs = +∞ limit of LSC theory. In other words, FKR theory corresponds to the double limits of = 2δ cs /L cs = 0 and κ = πL cs /l cs = +∞ in LSC theory.
Let us discuss about perspectives for FKR theory from LSC theory, moreover. In general, the perturbation theory must start from a rigorous equilibrium but the current sheet of finite thickness assumed in FKR theory cannot be equilibrium for the null flow field in the resistive MHD, i.e. = 0. There are the following strategies to avoid this problem.
First, if = 0 limit, i.e., ideal-MHD limit, is taken, the current sheet can be the rigorous equilibrium for the null flow field. In fact, with the exception of the δ cs = +∞ limit, the = 2δ cs /L cs = 0 limit is established for the L cs = +∞ limit which has been discussed above. In this case, the tearing instability in FKR theory will occur for κ (= 2πδ cs /l cs ) = 0 < 1.0915. However, as mentioned above, the growth rate will be close to zero. Note that the = 0 limit means when the thickness δ cs of the current sheet is extremely thin with respect to L cs .
Second, let us consider when the thickness δ cs of the current sheet is extremely thick, i.e., the case of δ cs = +∞ limit. Such a current sheet of infinite thickness seems to be a tricky consideration but ,in fact, it has been discussed (Cross & Hoven 1971) though potentially remaining a controversial topic. E.g., since Eq.(2.4) in the vicinity of the neutral sheet, i.e. ξ = 0, is almost a linear function with respect to ξ, if we focus only on the vicinity, the equilibrium of FKR theory can be rigorously established therein, in resistive MHD with the null flow field. As mentioned above, note that FKR theory   . The first tearing instability in 80 < t < 280. The figure format is the same as that in Fig.17. Lcs drastically changes between t = 80 (blue slope) and 200 (red slope), leading to the drastic change in the growth rate. requires L cs = +∞ limit in LSC theory. When the L cs = +∞ and δ cs = +∞ limits are simultaneously taken, = 2δ cs /L cs = 0 limit is not necessarily required to keep the rigorous equilibrium. Moreover, it leads to κ = πL cs /l cs = +∞ with the exception of l cs = +∞. It results in κ (= 2πδ cs /l cs ) > 1.0915. Hence, the tearing instability in FKR theory will be stabilized for the current sheet of infinite thickness.
Appendix D. The first tearing instabilities in MHD simulation Figure 27 shows the time variations of the x-profiles of V x along the x-axis, i.e., y = 0 in 80 < t < 280, in which the first tearing instability occurs, as shown in Figures 12(a) and (b). At this time, the X-point is constantly located at the origin. As shown in Figure  27, since the V x fluctuation observed during 0 < x < 30 at t = 80 is extremely weak, the ∂V x /∂x measured at the origin results in approximately L cs = 2000, which is much larger than that of the second tearing instability shown in Table 1. Moreover, since the weak V x fluctuation is strongly localized around the origin, l cs at t = 80 is approximately 30, which is not excessively large. Hence, κ reaches approximately 200, which is very large. Unfortunately, Figure 8 does not cover the κ range, but the predicted λ up will be extremely low. Hence, the first tearing instability develops slowly initially, i.e., t < 160. However, as the first tearing instability becomes fully developed at t > 160, the predicted growth rate accelerates because L cs drastically decreases to 200, resulting in a substantial decrease in κ. As a result, the first tearing instability slowly proceeds in t < 160 and gradually accelerates for 160 < t < 240. Then, the instability is terminated around t = 240, where the L cs switches from decreasing to increasing in time. This switching is the same as that observed at t = 356 in Table 1. Finally, this process leads to the formation of an elongated SP-like sheet, as shown in Figure 12(b).