{'fraction_non_alphanumeric': 0.054940074529157015, 'fraction_numerical': 0.018968006177191393, 'mean_word_length': 4.1477577119156654, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stability analysis of the numerical Method of characteristics applied to a class of energy-preserving systems. Part II: Nonreflecting boundary conditions

We show that imposition of non-periodic, in place of periodic, boundary conditions (BC) can alter stability of modes in the Method of characteristics (MoC) employing certain ordinary-differential equation (ODE) numerical solvers. Thus, using non-periodic BC may render some of the MoC schemes stable for most practical computations, even though they are unstable for periodic BC. This fact contradicts a statement, found in some literature, that an instability detected by the von Neumann analysis for a given numerical scheme implies an instability of that scheme with arbitrary (i.e., non-periodic) BC. We explain the mechanism behind this contradiction. We also show that, and explain why, for the MoC employing some other ODE solvers, stability of the modes may be unaffected by the BC.


Introduction
In Part I [1] of this study we considered the numerical stability of the Method of characteristics (MoC) applied to a class of hyperbolic partial differential equations (PDEs) with periodic boundary conditions (BC). More specifically, the ordinary differential equation (ODE) numerical solvers employed by the MoC in [1] along the characteristics were the simple Euler (SE), modified Euler (ME), and Leapfrog (LF) ones. The class of the PDEs considered in [1] has the linearized form: whereũ is a small perturbation on top of the background solution u (0) of the original nonlinear PDE, Σ = diag(I N , −I N ), with I N being the N -dimensional identity matrix, and P is a constant 2N ×2N matrix. Importantly, the considered class of the PDEs possesses a number of "conservation laws". In particular, the solution u ≡ u + , u − T of the nonlinear PDE satisfies: where | · · · | stands for the length of the corresponding N -dimensional vector. For periodic BC on hence the name "conservation law". (For future reference, let us note that the meaning of the integral here is the energy of the solution.) Therefore, we initially believed that the LF solver, which is well known to (almost) preserve conserved quantities of energy-preserving ODEs over indefinitely long integration times, would also perform the best among the three aforementioned solvers when employed by the MoC to integrate the above energy-preserving PDE.
In [1] we showed that for periodic BC, the MoC with any of those three ODE solvers exhibits numerical instability. For the SE solver, the strongest such an instability occurs in the "ODE" and "anti-ODE" limits, i.e. for the wavenumbers k = 0 and |k| = k max ≡ π/h, respectively. (Here h is the grid spacing in both space and time.) This mild instability of the SE applied to conservative ODEs (e.g., the harmonic oscillator model) is well known: see, e.g., [2]. The ME solver is known to exhibit a similar, although much weaker, instability when applied to conservative ODEs. We have found, however, that when the ME is used within the MoC framework, its most unstable modes occur in the "middle" 1 of the Fourier spectrum, i.e. for |k| ∼ k max /2, and that the instability's growth rate is much greater than that in the ODE (and anti-ODE) limit. Finally, the LF solver used along the characteristics was found to exhibit by far the strongest instability among those three solvers. In Fig. 1 we show the amplification factor, |λ| ≡ ũ n+1 / ũ n of the three numerical schemes: MoC-SE, MoC-ME, and MoC-LF; hereũ n is the numerical error at the nth time level.
Note the order-of-magnitude difference of the instability growth rate of the MoC-LF compared to that rate of the MoC-Euler methods. In this paper we extend the stability analysis of these three schemes to the case where the BC are not periodic. More specifically, we consider so-called nonreflecting BC, whose exact form will be stated in Section 2 and which are often relevant for hyperbolic PDEs whose left-hand side (lhs) is of the same form as that of (1). The motivation for our study came from the following numerical observation. When we simulated the given PDE system with nonreflecting BC by each of the three schemes, we found that the MoC-SE and MoC-LF exhibited the same growth rates of the numerical error as in the case of periodic BC. This is in agreement with the common knowledge that the instability of a problem with periodic BC implies an instability of a problem with any other BC (see below). However, the instability of the MoC-ME (see Fig. 1(b)) was suppressed when we used nonreflecting BC. 2 This surprising observation prompted two questions, which we address in this work.
First, as we have already mentioned, it is commonly stated in textbooks on numerical analysis of PDEs (see, e.g., [3]- [5]) that the numerical instability of a problem with periodic BC is sufficient for the same problem with some non-periodic BC to also be numerically unstable. This statement may be more familiar to the reader in its equivalent form: "The von Neumann stability analysis of a scheme (with constant coefficients) gives a necessary, but not always sufficient, condition for its stability". We emphasize here that this statement is not at all related to the fact that some of the eigenvalues found by the von Neumann analysis may be at the boundary of the stability region, with the non-periodicity "pushing" them into the instability region. Nor is it related here to the presence of eigenvalues with high algebraic but low geometric multiplicity (i.e., "repeated" eigenvalues). Rather, it is related to the fact, assumed in, e.g., [3]- [5], that the modes of the periodic problem (i.e., Fourier harmonics) approximate a proper subset of the modes of the non-periodic problem. We will discuss this in detail in Section 7.3.
Our finding about the MoC-ME clearly violates this common knowledge: the scheme is unstable for periodic BC but becomes stable (see footnote 2 above) for nonreflecting BC. Thus, the first question that we will answer is: (i) Why does this occur?
However, as we have mentioned, this suppression of numerical instability takes place for only one out of the three schemes considered. Thus, the second question is: (ii) Why does it occur for the MoC-ME but not for the MoC-SE and MoC-LF?
The main part of this paper is organized as follows. In Section 2 we give more details about the PDE system under study. In Section 3 we set up the framework for the stability analysis of the MoC schemes with non-periodic BC. Let us stress that this framework is different from the von Neumann stability analysis, which is used for problems with periodic BC. In fact, we have found only one paper, [6], where a similar (although less general) type of analysis for a MoC method had been used. In Section 4, we will apply our analysis to the MoC-SE with nonreflecting BC. Let us stress that this scheme is not at the focus of our study -because nonreflecting BC do not suppress instability for it, -and yet we will spend a substantial amount of effort on this case. The reason is that in this simplest case, we will be able to not only clearly outline the steps of the analysis, but also to quantitatively verify its predictions by direct numerical simulations. In Section 5 we will apply this analysis to the MoC-ME and show that for it, unlike for the MoC-SE, nonreflecting BC suppress instability. A qualitative reason for that is explained at the end of Section 4. This explanation will address the 'second question' listed above, but only partially: it will leave open the question as to why the same mechanism does not suppress an instability for the MoC-LF. We will address the MoC-LF case in Section 6, thereby completely answering the 'second question'.
In Section 7 we will summarize our conclusions and will present a qualitative explanation to the 'first question', i.e., why an instability of a problem with periodic BC does not always imply an instability of the same problem with non-periodic BC. As this explanation pertains to, perhaps, the most unexpected result of our study, we compartmentalized it in a separate subsection, 7.3. The reader who is not interested in the details of the analysis may proceed directly to that subsection.
In subsection 7.4 we complement the quantitative answer to the 'second question', given in Section 6, by a qualitative consideration. Appendices A and B contain technical derivations for Sections 3 and 4. In Appendix C we present a background on a system different from the constant-coefficient system considered in the main part of the paper. That other system is a soliton (i.e., localized and thus non-constant) solution of the Gross-Neveu model in the relativistic field theory, which has been actively studied in the past decade. In subsection 7.2 we illustrate with direct numerical simulations that suppression of numerical instability of the MoC-ME by nonreflecting BC, predicted by our analysis of the constant-coefficient system, also occurs for the Gross-Neveu soliton.

Physical model
While our study will focus on a linear problem of a rather general form (see Eq. (1) or (6a) below), we will begin by stating a specific nonlinear problem which had originally motivated our study and whose linearization leads to (6a). The vector form of the PDE system under consideration is: where S ± ≡ [S ± 1 , S ± 2 , S ± 3 ] T ,Ĵ = diag(1, −1, −2), and superscript 'T' denotes the transposition. This system is a representative of a class of models that arise in studying propagation of light in birefringent optical fibers with Kerr nonlinearity [7]- [10]. It should be noted that the specific numerical entries ofĴ arise from physical considerations and therefore should not be replaced with arbitrary numbers, as that would not correspond to a physical situation. The component form and a non-constant, soliton/kink solution of (3) can be found in, e.g., [1]. Here (as in [1]) we consider the numerical stability of a constant solution of (3): when it is simulated by the MoC. Considering only one representative, (3), of a broader class of models (see [1]) and its simplest solution, (4), will allow us to focus on the analysis of the numerical scheme without being distracted by complexities of the physical model. The relevance of the numerical (in)stability of the constant solution (4) to the numerical (in)stability of "more physically interesting", soliton/kink and related, solutions was discussed in [1] after Eq. (8).
We will linearize Eqs. (3) on the background of solution (4) using where S ± j 0 are the components of the exact solution (4) and s ± j are small perturbations. The linearized system (3) has the form: where s = [s + 1 , s + 3 , s − 1 , s − 3 ] T , the 4 × 4 matrices in (6a) are: and the 2 × 2 matrices in (7a) are: Given the trivial dynamics (6b) of s ± 2 , below we will consider only the dynamics of s ± 1,3 , given by (6a). As we pointed out in [1], systems of the latter form describe linear or linearized dynamics in a wide class of physical models in plasma physics, various areas of optics, and relativistic field theory.
Let us reiterate, however, that the focus of this work is not on a particular physical application but, rather, on understanding the behavior of the MoC schemes with non-periodic BC, as we explained in the Introduction. The reason 3 why we chose to consider the specific form of matrix P is that the forthcoming analysis simply cannot be performed for an unspecified form of P . (Note that P in (6a) cannot be diagonalized without affecting the matrix on the lhs.) While the results of the analysis, of course, depend on the specific form of P, its methodology does not. Moreover, the conclusion about suppression of the instability of the MoC-ME by nonreflecting BC holds for at least one other system of current research interest, as we demonstrate in Section 7.2 with direct numerical simulations. 4 In Eqs. (6) and (7) and in what follows we have adopted the following notations. Boldfaced quantities with an underline, S ± , will continue to denote 3×1 vectors, as in (3). Boldfaced quantities without an underline or a hat will denote 4 × 4 matrices or 4 × 1 vectors, as in (6a); the ambiguity of the same notations for matrices and vectors here will not cause any confusion. Finally, underlined letters in regular (not boldfaced) font will denote 2 × 1 vectors; e.g.: Clearly then, s ≡ (s + ) T , (s − ) T T .
The nonreflecting BC for system (3) are: where x ∈ [0, L] and F ± (t) are given. These BC specify the right (left)-propagating components at the left (right) boundary of the spatial domain. For the small perturbation in (5), BC (9) transform into: 3 Setup of the stability analysis of the MoC with BC (10) We will illustrate details of this setup using the simplest "flavor" of the MoC, the MoC-SE. This will allow us to skip most details in later sections, devoted to the more complex schemes, the MoC-ME and MoC-LF.
Scheme (15) with nonreflecting BC can be written in the form or, equivalently, Note that in (17a), O denotes the 4 × 4 zero matrix.
To analyze stability of the MoC-SE, we need to determine whether the largest-in-magnitude eigenvalue of N exceeds 1. Since N is block-tridiagonal and Toeplitz, we will find its eigenvalues using the method for non-block tridiagonal Toeplitz matrices (see, e.g., [11]); note that it is here where our approach differs from that in [6]. Namely, we first seek the 4 × 1 "components" of × in the form: In Appendix A we explain why ρ should be considered a scalar rather than a 4 × 4 matrix. Substituting (18) into (15a) and using where λ is the eigenvalue that we want to find, we obtain for the eigenvector ξ ≡ (s) m : We will now outline three steps of our analysis, which will be carried out in Sections 4 and 5 for the MoC-SE and MoC-ME, respectively.
Step 1: We will show that the characteristic equation for (20) yields four solutions ρ j (λ), j = 1, . . . , 4 for a given λ. Step 2: For the ρ j found in Step 1, we will find their respective eigenvectors ξ j from (20). These eigenvectors do not affect the spatial structure of the modes but simply reflect the distribution of "energy" among the components of s ≡ [s + 1 , s + 3 , s − 1 , s − 3 ] T at any given point in space.
Step 3: Using the form we will determine constants C j by requiring (21) to satisfy the BC (14). Namely, denoting and substituting (21) into (14), we have: These conditions can be rewritten as a homogeneous linear system: Note that the dependence of Φ on λ comes from such a dependence of ξ j and ρ j . Thus, the eigenvalue λ of the amplification matrix in (17) will be found from the characteristic equation Let us note that the above steps are standard in analyzing any intial-boundary value problem (see, e.g., [12]) with separable spatial and temporal variables. Importantly, they differ from the steps of the stability analysis of the initial-boundary value problems found in textbooks [3]- [5] in the following aspect. The textbooks postulate that all spatial modes of the problem fall into two groups: those localized near either of the boundaries and those resembling Fourier harmonics inside the spatial domain, sufficiently far away from the boundaries. In contrast, our Step 1 finds the spatial modes. Surprisingly, those modes turn out to violate the aforementioned categorization into the two groups. It is this circumstance that will be shown to be the reason behind the "disappearance" of the numerical instability of the MoC-ME, announced in the Introduction.

Stability analysis of the MoC-SE
In the first three subsections of this Section we will implement the three respective steps listed at the end of Section 3. The main result of our analysis will be obtained in subsection 4.3. In the fourth subsection, we will present a verification of this analysis by direct numerical simulations of scheme (11), (12). In the fifth subsection, we will discuss what bearing the results of this Section will have on those in Section 5.

4.1
Step 1: Finding ρ(λ) in (20) To obtain a perturbative (for h ≪ 1) solution of the characteristic polynomial of (20), note that it can be written as where A, B are defined in (7b). Consequently, in the main order, one has Since each of these roots is double, then for 0 < h ≪ 1 one will have four roots of (25). To find them, we substitute 5 into (25). The orders O(h 2 ) and O(h 3 ) of the resulting expression, found with software Mathematica, yield expressions for ρ (1) j and ρ (2) j , respectively: Let us note that the subscript notations '(±)' refer to distinct roots in (28). Thus, they are in no way related to the superscript notations '±', which were used in, e.g., (22) and will be used in similar context below. Those superscript notations denote the "forward"-and "backward"-propagating components of the numerical error, in analogy with that notation in (13). To emphasize this difference between the sub-and superscript ± notations, we will always use parentheses in subscript (±).

Formulae (28) complete
Step 1 of the analysis of the eigenvalue problem for the amplification matrix N in (17). A discussion about one of its assumptions will now be in order. 5 see a discussion in the paragraph that starts with Eqs. (29) It should be noted that seeking the solution of (25) in the form (27) is valid only as long as the two roots (26) are "sufficiently distinct", i.e.: or, equivalently, In other words, the analysis presented in this section will break down for where the implication is based on (26), (27). (Note that this does not preclude the possibility |λ| 2 ≈ 1; e.g., for λ 2 ≈ −1 our analysis will be valid.) The proper approach in the case (30a) is to for which a fourth-degree polynomial in α and β will result. The corresponding calculations are considerably more complex than those in Sections 4.1-4.3 and will not be carried out. 6 Therefore, it is now important to clarify the following issues: (i) What we will miss by not considering case (30); (ii) Why this missed information will not be crucial for our analysis; and (iii) How the curtailed analysis based on (27) will benefit us.
Before we address these issues, we need to establish a correspondence between modes of the problem with nonreflecting BC, characterized by parameter ρ, and Fourier harmonics of the problem with periodic BC. According to (18), ρ is the counterpart of e ikh in the von Neumann analysis (see the line before Eq. (18) in [1]). Therefore, the modes in (30a) correspond to the Fourier harmonics with the lowest and highest wavenumbers, kh = 0 and kh = π: We will now present answers to questions (i)-(iii) above.
(i) We will be unable to rigorously describe the instability of the MoC-SE with nonreflecting BC. Indeed, as Fig. 1(a) shows, the most unstable harmonics of the periodic problem are those with kh = 0 and kh = π, and our numerical simulations, described in Section 4.4, indicate that the corresponding modes (31) are the most unstable in the problem with nonreflecting BC. Thus, it may sound as if we will intentionally neglect the main goal of our study. However, in the next two paragraphs, we will explain why this is not so.
(ii) Let us focus on the mode with ρ ≈ 1; the case of the mode with ρ ≈ −1 can be considered similarly. While a nonzero mode with ρ = 1 will not satisfy the BC (10), one can still consider the 6 A related analysis, however, will be required for the MoC-LF and therefore will be presented in Section 6.
Then, by taking the limit β → 0 in the characteristic polynomial obtained from the substitution of (30b) into (20), one finds: The last two roots correspond exactly to the values (|λ|− 1)/h 2 found by the von Neumann analysis for the periodic BC [1]. Moreover, in our simulations of the MoC-SE, we found that the growth rates of the numerically unstable modes are the same for the periodic and nonreflecting BC. Thus, the simplified analysis above has been able to quantitatively predict the instability growth rate of the MoC-SE with periodic BC.
(iii) Then, what is the point of considering the less unstable (or, rather, as we will show, even stable) modes satisfying (27), (29)? The point is that we will thereby uncover a mechanism which suppresses the numerical instability of those modes, and it will turn out to be the same mechanism which suppresses the most unstable modes of the MoC-ME! We will preview this in Section 4.5.
Here we will only reiterate our earlier statement: Presenting details for the simplest scheme, the MoC-SE, will keep the analysis more transparent and will also allow us to proceed faster when considering the MoC-ME in Section 5.

Step 2: Finding ξ in (20)
We will now find the eigenvectors corresponding to the four roots (28). We will present details of seek the solution of (20) in the form where u, v = O(1) are to be determined. (It will become clear later that to explain the suppression of instability of the MoC-SE, we will not need higher-order terms in (33).) Using (16a), (28) with terms up to O(h), and (33), we obtain in the O(h) order of (20): The top 2 × 1 block of this equation determines u: here A (and B below) is defined in (7b). The bottom 2 × 1 block of (34) yields: Performing similar calculations for the roots ρ 1(−) and ρ 2(±) and using the explicit form of matrix B, one finds: Recall that, as stated after (28), the super-and subscript ± notations are unrelated to each other.
The superscripts '±' are used by analogy with that notation in (13) and thus denote the "forward"and "backward"-propagating components of the numerical error. On the other hand, the subscripts (±) simply refer to distinct roots in (28).

Step 3: Finding eigenvalues of N in (17)
Here we will obtain the main result of Section 4; it is given by in qualitative form by relation (45) and in quantitative form by Eq. (43).
As we explained in Section 3, the eigenvalues of the amplification matrix N are found from (24).
Recall that (M + 1) is the number of points in the spatial grid.
Before we proceed with finding λ from (24), let us demonstrate that a naive perturbation expansion of det Φ in powers of h will fail. Understanding the reason for that will help us guess the forthcoming main result of this section. Setting h → 0 in (37b) and using the asymptotic expression are nonsingular matrices. Hence in the limit h → 0, (24) formally yields lim h→0 λ 2M = 0. This, however, suggests that a more subtle calculation is required and, moreover, that one should expect a relation In what follows we will obtain a quantitative form of this result.

Let us rewrite (41) as
where (· · · ) 1/2 can denote either of the two values of the square root; thus, given (42b), z 1/2 can take on four distinct values. In our numerical simulations, we had Lh = O(1), and hence Remembering from (29b) that λ 2 − 1 = O(1), one has from (43) that Since M = L/h ≫ 1, this has two consequences. First, Second, after (t/h) time steps, the magnitude of the considered modes is proportional to which means that these modes decay in time, i.e., are stable.
To conclude this subsection and prepare the ground for the next one, let us show how (43) can be solved approximately. It follows from (46) that λ 2 lies inside, and very close to, the unit circle in the complex plane. Recall from (29b) that our analysis is valid away from the O(h)-vicinity of the point λ 2 = 1. The factor λ 2 − 1 can be written as and arg(λ 2 − 1) ∈ [π/2, 3π/2] (since |λ| 2 < 1). Then Note that the first term in the exponent is O(1/M ) and hence can be neglected. Since z 1/2 can take on four distinct values (see text after (43)), then (48) yields 4M eigenvalues λ l of the matrix N in (17). (This is consistent with the dimension, 4(M + 1) × 4(M + 1), of that matrix, given that four entries of × are fixed by the BC (14).) In the next subsection we will use direct simulations to verify expression (48) for l ≈ M/4, where arg(λ 2 ) ≈ π and hence |λ 2 − 1| ≈ 2. These eigenvalues correspond to the "middle" of the spectral domain, i.e., kh ≈ π/2 (see the text before (31)), in the periodic problem, shown in Fig. 1(a). This particular choice of arg(λ 2 ) is just the matter of convenience and does not limit our analysis, as long as condition (29b) holds.

Verification of relation (45)
More specifically, we will verify that in the "middle" of the spectrum of the eigenvalue problem (17), (14) one has where the denominator is the value of |λ 2 − 1| 2 for Verification of (49) requires several steps.
Step 1 We begin by extracting modes with λ 2 ≈ −1 from the numerical solution of Eqs. (3) by the MoC-SE. By (26), these modes have ρ ≈ ±i, which according to the text before (31) means that they correspond to harmonics with kh ≈ π/2 in the periodic problem. We first compute the numerical error by subtracting the exact solution (4) from the numerical solution obtained by the MoC-SE (11), (12). To force periodic BC on that error, so as to enable the application of Fourier transform, we multiply the error by a super-Gaussian "window" f (x) = exp −(x/(L/3)) 8 . In what follows we will refer to the so modified error as just "error" (without the quotation marks).
The Fourier spectrum of a representative error is shown in Fig. 2. It confirms that while the modes with kh ≈ 0 and π (i.e., ρ ≈ ±1) grow (see the end of Section 4.1), the rest of the modes decay, in qualitative agreement with (47). As mentioned above, we will focus on the harmonics with kh ≈ π/2.
Step 2 In order to smooth out a noisy profile of the error, we used its value averaged over (2m ave + 1) wavenumbers: where F [· · · ] denotes the Fourier transform, · · · is the Euclidean norm of a vector, and ∆k = Step 3 The "staircase" shape of the curves in this Figure is explained as follows. At point x m = mh of the spatial grid, the size of any any given mode of the error is proportional to |ρ m λ n |: see (18) and (19). Recalling from (26) that ρ ≈ λ −1 or λ, this gives |λ| n−m or |λ| n+m . These expressions correspond to the error's propagation, rather than mere decay, to the right and to the left, respectively. Due to the finite size of the domain, the error eventually exits it. (It is not reflected back into the domain because of the special kind -nonreflecting -BC that we consider here.) As a given "batch" of error leaves the domain, the magnitude of the error inside the domain drops (as in a jump); then some "leftover" error from the middle of the domain moves towards its boundaries, exits it, and so on.
Note that the temporal period of the error's evolution in Fig. 3(a) equals the length of the domain; this is consistent with the above scenario. Therefore, to find λ from the error's evolution plots, we "straighten" each "staircase" curve by taking the data points every t = L time units and perform a linear regression of these data. This is illustrated in Fig. 3(b). Then |λ| was computed from the slope r of the resulting line (the dotted line in Fig. 3(b)) as |λ| = 10 hr .
Step 4 As an initial confirmation of (49), we first note that the slope of ln |λ| 2M vs. ln h should be (2 + 3Lh), which for Lh ≪ 1 is approximately 2. This is verified in Fig. 4(a). Then in Fig. 4 we present a more detailed comparison between |λ| 2M predicted by (49) and that obtained from direct numerical simulations. The agreement between the two is seen to be good. Quantitative comparison (ratio) between |λ| 2M computed by direct numerics and by the theory, i.e., Eq. (49). In both panels the symbols represent measured data, while the lines connecting them are guides for the eye.

Insight into instability suppression in MoC-ME with nonreflecting BC
The following is a heuristic description of how the imposition of nonreflecting BC in the MoC-SE changes the stability of some of the modes compared to the case of periodic BC. After presenting that description, we will apply it to explain why we expect a suppression of the instability for the MoC-ME, even though the instability of the MoC-SE is not suppressed by nonreflecting BC.
The solid line in Fig. 5(a) shows the amplification factor of the MoC-SE with periodic BC (same as in Fig. 1(a)). While all the modes are unstable, the strongest instability occurs for kh = 0 and π. The dashed line shows schematically the locus of the modes for which the imposition of nonreflecting BC suppresses the instability via (43), (47). Recall here that the analysis that led to (43) if valid only for the modes sufficiently "away" from having ρ = ±1, or, equivalently, kh = 0 or π in Fig. 5(a); this is why the "box" there does not occupy the entire interval [0, π]. Figure   5(b) shows schematically the effect of the imposition of nonreflecting BC on the amplification factor: That factor becomes less than unity for the modes inside the "box" in panel (a). Thus, the instability of those modes is suppressed. However, this does not affect the most unstable modes with kh = 0 and π. Therefore, the observed instability of the MoC-SE is the same for the periodic and nonreflecting BC. Applying the same concept to the MoC-ME and using Fig. 1(b), we obtain the result shown in Fig. 6(a). Assuming that the nonreflecting BC suppress instability in the MoC-SE and MoC-ME in a similar way, we see in Fig. 6(b) that the instability of the most unstable modes of the MoC-ME is to be suppressed. The modes with arg(ρ) ≈ 0 and π may remain unstable with nonreflecting BC, but their instability is much weaker than that of modes with arg(ρ) ≈ π/2 (see Section 4 in [1]).
Such a weak instability can be ignored unless the simulation time is very long.

Stability analysis of the MoC-ME
Here we will follow the steps of Sections 3 and 4. We will skip most of the details which are similar for the MoC-ME and MoC-SE and will emphasize only the essential differences. To keep the numeration of subsections in Sections 5 and 4 the same, we will present the setup for the MoC-ME in the forthcoming preamble. The main result of our analysis, indicating suppression of the most strongly unstable modes, will be obtained in Section 5.3. It will then be verified numerically in Section 5.4.
The MoC-ME scheme for Eqs. (3) is: where f ± j are the nonlinear functions on the rhs of (3). The nonreflecting BC have the form (12). The numerical error satisfies the linearized form of (52): with the BC for it being given by (14). Equations (53a) can be cast in the form (15a), where now ; MoC-ME verbatim. We will now implement the three Steps listed in that Section.

5.1
Step 1: Finding ρ(λ) in (20) The characteristic equation (20) for the MoC-ME coincides with that for the MoC-SE, given by in Section 4, here we will not obtain explicit expressions for ρ (2) j as they will not be able to facilitate a quantitative comparison of the analytical and numerical results. However, in our analysis we will need to refer to the presence of such terms, and therefore below we will include them in some of the formulas. In the order O(h) in (27), the expressions for ρ (1) j are the same as for the MoC-SE. Thus, the counterpart of (28) for the MoC-ME is: 2(±) ≡ λ ρ 2(±) .
The reasons for this are as follows.
(i) The numerical instability for modes with ρ ≈ ±1 will still occur. However, it corresponds to the instability in the ODE-and anti-ODE limits of the von Neumann analysis of the periodic problem and thus will be much weaker than the instability (in the periodic problem) for modes with ρ = e ikh ≈ e iπ/2 : see Section 4.5 above and also Section 4 in [1]. As noted there, that weak instability will be inconsequential for moderately long simulation times considered here.
(ii) As we previewed in Section 4.5 and will demonstrate below, our forthcoming analysis, based on (55), will succeed in predicting the instability suppression for the modes with ρ ≈ e iπ/2 , which are the most unstable in the periodic problem. In other words, that analysis will describe the most important phenomenon that motivated this entire study.
(iii) The analysis for the case λ 2 ≈ 1 is much more technical than the one given below. Since its results will not affect the main conclusions of the analysis based on (55), then it seems appropriate not to carry it out here in order not to distract the reader's attention. Thus, we will proceed assuming the condition (29b).

5.2
Step 2: Finding ξ in (20) As in Section 4.2, we will outline the derivation of only ξ 1(+) and will simply state the results for the other three eigenvectors. Similarly to (33), we seek We will compute only u (0) and v (1) . The other terms in (56) will enter some of the expressions in Section 5.3. We will not need their precise values, but their form will require u (1), (2) , v (2) to be defined; this is why we listed them in (56).
Similarly to Section 4. instead of (16b). It is the difference in the structure of these expressions that will cause a substantial difference between the results for the MoC-ME and MoC-SE. The vector u (0) in (56) is found to still be given by (35a), but the counterpart of (35b) is now: The numerator in the fraction above, which is absent in (35b), will turn out to be "responsible" for the aforementioned substantial difference.

Step 3: Finding eigenvalues of N in (17)
Here we will obtain the main result of Section 5, which will be given by relation (63).
The equation from which the eigenvalues of the amplification matrix N are to be found is (24), where matrix Φ is defined by (23b). Substituting there expressions (55) and (58), we obtain a counterpart of (37). For reasons to be explained soon, in writing it, we will specify the first two h-orders of the entries (and will also slightly change the notations): 2(+) , u Note that the main-order entries above are the same as their counterparts in (37), except that Φ 12 , Φ (1) 21 have the extra factor (λ 2 + 1)/2. It is this factor that will substantially change the behavior of modes with ρ ≈ e ±iπ/2 in the MoC-ME compared to those modes in the MoC-SE. Therefore, we will now discuss what difference in the analysis this factor makes.
Most importantly, for modes with ρ ✚ ✚ ≈ ± i (and for those with ρ ✚ ✚ ≈ ± 1; see Section 5.1), the main result of Section 4, Eq. (45): |λ| M = O(h), remains valid. With time, these modes decay as O(h t/L ) (see (47)), just as they do for the MoC-SE. However, the modes in the narrow interval where will be shown to decay even faster, as O(h 2t/L ). This, of course, will not change the fact that these modes are stable in the MoC-ME, as they are in the MoC-SE, but it will change the Fourier spectrum of the numerical error, as we will show below.
Let us explain why condition (60) affects |λ| M . Under this condition, the two terms in each of the off-diagonal entries in (59a) both become O(h 2 ): see (59c) and (59d). Therefore, in this case it is more appropriate to rewrite (59a) as where by introducing the tilde-notations in the second row, we have also accounted for the λ Mdependence. The calculation of Φ 12 , Φ 21 would require finding the higher-order correction terms, denoted by ξ (2) j(±) in (58). This would be quite a tedious task, which would not provide any additional qualitative insight into the evolution of the modes (60). Therefore, we will not carry it out but will proceed with the general form (61). Furthermore, we will neglect the O(h)-terms in the diagonal entries of (61) because Φ  (59e)). Then, repeating the steps that led to (39), we obtain: While we do not compute the explicit form of two out of the four matrices in (62), we still know that their entries are of order one. Therefore, (62) can hold only when the coefficient in front of the four-matrix product is also of order one, or, equivalently, when This is the result for modes (60) of the MoC-ME, announced earlier. It predicts a faster temporal decay of these modes than of the rest of the modes, satisfying (45). We will now verify this by direct numerical simulations. Note that unlike in the MoC-SE case (see (49)), here we will not attempt to resolve the numeric coefficient in the O(h 4 ) notation.

Verification of (63) for modes (60)
We will perform this verification in two different ways. For both ways, we will make the numerical error periodic in space by multiplying it with a "window", as in Step 1 of Section 4.4. This will allow us to work with the Fourier spectrum of the so modified error.

Dependence of spectral components of error on time
The Fourier spectra of the error at equidistant times are shown in Fig. 7. As we announced after (63), modes near kh = π/2 decay faster than those farther away from kh = π/2. This faster decay of modes with kh ≈ π/2 is manifested by the deepening "dip" in the spectra. More specifically, a comparison of (63) and (45) shows that the modes near kh = π/2 decay twice as fast as the modes farther away from kh = π/2. Indeed, (47) (which in the MoC-ME holds for modes with kh ✚ ✚ ≈ π/2) and (63), for modes with kh ≈ π/2, yield: This implies that for any two moments of time, t 1,2 , one is to have: This observation can be employed to verify the theoretical predictions (64) against the numerical results presented in Fig. 7. Using the numbers listed in that Figure

Dependence of spectral components of error on h
We will now confirm (63) in yet another way. Recall that when we extract |λ| from numerically obtained spectra using formula (51), we perform averaging over m ave nodes around the wavenumber of interest. If we focus on such a narrow vicinity of kh = π/2 that all modes in the avegaging box satisfy condition (60) (which requires m ave ≪ M/4; see the narrower box in Fig. 7), then we expect that the numerically calculated |λ| will satisfy (63). That relation is equivalent to where the "const" may actually be a slow function of ln h, as in (49), and tend to a "true" constant for h → 0. Thus, the slope of the curve ln |λ| 2M vs. ln h should be (approximately) 4.
On the other hand, when m ave in (51) becomes comparable to M/4, the averaging interval, represented by the wider box in Fig. 7, includes mostly the modes satisfying (45). Then the slope of ln |λ| 2M vs. ln h should be (approximately) 2. This is confirmed by Fig. 8,   6 Instability of modes with ρ ≈ ±i for the MoC-LF In Sections 4 and 5 we have seen that modes in the "middle" of the spectrum (i.e. having ρ ≈ ±i), which are unstable in the MoC-SE and MoC-ME schemes with periodic BC (see Fig. 1), are "made" stable by the nonreflecting BC. In contrast, as we announced in the Introduction, direct numerical simulations of the MoC-LF show that the scheme exhibits very similar instability (see Fig. 1(c)) for both periodic and nonreflecting BC. In this section we will explain why this is so. Since the instability in the periodic case exists only for the modes with kh ≈ π/2, or, equivalently, satisfying (60), we will proceed to find the amplification factor, λ, only for these modes.
Since the numerical error satisfies the linearized equations, we will state the BC only for them.
The MoC-LF scheme involves three time levels, and therefore the BC (14) In writing the last step in each of (69a) and (69b), we have used (14). Thus, (14) and (69) form the BC for the MoC-LF.
We will now follow the steps of Section 3 and arrive at a counterpart of the characteristic polynomial (24), whose roots will give the amplification factor λ of the MoC-LF. In Sections 6.1-6.3 we will follow the steps of finding λ as outlined at the end of Section 3.
Since we are only interested in the modes satisfying (60), these roots have the form: with β j = O(1), which will be found in Section 6.1. The eigenvectors corresponding to ρ (±, j) will be denoted by ξ (±, j) . Seeking, similarly to (21), where C (r,j) are constants, one rewrites the BC (14) as: and the BC (69) as: In principle, arranging (73) and (74) into a matrix equation like (23b) will then allow one to determine the amplification factor |λ| via a counterpart of (24). However, before we proceed along these lines, we will simplify the form of (74) in order to simplify subsequent calculations.
In Section 6.1 we will show that there is a set of modes for which β j ∈ R, as this will suffice for our purpose of explaining the instability of the MoC-LF. In what follows we will consider only such modes. Then, for h ≪ 1, where we have used M = L/h. Next, later on we will show that all entries of ξ (r,j) are O(1).
Therefore, coefficients multiplying C (r,j) on the lhs of (74a where C (r) = [C (r,1) , . . . , C (r,4) ] T and Φ (r) are defined as in (23b) using the components of the eigenvectors ξ (r,j) . Similarly, (76) can be written as Together, (77) imply that which seems to impose two equations on λ. However, in Sections 6.2 and 6.3 we will show that ξ (−,j) are related to ξ (+,j) in such a way that the two equations in (78) are, in fact, the same. Therefore, solving, say, the first of them will yield the values of λ. Our goal, to be achieved in Section 6.3, will be to show that these values are very close to the eigenvalues of the unstable modes in the periodic problem.
Let us now note that the discarding of certain O(h) terms, which has reduced (74) to (76), implies that the order of approximation with which one computes the boundary values (s + ) 1 and (s − ) M −1 is inconsequential for the instability of the MoC-LF. This justifies our earlier corresponding statement, found before (69).
We will now follow the three steps, listed at the end of Section 3, of solving the characteristic equation (78). In our presentations we will focus on their differences from the corresponding calculations in Sections 4 and 5. 6.1 Step 1: Finding ρ(λ) in (70) The characteristic polynomial for (70) is: In the order O(1), it has four double roots: 8 Recall that we are specifically looking for ρ ≈ ±i (see (71)), in which case (80a) implies λ ≈ ±i and then Due to this degeneracy, a perturbation expansion in h should proceed not similarly to (27), but instead similarly to (30b). Selecting the case λ ≈ +i, we seek solutions of (79) in the form here α and β are O(1) quantities to be determined below. A similar ansatz can be used for λ ≈ −i.
Below we will present details for the '+' sign in (81) and will only state the final result for the '−' sign. Substituting (81) into (79), in the lowest nontrivial order, O(h 4 ), one finds: whence For the '−' sign in (81), one obtains the same relation (82a). Note that (81) and (82) are to be interpreted as follows: For a given λ (or, equivalently, α), which will be found in Section 6.3, one has to find eight modes (i.e., eight values of ρ, with four corresponding to each of the two signs in (81)), which satisfy the nonreflecting BC (14), (69). This appears to be converse to the procedure that one follows in the case of periodic BC. Namely, there, one first finds the modes (Fourier harmonics) exp [ikhm] with k that are consistent with the periodic BC, and then for those modes finds λ; see, e.g., Section 5 in [1]. Now recall our goal: We want to explain why the instability of the MoC-LF is (almost) the same for the periodic and nonreflecting BC. This means that for the nonreflecting BC, we need to focus on the modes that are counterparts of the Fourier harmonics exp[ikhm], m = 0, 1, . . . , M in the periodic case. For the ρ in (81) to mimic exp[ikh] (with kh ≈ π/2), one needs to have β ∈ R.
For α ∈ R (which will be justified by our results in Section 6.3), this occurs for see (82b). Therefore, in what follows we will consider only the range (83) of values of α.

6.2
Step 2: Finding ξ in (70) Substituting (81) in (70), one finds that all terms in the order O(1) cancel out, and then in the order O(h) one finds: Using the explicit form of Γ 0,1 and Ω 0,1 from (16b) and (7) and relation (82a) between α and β, one finds the eigenvector ξ. Below we state the result for both signs in (81): As explained before (83), we consider the case α, β ∈ R; then (85) yields: where σ 3 = diag(1, −1) is a Pauli matrix. Also, for future use, note that from (75) it follows that where we have used the fact that all β j that we consider are real.

Step 3: Finding amplification factor of MoC-LF from (78)
Here we will obtain the main result of Section 6. Unlike in the previous two sections, here λ does not have an analytical expression, but can be easily found by graphing an analytically defined function.
Recall that our goal is to find those λ, or, equivalently, α in (81), which make (78) hold true.
Since the eigenvalue problem (70), (73), (76) has a finite dimension, it is to be expected that the set of its eigenvalues (or, α), is discrete. 9 Recall from the discussion after (82) that it is the value of λ that determines the allowed modes. Thus, we are to compute det Φ (+) while considering β j to be functions of α as per (82b)  7 Summary and discussion

Summary and an example of an unusual numerical instability
The main contributions of this study are as follows. First, we have presented an approach to carry out a stability analysis of the MoC for a PDE system with constant coefficients while placing emphasis on its having non-periodic BC. The two ingredients of this approach have been known before. One ingredient was Zió lko's paper [6] presenting the problem in terms of a block-tridiagonal  our analysis was the well-known method of finding eigenvalues of a band Toeplitz matrix; see, e.g., [11,13]. In Section 3 we outlined how these two ingredients can be combined to enable a stability analysis of the MoC applied to a constant-coefficient PDE with non-periodic BC. While we have considered a special kind of such BC -nonreflecting, -the approach can be straightforwardly generalized to any other kind of BC.
Second, we have applied our analysis to three "flavors" of the MoC, where the ODEs along the characteristics were solved by the SE, ME, and LF methods. For all three of them, a previous study [1] has found, for periodic BC, a conspicuous instability for modes in the "middle" of the Fourier spectrum, i.e., for kh ≈ π/2. Here, we have shown that this instability is suppressed for the MoC-SE and MoC-ME schemes when nonreflecting BC are applied. As we announced in the Introduction, this finding contradicts a statement made in textbooks [3]- [5] that an instability of a scheme with periodic BC always implies its instability for other types of BC. We will provide a conceptual explanation of this contradiction in Section 7.3.
see (47). This holds for all modes with |kh| ✚ ✚ ≈ 0, π in the MoC-SE and for modes with |kh| ✚ ✚ ≈ 0, π, π/2 10 after "straightening out" the staircase shape; see Step 3 in Sec. 4.4 in the MoC-ME. The modes in the MoC-ME with |kh| ≈ π/2 decay in time even faster: as follows from comparison of (45) with (63). Thus, the modes' decay rate depends both on the discretization step h and on the domain length L in an unusual way. To our knowledge, neither of these features has been previously reported for other numerical schemes applied to PDEs with spatially-constant coefficients. 11 Moreover, our analysis has predicted yet another unusual phenomenon: Merely by increasing the product Lh, one can turn stable modes of the MoC-SE with nonreflecting BC into unstable ones; see (43) or (49). We illustrate this via slightly altering the case shown in Fig. 2, which has stable modes in the "middle" of the spectrum. The only difference between that case and that shown in Fig. 10(a) is that in the latter L is four times greater. According to (49), the magnitude of the modes in the "middle" of the spectrum is to increase by a factor of approximately 8 over time t = L, and this agrees well with the results seen in Fig. 10 to expect that such a factor is indeed present in (63) (for modes with ρ ≈ ±i) and in (45) (for modes with ρ ✚ ✚ ≈ ± i). This expectation is borne out by the result shown in Fig. 10(b), which should be compared with the stable evolution of the MoC-ME error in Fig. 7. By trial-and-error we have determined that for L = 300 the scheme is still stable (as it is also in the case shown in Fig. 7), while for L = 400 the growth rate, while positive, is quite small. Therefore, we showed the case with L = 600 and hence a higher growth rate. Let us note that we are not aware of any previous reports where merely changing the length of the spatial domain would alter stability of some of the modes and, as in the MoC-ME example above, even of the entire scheme.

On validity of our analysis for other models
Let us emphasize that the second and third contributions stated above are specific to the energypreserving PDEs whose linearized form is (1). Such PDEs arise not only in the theory of birefringent optical fibers, but, more generally, in any theory that involves two coupled waves (or two distinct, e.g., counter-propagating groups of waves); see Section 2 of Ref. [1] for a more detailed discussion and references. Thus, these hyperbolic PDEs are fairly common. In Appendix C we present a setup of another model -the Gross-Neveu model from the relativistic field theory -whose linearized  (The peak near kh = 0 in Fig. 11(a) is not related to numerical instability; its origin is explained in [1].) If we now simulate the soliton of the Gross-Neveu model with the same parameters as for Fig. 11(a) but with nonreflecting BC (101), the Fourier spectrum of the numerical error acquires the shape shown in Fig. 11(b). The error is seen to be several orders of magnitude higher than the initial error (which has the order 10 −10 in the spectral domain). However, this occurs not because of its exponential growth, but because a non-periodicity at the boundaries, created by the nonreflecting BC (see Fig. 11(c)), led to a spectrum that decays with k very slowly. To confirm that the error does not grow in time (but, in fact, slightly decays), we presented the spectra for t = 1000 and t = 5000 in Fig. 11(b). Thus, as predicted by our theory for the simpler model based on Eqs. (3), (4), the MoC-ME scheme for the soliton of the Gross-Neveu model is also stabilized by nonreflecting BC.
The reader may notice that the spectral shape of the error in Fig. 11(b) is missing the "dip" around kh = π/2, which is seen in Fig. 7 and which was predicted by our analysis of the constantcoefficient model in Section 5. This occurs due to the same reason, related to the non-periodicity of the numerical solution at the boundaries, as illustrated in Fig. 11(c). In order to recover the shape of the spectrum as seen in Fig. 7, we have made the numerical error periodic by multiplying 12 which is equivalent to Fig. 5(c)   it (but not the numerical solution) by a spatial super-Gaussian "window", as described in Step 1 of Section 4.4. 13 In Fig. 12 we show the evolution of the so modified, periodic error. Figure 12(a) shows the case for L = 128 and h = L/2 12 ≈ 0.031. Note that we had to increase both L and h compared to the case shown in Fig. 11 since otherwise, in agreement with estimates (89), the error would reach the size of the computer round-off error very quickly. For the same reason, in the initial condition we added to the soliton a white noise of the size 10 −6 as opposed to 10 −12 .
The profile with a "dip", seen in Fig. 7, is also observed in Fig. 12(a). Also in agreement with our analysis for the simpler model, the "dip" decreases approximately twice as fast as the "wings" of the spectrum: see Fig. 12(b).
Finally, as we pointed out at the end of Section 7.1, it should be possible to make the MoC-ME unstable again by simply increasing the length L of the computational domain. This, indeed, also holds true for the soliton of the Gross-Neveu model. We chose the same L = 600 as in Fig. 10(b), kept the same h as in Fig. 12(a), and added to the initial soliton a white noise of magnitude 10 −12 (since the error will now grow, not decay). The corresponding Fourier spectra at different times are shown in Fig. 12(c). The shapes of the spectra for our simpler model (3), (4) (Fig. 10(b)) and for the soliton of the Gross-Neveu model (Fig. 12(c)) are seen to be very similar.
Let us point out that this dependence of the stability of the MoC-ME scheme on L for the spatially localized solution of (99) provides an indirect evidence that conclusions of our analysis are valid for an even broader class of models. More specifically, a sufficient condition for the same (in)stability behavior of a numerical scheme appears to be the presence of the constant, linear coupling between the forward-and backward-propagating components of the solution, given by 13 In fact, the error used in plotting Fig. 7 Figure 12: (a) Fourier spectra of the periodic numerical error of the MoC-ME with nonreflecting BC for the same soliton solution as in Fig. 11, but for L = 128; see text for details. The times increase from top to bottom with the increment t incr = L = 128, as discussed in Section 4.4. Let us now explain the statement made two sentences above. From the analysis in Section 4.3 (and a similar analysis in Section 5.3) it follows that the dependence of λ on h occurs due to the linearization matrix P in (7) (and, more generally, in (1)) being nonzero. For Eqs. (99), entries of the corresponding P have two contributions: from the nonlinear terms and from the linear terms on the rhs. Now note that the nonlinear terms are essentially nonzero only in an interval of width O(1) where the soliton is essentially nonzero. Thus, their contributions to entries of P are nonzero only in that O(1)-long interval, and hence they cannot affect λ by a factor that is related to the length L ≫ 1 of the entire spatial domain. On the other hand, the linear terms on the rhs of (99) make spatially constant contributions to entries of P and thus must be solely responsible for the changes of the (in)stability behavior of the numerical scheme with L.
Thus, we have shown that our analysis for a constant-coefficient, non-dissipative hyperbolic system whose linearization is given by (1) has predictive value for another non-dissipative hyperbolic model with non-constant coefficients. An extension of our stability result to dissipative hyperbolic PDEs, or non-dissipative ones whose linearization in some way substantially differs from (1) 14 or to more general BC (e.g., those considered in [6]), remains an open problem. However, the approach will remain the same as that outlined in Section 3. 14 E.g., there are more than two groups of waves propagating with distinctly different group velocities. 7.3 Why von Neumann analysis may incorrectly "predict" instability for nonperiodic BC We will now explain why the instability of the modes with kh ✚ ✚ ≈ 0, π in the MoC-SE and MoC-ME with periodic BC did not imply an instability of their counterparts when the BC are changed to nonreflecting. Let us recall that the opposite statement was made in textbooks [3]- [5], and thus, our contradicting that common knowledge may be the most unexpected conclusion of this study.
We will begin by stating the reason behind this contradiction in general terms. Our system (1) (or (6a)) has the same form as that considered in [3]- [5]. However, one of the underlying assumptions made in [3]- [5] does not actually hold for the spatial modes of the MoC-SE and MoC-ME.
We will now explain this in detail. In [3]- [5], it was assumed that all modes in a problem with non-periodic BC fell into two groups. Modes in one group are localized near the boundary (or either of the boundaries). These modes become exponentially small in the "bulk" of the domain.
In particular, the mode localized near the left boundary does not "feel" the right boundary, and vice versa; see Fig. 13(a). Modes in the other group do not exhibit a monotonic decay away from the boundary(ies) of the spatial domain. In sufficiently large domains, such modes resemble Fourier harmonics (i.e., modes of the periodic problem) in the "bulk" of the domain; see Fig. 13(a).
(Near the boundary(ies), their profile would, of course, be modified to satisfy the BC.) Given this resemblance, one can reasonably expect the (in)stability of modes in this group to be similar to that of the corresponding Fourier harmonics. Therefore, if any of the modes in the periodic problem is unstable, so should be the corresponding mode from this "second group" of the non-periodic problem. It does not matter in this case whether modes from the "first group" are stable or not, for the numerical scheme is already unstable due to modes from the "second group". This is the reason why it was stated in [3]- [5] that an instability of a scheme predicted by the von Neumann analysis implies an instability of that scheme with arbitrary BC. The reason why the above statement does not hold for the MoC-SE and MoC-ME is schemat-ically shown in Fig. 13(b). Namely, the modes in the problem with nonreflecting BC do not split into the two groups described in the previous paragraph. Thus, the underlying assumption that led to the conclusion in textbooks [3]- [5] does not hold for the MoC-SE and MoC-ME.
Indeed, the spatial dependence of the eigenmodes in these methods is given by (18) as depicted in Fig. 13(b). In particular, the mode decaying away from the left boundary is not negligibly small at the right boundary no matter how large the spatial domain is; and vice versa. Even more importantly, there are no modes in the problem with nonreflecting BC that resemble Fourier harmonics in the "bulk" of the domain. Indeed, all modes either decay or grow monotonically in magnitude. This exponential decay or growth 15 appears to make these modes sufficiently different from Fourier harmonics, so that a prediction of the von Neumann analysis about instability of Fourier harmonics has no bearing for the modes of the problem with nonreflecting BC.
As a minor clarification about Fig. 13(b), let us note that for a given λ, there are two spatially growing and two decaying modes, according to the number of roots ρ in (28). We stress that, according to (23a), all four of those modes are required to satisfy the nonreflecting BC at each boundary.

When can von Neumann analysis correctly predict instability for a MoC scheme with non-periodic BC
Finally, let us give a qualitative reason why the above explanation did not apply to the MoC-LF scheme, for which, as we showed in Section 6, the nonreflecting BC did not suppress the instability that existed for periodic BC. The reason is precisely that (unstable) modes of the "second group", which resemble Fourier harmonics in the "bulk" of the domain, do exist in the MoC-LF scheme.
This, in turn, occurs due to the following difference between the LF and Euler ODE solvers. Unlike the latter, the LF solver involves more than two time levels. This leads, in the limit h → 0, to the existence of "parasitic" amplification factors, ρ parasitic = −λ, −λ −1 , in addition to the true ones, ρ true = λ, λ −1 . For ρ = ±i, which is where the modes of the periodic problem are found to be unstable, each of the "parasitic" factors coincides with a true one: −λ = λ −1 , −λ −1 = λ. In this situation and for 0 < h ≪ 1, it is possible to have |λ| > 1 but |ρ| = 1; (92) see (81) and the paragraph before (83), where we specified that α, β ∈ R. Condition (92), which is possible for the MoC-LF but impossible for the MoC-SE and MoC-ME (see (28)), means that for the MoC-LF, there do exist modes which resemble Fourier harmonics in the "bulk" of the spatial domain (|ρ| = 1) and are unstable (|λ| > 1). This is why the von Neumann analysis for the MoC-LF predicts nearly the same instability that exists for that scheme with non-periodic BC. Thus, a necessary ingredient for this to occur for an arbitrary MoC scheme is that its ODE solver must have some ρ parasitic that coincides with a ρ true for h → 0.