1 Introduction
In magnetohydrodynamics (MHD), the tearing instability (TI) occurs in highly sheared magnetic field configurations called current sheets. In a current sheet there is a thin (compared to larger length scales outside the current sheet) layer of intense current density where the magnetic field changes direction rapidly. If the conditions of the TI are met, the current sheet begins to ‘tear’ or, to be more precise, the topology of the magnetic field changes to form multiple islands (or plasmoids in three dimensions) of magnetic flux. Since the seminal work of Furth, Kileen & Rosenbluth (Reference Furth, Kileen and Rosenbluth1963), the onset of the TI has been traditionally studied using normal mode analysis, to the extent that the terms ‘tearing instability’ and ‘tearing mode’ are often used synonymously.
Recent studies that address the linear onset of TI in high aspect ratio current sheets, also known as the plasmoid instability (PI), (e.g. Loureiro, Schekochihin & Cowley Reference Loureiro, Schekochihin and Cowley2007; Bhattacharjee et al.
Reference Bhattacharjee, Huang, Yang and Rogers2009; Pucci & Velli Reference Pucci and Velli2014; Tenerani et al.
Reference Tenerani, Velli, Pucci, Landi and Rappazzo2016; Uzdensky & Loureiro Reference Uzdensky and Loureiro2016) do so from the point of view of normal mode analysis. For the TI (and the PI), however, normal mode analysis cannot give a complete picture of its linear onset. The operators in the equations describing the onset of the TI are non-normal. This means that the eigenmodes are not orthogonal and for the application in hand are heavily ill conditioned (Borba et al.
Reference Borba, Riedel, Kerner, Huysmans, Ottaviani and Schmid1994). Therefore, although eigenmodes may be damped as
$t\rightarrow \infty$
, they can result in significant transient (or algebraic) growth within a finite time. Performing normal mode analysis on equations with non-normal operators results in the translation to a later time when the transient growth has been damped away. Therefore, if significant transient growth is possible, it is ignored in normal mode analysis.
Although stability theory in plasma physics is dominated by normal mode (eigenvalue) analysis, studies of non-modal behaviour are on the increase. In the MHD literature, one early suggestion that subcritical behaviour may be important for the tearing instability was made by Dahlburg et al. (Reference Dahlburg, Zang, Montgomery and Hussaini1983), although the mechanism was thought to be nonlinear rather than linear. Later, Dahlburg (Reference Dahlburg1994) studied the algebraic growth of current sheets in ideal MHD as a possible route to turbulent reconnection through the creation of smaller length scales. Borba et al. (Reference Borba, Riedel, Kerner, Huysmans, Ottaviani and Schmid1994) investigated the eigenmodes of resistive MHD using pseudospectra but did not focus on the TI. They argue that the non-orthogonality of the eigenmodes implies that normal mode analysis can only describe instability growth on a long time scale (of
$O(S^{1/2})$
, where
$S$
is the Lundquist number that will be defined later). Other researchers have recognized the importance of non-modal growth in other MHD applications, including kinematic dynamo theory (e.g. Farrell & Ioannou Reference Farrell and Ioannou1999a
,Reference Farrell and Ioannou
b
; Livermore & Jackson Reference Livermore and Jackson2006; Chen et al.
Reference Chen, Herreman, Li, Livermore, Luo and Jackson2018), the magnetorotational instability (e.g. Squire & Bhattacharjee Reference Squire and Bhattacharjee2014a
,Reference Squire and Bhattacharjee
b
) and the tearing instability (MacTaggart & Stewart Reference MacTaggart and Stewart2017). There is also a growing interest in the subcritical transition to turbulence in tokamak plasmas (e.g. Landremann, Plunk & Dorland Reference Landremann, Plunk and Dorland2015; van Wyk et al.
Reference van Wyk, Highcock, Schekochihin, Roach, Field and Dorland2016) and the non-modal consequences of shearing on microinstabilities (e.g. Newton, Cowley & Loureiro Reference Newton, Cowley and Loureiro2010).
The purpose of this article is to investigate the non-modal transient growth at the onset of the classical TI. In particular, our aim is to determine the dependence of the maximum possible transient growth on the Lundquist number and to understand the relationship between the eigenvalue spectrum and the underlying physics. We solve the linearized equations numerically and use pseudospectral theory to help us understand how the spectrum relates to (i) the transient growth and (ii) the optimal initial conditions that give rise to the maximum possible transient growth.
2 Model description
To study the TI, we consider the non-dimensional, incompressible, visco-resistive MHD equations
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}=-\unicode[STIX]{x1D735}p+(\unicode[STIX]{x1D735}\times \boldsymbol{B})\times \boldsymbol{B}+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{B}}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D735}\times (\boldsymbol{u}\times \boldsymbol{B})+\frac{1}{S}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{B}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{B}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{B}$
is the magnetic field,
$\boldsymbol{u}$
is the velocity,
$p$
is the plasma pressure,
$Re$
is the Reynolds number and
$S$
is the Lundquist number.
For our background equilibrium,
$$\begin{eqnarray}\displaystyle p_{0}=p_{0}(x),\quad \boldsymbol{B}_{0}=B_{0}(x)\boldsymbol{e}_{z},\quad \boldsymbol{u}_{0}=U_{0}(x)\boldsymbol{e}_{z}, & & \displaystyle\end{eqnarray}$$
where the subscript 0 corresponds to the equilibrium and
$$\begin{eqnarray}\displaystyle & \displaystyle \mathbf{0}=-\unicode[STIX]{x1D735}p_{0}+(\unicode[STIX]{x1D735}\times \boldsymbol{B}_{0})\times \boldsymbol{B}_{0}+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{0}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \mathbf{0}=\unicode[STIX]{x1D735}\times (\boldsymbol{u}_{0}\times \boldsymbol{B}_{0}). & \displaystyle\end{eqnarray}$$
Magnetic diffusion is not included in the background equilibrium (2.6) as we are only interested in time scales much shorter than the global magnetic diffusion time scale. Clearly, for the assumed forms of the equilibrium magnetic and velocity fields (2.4)
$_{2}$
and (2.4)
$_{3}$
, equation (2.6) is satisfied. Once
$\boldsymbol{u}_{0}$
and
$\boldsymbol{B}_{0}$
are chosen, the background pressure
$p_{0}$
is determined from (2.5).
We will now linearize equations (2.1) and (2.2) about a background equilibrium by setting
$(\boldsymbol{u},\boldsymbol{B},p)=(\boldsymbol{u}_{0},\boldsymbol{B}_{0},p_{0})+(\boldsymbol{u}_{1},\boldsymbol{B}_{1},p_{1})$
and focus on the two-dimensional version of the equations. Assuming perturbations of the form
$$\begin{eqnarray}\displaystyle \boldsymbol{u}_{1}=[u(x,t),0,u_{z}(x,t)]^{\text{T}}\exp (\text{i}kz),\quad \boldsymbol{B}_{1}=[b(x,t),0,b_{z}(x,t)]^{\text{T}}\exp (\text{i}kz), & & \displaystyle\end{eqnarray}$$
the linearized form of (2.1) and (2.2) can be written as
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(D^{2}-k^{2})u=L_{B_{0}}b-L_{U_{0}}u+\frac{1}{Re}(D^{2}-k^{2})^{2}u, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}b}{\unicode[STIX]{x2202}t}=\text{i}k(B_{0}u+U_{0}b)+\frac{1}{S}(D^{2}-k^{2})b, & \displaystyle\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle L_{U_{0}}=\text{i}k[U_{0}(D^{2}-k^{2})-U_{0}^{\prime \prime }],\quad L_{B_{0}}=\text{i}k[B_{0}(D^{2}-k^{2})-B_{0}^{\prime \prime }],\quad D=\unicode[STIX]{x2202}/\unicode[STIX]{x2202}x, & & \displaystyle\end{eqnarray}$$
and the prime refers to differentiation with respect to
$x$
in the background equilibrium fields. Equations (2.8) and (2.9) have essentially the same form as those obtained with reduced MHD in the presence of a large guide field (e.g. Loureiro et al.
Reference Loureiro, Schekochihin and Cowley2007).
To complete the set-up of the model, we require boundary conditions for (2.8) and (2.9). In this paper we will consider no-slip and perfectly conducting boundary conditions,
$$\begin{eqnarray}\displaystyle u=Du=b=0\quad \text{at}~x=\pm d, & & \displaystyle\end{eqnarray}$$
where
$d$
is a non-dimensional distance. Since the tearing instability grows in a thin boundary layer at
$x=0$
, the choice of boundary conditions should not have a large effect on the initial development of the instability if
$d$
is sufficiently large.
To facilitate the discussion of our analysis later, we rewrite (2.8) and (2.9) in the form
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}M\boldsymbol{v}=L\boldsymbol{v}, & & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{v}=(u,b)^{\text{T}}$
,
$$\begin{eqnarray}\displaystyle M=\left(\begin{array}{@{}cc@{}}D^{2}-k^{2} & 0\\ 0 & I\end{array}\right),\quad L=\left(\begin{array}{@{}cc@{}}\displaystyle \frac{1}{Re}(D^{2}-k^{2})^{2}-L_{U_{0}} & L_{B_{0}}\\ \text{i}kB_{0} & \text{i}kU_{0}+\displaystyle \frac{1}{S}(D^{2}-k^{2})\end{array}\right), & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
and
$I$
represents the identity operator.
3 Numerical implementation
In this section we describe briefly the numerical techniques used to discretize and study the behaviour of (2.12). Henceforth, we will describe matrices rather than operators and eigenvectors rather than eigenmodes as the equations will be discretized. To signify this change, the notation for a matrix will have the same letter as the operator it represents but it will now be in bold, e.g.
$A$
is an operator and
$\unicode[STIX]{x1D63C}$
is the finite matrix discretization of that operator. In order not to introduce too much extra notation, eigenvectors will have the same notation as the eigenmodes.
3.1 Solving the eigenvalue problem
Assuming a time dependence of
$\exp (\unicode[STIX]{x1D70E}t)$
, the initial value problem (2.12) becomes the generalized eigenvalue problem
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}\unicode[STIX]{x1D648}\boldsymbol{v}=\unicode[STIX]{x1D647}\boldsymbol{v}. & & \displaystyle\end{eqnarray}$$
The matrix
$\unicode[STIX]{x1D648}^{-1}\unicode[STIX]{x1D647}$
involves (discrete) derivatives (up to fourth order) in
$x\in [-d,d]$
. It is trivial to move from this domain to
$[-1,1]$
and so we discretize the equations using a Chebyshev pseudospectral method (Trefethen Reference Trefethen1999, Reference Trefethen2000). If
$N$
is a positive integer, the
$N+1$
Chebyshev points are given by
$$\begin{eqnarray}\displaystyle x_{i}=\cos \left(\frac{i\unicode[STIX]{x03C0}}{N}\right),\quad i=0,\ldots ,N. & & \displaystyle\end{eqnarray}$$
Note that we could define the
$x_{i}$
with a minus sign in front so as to move from
$-1$
to 1. However, due to the symmetry of our problem, there is no need to make this step. On the domain
$[-1,1]$
the first order spectral differentiation matrix
$\unicode[STIX]{x1D63F}_{N}$
is given by
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle (\unicode[STIX]{x1D63F}_{N})_{00}=\frac{2N^{2}+1}{6},\quad (\unicode[STIX]{x1D63F}_{N})_{NN}=-\frac{2N^{2}+1}{6},\\ \displaystyle (\unicode[STIX]{x1D63F}_{N})_{jj}=-\frac{x_{j}}{2(1-x_{j}^{2})}\quad \text{for}~1\leqslant j\leqslant N-1,\\ \displaystyle (\unicode[STIX]{x1D63F}_{N})_{ij}=\frac{c_{i}}{c_{j}}\frac{(-1)^{i+j}}{x_{i}-x_{j}}\quad \text{for}~i\neq j.\end{array}\right\} & & \displaystyle\end{eqnarray}$$
The above coefficients are defined as
$$\begin{eqnarray}\displaystyle c_{i}=\left\{\begin{array}{@{}ll@{}}2 & \text{for}~i=0~\text{or}~N,\\ 1 & \text{for}~1\leqslant i\leqslant N-1.\end{array}\right. & & \displaystyle\end{eqnarray}$$
The second order differentiation matrix is given simply by
$\unicode[STIX]{x1D63F}_{N}^{2}$
. Due to the boundary conditions (2.11), we strip the first and last rows of
$\unicode[STIX]{x1D63F}_{N}$
and
$\unicode[STIX]{x1D63F}_{N}^{2}$
.
The above matrices are constructed via polynomial interpolation (Trefethen Reference Trefethen2000). To define a fourth order differentiation matrix, we require a polynomial interpolant that satisfies two more boundary conditions than that for the second order differentiation matrix. Again following Trefethen (Reference Trefethen2000), the resulting fourth order differentiation matrix is
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D64E}_{N}=[\text{diag}(1-x_{j}^{2})\unicode[STIX]{x1D63F}_{N}^{4}-8\text{diag}(x_{j})\unicode[STIX]{x1D63F}_{N}^{3}-12\unicode[STIX]{x1D63F}_{N}^{2}]\text{diag}\left(\frac{1}{1-x_{j}^{2}}\right), & & \displaystyle\end{eqnarray}$$
where
$j=1,\ldots ,N-1$
after stripping the first and last rows.
With all the differentiation matrices defined, the matrices of the eigenvalue problem are constructed as they are displayed in (2.13). We then solve the eigenvalue problem (3.1) using the QZ algorithm (Golub & Van Loan Reference Golub and Van Loan1996).
3.2 Quadrature
Once the eigenvalue spectrum is obtained, this describes the asymptotic phase of the linear onset of the TI. To aid our understanding of how the transient growth relates to the spectrum, we require a generalization of the eigenvalue spectrum known as the pseudospectrum (Borba et al. Reference Borba, Riedel, Kerner, Huysmans, Ottaviani and Schmid1994; Trefethen & Embree Reference Trefethen and Embree2005). The calculation of pseudospectra (described in more detail later) and other useful quantities will require the evaluation of norms. To evaluate norms accurately, we must take into account appropriate quadrature weights for the spectral discretizations on an irregular grid. Following the description given in Trefethen (Reference Trefethen1999), we consider a weight matrix of the form
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D652}=\text{diag}(w_{1},\ldots ,w_{N},w_{1},\ldots ,w_{N}), & & \displaystyle\end{eqnarray}$$
with Gauss–Lobatto weights defined by
$$\begin{eqnarray}\displaystyle w_{j}^{2}=\frac{\unicode[STIX]{x03C0}\sqrt{d^{2}-x_{j}^{2}}}{2(N+1)}. & & \displaystyle\end{eqnarray}$$
If
$\Vert \cdot \Vert$
denotes the weighted vector norm that approximates the continuous
$L_{2}$
-norm, then
$$\begin{eqnarray}\displaystyle \Vert \boldsymbol{u}\Vert =\Vert \unicode[STIX]{x1D652}\boldsymbol{u}\Vert _{2}, & & \displaystyle\end{eqnarray}$$
for some vector
$\boldsymbol{u}$
of length
$2N$
. The corresponding matrix norm is
$$\begin{eqnarray}\displaystyle \Vert \unicode[STIX]{x1D63C}\Vert =\Vert \unicode[STIX]{x1D652}\unicode[STIX]{x1D63C}\unicode[STIX]{x1D652}^{-1}\Vert _{2}, & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D63C}$
is a
$2N\times 2N$
matrix. More detailed descriptions of these and related results can be found in Reddy, Schmid & Henningson (Reference Reddy, Schmid and Henningson1993), Trefethen (Reference Trefethen1999), Trefethen & Embree (Reference Trefethen and Embree2005). Throughout the rest of this article, the vector and matrix norms that we will use are those defined in (3.8) and (3.9).
3.3 Matrix projection
One practical issue related to the numerical solution of the eigenvalue problem is that ‘spurious eigenvalues’ (e.g. Bourne Reference Bourne2003) are generated by the numerical scheme which do not have any physical interpretation related to the TI. There are various methods of removing these values, such as solving the equivalent adjoint eigenvalue problem and removing any eigenvalues that do not appear in both the original and adjoint calculations (e.g. Stewart et al.
Reference Stewart, Waters, Billingham and Jensen2009). In this article, we bypass the problem of spurious eigenvalues by projecting
$\unicode[STIX]{x1D648}^{-1}\unicode[STIX]{x1D647}$
onto a lower-dimensional subspace. That is, we focus on a (physically interesting) part of the complex plane and only consider the eigenvalues in this region, cutting out the spurious eigenvalues. This approach also aids to accelerate calculations. There are many ways to achieve projection. In this article we make use of the QR algorithm (Golub & Van Loan Reference Golub and Van Loan1996). Let
$\unicode[STIX]{x1D63D}=\unicode[STIX]{x1D652}\unicode[STIX]{x1D648}^{-1}\unicode[STIX]{x1D647}\unicode[STIX]{x1D652}^{-1}$
and
$\unicode[STIX]{x1D651}$
be an
$2N\times n$
matrix whose columns are selected linearly independent eigenvectors of
$\unicode[STIX]{x1D63D}$
satisfying
$\unicode[STIX]{x1D63D}\unicode[STIX]{x1D651}=\unicode[STIX]{x1D651}\unicode[STIX]{x1D63F}$
for some
$n\times n$
diagonal matrix
$\unicode[STIX]{x1D63F}$
of corresponding eigenvalues. If
$\unicode[STIX]{x1D651}=\unicode[STIX]{x1D64C}\unicode[STIX]{x1D64D}$
is a QR decomposition of
$\unicode[STIX]{x1D651}$
, we can find an upper-triangular
$n\times n$
matrix
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D64F}=\unicode[STIX]{x1D64D}\unicode[STIX]{x1D63F}\unicode[STIX]{x1D64D}^{-1}, & & \displaystyle\end{eqnarray}$$
which is the matrix representation of the projection of
$\unicode[STIX]{x1D63D}$
onto a space of the selected eigenvectors.
Later, we will investigate what parts of the eigenvalue spectrum are ‘physically interesting’, i.e. what parts contribute the largest transient growth. Knowing these locations will allow us to select a suitable projection using the above method.
4 Transient growth analysis
4.1 Parameter selection
We will consider a domain size given by
$d=10$
, which allows us to make comparisons with previous work (MacTaggart & Stewart Reference MacTaggart and Stewart2017) and is a value large enough to not allow the boundary conditions to interfere strongly with the onset of the TI. Here we will only focus on the classical TI and set
$U_{0}=0$
and
$Re=10^{6}$
. In future work we will consider the effects of different background flows and current sheet thicknesses. The background magnetic field equilibrium is given by
$$\begin{eqnarray}\displaystyle B_{0}(x)=\tanh (x), & & \displaystyle\end{eqnarray}$$
corresponding to the Harris sheet (Harris Reference Harris1962).
4.2 Spectra
The eigenvalue spectrum for the onset of the TI consists of a branched structure below
$\text{Re}(\unicode[STIX]{x1D70E})=0$
and a unique eigenvalue on the positive side of this line corresponding to the tearing mode. An example of the spectrum is given in figure 1. This spectrum has been produced with
$S=1000$
,
$k=0.5$
and
$N=700$
. The spectrum consists of three main branches (labelled
$b_{1}$
,
$b_{2}$
and
$c$
) and two (sub)branches connecting
$b_{1}$
and
$b_{2}$
to
$c$
, labelled
$s_{1}$
and
$s_{2}$
. This branching structure is qualitatively similar to that found by Riedel (Reference Riedel1986) by means of WKB analysis and also matches other numerical studies (Goedbloed, Keppens & Poedts Reference Goedbloed, Keppens and Poedts2010; MacTaggart & Stewart Reference MacTaggart and Stewart2017). The spectrum is also symmetric about
$\text{Im}(\unicode[STIX]{x1D70E})=0$
. It can easily be demonstrated that (3.1) possesses the symmetry
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}\rightarrow \unicode[STIX]{x1D70E}^{\ast },\quad u\rightarrow -u^{\ast },\quad b\rightarrow b^{\ast }, & & \displaystyle\end{eqnarray}$$
for the background equilibria we have selected. The asterisk denotes the complex-conjugate.
Eigenvalue spectrum with
$\text{Re}(\unicode[STIX]{x1D70E})\geqslant -0.5$
for
$S=10^{3}$
and
$k=0.5$
. Eigenvalues are shown as solid dots apart from the unique eigenvalue corresponding to the tearing mode, which is shown as a hollow circle.
The general branching structure outlined above is maintained as
$S$
is increased, however, the intersections of the
$b$
and
$s$
branches become more compressed. To illustrate this, figure 2 displays the part of the spectrum, for the case
$S=10^{6}$
,
$k=0.5$
and
$N=2000$
, where the branches
$b_{2}$
and
$s_{2}$
meet. By the symmetry of (4.2), there is an equivalent structure at the intersection of
$b_{1}$
and
$s_{1}$
.
Eigenvalue spectrum for
$S=10^{6}$
and
$k=0.5$
at the intersection of branches
$b_{2}$
and
$s_{2}$
.
The
$b_{2}$
branch has been pushed closer to the
$\text{Re}(\unicode[STIX]{x1D70E})=0$
line and the eigenvalues have become densely packed into a small triangular region. This mirrors the behaviour of the spectrum of the Orr–Sommerfield operator at the intersection of its eigenvalue branches for large
$Re$
(e.g. Schmid & Henningson Reference Schmid and Henningson2001). Eigenvalues in such intersection regions are highly sensitive and full numerical convergence is difficult to achieve (Kerner Reference Kerner1998). However, it is this sensitivity that makes these regions important for transient growth, as will be demonstrated later. Despite the numerical difficulties associated with the calculation of these spectra, our transient growth calculations, which depend on the spectra, are converged and will be shown to follow a clear scaling law (see Hanifi, Schmid & Henningson (Reference Hanifi, Schmid and Henningson1996) for a similar situation).
4.3 Pseudospectra
The individual eigenvalues and eigenvectors only describe the behaviour of a linear instability on a large time scale. They give no information about transient growth that can occur much sooner. One mathematical structure which can provide information about transient growth is a generlization of the spectrum known as the pseudospectrum (Trefethen & Embree Reference Trefethen and Embree2005).
Definition 1. Let
$\unicode[STIX]{x1D63C}=\unicode[STIX]{x1D648}^{-1}\unicode[STIX]{x1D647}$
and
$\unicode[STIX]{x1D716}>0$
be arbitrary. The
$\unicode[STIX]{x1D716}$
-pseudospectrum
$\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D716}}(\unicode[STIX]{x1D63C})$
of
$\unicode[STIX]{x1D63C}$
is the set of
$z\in \mathbb{C}$
such that
$$\begin{eqnarray}\displaystyle \Vert (z\unicode[STIX]{x1D644}-\unicode[STIX]{x1D63C})^{-1}\Vert >\unicode[STIX]{x1D716}^{-1}, & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D644}$
is the identity matrix.
Note that in this definition we could replace
$\unicode[STIX]{x1D63C}$
by
$\unicode[STIX]{x1D63D}$
, as defined in § 3.2. All numerical calculations involving norms will use
$\unicode[STIX]{x1D63D}$
in this article due to the spectral discretization. Note also that pseudospectra are not related to the pseudospectral discretization of the differential equations described in the previous section, both topics just share the same name.
To quote the monograph on pseudospectra, Trefethen & Embree (Reference Trefethen and Embree2005) (to which the reader is directed for a comprehensive account of the subject), ‘the
$\unicode[STIX]{x1D716}$
-pseudospectrum is the open subset of the complex plane bounded by the
$\unicode[STIX]{x1D716}^{-1}$
level curve of the norm of the resolvent’. For non-normal matrices, these level curves can extend
$O(1)$
distances from the position of the eigenvalues. For normal matrices, the curves form
$O(\unicode[STIX]{x1D716})$
balls around the eigenvalues. The
$\unicode[STIX]{x1D716}$
-pseudospectra for the projection of
$\unicode[STIX]{x1D63C}$
, where the eigenvalues satisfying
$-0.5<\text{Re}(\unicode[STIX]{x1D70E})<0$
are kept, are shown in figure 3 for a range of
$\unicode[STIX]{x1D716}$
. This figure has been produced using EigTool (Wright Reference Wright2002).
The eigenvalue spectrum and pseudospectra for the projection covering
$-0.5<\text{Re}(\unicode[STIX]{x1D70E})<0$
. Eigenvalues are shown as solid dots. The boundaries of
$\unicode[STIX]{x1D716}$
-pseudospectra are for
$\unicode[STIX]{x1D716}=10^{-8},\ldots ,10^{-2}$
.
The extension of the level curves far from the position of the eigenvalues, particularly at the intersection points of the
$b$
and
$s$
branches, is a clear sign of the non-normality of
$\unicode[STIX]{x1D63C}$
. Further, a level curve is displayed crossing
$\text{Re}(\unicode[STIX]{x1D70E})=0$
into the positive half-plane. This fact gives important information on the behaviour of transient growth. To see why, first note that the formal solution of the discretized version of (2.12) can be written as
$$\begin{eqnarray}\displaystyle \boldsymbol{v}(t)=\exp (t\unicode[STIX]{x1D63C})\boldsymbol{v}(0). & & \displaystyle\end{eqnarray}$$
In normal mode analysis, the growth rate of the linear onset of an instability is given by the rightmost eigenvalue,
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D63C})=\sup _{z\in \unicode[STIX]{x1D70E}(\unicode[STIX]{x1D63C})}\text{Re}(z). & & \displaystyle\end{eqnarray}$$
The quantity
$\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D63C})$
is also known as the spectral abscissa of
$\unicode[STIX]{x1D63C}$
. An analogous definition for pseudospectra is
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}_{\unicode[STIX]{x1D716}}(\unicode[STIX]{x1D63C})=\sup _{z\in \unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D716}}(\unicode[STIX]{x1D63C})}\text{Re}(z). & & \displaystyle\end{eqnarray}$$
The envelope of transient growth is given by
$\Vert \text{exp}(t\unicode[STIX]{x1D63C})\Vert$
for
$t\geqslant 0$
. This envelope gives the maximum possible transient growth that can be achieved at a time
$t$
optimized over all normalized initial conditions (Schmid & Henningson Reference Schmid and Henningson2001). This fact can found easily from (4.4), i.e.
$$\begin{eqnarray}\displaystyle \sup _{\boldsymbol{v}(0)}\frac{\Vert \boldsymbol{v}(t)\Vert _{2}^{2}}{\Vert \boldsymbol{v}(0)\Vert _{2}^{2}}=\sup _{\boldsymbol{v}(0)}\frac{\Vert \text{exp}(t\unicode[STIX]{x1D63C})\boldsymbol{v}(0)\Vert _{2}^{2}}{\Vert \boldsymbol{v}(0)\Vert _{2}^{2}}=\Vert \text{exp}(t\unicode[STIX]{x1D63C})\Vert _{2}^{2}. & & \displaystyle\end{eqnarray}$$
A simple and practical lower bound on the envelope height is given by
$$\begin{eqnarray}\displaystyle \sup _{t\geqslant 0}\Vert \text{exp}(t\unicode[STIX]{x1D63C})\Vert \geqslant \frac{\unicode[STIX]{x1D6FC}_{\unicode[STIX]{x1D716}}(\unicode[STIX]{x1D63C})}{\unicode[STIX]{x1D716}}\quad \forall \unicode[STIX]{x1D716}>0. & & \displaystyle\end{eqnarray}$$
This result is related to the Kreiss matrix theorem (Trefethen & Embree Reference Trefethen and Embree2005).
For our choice of
$\unicode[STIX]{x1D63C}$
,
$\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D63C})<0$
describes only an asymptotically decaying solution. Since the pseudospectra in figure 3 pass beyond
$\text{Re}(\unicode[STIX]{x1D70E})=0$
, the bound in (4.8) states that transient growth can be expected. Looking at the level curve for
$\unicode[STIX]{x1D716}=10^{-2}$
in figure 3, the maximum value of
$\text{Re}(z)\approx 0.023$
. Hence, it follows from (4.8) that
$\sup _{t\geqslant 0}\Vert \text{exp}(t\unicode[STIX]{x1D63C})\Vert \gtrsim 2.23$
. Later we will demonstrate that the maximum of transient growth follows a precise scaling as a function of
$S$
.
4.4 Transient growth dependence on the spectrum
In order to understand how the transient growth depends on different parts of the spectrum we will consider four different projections of
$\unicode[STIX]{x1D63C}$
(using the parameters
$S=10^{3}$
,
$k=0.5$
and
$N=700$
) that focus on different parts of
$\mathbb{C}$
, i.e. on different groups of eigenvalues. Figure 4 displays the eigenvalue spectra of the four cases that we consider. Case (a) contains only the
$c$
branch with eigenvalues satisfying
$-0.5<\text{Re}(\unicode[STIX]{x1D70E})<0$
. Case (b) includes parts of all three main branches (
$b_{1}$
,
$b_{2}$
and
$c$
) and branches
$s_{1}$
and
$s_{2}$
. These eigenvalues are in the range
$-0.25<\text{Re}(\unicode[STIX]{x1D70E})<0$
. Case (c) omits the branch connections of
$b_{1}$
with
$s_{1}$
and
$b_{2}$
with
$s_{2}$
. These eigenvalues are in the range
$-0.1<\text{Re}(\unicode[STIX]{x1D70E})<0$
. Case (d) considers a small selection of eigenvalues in the range
$-0.01<\text{Re}(\unicode[STIX]{x1D70E})<0$
, which lie on the three main branches.
The eigenvalue spectra of projections of
$\unicode[STIX]{x1D63C}$
focussing on different parts of
$\mathbb{C}$
. Details are given in the main text.
The transient growth envelopes for the four projections whose spectra are displayed in figure 4.
Looking at the overall behaviour of the transient growth displayed in figure 5, cases (a), (c) and (d) are similar in that all reach maximum growth at
$t\approx 100{-}200$
and have similar maxima in the range 3.5–5. The main visual difference between these cases is that the transient growth in case (a) is smooth, whereas those in cases (c) and (d) possess many oscillations. This difference is down to only the central branch of eigenvalues being included in (a), whereas portions of the three main branches are included in (c) and (d). Case (c) also has a faster initial growth rate compared to cases (a) and (d) and attains the highest transient growth out of the three cases. The spectrum of case (c) contains the upper parts of the
$s$
branches but not the connection points with the
$b$
branches.
Case (b) is strikingly different to the rest. Its transient growth exhibits a sharp rise to a maximum that is approximately double that of the other cases. The eigenvalues considered for case (b) include all the branches and the intersection points. The transient growth curve in figure 5(b) and the behaviour of the
$\unicode[STIX]{x1D716}$
-pseudospectra in figure 3 indicate that the intersection points of the
$b$
and
$s$
branches in the spectrum are important for strong transient growth.
4.5 The relationship of physical quantities to the spectrum and transient growth
The previous subsection demonstrated that including the branching structure of the eigenvalue spectrum in the projection is important for fast and strong transient growth. We will now focus on how the branches depend on physical quantities, namely forces and energy balance.
Let us first consider forces in the momentum balance equation (2.1). In the present set-up of the TI with no equilibrium flow, the viscosity plays little role. This can be seen by comparing figure 1 with figure 1(a) of MacTaggart & Stewart (Reference MacTaggart and Stewart2017), which plots the spectrum for the same parameters but for the inviscid MHD equations. Note that in MacTaggart & Stewart (Reference MacTaggart and Stewart2017), time dependence is based on
$\exp (-\text{i}\unicode[STIX]{x1D70E}t)$
rather than
$\exp (\unicode[STIX]{x1D70E}t)$
, so the plot of the spectrum is rotated.
In order to assess the effect of the Lorentz force, we can compare the spectrum of resistive MHD to that of kinematic MHD. Kinematic MHD is concerned with the solution of the induction equation and normally does not consider the momentum equation (the velocity field is prescribed). Here, we can solve a form of kinematic MHD where the momentum equation is solved but the Lorentz force is removed. In the linearized equations, this is equivalent to setting
$L_{B_{0}}=0$
. Figure 6 displays the spectrum of kinematic MHD for the same parameters used for figure 1.
The spectrum of kinematic MHD with the same parameters used for figure 1.
Notice that the
$c$
branch is the only one remaining and that this system is asymptotically stable to linear perturbations. From the results of the previous subsection, this means that the possible transient growth is much more limited. This result makes sense physically since the energy source of the TI is the equilibrium magnetic field and the bending of field lines via the Lorentz force can lead to larger perturbations.
Another way to consider what physical quantities play a role in transient growth is to determine the energy balance. Borba et al. (Reference Borba, Riedel, Kerner, Huysmans, Ottaviani and Schmid1994) did this for inviscid MHD and we will extend the analysis for the present visco-resistive case. Using the notation of § 2, the application of vector calculus leads to
$$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}\frac{1}{2}\int (|\boldsymbol{u}_{1}|^{2}+|\boldsymbol{B}_{1}|^{2})\,\text{d}V & = & \displaystyle \int [\boldsymbol{j}_{1}\boldsymbol{\cdot }(\boldsymbol{u}_{0}\times \boldsymbol{B}_{1})-\boldsymbol{j}_{0}\boldsymbol{\cdot }(\boldsymbol{u}_{1}\times \boldsymbol{B}_{1})]\,\text{d}V-\frac{1}{S}\int |\unicode[STIX]{x1D735}\times \boldsymbol{B}_{1}|^{2}\,\text{d}V\nonumber\\ \displaystyle & & \displaystyle -\,\int (\boldsymbol{u}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}_{0}\boldsymbol{\cdot }\boldsymbol{u}_{1}\,\text{d}V-\frac{1}{Re}\int |\unicode[STIX]{x1D735}\times \boldsymbol{u}_{1}|^{2}\,\text{d}V,\end{eqnarray}$$
where
$\boldsymbol{j}_{0}=\unicode[STIX]{x1D735}\times \boldsymbol{B}_{0}$
and
$\boldsymbol{j}_{1}=\unicode[STIX]{x1D735}\times \boldsymbol{B}_{1}$
. The integrals associated with the diffusion terms are either zero or negative (including the negative sign), so these terms can only describe decay and not transient growth. The first integral on the right-hand side describes how energy can be taken from the equilibrium magnetic field and background flow due to the effect of the Lorentz force, leading to transient growth. The third integral on the right-hand side describes how transient growth is possible, without the aid of the Lorentz force, by the perturbed flow extracting energy from the equilibrium flow. In the examples considered so far, the terms involving
$\boldsymbol{u}_{0}$
are zero. This fact, combined with the large Reynolds number we have selected (
$Re=10^{6}$
) explains the similarity of the spectra of the viscous and inviscid cases.
4.6 Scaling of maximum transient growth
We now return to the projection of
$\unicode[STIX]{x1D63C}$
with eigenvalues in the range
$-0.5<\text{Re}(\unicode[STIX]{x1D70E})<0$
. Since we are considering the case
$k=0.5$
, we take
$\text{Re}(\unicode[STIX]{x1D70E})=-0.5$
as the ‘lower boundary’ of our projection. Beyond this point, only the
$c$
branch continues and this part of the spectrum does not contribute significantly to transient growth, as would be expected from the analysis above. In order to determine how the transient growth scales with the Lundquist number
$S$
, we explicitly calculate
$\Vert \text{exp}(t\unicode[STIX]{x1D63C})\Vert$
for
$t\in [0,300]$
and
$S=10^{n}$
,
$n=3,4,5,6$
. The transient growth envelopes are shown in figure 7 and are calculated with
$N=2000$
.
The transient growth envelopes for
$k=0.5$
and
$S=10^{3},\ldots ,10^{6}$
.
After an initial perturbation, all the curves follow the same gradient before they each turn toward their maxima. Each curve contains a ‘hump’ where it reaches its maximum and this is due to the inclusion of the branch intersection points of
$b$
and
$s$
in the spectrum, as in case (b) from § 4.4. If we consider the maximum transient growth of each curve and the times when the maxima occur, we can find a simple scaling law. These quantities are displayed in figure 8.
The maximum transient growth,
$\max \Vert \text{exp}(t\unicode[STIX]{x1D63C})\Vert$
, and the time when it occurs,
$t_{\max }$
, as a function of
$S^{1/4}$
. The dashed lines are lines of best fit.
Figure 8 demonstrates that the maximum transient growth,
$\max \Vert \text{exp}(t\unicode[STIX]{x1D63C})\Vert$
, and the time at which it occurs,
$t_{max}$
, both depend linearly on
$S^{1/4}$
. This simple scaling relation is robust for ‘tearing-unstable’ wavenumbers, i.e.
$0<k<1$
. Figure 9 demonstrates the same scaling profiles for
$k=0.2$
and
$k=0.8$
.
The maximum transient growth,
$\max \Vert \text{exp}(t\unicode[STIX]{x1D63C})\Vert$
, and the time when it occurs,
$t_{\max }$
, as a function of
$S^{1/4}$
. The dashed lines are lines of best fit. (a) and (b) refer to
$k=0.2$
, (c) and (d) refer to
$k=0.8$
.
The projection we have used for the cases
$k=0.2$
and
$k=0.8$
is the same as that for the
$k=0.5$
case. This means that for the
$k=0.2$
case, more eigenvalues on the
$c$
branch are used beyond the point where the
$c$
and
$b$
branches meet. For the
$k=0.8$
case, however, the spectrum is truncated before the
$b$
branches meet the
$c$
branch. Despite these different truncations of the spectra, since both cases include the vital (as emphasized by the previous analysis) branch points, the optimal transient growth is insensitive to the different truncation locations and follows the same scaling law as the
$k=0.5$
case. Our numerical results suggest the scaling that a maximum transient growth of
$O(S^{1/4})$
, optimized using the
$L_{2}$
-norm, is possible in a time of
$O(S^{1/4})$
for
$S\gg 1$
.
4.7 Optimal initial perturbations
As well as determining the maximum possible transient growth at any time
$t$
, we can also determine the initial condition which produces that growth at time
$t$
. This can be achieved via a singular value decomposition (SVD),
$$\begin{eqnarray}\displaystyle \exp (t\unicode[STIX]{x1D63C})=\unicode[STIX]{x1D650}\unicode[STIX]{x1D72E}\unicode[STIX]{x1D651}^{\ast }, & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D650}$
and
$\unicode[STIX]{x1D651}$
are unitary matrices and
$\unicode[STIX]{x1D72E}$
is a matrix containing the singular values ordered by size. The asterisk denotes the complex-conjugate transpose. The first column of
$\unicode[STIX]{x1D651}$
corresponds to the optimal initial condition (Schmid & Henningson Reference Schmid and Henningson2001; Trefethen & Embree Reference Trefethen and Embree2005). To illustrate the form of such initial conditions, figure 10 displays the non-zero parts of
$u$
and
$b$
at
$t=0$
which produce the maximum transient growth at time
$t=50$
, for
$S=10^{4}$
and
$k=0.5$
. The forms shown in figure 10 have a striking resemblance to ‘wave packets’ located at the current sheet. Such wave packet solutions can be understood with the help of pseudospectral theory. A definition of pseudospectra equivalent to Definition 1 is
The optimal initial conditions that produce the maximum transient growth at
$t=50$
for
$S=10^{4}$
and
$k=0.5$
.
Definition 2.
$\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D716}}(\unicode[STIX]{x1D63C})$
is the set of
$z\in \mathbb{C}$
such that
$$\begin{eqnarray}\displaystyle \Vert (z\unicode[STIX]{x1D644}-\unicode[STIX]{x1D63C})\boldsymbol{p}\Vert <\unicode[STIX]{x1D716}, & & \displaystyle\end{eqnarray}$$
for some vector
$\boldsymbol{p}$
with
$\Vert \boldsymbol{p}\Vert =1$
.
The vector
$\boldsymbol{p}$
is known as a pseudomode and can be thought of as a generalization of an eigenmode in the same way that a pseudospectrum is a generalization of a spectrum. That is, a pseudomode grows algebraically and an eigenmode grows exponentially. For the onset of the TI, this algebraic growth can be much faster and greater, within a fixed time, than the exponential growth of the tearing mode. Pseudomodes have a close relationship to the WKB approximation of eigenmodes (Obrist & Schmid Reference Obrist and Schmid2010), resulting in their wave packet form. They can also be interpreted as ‘noise’ in the system (Vanneste & Byatt-Smith Reference Vanneste and Byatt-Smith2007). Further details of the theory of pseudomodes can be found in the monograph of Trefethen & Embree (Reference Trefethen and Embree2005).
5 Conclusions
5.1 Summary
In this article we have studied the non-modal onset of the classical tearing instability for visco-resistive MHD. We have paid particular attention to the eigenvalue spectrum of the linearized MHD equations. We demonstrate that the branching structure of the spectrum, which exists in the ‘damped half-plane’ of
$\mathbb{C}$
, is important for transient growth. We reveal this behaviour through the calculation of pseudospectra and by finding the maximum possible transient growth due to subsets of the spectrum. The spectrum branches are also closely linked to the Lorentz force, which is needed for strong transient growth. The importance of the Lorentz force in driving transient growth is also found from considering the energy balance of the system. A simple scaling law is determined for tearing-unstable wavenumbers, revealing that the maximum possible transient growth, measured in the
$L_{2}$
-norm, can grow to
$O(S^{1/4})$
in a time of
$O(S^{1/4})$
. Optimal initial conditions which produce the maximum transient growth are shown to take the shape of wave packets and can be interpreted as noise in the system. Although significant transient growth is possible during the linear onset of the TI for
$S\gg 1$
, it will only occur if the form of the initial perturbation allows it. Both non-modal and modal growth are required to give a complete description of the linear onset of the TI.
5.2 Discussion
5.2.1 Tearing-stable cases
In this article we have focussed on wavenumbers for which the current sheet is unstable to the TI
$(0<k<1)$
. For wavenumbers
$k>1$
, transient growth is also possible and its maximum possible size increases with increasing
$S$
. There does not, however, appear to be a simple scaling law as we derived for tearing-unstable values of
$k$
. Figure 11 displays the optimal initial conditions that produce the maximum size of transient growth at
$t=50$
for
$S=10^{4}$
and
$k=1.01$
(the size of this optimal transient growth is 11.94).
The optimal initial conditions that produce the maximum transient growth at
$t=50$
for
$S=10^{4}$
and
$k=1.01$
.
In this example, the optimal initial conditions take the form of two wave packets on either side of the current sheet. Perturbations such as those shown in figure 11 should be used as initial conditions in nonlinear MHD simulations in order to determine the nonlinear consequences of linear transient growth. It may be the case that optimal initial conditions, through transient growth, can excite the tearing instability for values of
$k$
which are linearly stable in normal mode analysis.
5.2.2 Choice of norm
Unlike asymptotic stability, the size of transient growth is dependent on the norm used to measure it. Optimizing with respect to different norms will give different results. Therefore, in calculations of transient growth, it is important to choose a norm with a clear physical meaning. In this article, we have focussed on the
$L_{2}$
-norm, which can be thought of as the ‘root mean square’ of the variables and measures typical size that the variables can be amplified to. Other useful norms are the infinity norm, which measures the maximum amplitude of the perturbation and the energy norm. For a given
$k$
, the energy norm (disturbance kinetic plus magnetic energies) can be written as (e.g. MacTaggart & Stewart Reference MacTaggart and Stewart2017)
$$\begin{eqnarray}\displaystyle \Vert \boldsymbol{v}\Vert _{E}^{2}=\frac{1}{2k^{2}}\int _{-d}^{d}(|Du|^{2}+k^{2}|u|^{2}+|Db|^{2}+k^{2}|b|^{2})\,\text{d}x=\frac{1}{2k^{2}}\Vert D\boldsymbol{v}\Vert _{2}^{2}+\frac{1}{2}\Vert \boldsymbol{v}\Vert _{2}^{2}, & & \displaystyle\end{eqnarray}$$
where use has been made of (2.3)
$_{1}$
and (2.3)
$_{2}$
. Notice from (5.1) that if a certain value of transient growth is found in the
$L_{2}$
-norm, this is a lower bound of the resulting energy measure. The converse is not generally true, however, as a large measure in the energy norm need not imply a large
$L_{2}$
-norm. With regard to transient growth in the TI, however, there is still the possibility of increasing transient growth with increasing
$S$
, optimized with respect to the energy norm. This result was first looked at in MacTaggart & Stewart (Reference MacTaggart and Stewart2017). Using the technique described in MacTaggart & Stewart (Reference MacTaggart and Stewart2017) we present, in figure 12, numerical estimates of the maximum growth envelopes for different values of
$S$
. As in our previous analysis using the
$L_{2}$
-norm, we take
$k=0.5$
and consider the contribution from asymptotically stable modes with eigenvalues in the range
$-0.5<\text{Im}(\unicode[STIX]{x1D70E})<0$
. The resolution is
$N=1600$
.
The maximum energy growth envelopes for different
$S$
.
Although the shapes and maxima of the transient growth envelopes in figure 12 are different compared to those in figure 7, there is still an increasing possibility of significant transient growth as
$S$
increases. Further work is required, using both the
$L_{2}$
-norm and the energy norm, to investigate transient growth at very high (astrophysical) values of
$S$
.
5.2.3 Transient growth in other TI research?
Simulations of the TI have generally skipped directly to the nonlinear phase of the instability, e.g. the GEM Magnetic Reconnection Challenge (Birn et al.
Reference Birn, Drake, Shay, Rogers, Denton, Hesse, Kuznetsova, Ma, Bhattacharjee and Otto2001a
; Birn & Hesse Reference Birn and Hesse2001b
), or have been interpreted in terms of the tearing mode. In some high Lundquist number simulations, such as Samtaney et al. (Reference Samtaney, Loureiro, Uzdensky, Schekochihin and Cowley2009) with
$S=10^{4}$
, transient growth is not reported. However, significant transient amplification will only occur if the initial condition (perturbation) is of a suitable form. It is therefore possible for transient growth not to be detected in simulations if the initial condition is not one that leads to significant transient growth. We have shown that optimal initial conditions take the form of wave packets, which can be interpreted as noise. Interestingly, recent simulations by Huang, Comisso & Bhattacharjee (Reference Huang, Comisso and Bhattacharjee2017) clearly show that different levels of noise in the initial condition affect when current sheet disruption occurs (see their figure 12). It is not unreasonable to suggest that different noise patterns could lead to different transient growth, affecting when current sheet disruption takes place. A systematic study of how initial conditions, transient growth and nonlinear consequences are related in MHD simulations will be carried out in future work.
Acknowledgements
The author would like to thank Professor A. Valli and Drs J. Pestana and P. Stewart for stimulating conversations on various topics related to this article. The author would also like to thank the anonymous referees for helping to improve this article.
