Instabilities in a current sheet with plasma jet

For nearly one century, scientists have been investigating the mechanisms that cause space weather events such as the magnetic storms (e.g. Chapman & Ferraro, 1940; Ferraro, 1952). Magnetic reconnection is identified as one of the most important processes in space plasmas that drives various explosive phenomena, such as the solar flares (e.g. Masuda et al., 1994), coronal mass ejections (e.g. Gosling et al., 1995), and the magnetic substorms in the Earth’s magnetotail (e.g. Angelopoulos et al., 2008). It efficiently converts the magnetic energy in a current sheet to the kinetic and internal energies of the plasma through re-configuration of the topology of the magnetic field and connectivity of the magnetic field lines. In the laboratory, magnetic reconnection destabilizes the plasma (e.g. Yamada et al., 1994; Dorfman et al., 2013) and is fatal to stable controlled fusion. Understanding of how the reconnection triggers and evolves is crucial for a complete description of how energy is stored and released in different plasma environments.

As the plasma is a multi-scale system, various types of waves and instabilities exist on the largest magnetohydrodynamic (MHD) scale, the intermediate ion-kinetic scale, and the smallest electron-kinetic scale. Consequently, magnetic reconnection is a multi-scale process. In the past two decades, numerical simulations (e.g. Daughton et al., 2006; Guo et al., 2015; Cassak et al., 2017; Lu et al., 2019) as well as in-situ measurements of the space plasma (e.g. Burch et al., 2016; Torbert et al., 2018) have greatly enhanced our understanding of the microscopic kinetic physics of reconnection in the ion diffusion region and electron diffusion region. But it still remains unclear as how fast reconnection triggers in the macroscopic current sheets whose dimensions are much larger than any ion kinetic scales, e.g. current sheets in the pre-flare configurations in the solar corona.

In ideal-MHD regime, the magnetic field lines are “frozen in” the plasma and change of their connectivity is prohibited. Thus, a certain mechanism that breaks the ideal-MHD condition must play a role for the reconnection to happen. In macroscopic current sheets, this mechanism is either collision-induced resistivity, or some kind of effective resistivity caused by microscopic wave-particle interactions (e.g. Büchner & Elkina, 2006; Ma et al., 2018). Hence, the triggering problem of reconnection at MHD scales is essentially the stability problem of the resistive current sheet. Since 1960s, many works have been conducted on the resistive instability, i.e. the so-called “tearing instability”, of the current sheet (Furth et al., 1963; Coppi et al., 1966). The tearing mode grows with the help of resistivity that transfers the magnetic energy stored in the shear magnetic field to the growing perturbations, leading to the formation of a chain of plasmoids. Considering an infinitely long current sheet with thickness $a$, one can define the dimensionless Lundquist number $S=aV_{A}/\eta$, where $V_{A}=B/\sqrt{\mu_{0}\rho}$ is the characteristic Alfvén speed with $\rho$ and $B$ being the plasma density and asymptotic magnetic field strength, and $\eta$ is the magnetic diffusivity, which is resistivity divided by the permeability $\mu_{0}$. For simplicity, we will refer to $\eta$ as “resistivity” hereinafter. Linear theory predicts that the most unstable tearing mode has a growth rate $\gamma\tau_{a}\sim S^{-1/2}$ where $\tau_{a}=a/V_{A}$ is the Alfvén crossing time. This relation implies a faster growth of the instability with larger resistivity. As plasma in most of the space environments and laboratories is weakly collisional (Ji & Daughton, 2011; Pucci et al., 2017), the tearing mode seems to grow at very slow speed.

However, considering a two-dimensional current sheet whose length is $L$, in most pre-reconnection configurations, its aspect ratio can be very large ($L\gg a$). In this case, we should use $L$ instead of $a$ to measure the growth rate of tearing instability. After re-defining the Alfvén crossing time and Lundquist number such that $\tau_{L}=L/V_{A}$ and $S_{L}=LV_{A}/\eta$, it can be shown that the maximum growth rate of tearing mode is $\gamma\tau_{L}\sim S_{L}^{-1/2}\times(a/L)^{-3/2}$ (Pucci & Velli, 2013), implying that the aspect ratio of the current sheet is an important factor in determining how fast the mode grows. The growth rate can be extremely large at low-resistivity limit ($S_{L}\rightarrow\infty$) if the current sheet is thinner than a critical value $a/L\sim S_{L}^{1/3}$. Especially, in the classic model for steady reconnection with resistivity, i.e., the Sweet-Parker type current sheet, whose aspect ratio is $a/L\sim S_{L}^{-1/2}$ (Sweet, 1958; Parker, 1957), the maximum growth rate of tearing is $\gamma\tau_{L}\sim S_{L}^{1/4}$ (Tajima & Shibata, 2018; Loureiro et al., 2007). The positive power-law index means that the growth rate can be extremely large in the limit $S_{L}\rightarrow\infty$ (Bhattacharjee et al., 2009; Huang & Bhattacharjee, 2013).

The above analysis leads to a plausible scenario of the explosive energy release of the macroscopic current sheet. Initially, the current sheet is thick with $(a/L)>S_{L}^{-1/3}$, and thus is stable to tearing mode. Then some external forces gradually build up magnetic energy and result in thinning of the current sheet. Once the current sheet thins to the critical aspect ratio $(a/L)\sim S_{L}^{-1/3}$, the growth of the tearing mode suddenly becomes very fast, and the current sheet breaks up into many plasmoids and smaller-scale current sheets. This process can happen recursively in the newly formed current sheets and dissipates the magnetic energy rapidly (Shibata & Tanuma, 2001; Tenerani et al., 2015b; Landi et al., 2015; Papini et al., 2019), until it is terminated due to the decrease of Lundquist number (Shi et al., 2018) or the ion kinetic effect (Shi et al., 2019; Bora et al., 2021). Thus, tearing instability is an important and fundamental mechanism that facilitates fast reconnection in the large-scale current sheets. Consequently, it is important to thoroughly study it under different configurations. Recent progresses on this topic include the calculations of its linear growth rate with viscosity (Tenerani et al., 2015a), different background magnetic field profiles (Pucci et al., 2018), Hall effect (Pucci et al., 2017), guide field (Shi et al., 2020), ion-neutral collisions (Pucci et al., 2020), and normal component of magnetic field (Shi et al., 2021).

In space plasma, current sheets are frequently accompanied by plasma flows. In the dayside magnetosheath, reconnection events are often observed within highly turbulent plasma (e.g. Huang et al., 2016), and so for the solar wind (e.g. Osman et al., 2014). As a result, the reconnecting current sheets are likely to be affected by plasma flows of all directions. Plasma flows are also detected in the pre-flare corona (e.g. Wallace et al., 2010) and nightside magnetotail current sheet (e.g. Lane et al., 2021). At the tip of the helmet streamer where the heliospheric current sheet forms, growth of tearing instability accompanied by an outward propagating solar wind stream is observed in MHD simulations (Réville et al., 2020, 2022). Thus, study of how the tearing mode instability is modified by plasma flows is necessary. Many works have been conducted on the effect of a shear flow, i.e., flow parallel to the shear magnetic field (e.g. Hofman, 1975; Paris & Sy, 1983; Einaudi & Rubini, 1986; Chen & Morrison, 1990; Ofman et al., 1991; Paris et al., 1993; Dahlburg et al., 1997; Chen et al., 1997; Faganello et al., 2010). For example, Chen et al. (1997), through 2D MHD simulations, show that a sub-Alfvénic shear flow stabilizes tearing mode, and with a super-Alfvénic shear flow the instability is dominated by Kelvin-Helmholtz mode. This result is confirmed by a recent work (Shi et al., 2021) that calculates the linear instability growth rate by a boundary-value-problem approach. They also show that when the flow is exactly Alfvénic, the current sheet is extremely stable and the perturbation only grows at the rate of diffusion. Compared with the shear flow, the case where a plasma jet exists at the center of the current sheet is more complicated, because the jet itself is susceptible to two types of streaming instabilities, i.e., the sausage (varicose) mode and the kink (sinuous) mode. Wang et al. (1988a), using an initial value solver of compressible MHD equations, show that a super-Alfvénic plasma jet can increase the growth rate of the tearing mode. Subsequent works (Wang et al., 1988b; Lee et al., 1988) show that under this type of configuration both the kink mode and sausage mode exist and the sausage mode mixes with the tearing mode in the presence of resistivity. 2D MHD (Bettarini et al., 2006) and Hall-MHD simulations (Hoshino & Higashimori, 2015) confirm these early results.

In this study, we carry out a comprehensive investigation of the stability problem of the current sheet with a plasma jet in the framework of linear incompressible MHD, using an eigenvalue problem solver. We examine both the streaming sausage mode and streaming kink mode with and without resistivity. We derive the controlling equation set under a generalized configuration such that the plasma jet can have arbitrary angle with respect to the magnetic field and a finite guide field is allowed. The paper is organized as follows. In section 2, we describe the background fields used in this study. In section 3, we derive the equation set for the perturbation field. In section 4 we present the detailed results of our calculation. In section 5 we summarize the results and discuss the possible applications of the results to space plasma.

We start from the resistive-MHD equation set

	$$\frac{\partial\rho}{\partial t}+\nabla\cdot\left(\rho\mathbf{V}\right)=0$$		(1a)
	$$\rho\frac{\partial\mathbf{V}}{\partial t}+\rho\mathbf{V}\cdot\nabla\mathbf{V}=-\nabla P+\mathbf{J}\times\mathbf{B}$$		(1b)
	$$\frac{\partial\mathbf{B}}{\partial t}=\nabla\times\left(\mathbf{V}\times\mathbf{B}\right)+\frac{1}{S}\nabla^{2}\mathbf{B}$$		(1c)
	$$\frac{\partial P}{\partial t}+\mathbf{V}\cdot\nabla P+\kappa\left(\nabla\cdot\mathbf{V}\right)P=0$$		(1d)

where $\rho,\mathbf{V},\mathbf{B},P$ are density, velocity, magnetic field, and pressure, $\kappa$ is the adiabatic index, and $S$ is the Lundquist number. Incompressiblility is assumed throughout the study, i.e., $\rho(t,\mathbf{x})\equiv\rho_{0}$ is a constant. In a generalized configuration, the background magnetic field and velocity consist of both $x$ and $z$ components but are functions of $y$ only:

$$\mathbf{B}=B_{x}(y)\hat{e}_{x}+B_{z}(y)\hat{e}_{z},\quad\mathbf{V}=V_{x}(y)\hat{e}_{x}+V_{z}(y)\hat{e}_{z}.$$

(2)

Consequently, the momentum convection term and the magnetic tension force are both zero: $\mathbf{V}\cdot\nabla\mathbf{V}\equiv 0,\,\mathbf{B}\cdot\nabla\mathbf{B}\equiv 0$, and the divergences of $\mathbf{V}$ and $\mathbf{B}$ are also zero. The zeroth-order scalar pressure ensures the pressure balance

$$P(y)=P^{T}-\frac{B^{2}(y)}{2\mu_{0}},$$

(3)

where $P^{T}$ is the uniform total pressure. With the above configuration, the background field is in equilibrium without resistivity, but will diffuse with a finite resistivity. However, in most of the space and laboratory plasmas, the resistivity is extremely small, thus the diffusion time is much longer than the growth time of instabilities of interest. Therefore, we are able to neglect the diffusion of background fields.

In this study, we adopt the Harris type current sheet model for the magnetic field such that $B_{x}(y)=B_{0}\tanh(y/a)$. We also allow a uniform guide field $B_{z}(y)=B_{g}$. We assume the velocity is of the following form:

$$\mathbf{V}=V(y)\left(\cos(\alpha)\,\hat{e}_{x}+\sin(\alpha)\,\hat{e}_{z}\right)$$

(4)

i.e., the jet rotates from the $x$ direction by an angle $\alpha$ (Figure 1a). We adopt a flow function:

$$V(y)=V_{0}\,\mathrm{sech}^{2}\left(\frac{y}{d}\right),$$

(5)

where $d$ is the half-thickness of the jet, and $V_{0}$ is the flow speed at the center of the current sheet ($y=0$) and is a variable parameter. Throughout the study, we fix $d=a$, i.e., the width of the jet is the same as the width of the current sheet. The profile of $V(y)$ is plotted in Figure (1b) together with the profile of the $x$-component of the magnetic field $B(y)$.

For the perturbations, we consider a Fourier mode whose growth rate $\gamma$ is a complex number and wave vector $\mathbf{k}$ is in the $x-z$ plane with an arbitrary angle $\theta$ with respect to the $x$ direction (Figure 1a):

$$\mathbf{k}=k\cos(\theta)\,\hat{e}_{x}+k\sin(\theta)\,\hat{e}_{z}.$$

(6)

Hence, the perturbation fields have the form:

$$\left(\begin{array}[]{c}\mathbf{u}(t,\mathbf{x})\\ \mathbf{b}(t,\mathbf{x})\\ p(t,\mathbf{x})\end{array}\right)=\left(\begin{array}[]{c}\mathbf{u}(y)\\ \mathbf{b}(y)\\ p(y)\end{array}\right)\exp(\gamma t+i\mathbf{k}\cdot\mathbf{x})$$

(7)

For simplicity, one can rotate the $x-z$ coordinates with respect to the $y$-axis by angle $\theta$ and get a new coordinate system $\tilde{x}-\tilde{z}$ (Figure 1a) such that $\mathbf{k}=k\hat{e}_{\tilde{x}}$, and thus there is $\partial_{\tilde{z}}\equiv 0$. The background magnetic field and velocity can then be written as $\mathbf{B}=B_{\tilde{x}}(y)\,\hat{e}_{\tilde{x}}+B_{\tilde{z}}(y)\,\hat{e}_{\tilde{z}}$ and $\mathbf{V}=V_{\tilde{x}}(y)\,\hat{e}_{\tilde{x}}+V_{\tilde{z}}(y)\,\hat{e}_{\tilde{z}}$ after projection to the new coordinate system.

The next step is to derive a closed linear equation set for the eigenvalue problem. We start from the linearized momentum equation (with uniform density $\rho$)

$$\gamma\mathbf{u}+\mathbf{V}\cdot\nabla\mathbf{u}+\mathbf{u}\cdot\nabla\mathbf{V}=-\frac{1}{\rho}\nabla p_{1}^{T}+\left(\mathbf{B}\cdot\nabla\mathbf{b}+\mathbf{b}\cdot\nabla\mathbf{B}\right),$$

(8)

where we have normalized the magnetic field by $\sqrt{\mu_{0}\rho}$ so that it is in the unit of speed. To get rid of the 1st-order pressure, we can take the curl of the equation and get

$$\gamma\nabla\times\mathbf{u}+\nabla\times\left(\mathbf{V}\cdot\nabla\mathbf{u}+\mathbf{u}\cdot\nabla\mathbf{V}\right)=\nabla\times\left(\mathbf{B}\cdot\nabla\mathbf{b}+\mathbf{b}\cdot\nabla\mathbf{B}\right).$$

(9)

Using $\nabla\cdot\mathbf{u}=0$ and $\nabla\cdot\mathbf{b}=0$ to eliminate $u_{\tilde{x}}$ and $b_{\tilde{x}}$, the $\tilde{z}$ component of the above equation can be re-arranged in the following form:

$$\gamma\left(u_{y}^{\prime\prime}-k^{2}u_{y}\right)+ik\left[V_{\tilde{x}}\left(u_{y}^{\prime\prime}-k^{2}u_{y}\right)-V_{\tilde{x}}^{\prime\prime}u_{y}\right]=ik\left[B_{\tilde{x}}\left(b_{y}^{\prime\prime}-k^{2}b_{y}\right)-B_{\tilde{x}}^{\prime\prime}b_{y}\right]$$

(10)

where the prime indicates $\partial_{y}$. We note that the above equation contains only the $y$-component of $\mathbf{u}$ and $\mathbf{b}$. The linearized induction equation is

$$\gamma\mathbf{b}=\mathbf{b}\cdot\nabla\mathbf{V}-\mathbf{V}\cdot\nabla\mathbf{b}+\mathbf{B}\cdot\nabla\mathbf{u}-\mathbf{u}\cdot\nabla\mathbf{B}+\frac{1}{S}\nabla^{2}\mathbf{b}$$

(11)

where $S=aV_{A}/\eta$ is defined with the half-thickness of the current sheet $a$ and the upstream Alfvén speed $V_{A}=B_{0}/\sqrt{\mu_{0}\rho}$. The $y$-component of the linearized induction equation is

$$\gamma b_{y}=-ikV_{\tilde{x}}b_{y}+ikB_{\tilde{x}}u_{y}+\frac{1}{S}\left(b_{y}^{\prime\prime}-k^{2}b_{y}\right)$$

(12)

which also contains only the $y$-component of $\mathbf{u}$ and $\mathbf{b}$. Thus, equations (10) and (12) form a closed equation set for $u_{y}$ and $b_{y}$:

	$$\gamma\left(u_{y}^{\prime\prime}-k^{2}u_{y}\right)+ik\left[V_{\tilde{x}}\left(u_{y}^{\prime\prime}-k^{2}u_{y}\right)-V_{\tilde{x}}^{\prime\prime}u_{y}\right]=k\left[B_{\tilde{x}}\left(b_{y}^{\prime\prime}-k^{2}b_{y}\right)-B_{\tilde{x}}^{\prime\prime}b_{y}\right]$$		(13a)
	$$\gamma b_{y}=-ikV_{\tilde{x}}b_{y}-kB_{\tilde{x}}u_{y}+\frac{1}{S}\left(b_{y}^{\prime\prime}-k^{2}b_{y}\right)$$		(13b)

Here we have assimilated a $\pi/2$ phase difference between $u_{y}$ and $b_{y}$, i.e., we have replaced $ib_{y}$ with $b_{y}$.

One can immediately find that, if $V_{\tilde{x}}=0$, i.e., if $\mathbf{k}\cdot\mathbf{V}=0$, the system is purely determined by $B_{\tilde{x}}$, similar to the classic tearing case, and the growth rate is purely real. As an example, consider an anti-parallel magnetic field $\mathbf{B}=B_{x}(y)\hat{e}_{x}$ and an out-of-plane flow $\mathbf{V}=V_{z}(y)\hat{e}_{z}$. In this case, if $\mathbf{k}=k\hat{e}_{x}$, the system reduces to the classic tearing, i.e., the flow has no effect on the solution. But if the wave vector is not along $x$ and has a finite $z$-component, the flow will alter the growth rate and introduce an oscillation (a non-zero frequency) to the solution. Similarly, if $\mathbf{k}\cdot\mathbf{B}=0$, the system is determined purely by the flow $V_{\tilde{x}}$ and there are only stream-induced instabilities.

As a final remark, we note that Equation (13) is in generalized form, and works for any functions $\mathbf{V}(y)$ and $\mathbf{B}(y)$ once they have the form of equation (2) and the system is incompressible.

Equation (13) is a boundary-value-problem (BVP) with the boundary conditions $u_{y}(y\rightarrow\pm\infty)=0,b_{y}(y\rightarrow\pm\infty)=0$. Far from the center of the current sheet ($y\rightarrow\pm\infty$), the derivatives of the background fields reduce to zero. Consequently, one can see that $u_{y},b_{y}\propto\exp(-k|y|)$ satisfy equation (13a), and equation (13b) just gives the ratio $u_{y}/b_{y}$ at $y\rightarrow\pm\infty$. In this study, we use the numerical BVP solver implemented in the Python package SciPy (Virtanen et al., 2020) to solve equation (13). In practice, we set the boundaries at $y=\pm 15a$, which are large enough to acquire accurate solutions.

Before solving the equation set, we need to define the parities of $u_{y}$ and $b_{y}$ first. Given that $V_{\tilde{x}}(y)$ is an even function and $B_{\tilde{x}}(y)$ is an odd function, one can see from equation (13) that there are two possible combinations of the parities of $u_{y}$ and $b_{y}$: (1) $u_{y}(y)$ is odd and $b_{y}(y)$ is even, (2) $u_{y}(y)$ is even and $b_{y}(y)$ is odd. The first case is the so-called “sausage” (varicose) mode, which leads to the formation of a chain of blobs and plasmoids (if resistivity is non-zero). Tearing instability is categorized to the sausage mode. The second case is the “kink” (sinuous) mode, which results in wavy distortion of the current sheet and stream. In solving the problem, we carefully search for both of the two modes.

In this study, we adopt a boundary-value-problem solver to study the instabilities inside a current sheet with the presence of a plasma jet. When the jet is collimated with the anti-parallel component of the magnetic field, both of the sausage mode and kink mode can be stabilized by the magnetic field. Without resistivity, the stability thresholds for the kink mode and sausage mode are $V_{0}/B_{0}\approx 2$ and $V_{0}/B_{0}\approx 1$ respectively (Figure 3). With a finite resistivity, the streaming sausage mode couples with the tearing mode, but the streaming kink mode is not modified by the resistivity significantly unless the resistivity is very large ($S<100$). Thus, in most of space and laboratory current sheets where $S$ is extremely large, the kink mode can be excited only if the jet speed is large ($V_{0}/B_{0}\gtrsim 2$). For $V_{0}/B_{0}\lesssim 1$, the sausage mode is tearing-like, with a power-law relation between the maximum growth rate and the Lundquist number max$(\gamma)\propto S^{-1/2}$ in the large $S$ limit, same as the tearing mode without flow, while the values of the maximum growth rate increase with the ratio $V_{0}/B_{0}$. For $V_{0}/B_{0}\gtrsim 1$, the sausage mode gradually transits to more streaming-like, and the maximum growth rate becomes less dependent on $S$ (Figure 4). In the case of a jet flowing along the direction perpendicular to the anti-parallel component of the magnetic field, our result reveals that, once the jet speed exceeds a threshold which is determined by the Lundquist number, the maximum growth rate of the sausage mode may increase with the angle between the wave vector and the reconnecting magnetic field component (Figure 6). This is because the mode transits from the pure tearing to pure streaming as the wave vector rotates from the anti-parallel magnetic field direction to the jet direction. Last, the out-of-plane jet combined with a finite guide field leads to a complex behavior of the maximum growth rate of the sausage mode. With certain $V_{0}/B_{0}$ and $\theta$ values (panel (c) of Figure 7 and Figure 8), the maximum growth rate increases with the guide field strength. But the increase is not very large and is non-monotonically dependent on $\theta$. For example, panel (c) of Figure 7 shows that the increase in max$(\gamma)$ with $B_{g}$ from $B_{g}/B_{0}=0$ to $B_{g}/B_{0}=0.5$ is larger for $\theta=30^{\circ}$ and 45^∘ than that for $\theta=40^{\circ}$. More importantly, overall the guide field quenches both the oblique tearing mode and the streaming sausage mode. As a result, increasing $B_{g}$ will gradually turn the monotonically increasing max$(\gamma)$-$\theta$ curve to a curve that increases at first and then drops (Figure 8).

These results indicate that plasma flow plays an important role in destabilizing the current sheets in space and laboratory plasma. A jet whose width is comparable to that of the current sheet and peak speed similar to the upstream Alfvén speed can enhance the maximum growth rate of the tearing mode to more than twice of that in the no-flow case (Figure 3). When the jet has a finite component along the direction perpendicular to the anti-parallel component of the magnetic field, even if the component is much smaller than the upstream Alfvén speed, the oblique sausage mode ($\theta>0$) may have comparable or even larger growth rate than the parallel sausage mode ($\theta=0$), and the most unstable mode may be perpendicular ($\mathbf{k}=k\hat{e}_{z}$) (Figure 6). The reason is that the out-of-plane flow (along $z$ direction) does not feel the stabilization effect by the magnetic field along $x$, and the growth rate of the pure streaming sausage mode is usually much larger than the pure tearing mode in the large $S$ limit. When the out-of-plane jet and guide field coexist, the most unstable mode may be oblique rather than parallel or perpendicular (Figure 8).

We note that several factors which are absent in this study may have non-negligible effects on the analyzed instabilities. Here we assume a uniform density profile and incompressibility. However, compressible MHD simulations show that a nonuniform background plasma density such as in the magnetotail can modify the growth rate of both tearing and streaming modes (Hoshino & Higashimori, 2015). In addition, if Hall effect is included, out-of-plane components of the magnetic field and velocity perturbations are generated even for the parallel mode ($\mathbf{k}=k\hat{e}_{x}$). Therefore, the out-of-plane jet will modify both the oblique and the parallel modes. Moreover, different widths of the jet and current sheet will change the results (Einaudi & Rubini, 1986; Hoshino & Higashimori, 2015). As a final remark, it is worth noting that in collisionless regime where the electron inertia is the only mechanism that breaks the frozen-in condition, an out-of-plane plasma jet plays a stabilizing role of the tearing mode even if the mode is parallel ($\mathbf{k}=k\hat{e}_{x}$) (Tassi et al., 2014). This is very different from the resistive-MHD regime where the out-of-plane jet only modifies the oblique tearing mode.

Acknowledgements: The author thanks Prof. Marco Velli and Dr. Kun Zhang for very helpful suggestions and comments, and the SciPy team for implementing the boundary-value-solver in Python (Virtanen et al., 2020). The work was supported by NASA HERMES DRIVE Science Center grant No. 80NSSC20K0604.